001:       SUBROUTINE ZHEEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
002:      $                   LRWORK, IWORK, LIWORK, INFO )
003: *
004: *  -- LAPACK driver routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       CHARACTER          JOBZ, UPLO
010:       INTEGER            INFO, LDA, LIWORK, LRWORK, LWORK, N
011: *     ..
012: *     .. Array Arguments ..
013:       INTEGER            IWORK( * )
014:       DOUBLE PRECISION   RWORK( * ), W( * )
015:       COMPLEX*16         A( LDA, * ), WORK( * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  ZHEEVD computes all eigenvalues and, optionally, eigenvectors of a
022: *  complex Hermitian 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: *  Arguments
033: *  =========
034: *
035: *  JOBZ    (input) CHARACTER*1
036: *          = 'N':  Compute eigenvalues only;
037: *          = 'V':  Compute eigenvalues and eigenvectors.
038: *
039: *  UPLO    (input) CHARACTER*1
040: *          = 'U':  Upper triangle of A is stored;
041: *          = 'L':  Lower triangle of A is stored.
042: *
043: *  N       (input) INTEGER
044: *          The order of the matrix A.  N >= 0.
045: *
046: *  A       (input/output) COMPLEX*16 array, dimension (LDA, N)
047: *          On entry, the Hermitian matrix A.  If UPLO = 'U', the
048: *          leading N-by-N upper triangular part of A contains the
049: *          upper triangular part of the matrix A.  If UPLO = 'L',
050: *          the leading N-by-N lower triangular part of A contains
051: *          the lower triangular part of the matrix A.
052: *          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
053: *          orthonormal eigenvectors of the matrix A.
054: *          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
055: *          or the upper triangle (if UPLO='U') of A, including the
056: *          diagonal, is destroyed.
057: *
058: *  LDA     (input) INTEGER
059: *          The leading dimension of the array A.  LDA >= max(1,N).
060: *
061: *  W       (output) DOUBLE PRECISION array, dimension (N)
062: *          If INFO = 0, the eigenvalues in ascending order.
063: *
064: *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
065: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
066: *
067: *  LWORK   (input) INTEGER
068: *          The length of the array WORK.
069: *          If N <= 1,                LWORK must be at least 1.
070: *          If JOBZ  = 'N' and N > 1, LWORK must be at least N + 1.
071: *          If JOBZ  = 'V' and N > 1, LWORK must be at least 2*N + N**2.
072: *
073: *          If LWORK = -1, then a workspace query is assumed; the routine
074: *          only calculates the optimal sizes of the WORK, RWORK and
075: *          IWORK arrays, returns these values as the first entries of
076: *          the WORK, RWORK and IWORK arrays, and no error message
077: *          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
078: *
079: *  RWORK   (workspace/output) DOUBLE PRECISION array,
080: *                                         dimension (LRWORK)
081: *          On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
082: *
083: *  LRWORK  (input) INTEGER
084: *          The dimension of the array RWORK.
085: *          If N <= 1,                LRWORK must be at least 1.
086: *          If JOBZ  = 'N' and N > 1, LRWORK must be at least N.
087: *          If JOBZ  = 'V' and N > 1, LRWORK must be at least
088: *                         1 + 5*N + 2*N**2.
089: *
090: *          If LRWORK = -1, then a workspace query is assumed; the
091: *          routine only calculates the optimal sizes of the WORK, RWORK
092: *          and IWORK arrays, returns these values as the first entries
093: *          of the WORK, RWORK and IWORK arrays, and no error message
094: *          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
095: *
096: *  IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
097: *          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
098: *
099: *  LIWORK  (input) INTEGER
100: *          The dimension of the array IWORK.
101: *          If N <= 1,                LIWORK must be at least 1.
102: *          If JOBZ  = 'N' and N > 1, LIWORK must be at least 1.
103: *          If JOBZ  = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
104: *
105: *          If LIWORK = -1, then a workspace query is assumed; the
106: *          routine only calculates the optimal sizes of the WORK, RWORK
107: *          and IWORK arrays, returns these values as the first entries
108: *          of the WORK, RWORK and IWORK arrays, and no error message
109: *          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
110: *
111: *  INFO    (output) INTEGER
112: *          = 0:  successful exit
113: *          < 0:  if INFO = -i, the i-th argument had an illegal value
114: *          > 0:  if INFO = i and JOBZ = 'N', then the algorithm failed
115: *                to converge; i off-diagonal elements of an intermediate
116: *                tridiagonal form did not converge to zero;
117: *                if INFO = i and JOBZ = 'V', then the algorithm failed
118: *                to compute an eigenvalue while working on the submatrix
119: *                lying in rows and columns INFO/(N+1) through
120: *                mod(INFO,N+1).
121: *
122: *  Further Details
123: *  ===============
124: *
125: *  Based on contributions by
126: *     Jeff Rutter, Computer Science Division, University of California
127: *     at Berkeley, USA
128: *
129: *  Modified description of INFO. Sven, 16 Feb 05.
130: *  =====================================================================
131: *
132: *     .. Parameters ..
133:       DOUBLE PRECISION   ZERO, ONE
134:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
135:       COMPLEX*16         CONE
136:       PARAMETER          ( CONE = ( 1.0D0, 0.0D0 ) )
137: *     ..
138: *     .. Local Scalars ..
139:       LOGICAL            LOWER, LQUERY, WANTZ
140:       INTEGER            IINFO, IMAX, INDE, INDRWK, INDTAU, INDWK2,
141:      $                   INDWRK, ISCALE, LIOPT, LIWMIN, LLRWK, LLWORK,
142:      $                   LLWRK2, LOPT, LROPT, LRWMIN, LWMIN
143:       DOUBLE PRECISION   ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
144:      $                   SMLNUM
145: *     ..
146: *     .. External Functions ..
147:       LOGICAL            LSAME
148:       INTEGER            ILAENV
149:       DOUBLE PRECISION   DLAMCH, ZLANHE
150:       EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANHE
151: *     ..
152: *     .. External Subroutines ..
153:       EXTERNAL           DSCAL, DSTERF, XERBLA, ZHETRD, ZLACPY, ZLASCL,
154:      $                   ZSTEDC, ZUNMTR
155: *     ..
156: *     .. Intrinsic Functions ..
157:       INTRINSIC          MAX, SQRT
158: *     ..
159: *     .. Executable Statements ..
160: *
161: *     Test the input parameters.
162: *
163:       WANTZ = LSAME( JOBZ, 'V' )
164:       LOWER = LSAME( UPLO, 'L' )
165:       LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
166: *
167:       INFO = 0
168:       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
169:          INFO = -1
170:       ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
171:          INFO = -2
172:       ELSE IF( N.LT.0 ) THEN
173:          INFO = -3
174:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
175:          INFO = -5
176:       END IF
177: *
178:       IF( INFO.EQ.0 ) THEN
179:          IF( N.LE.1 ) THEN
180:             LWMIN = 1
181:             LRWMIN = 1
182:             LIWMIN = 1
183:             LOPT = LWMIN
184:             LROPT = LRWMIN
185:             LIOPT = LIWMIN
186:          ELSE
187:             IF( WANTZ ) THEN
188:                LWMIN = 2*N + N*N
189:                LRWMIN = 1 + 5*N + 2*N**2
190:                LIWMIN = 3 + 5*N
191:             ELSE
192:                LWMIN = N + 1
193:                LRWMIN = N
194:                LIWMIN = 1
195:             END IF
196:             LOPT = MAX( LWMIN, N +
197:      $                  ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 ) )
198:             LROPT = LRWMIN
199:             LIOPT = LIWMIN
200:          END IF
201:          WORK( 1 ) = LOPT
202:          RWORK( 1 ) = LROPT
203:          IWORK( 1 ) = LIOPT
204: *
205:          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
206:             INFO = -8
207:          ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
208:             INFO = -10
209:          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
210:             INFO = -12
211:          END IF
212:       END IF
213: *
214:       IF( INFO.NE.0 ) THEN
215:          CALL XERBLA( 'ZHEEVD', -INFO )
216:          RETURN
217:       ELSE IF( LQUERY ) THEN
218:          RETURN
219:       END IF
220: *
221: *     Quick return if possible
222: *
223:       IF( N.EQ.0 )
224:      $   RETURN
225: *
226:       IF( N.EQ.1 ) THEN
227:          W( 1 ) = A( 1, 1 )
228:          IF( WANTZ )
229:      $      A( 1, 1 ) = CONE
230:          RETURN
231:       END IF
232: *
233: *     Get machine constants.
234: *
235:       SAFMIN = DLAMCH( 'Safe minimum' )
236:       EPS = DLAMCH( 'Precision' )
237:       SMLNUM = SAFMIN / EPS
238:       BIGNUM = ONE / SMLNUM
239:       RMIN = SQRT( SMLNUM )
240:       RMAX = SQRT( BIGNUM )
241: *
242: *     Scale matrix to allowable range, if necessary.
243: *
244:       ANRM = ZLANHE( 'M', UPLO, N, A, LDA, RWORK )
245:       ISCALE = 0
246:       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
247:          ISCALE = 1
248:          SIGMA = RMIN / ANRM
249:       ELSE IF( ANRM.GT.RMAX ) THEN
250:          ISCALE = 1
251:          SIGMA = RMAX / ANRM
252:       END IF
253:       IF( ISCALE.EQ.1 )
254:      $   CALL ZLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO )
255: *
256: *     Call ZHETRD to reduce Hermitian matrix to tridiagonal form.
257: *
258:       INDE = 1
259:       INDTAU = 1
260:       INDWRK = INDTAU + N
261:       INDRWK = INDE + N
262:       INDWK2 = INDWRK + N*N
263:       LLWORK = LWORK - INDWRK + 1
264:       LLWRK2 = LWORK - INDWK2 + 1
265:       LLRWK = LRWORK - INDRWK + 1
266:       CALL ZHETRD( UPLO, N, A, LDA, W, RWORK( INDE ), WORK( INDTAU ),
267:      $             WORK( INDWRK ), LLWORK, IINFO )
268: *
269: *     For eigenvalues only, call DSTERF.  For eigenvectors, first call
270: *     ZSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
271: *     tridiagonal matrix, then call ZUNMTR to multiply it to the
272: *     Householder transformations represented as Householder vectors in
273: *     A.
274: *
275:       IF( .NOT.WANTZ ) THEN
276:          CALL DSTERF( N, W, RWORK( INDE ), INFO )
277:       ELSE
278:          CALL ZSTEDC( 'I', N, W, RWORK( INDE ), WORK( INDWRK ), N,
279:      $                WORK( INDWK2 ), LLWRK2, RWORK( INDRWK ), LLRWK,
280:      $                IWORK, LIWORK, INFO )
281:          CALL ZUNMTR( 'L', UPLO, 'N', N, N, A, LDA, WORK( INDTAU ),
282:      $                WORK( INDWRK ), N, WORK( INDWK2 ), LLWRK2, IINFO )
283:          CALL ZLACPY( 'A', N, N, WORK( INDWRK ), N, A, LDA )
284:       END IF
285: *
286: *     If matrix was scaled, then rescale eigenvalues appropriately.
287: *
288:       IF( ISCALE.EQ.1 ) THEN
289:          IF( INFO.EQ.0 ) THEN
290:             IMAX = N
291:          ELSE
292:             IMAX = INFO - 1
293:          END IF
294:          CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
295:       END IF
296: *
297:       WORK( 1 ) = LOPT
298:       RWORK( 1 ) = LROPT
299:       IWORK( 1 ) = LIOPT
300: *
301:       RETURN
302: *
303: *     End of ZHEEVD
304: *
305:       END
306: