001:       SUBROUTINE DSYEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, IWORK,
002:      $                   LIWORK, INFO )
003: *
004: *  -- LAPACK driver routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       CHARACTER          JOBZ, UPLO
011:       INTEGER            INFO, LDA, LIWORK, LWORK, N
012: *     ..
013: *     .. Array Arguments ..
014:       INTEGER            IWORK( * )
015:       DOUBLE PRECISION   A( LDA, * ), W( * ), WORK( * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  DSYEVD computes all eigenvalues and, optionally, eigenvectors of a
022: *  real symmetric matrix A. If eigenvectors are desired, it uses a
023: *  divide and conquer algorithm.
024: *
025: *  The divide and conquer algorithm makes very mild assumptions about
026: *  floating point arithmetic. It will work on machines with a guard
027: *  digit in add/subtract, or on those binary machines without guard
028: *  digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
029: *  Cray-2. It could conceivably fail on hexadecimal or decimal machines
030: *  without guard digits, but we know of none.
031: *
032: *  Because of large use of BLAS of level 3, DSYEVD needs N**2 more
033: *  workspace than DSYEVX.
034: *
035: *  Arguments
036: *  =========
037: *
038: *  JOBZ    (input) CHARACTER*1
039: *          = 'N':  Compute eigenvalues only;
040: *          = 'V':  Compute eigenvalues and eigenvectors.
041: *
042: *  UPLO    (input) CHARACTER*1
043: *          = 'U':  Upper triangle of A is stored;
044: *          = 'L':  Lower triangle of A is stored.
045: *
046: *  N       (input) INTEGER
047: *          The order of the matrix A.  N >= 0.
048: *
049: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA, N)
050: *          On entry, the symmetric matrix A.  If UPLO = 'U', the
051: *          leading N-by-N upper triangular part of A contains the
052: *          upper triangular part of the matrix A.  If UPLO = 'L',
053: *          the leading N-by-N lower triangular part of A contains
054: *          the lower triangular part of the matrix A.
055: *          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
056: *          orthonormal eigenvectors of the matrix A.
057: *          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
058: *          or the upper triangle (if UPLO='U') of A, including the
059: *          diagonal, is destroyed.
060: *
061: *  LDA     (input) INTEGER
062: *          The leading dimension of the array A.  LDA >= max(1,N).
063: *
064: *  W       (output) DOUBLE PRECISION array, dimension (N)
065: *          If INFO = 0, the eigenvalues in ascending order.
066: *
067: *  WORK    (workspace/output) DOUBLE PRECISION array,
068: *                                         dimension (LWORK)
069: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
070: *
071: *  LWORK   (input) INTEGER
072: *          The dimension of the array WORK.
073: *          If N <= 1,               LWORK must be at least 1.
074: *          If JOBZ = 'N' and N > 1, LWORK must be at least 2*N+1.
075: *          If JOBZ = 'V' and N > 1, LWORK must be at least
076: *                                                1 + 6*N + 2*N**2.
077: *
078: *          If LWORK = -1, then a workspace query is assumed; the routine
079: *          only calculates the optimal sizes of the WORK and IWORK
080: *          arrays, returns these values as the first entries of the WORK
081: *          and IWORK arrays, and no error message related to LWORK or
082: *          LIWORK is issued by XERBLA.
083: *
084: *  IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
085: *          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
086: *
087: *  LIWORK  (input) INTEGER
088: *          The dimension of the array IWORK.
089: *          If N <= 1,                LIWORK must be at least 1.
090: *          If JOBZ  = 'N' and N > 1, LIWORK must be at least 1.
091: *          If JOBZ  = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
092: *
093: *          If LIWORK = -1, then a workspace query is assumed; the
094: *          routine only calculates the optimal sizes of the WORK and
095: *          IWORK arrays, returns these values as the first entries of
096: *          the WORK and IWORK arrays, and no error message related to
097: *          LWORK or LIWORK is issued by XERBLA.
098: *
099: *  INFO    (output) INTEGER
100: *          = 0:  successful exit
101: *          < 0:  if INFO = -i, the i-th argument had an illegal value
102: *          > 0:  if INFO = i and JOBZ = 'N', then the algorithm failed
103: *                to converge; i off-diagonal elements of an intermediate
104: *                tridiagonal form did not converge to zero;
105: *                if INFO = i and JOBZ = 'V', then the algorithm failed
106: *                to compute an eigenvalue while working on the submatrix
107: *                lying in rows and columns INFO/(N+1) through
108: *                mod(INFO,N+1).
109: *
110: *  Further Details
111: *  ===============
112: *
113: *  Based on contributions by
114: *     Jeff Rutter, Computer Science Division, University of California
115: *     at Berkeley, USA
116: *  Modified by Francoise Tisseur, University of Tennessee.
117: *
118: *  Modified description of INFO. Sven, 16 Feb 05.
119: *  =====================================================================
120: *
121: *     .. Parameters ..
122:       DOUBLE PRECISION   ZERO, ONE
123:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
124: *     ..
125: *     .. Local Scalars ..
126: *
127:       LOGICAL            LOWER, LQUERY, WANTZ
128:       INTEGER            IINFO, INDE, INDTAU, INDWK2, INDWRK, ISCALE,
129:      $                   LIOPT, LIWMIN, LLWORK, LLWRK2, LOPT, LWMIN
130:       DOUBLE PRECISION   ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
131:      $                   SMLNUM
132: *     ..
133: *     .. External Functions ..
134:       LOGICAL            LSAME
135:       INTEGER            ILAENV
136:       DOUBLE PRECISION   DLAMCH, DLANSY
137:       EXTERNAL           LSAME, DLAMCH, DLANSY, ILAENV
138: *     ..
139: *     .. External Subroutines ..
140:       EXTERNAL           DLACPY, DLASCL, DORMTR, DSCAL, DSTEDC, DSTERF,
141:      $                   DSYTRD, XERBLA
142: *     ..
143: *     .. Intrinsic Functions ..
144:       INTRINSIC          MAX, SQRT
145: *     ..
146: *     .. Executable Statements ..
147: *
148: *     Test the input parameters.
149: *
150:       WANTZ = LSAME( JOBZ, 'V' )
151:       LOWER = LSAME( UPLO, 'L' )
152:       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
153: *
154:       INFO = 0
155:       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
156:          INFO = -1
157:       ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
158:          INFO = -2
159:       ELSE IF( N.LT.0 ) THEN
160:          INFO = -3
161:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
162:          INFO = -5
163:       END IF
164: *
165:       IF( INFO.EQ.0 ) THEN
166:          IF( N.LE.1 ) THEN
167:             LIWMIN = 1
168:             LWMIN = 1
169:             LOPT = LWMIN
170:             LIOPT = LIWMIN
171:          ELSE
172:             IF( WANTZ ) THEN
173:                LIWMIN = 3 + 5*N
174:                LWMIN = 1 + 6*N + 2*N**2
175:             ELSE
176:                LIWMIN = 1
177:                LWMIN = 2*N + 1
178:             END IF
179:             LOPT = MAX( LWMIN, 2*N +
180:      $                  ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 ) )
181:             LIOPT = LIWMIN
182:          END IF
183:          WORK( 1 ) = LOPT
184:          IWORK( 1 ) = LIOPT
185: *
186:          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
187:             INFO = -8
188:          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
189:             INFO = -10
190:          END IF
191:       END IF
192: *
193:       IF( INFO.NE.0 ) THEN
194:          CALL XERBLA( 'DSYEVD', -INFO )
195:          RETURN
196:       ELSE IF( LQUERY ) THEN
197:          RETURN
198:       END IF
199: *
200: *     Quick return if possible
201: *
202:       IF( N.EQ.0 )
203:      $   RETURN
204: *
205:       IF( N.EQ.1 ) THEN
206:          W( 1 ) = A( 1, 1 )
207:          IF( WANTZ )
208:      $      A( 1, 1 ) = ONE
209:          RETURN
210:       END IF
211: *
212: *     Get machine constants.
213: *
214:       SAFMIN = DLAMCH( 'Safe minimum' )
215:       EPS = DLAMCH( 'Precision' )
216:       SMLNUM = SAFMIN / EPS
217:       BIGNUM = ONE / SMLNUM
218:       RMIN = SQRT( SMLNUM )
219:       RMAX = SQRT( BIGNUM )
220: *
221: *     Scale matrix to allowable range, if necessary.
222: *
223:       ANRM = DLANSY( 'M', UPLO, N, A, LDA, WORK )
224:       ISCALE = 0
225:       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
226:          ISCALE = 1
227:          SIGMA = RMIN / ANRM
228:       ELSE IF( ANRM.GT.RMAX ) THEN
229:          ISCALE = 1
230:          SIGMA = RMAX / ANRM
231:       END IF
232:       IF( ISCALE.EQ.1 )
233:      $   CALL DLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO )
234: *
235: *     Call DSYTRD to reduce symmetric matrix to tridiagonal form.
236: *
237:       INDE = 1
238:       INDTAU = INDE + N
239:       INDWRK = INDTAU + N
240:       LLWORK = LWORK - INDWRK + 1
241:       INDWK2 = INDWRK + N*N
242:       LLWRK2 = LWORK - INDWK2 + 1
243: *
244:       CALL DSYTRD( UPLO, N, A, LDA, W, WORK( INDE ), WORK( INDTAU ),
245:      $             WORK( INDWRK ), LLWORK, IINFO )
246:       LOPT = 2*N + WORK( INDWRK )
247: *
248: *     For eigenvalues only, call DSTERF.  For eigenvectors, first call
249: *     DSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
250: *     tridiagonal matrix, then call DORMTR to multiply it by the
251: *     Householder transformations stored in A.
252: *
253:       IF( .NOT.WANTZ ) THEN
254:          CALL DSTERF( N, W, WORK( INDE ), INFO )
255:       ELSE
256:          CALL DSTEDC( 'I', N, W, WORK( INDE ), WORK( INDWRK ), N,
257:      $                WORK( INDWK2 ), LLWRK2, IWORK, LIWORK, INFO )
258:          CALL DORMTR( 'L', UPLO, 'N', N, N, A, LDA, WORK( INDTAU ),
259:      $                WORK( INDWRK ), N, WORK( INDWK2 ), LLWRK2, IINFO )
260:          CALL DLACPY( 'A', N, N, WORK( INDWRK ), N, A, LDA )
261:          LOPT = MAX( LOPT, 1+6*N+2*N**2 )
262:       END IF
263: *
264: *     If matrix was scaled, then rescale eigenvalues appropriately.
265: *
266:       IF( ISCALE.EQ.1 )
267:      $   CALL DSCAL( N, ONE / SIGMA, W, 1 )
268: *
269:       WORK( 1 ) = LOPT
270:       IWORK( 1 ) = LIOPT
271: *
272:       RETURN
273: *
274: *     End of DSYEVD
275: *
276:       END
277: