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