ScaLAPACK  2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
pdsyev.f
Go to the documentation of this file.
00001       SUBROUTINE PDSYEV( JOBZ, UPLO, N, A, IA, JA, DESCA, W,
00002      $                   Z, IZ, JZ, DESCZ, WORK, LWORK, INFO )
00003 *
00004 *  -- ScaLAPACK routine (version 1.7) --
00005 *     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
00006 *     and University of California, Berkeley.
00007 *     May 25, 2001
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          JOBZ, UPLO
00011       INTEGER            IA, INFO, IZ, JA, JZ, LWORK, N
00012 *     ..
00013 *     .. Array Arguments ..
00014       INTEGER            DESCA( * ), DESCZ( * )
00015       DOUBLE PRECISION   A( * ), W( * ), WORK( * ), Z( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  PDSYEV computes all eigenvalues and, optionally, eigenvectors
00022 *  of a real symmetric matrix A by calling the recommended sequence
00023 *  of ScaLAPACK routines.
00024 *
00025 *  In its present form, PDSYEV assumes a homogeneous system and makes
00026 *  no checks for consistency of the eigenvalues or eigenvectors across
00027 *  the different processes.  Because of this, it is possible that a
00028 *  heterogeneous system may return incorrect results without any error
00029 *  messages.
00030 *
00031 *  Notes
00032 *  =====
00033 *  A description vector is associated with each 2D block-cyclicly dis-
00034 *  tributed matrix.  This vector stores the information required to
00035 *  establish the mapping between a matrix entry and its corresponding
00036 *  process and memory location.
00037 *
00038 *  In the following comments, the character _ should be read as
00039 *  "of the distributed matrix".  Let A be a generic term for any 2D
00040 *  block cyclicly distributed matrix.  Its description vector is DESCA:
00041 *
00042 *  NOTATION        STORED IN      EXPLANATION
00043 *  --------------- -------------- --------------------------------------
00044 *  DTYPE_A(global) DESCA( DTYPE_) The descriptor type.
00045 *  CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
00046 *                                 the BLACS process grid A is distribu-
00047 *                                 ted over. The context itself is glo-
00048 *                                 bal, but the handle (the integer
00049 *                                 value) may vary.
00050 *  M_A    (global) DESCA( M_ )    The number of rows in the distributed
00051 *                                 matrix A.
00052 *  N_A    (global) DESCA( N_ )    The number of columns in the distri-
00053 *                                 buted matrix A.
00054 *  MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
00055 *                                 the rows of A.
00056 *  NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
00057 *                                 the columns of A.
00058 *  RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
00059 *                                 row of the matrix A is distributed.
00060 *  CSRC_A (global) DESCA( CSRC_ ) The process column over which the
00061 *                                 first column of A is distributed.
00062 *  LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
00063 *                                 array storing the local blocks of the
00064 *                                 distributed matrix A.
00065 *                                 LLD_A >= MAX(1,LOCr(M_A)).
00066 *
00067 *  Let K be the number of rows or columns of a distributed matrix,
00068 *  and assume that its process grid has dimension p x q.
00069 *  LOCr( K ) denotes the number of elements of K that a process
00070 *  would receive if K were distributed over the p processes of its
00071 *  process column.
00072 *  Similarly, LOCc( K ) denotes the number of elements of K that a
00073 *  process would receive if K were distributed over the q processes of
00074 *  its process row.
00075 *  The values of LOCr() and LOCc() may be determined via a call to the
00076 *  ScaLAPACK tool function, NUMROC:
00077 *          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
00078 *          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
00079 *
00080 *  Arguments
00081 *  =========
00082 *
00083 *     NP = the number of rows local to a given process.
00084 *     NQ = the number of columns local to a given process.
00085 *
00086 *  JOBZ    (global input) CHARACTER*1
00087 *          Specifies whether or not to compute the eigenvectors:
00088 *          = 'N':  Compute eigenvalues only.
00089 *          = 'V':  Compute eigenvalues and eigenvectors.
00090 *
00091 *  UPLO    (global input) CHARACTER*1
00092 *          Specifies whether the upper or lower triangular part of the
00093 *          symmetric matrix A is stored:
00094 *          = 'U':  Upper triangular
00095 *          = 'L':  Lower triangular
00096 *
00097 *  N       (global input) INTEGER
00098 *          The number of rows and columns of the matrix A.  N >= 0.
00099 *
00100 *  A       (local input/workspace) block cyclic DOUBLE PRECISION array,
00101 *          global dimension (N, N), local dimension ( LLD_A,
00102 *          LOCc(JA+N-1) )
00103 *
00104 *          On entry, the symmetric matrix A.  If UPLO = 'U', only the
00105 *          upper triangular part of A is used to define the elements of
00106 *          the symmetric matrix.  If UPLO = 'L', only the lower
00107 *          triangular part of A is used to define the elements of the
00108 *          symmetric matrix.
00109 *
00110 *          On exit, the lower triangle (if UPLO='L') or the upper
00111 *          triangle (if UPLO='U') of A, including the diagonal, is
00112 *          destroyed.
00113 *
00114 *  IA      (global input) INTEGER
00115 *          A's global row index, which points to the beginning of the
00116 *          submatrix which is to be operated on.
00117 *
00118 *  JA      (global input) INTEGER
00119 *          A's global column index, which points to the beginning of
00120 *          the submatrix which is to be operated on.
00121 *
00122 *  DESCA   (global and local input) INTEGER array of dimension DLEN_.
00123 *          The array descriptor for the distributed matrix A.
00124 *          If DESCA( CTXT_ ) is incorrect, PDSYEV cannot guarantee
00125 *          correct error reporting.
00126 *
00127 *  W       (global output) DOUBLE PRECISION array, dimension (N)
00128 *          If INFO=0, the eigenvalues in ascending order.
00129 *
00130 *  Z       (local output) DOUBLE PRECISION array,
00131 *          global dimension (N, N),
00132 *          local dimension ( LLD_Z, LOCc(JZ+N-1) )
00133 *          If JOBZ = 'V', then on normal exit the first M columns of Z
00134 *          contain the orthonormal eigenvectors of the matrix
00135 *          corresponding to the selected eigenvalues.
00136 *          If JOBZ = 'N', then Z is not referenced.
00137 *
00138 *  IZ      (global input) INTEGER
00139 *          Z's global row index, which points to the beginning of the
00140 *          submatrix which is to be operated on.
00141 *
00142 *  JZ      (global input) INTEGER
00143 *          Z's global column index, which points to the beginning of
00144 *          the submatrix which is to be operated on.
00145 *
00146 *  DESCZ   (global and local input) INTEGER array of dimension DLEN_.
00147 *          The array descriptor for the distributed matrix Z.
00148 *          DESCZ( CTXT_ ) must equal DESCA( CTXT_ )
00149 *
00150 *  WORK    (local workspace/output) DOUBLE PRECISION array,
00151 *          dimension (LWORK)
00152 *          Version 1.0:  on output, WORK(1) returns the workspace
00153 *          needed to guarantee completion.
00154 *          If the input parameters are incorrect, WORK(1) may also be
00155 *          incorrect.
00156 *
00157 *          If JOBZ='N' WORK(1) = minimal=optimal amount of workspace
00158 *          If JOBZ='V' WORK(1) = minimal workspace required to
00159 *             generate all the eigenvectors.
00160 *
00161 *
00162 *  LWORK   (local input) INTEGER
00163 *          See below for definitions of variables used to define LWORK.
00164 *          If no eigenvectors are requested (JOBZ = 'N') then
00165 *             LWORK >= 5*N + SIZESYTRD + 1
00166 *          where
00167 *             SIZESYTRD = The workspace requirement for PDSYTRD
00168 *                         and is MAX( NB * ( NP +1 ), 3 * NB )
00169 *          If eigenvectors are requested (JOBZ = 'V' ) then
00170 *             the amount of workspace required to guarantee that all
00171 *             eigenvectors are computed is:
00172 *
00173 *             QRMEM = 2*N-2
00174 *             LWMIN = 5*N + N*LDC + MAX( SIZEMQRLEFT, QRMEM ) + 1
00175 *
00176 *          Variable definitions:
00177 *             NB = DESCA( MB_ ) = DESCA( NB_ ) =
00178 *                  DESCZ( MB_ ) = DESCZ( NB_ )
00179 *             NN = MAX( N, NB, 2 )
00180 *             DESCA( RSRC_ ) = DESCA( RSRC_ ) = DESCZ( RSRC_ ) =
00181 *                              DESCZ( CSRC_ ) = 0
00182 *             NP = NUMROC( NN, NB, 0, 0, NPROW )
00183 *             NQ = NUMROC( MAX( N, NB, 2 ), NB, 0, 0, NPCOL )
00184 *             NRC = NUMROC( N, NB, MYPROWC, 0, NPROCS)
00185 *             LDC = MAX( 1, NRC )
00186 *             SIZEMQRLEFT = The workspace requirement for PDORMTR
00187 *                           when it's SIDE argument is 'L'.
00188 *
00189 *          With MYPROWC defined when a new context is created as:
00190 *             CALL BLACS_GET( DESCA( CTXT_ ), 0, CONTEXTC )
00191 *             CALL BLACS_GRIDINIT( CONTEXTC, 'R', NPROCS, 1 )
00192 *             CALL BLACS_GRIDINFO( CONTEXTC, NPROWC, NPCOLC, MYPROWC,
00193 *                                  MYPCOLC )
00194 *
00195 *          If LWORK = -1, the LWORK is global input and a workspace
00196 *          query is assumed; the routine only calculates the minimum
00197 *          size for the WORK array.  The required workspace is returned
00198 *          as the first element of WORK and no error message is issued
00199 *          by PXERBLA.
00200 *
00201 *  INFO    (global output) INTEGER
00202 *          = 0:  successful exit
00203 *          < 0:  If the i-th argument is an array and the j-entry had
00204 *                an illegal value, then INFO = -(i*100+j), if the i-th
00205 *                argument is a scalar and had an illegal value, then
00206 *                INFO = -i.
00207 *          > 0:  If INFO = 1 through N, the i(th) eigenvalue did not
00208 *                converge in DSTEQR2 after a total of 30*N iterations.
00209 *                If INFO = N+1, then PDSYEV has detected heterogeneity
00210 *                by finding that eigenvalues were not identical across
00211 *                the process grid.  In this case, the accuracy of
00212 *                the results from PDSYEV cannot be guaranteed.
00213 *
00214 *  Alignment requirements
00215 *  ======================
00216 *
00217 *  The distributed submatrices A(IA:*, JA:*) and Z(IZ:IZ+M-1,JZ:JZ+N-1)
00218 *  must verify some alignment properties, namely the following
00219 *  expressions should be true:
00220 *
00221 *  ( MB_A.EQ.NB_A.EQ.MB_Z .AND. IROFFA.EQ.IROFFZ .AND. IROFFA.EQ.0 .AND.
00222 *    IAROW.EQ.IZROW )
00223 *  where
00224 *  IROFFA = MOD( IA-1, MB_A ) and ICOFFA = MOD( JA-1, NB_A ).
00225 *
00226 *  =====================================================================
00227 *
00228 *  Version 1.4 limitations:
00229 *     DESCA(MB_) = DESCA(NB_)
00230 *     DESCA(M_) = DESCZ(M_)
00231 *     DESCA(N_) = DESCZ(N_)
00232 *     DESCA(MB_) = DESCZ(MB_)
00233 *     DESCA(NB_) = DESCZ(NB_)
00234 *     DESCA(RSRC_) = DESCZ(RSRC_)
00235 *
00236 *     .. Parameters ..
00237       INTEGER            BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
00238      $                   MB_, NB_, RSRC_, CSRC_, LLD_
00239       PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
00240      $                   CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
00241      $                   RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
00242       DOUBLE PRECISION   FIVE, ONE, TEN, ZERO
00243       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0,
00244      $                     TEN = 10.0D+0, FIVE = 5.0D+0 )
00245       INTEGER            IERREIN, IERRCLS, IERRSPC, IERREBZ, ITHVAL
00246       PARAMETER          ( IERREIN = 1, IERRCLS = 2, IERRSPC = 4,
00247      $                   IERREBZ = 8, ITHVAL = 10 )
00248 *     ..
00249 *     .. Local Scalars ..
00250       LOGICAL            LOWER, WANTZ
00251       INTEGER            CONTEXTC, CSRC_A, I, IACOL, IAROW, ICOFFA, 
00252      $                   IINFO, INDD, INDD2, INDE, INDE2, INDTAU, 
00253      $                   INDWORK, INDWORK2, IROFFA, IROFFZ, ISCALE, 
00254      $                   IZROW, J, K, LDC, LLWORK, LWMIN, MB_A, MB_Z, 
00255      $                   MYCOL, MYPCOLC, MYPROWC, MYROW, NB, NB_A, NB_Z,
00256      $                   NP, NPCOL, NPCOLC, NPROCS, NPROW, NPROWC, NQ, 
00257      $                   NRC, QRMEM, RSRC_A, RSRC_Z, SIZEMQRLEFT, 
00258      $                   SIZESYTRD
00259       DOUBLE PRECISION   ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, 
00260      $                   SMLNUM
00261 *     ..
00262 *     .. Local Arrays ..
00263       INTEGER            DESCQR( 9 ), IDUM1( 3 ), IDUM2( 3 )
00264 *     ..
00265 *     .. External Functions ..
00266       LOGICAL            LSAME
00267       INTEGER            INDXG2P, NUMROC, SL_GRIDRESHAPE
00268       DOUBLE PRECISION   PDLAMCH, PDLANSY
00269       EXTERNAL           LSAME, NUMROC, PDLAMCH, PDLANSY,
00270      $                   SL_GRIDRESHAPE
00271 *     ..
00272 *     .. External Subroutines ..
00273       EXTERNAL           BLACS_GRIDEXIT, BLACS_GRIDINFO, CHK1MAT, DCOPY,
00274      $                   DESCINIT, DSCAL, DSTEQR2, PCHK1MAT, PCHK2MAT,
00275      $                   PDELGET, PDGEMR2D, PDLASCL, PDLASET, PDORMTR,
00276      $                   PDSYTRD, PXERBLA
00277 *     ..
00278 *     .. Intrinsic Functions ..
00279       INTRINSIC          ABS, DBLE, ICHAR, MAX, MIN, MOD, SQRT, INT
00280 *     ..
00281 *     .. Executable Statements ..
00282 *       This is just to keep ftnchek and toolpack/1 happy
00283       IF( BLOCK_CYCLIC_2D*CSRC_*CTXT_*DLEN_*DTYPE_*LLD_*MB_*M_*NB_*N_*
00284      $    RSRC_.LT.0 )RETURN
00285 *
00286 *     Quick return
00287 *
00288       IF( N.EQ.0 ) RETURN
00289 *
00290 *     Test the input arguments.
00291 *
00292       CALL BLACS_GRIDINFO( DESCA( CTXT_ ), NPROW, NPCOL, MYROW, MYCOL )
00293       INFO = 0
00294 *
00295       WANTZ = LSAME( JOBZ, 'V' )
00296       IF( NPROW.EQ.-1 ) THEN
00297          INFO = -( 700+CTXT_ )
00298       ELSE IF( WANTZ ) THEN
00299          IF( DESCA( CTXT_ ).NE.DESCZ( CTXT_ ) ) THEN
00300             INFO = -( 1200+CTXT_ )
00301          END IF
00302       END IF
00303       IF( INFO .EQ. 0 ) THEN
00304          CALL CHK1MAT( N, 3, N, 3, IA, JA, DESCA, 7, INFO )
00305          IF( WANTZ )
00306      $      CALL CHK1MAT( N, 3, N, 3, IZ, JZ, DESCZ, 12, INFO )
00307 *
00308          IF( INFO.EQ.0 ) THEN
00309 *
00310 *           Get machine constants.
00311 *
00312             SAFMIN = PDLAMCH( DESCA( CTXT_ ), 'Safe minimum' )
00313             EPS = PDLAMCH( DESCA( CTXT_ ), 'Precision' )
00314             SMLNUM = SAFMIN / EPS
00315             BIGNUM = ONE / SMLNUM
00316             RMIN = SQRT( SMLNUM )
00317             RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
00318 *
00319             NPROCS = NPROW*NPCOL
00320             NB_A = DESCA( NB_ )
00321             MB_A = DESCA( MB_ )
00322             NB = NB_A
00323             LOWER = LSAME( UPLO, 'L' )
00324 *
00325             RSRC_A = DESCA( RSRC_ )
00326             CSRC_A = DESCA( CSRC_ )
00327             IROFFA = MOD( IA-1, MB_A )
00328             ICOFFA = MOD( JA-1, NB_A )
00329             IAROW = INDXG2P( 1, NB_A, MYROW, RSRC_A, NPROW )
00330             IACOL = INDXG2P( 1, MB_A, MYCOL, CSRC_A, NPCOL )
00331             NP = NUMROC( N+IROFFA, NB, MYROW, IAROW, NPROW )
00332             NQ = NUMROC( N+ICOFFA, NB, MYCOL, IACOL, NPCOL )
00333 
00334             IF( WANTZ ) THEN
00335                NB_Z = DESCZ( NB_ )
00336                MB_Z = DESCZ( MB_ )
00337                RSRC_Z = DESCZ( RSRC_ )
00338                IROFFZ = MOD( IZ-1, MB_A )
00339                IZROW = INDXG2P( 1, NB_A, MYROW, RSRC_Z, NPROW )
00340                SIZEMQRLEFT = MAX( ( NB_A*( NB_A-1 ) ) / 2, ( NP+NQ )*
00341      $                            NB_A ) + NB_A*NB_A
00342             ELSE
00343                SIZEMQRLEFT = 0
00344                IROFFZ = 0
00345                IZROW = 0
00346             END IF
00347             SIZESYTRD = MAX( NB * ( NP +1 ), 3 * NB )
00348 *
00349 *           Initialize the context of the single column distributed
00350 *           matrix required by DSTEQR2. This specific distribution
00351 *           allows each process to do 1/pth of the work updating matrix
00352 *           Q during DSTEQR2 and achieve some parallelization to an
00353 *           otherwise serial subroutine.
00354 *
00355             LDC = 0
00356             IF( WANTZ ) THEN
00357                CONTEXTC = SL_GRIDRESHAPE( DESCA( CTXT_ ), 0, 1, 1,
00358      $                                    NPROCS, 1 )
00359                CALL BLACS_GRIDINFO( CONTEXTC, NPROWC, NPCOLC, MYPROWC,
00360      $                              MYPCOLC )
00361                NRC = NUMROC( N, NB_A, MYPROWC, 0, NPROCS)
00362                LDC = MAX( 1, NRC )
00363                CALL DESCINIT( DESCQR, N, N, NB, NB, 0, 0, CONTEXTC,
00364      $                        LDC, INFO )
00365             END IF
00366 
00367 *
00368 *           Set up pointers into the WORK array
00369 *
00370             INDTAU = 1
00371             INDE = INDTAU + N
00372             INDD = INDE + N
00373             INDD2 = INDD + N
00374             INDE2 = INDD2 + N
00375             INDWORK = INDE2 + N
00376             INDWORK2 = INDWORK + N*LDC
00377             LLWORK = LWORK - INDWORK + 1
00378 *
00379 *           Compute the total amount of space needed
00380 *
00381             QRMEM = 2*N-2
00382             IF( WANTZ ) THEN
00383                LWMIN = 5*N + N*LDC + MAX( SIZEMQRLEFT, QRMEM ) + 1
00384             ELSE
00385                LWMIN = 5*N + SIZESYTRD + 1
00386             END IF
00387 *
00388          END IF
00389          IF( INFO.EQ.0 ) THEN
00390             IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
00391                INFO = -1
00392             ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
00393                INFO = -2
00394             ELSE IF( LWORK.LT.LWMIN .AND. LWORK.NE.-1 ) THEN
00395                INFO = -14
00396             ELSE IF( IROFFA.NE.0 ) THEN
00397                INFO = -5
00398             ELSE IF( DESCA( MB_ ).NE.DESCA( NB_ ) ) THEN
00399                INFO = -( 700+NB_ )
00400             END IF
00401             IF( WANTZ ) THEN
00402                IF( IROFFA.NE.IROFFZ ) THEN
00403                   INFO = -10
00404                ELSE IF( IAROW.NE.IZROW ) THEN
00405                   INFO = -10
00406                ELSE IF( DESCA( M_ ).NE.DESCZ( M_ ) ) THEN
00407                   INFO = -( 1200+M_ )
00408                ELSE IF( DESCA( N_ ).NE.DESCZ( N_ ) ) THEN
00409                   INFO = -( 1200+N_ )
00410                ELSE IF( DESCA( MB_ ).NE.DESCZ( MB_ ) ) THEN
00411                   INFO = -( 1200+MB_ )
00412                ELSE IF( DESCA( NB_ ).NE.DESCZ( NB_ ) ) THEN
00413                   INFO = -( 1200+NB_ )
00414                ELSE IF( DESCA( RSRC_ ).NE.DESCZ( RSRC_ ) ) THEN
00415                   INFO = -( 1200+RSRC_ )
00416                ELSE IF( DESCA( CTXT_ ).NE.DESCZ( CTXT_ ) ) THEN
00417                   INFO = -( 1200+CTXT_ )
00418                ENDIF
00419             END IF
00420          END IF
00421          IF( WANTZ ) THEN
00422             IDUM1( 1 ) = ICHAR( 'V' )
00423          ELSE
00424             IDUM1( 1 ) = ICHAR( 'N' )
00425          END IF
00426          IDUM2( 1 ) = 1
00427          IF( LOWER ) THEN
00428             IDUM1( 2 ) = ICHAR( 'L' )
00429          ELSE
00430             IDUM1( 2 ) = ICHAR( 'U' )
00431          END IF
00432          IDUM2( 2 ) = 2
00433          IF( LWORK.EQ.-1 ) THEN
00434             IDUM1( 3 ) = -1
00435          ELSE
00436             IDUM1( 3 ) = 1
00437          END IF
00438          IDUM2( 3 ) = 3
00439          IF( LSAME( JOBZ, 'V' ) ) THEN
00440             CALL PCHK2MAT( N, 3, N, 3, IA, JA, DESCA, 7, N, 3, N, 3,
00441      $                     IZ, JZ, DESCZ, 12, 3, IDUM1, IDUM2, INFO )
00442          ELSE
00443             CALL PCHK1MAT( N, 3, N, 3, IA, JA, DESCA, 7, 3, IDUM1,
00444      $                     IDUM2, INFO )
00445          END IF
00446 *
00447 *     Write the required workspace for lwork queries.
00448 *
00449          WORK( 1 ) = DBLE( LWMIN )
00450       END IF
00451 *
00452       IF( INFO.NE.0 ) THEN
00453          CALL PXERBLA( DESCA( CTXT_ ), 'PDSYEV', -INFO )
00454          IF( WANTZ ) CALL BLACS_GRIDEXIT( CONTEXTC )
00455          RETURN
00456       ELSE IF( LWORK .EQ. -1 ) THEN
00457          IF( WANTZ ) CALL BLACS_GRIDEXIT( CONTEXTC )
00458          RETURN
00459       END IF
00460 *
00461 *     Scale matrix to allowable range, if necessary.
00462 *
00463       ISCALE = 0
00464 *
00465       ANRM = PDLANSY( 'M', UPLO, N, A, IA, JA, DESCA, WORK( INDWORK ) )
00466 *
00467 
00468       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
00469          ISCALE = 1
00470          SIGMA = RMIN / ANRM
00471       ELSE IF( ANRM.GT.RMAX ) THEN
00472          ISCALE = 1
00473          SIGMA = RMAX / ANRM
00474       END IF
00475 *
00476       IF( ISCALE.EQ.1 ) THEN
00477          CALL PDLASCL( UPLO, ONE, SIGMA, N, N, A, IA, JA, DESCA, IINFO )
00478       END IF
00479 *
00480 *     Reduce symmetric matrix to tridiagonal form.
00481 *
00482       CALL PDSYTRD( UPLO, N, A, IA, JA, DESCA, WORK( INDD ),
00483      $              WORK( INDE ), WORK( INDTAU ), WORK( INDWORK ),
00484      $              LLWORK, IINFO )
00485 *
00486 *     Copy the values of D, E to all processes.
00487 *
00488       DO 10 I=1,N
00489          CALL PDELGET( 'A', ' ', WORK(INDD2+I-1), A,
00490      $                 I+IA-1, I+JA-1, DESCA )
00491  10   CONTINUE
00492       IF( LSAME( UPLO, 'U') ) THEN
00493           DO 20 I=1,N-1
00494              CALL PDELGET( 'A', ' ', WORK(INDE2+I-1), A, 
00495      $                     I+IA-1, I+JA, DESCA )
00496  20       CONTINUE
00497       ELSE
00498           DO 30 I=1,N-1
00499              CALL PDELGET( 'A', ' ', WORK(INDE2+I-1), A,
00500      $                     I+IA, I+JA-1, DESCA )
00501  30       CONTINUE
00502       ENDIF
00503 *
00504       IF( WANTZ ) THEN
00505 *
00506          CALL PDLASET( 'Full', N, N, ZERO, ONE, WORK( INDWORK ), 1, 1,
00507      $                 DESCQR )
00508 *
00509 *        DSTEQR2 is a modified version of LAPACK's DSTEQR.  The
00510 *        modifications allow each process to perform partial updates
00511 *        to matrix Q.
00512 *
00513          CALL DSTEQR2( 'I', N, WORK( INDD2 ), WORK( INDE2 ),
00514      $                 WORK( INDWORK ), LDC, NRC, WORK( INDWORK2 ), 
00515      $                 INFO )
00516 *
00517          CALL PDGEMR2D( N, N, WORK( INDWORK ), 1, 1, DESCQR, Z, IA, JA,
00518      $                  DESCZ, CONTEXTC )
00519 *
00520          CALL PDORMTR( 'L', UPLO, 'N', N, N, A, IA, JA, DESCA,
00521      $                 WORK( INDTAU ), Z, IZ, JZ, DESCZ,
00522      $                 WORK( INDWORK ), LLWORK, IINFO )
00523 *
00524       ELSE
00525 *
00526          CALL DSTEQR2( 'N', N, WORK( INDD2 ), WORK( INDE2 ),
00527      $                 WORK( INDWORK ), 1, 1, WORK( INDWORK2 ),
00528      $                 INFO )
00529       ENDIF
00530 *
00531 *     Copy eigenvalues from workspace to output array
00532 *
00533       CALL DCOPY( N, WORK( INDD2 ), 1, W, 1 )
00534 *
00535 *     If matrix was scaled, then rescale eigenvalues appropriately.
00536 *
00537       IF( ISCALE .EQ. 1 ) THEN
00538          CALL DSCAL( N, ONE / SIGMA, W, 1 )
00539       END IF
00540 *
00541 *     Free up resources
00542 *
00543       IF( WANTZ ) THEN
00544            CALL BLACS_GRIDEXIT( CONTEXTC )
00545       END IF
00546 *
00547 *     Compare every ith eigenvalue, or all if there are only a few,
00548 *     across the process grid to check for heterogeneity.
00549 *
00550       IF( N.LE.ITHVAL ) THEN
00551          J = N
00552          K = 1
00553       ELSE
00554          J = N/ITHVAL
00555          K = ITHVAL
00556       END IF
00557 *
00558       DO 40 I = 1, J
00559          WORK( I+INDTAU ) = W( (I-1)*K+1 )
00560          WORK( I+INDE ) = W( (I-1)*K+1 )
00561  40   CONTINUE
00562 *
00563       CALL DGAMN2D( DESCA( CTXT_ ), 'a', ' ', J, 1, WORK( 1+INDTAU ),
00564      $              J, 1, 1, -1, -1, 0 )
00565       CALL DGAMX2D( DESCA( CTXT_ ), 'a', ' ', J, 1, WORK( 1+INDE ),
00566      $              J, 1, 1, -1, -1, 0 )
00567 *
00568       DO 50 I = 1, J
00569          IF( INFO.EQ.0 .AND. ( WORK( I+INDTAU )-WORK( I+INDE )
00570      $        .NE. ZERO ) )THEN 
00571             INFO = N+1
00572          END IF
00573  50   CONTINUE
00574 *
00575       RETURN
00576 *
00577 *     End of PDSYEV
00578 *
00579       END