001:       SUBROUTINE CHEEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU,
002:      $                   ABSTOL, M, W, Z, LDZ, WORK, LWORK, RWORK,
003:      $                   IWORK, IFAIL, INFO )
004: *
005: *  -- LAPACK driver routine (version 3.2) --
006: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       CHARACTER          JOBZ, RANGE, UPLO
011:       INTEGER            IL, INFO, IU, LDA, LDZ, LWORK, M, N
012:       REAL               ABSTOL, VL, VU
013: *     ..
014: *     .. Array Arguments ..
015:       INTEGER            IFAIL( * ), IWORK( * )
016:       REAL               RWORK( * ), W( * )
017:       COMPLEX            A( LDA, * ), WORK( * ), Z( LDZ, * )
018: *     ..
019: *
020: *  Purpose
021: *  =======
022: *
023: *  CHEEVX computes selected eigenvalues and, optionally, eigenvectors
024: *  of a complex Hermitian matrix A.  Eigenvalues and eigenvectors can
025: *  be selected by specifying either a range of values or a range of
026: *  indices for the desired eigenvalues.
027: *
028: *  Arguments
029: *  =========
030: *
031: *  JOBZ    (input) CHARACTER*1
032: *          = 'N':  Compute eigenvalues only;
033: *          = 'V':  Compute eigenvalues and eigenvectors.
034: *
035: *  RANGE   (input) CHARACTER*1
036: *          = 'A': all eigenvalues will be found.
037: *          = 'V': all eigenvalues in the half-open interval (VL,VU]
038: *                 will be found.
039: *          = 'I': the IL-th through IU-th eigenvalues will be found.
040: *
041: *  UPLO    (input) CHARACTER*1
042: *          = 'U':  Upper triangle of A is stored;
043: *          = 'L':  Lower triangle of A is stored.
044: *
045: *  N       (input) INTEGER
046: *          The order of the matrix A.  N >= 0.
047: *
048: *  A       (input/output) COMPLEX array, dimension (LDA, N)
049: *          On entry, the Hermitian matrix A.  If UPLO = 'U', the
050: *          leading N-by-N upper triangular part of A contains the
051: *          upper triangular part of the matrix A.  If UPLO = 'L',
052: *          the leading N-by-N lower triangular part of A contains
053: *          the lower triangular part of the matrix A.
054: *          On exit, the lower triangle (if UPLO='L') or the upper
055: *          triangle (if UPLO='U') of A, including the diagonal, is
056: *          destroyed.
057: *
058: *  LDA     (input) INTEGER
059: *          The leading dimension of the array A.  LDA >= max(1,N).
060: *
061: *  VL      (input) REAL
062: *  VU      (input) REAL
063: *          If RANGE='V', the lower and upper bounds of the interval to
064: *          be searched for eigenvalues. VL < VU.
065: *          Not referenced if RANGE = 'A' or 'I'.
066: *
067: *  IL      (input) INTEGER
068: *  IU      (input) INTEGER
069: *          If RANGE='I', the indices (in ascending order) of the
070: *          smallest and largest eigenvalues to be returned.
071: *          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
072: *          Not referenced if RANGE = 'A' or 'V'.
073: *
074: *  ABSTOL  (input) REAL
075: *          The absolute error tolerance for the eigenvalues.
076: *          An approximate eigenvalue is accepted as converged
077: *          when it is determined to lie in an interval [a,b]
078: *          of width less than or equal to
079: *
080: *                  ABSTOL + EPS *   max( |a|,|b| ) ,
081: *
082: *          where EPS is the machine precision.  If ABSTOL is less than
083: *          or equal to zero, then  EPS*|T|  will be used in its place,
084: *          where |T| is the 1-norm of the tridiagonal matrix obtained
085: *          by reducing A to tridiagonal form.
086: *
087: *          Eigenvalues will be computed most accurately when ABSTOL is
088: *          set to twice the underflow threshold 2*SLAMCH('S'), not zero.
089: *          If this routine returns with INFO>0, indicating that some
090: *          eigenvectors did not converge, try setting ABSTOL to
091: *          2*SLAMCH('S').
092: *
093: *          See "Computing Small Singular Values of Bidiagonal Matrices
094: *          with Guaranteed High Relative Accuracy," by Demmel and
095: *          Kahan, LAPACK Working Note #3.
096: *
097: *  M       (output) INTEGER
098: *          The total number of eigenvalues found.  0 <= M <= N.
099: *          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
100: *
101: *  W       (output) REAL array, dimension (N)
102: *          On normal exit, the first M elements contain the selected
103: *          eigenvalues in ascending order.
104: *
105: *  Z       (output) COMPLEX array, dimension (LDZ, max(1,M))
106: *          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
107: *          contain the orthonormal eigenvectors of the matrix A
108: *          corresponding to the selected eigenvalues, with the i-th
109: *          column of Z holding the eigenvector associated with W(i).
110: *          If an eigenvector fails to converge, then that column of Z
111: *          contains the latest approximation to the eigenvector, and the
112: *          index of the eigenvector is returned in IFAIL.
113: *          If JOBZ = 'N', then Z is not referenced.
114: *          Note: the user must ensure that at least max(1,M) columns are
115: *          supplied in the array Z; if RANGE = 'V', the exact value of M
116: *          is not known in advance and an upper bound must be used.
117: *
118: *  LDZ     (input) INTEGER
119: *          The leading dimension of the array Z.  LDZ >= 1, and if
120: *          JOBZ = 'V', LDZ >= max(1,N).
121: *
122: *  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
123: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
124: *
125: *  LWORK   (input) INTEGER
126: *          The length of the array WORK.  LWORK >= 1, when N <= 1;
127: *          otherwise 2*N.
128: *          For optimal efficiency, LWORK >= (NB+1)*N,
129: *          where NB is the max of the blocksize for CHETRD and for
130: *          CUNMTR as returned by ILAENV.
131: *
132: *          If LWORK = -1, then a workspace query is assumed; the routine
133: *          only calculates the optimal size of the WORK array, returns
134: *          this value as the first entry of the WORK array, and no error
135: *          message related to LWORK is issued by XERBLA.
136: *
137: *  RWORK   (workspace) REAL array, dimension (7*N)
138: *
139: *  IWORK   (workspace) INTEGER array, dimension (5*N)
140: *
141: *  IFAIL   (output) INTEGER array, dimension (N)
142: *          If JOBZ = 'V', then if INFO = 0, the first M elements of
143: *          IFAIL are zero.  If INFO > 0, then IFAIL contains the
144: *          indices of the eigenvectors that failed to converge.
145: *          If JOBZ = 'N', then IFAIL is not referenced.
146: *
147: *  INFO    (output) INTEGER
148: *          = 0:  successful exit
149: *          < 0:  if INFO = -i, the i-th argument had an illegal value
150: *          > 0:  if INFO = i, then i eigenvectors failed to converge.
151: *                Their indices are stored in array IFAIL.
152: *
153: *  =====================================================================
154: *
155: *     .. Parameters ..
156:       REAL               ZERO, ONE
157:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
158:       COMPLEX            CONE
159:       PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
160: *     ..
161: *     .. Local Scalars ..
162:       LOGICAL            ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
163:      $                   WANTZ
164:       CHARACTER          ORDER
165:       INTEGER            I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
166:      $                   INDISP, INDIWK, INDRWK, INDTAU, INDWRK, ISCALE,
167:      $                   ITMP1, J, JJ, LLWORK, LWKMIN, LWKOPT, NB,
168:      $                   NSPLIT
169:       REAL               ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
170:      $                   SIGMA, SMLNUM, TMP1, VLL, VUU
171: *     ..
172: *     .. External Functions ..
173:       LOGICAL            LSAME
174:       INTEGER            ILAENV
175:       REAL               CLANHE, SLAMCH
176:       EXTERNAL           LSAME, ILAENV, CLANHE, SLAMCH
177: *     ..
178: *     .. External Subroutines ..
179:       EXTERNAL           CHETRD, CLACPY, CSSCAL, CSTEIN, CSTEQR, CSWAP,
180:      $                   CUNGTR, CUNMTR, SCOPY, SSCAL, SSTEBZ, SSTERF,
181:      $                   XERBLA
182: *     ..
183: *     .. Intrinsic Functions ..
184:       INTRINSIC          MAX, MIN, REAL, SQRT
185: *     ..
186: *     .. Executable Statements ..
187: *
188: *     Test the input parameters.
189: *
190:       LOWER = LSAME( UPLO, 'L' )
191:       WANTZ = LSAME( JOBZ, 'V' )
192:       ALLEIG = LSAME( RANGE, 'A' )
193:       VALEIG = LSAME( RANGE, 'V' )
194:       INDEIG = LSAME( RANGE, 'I' )
195:       LQUERY = ( LWORK.EQ.-1 )
196: *
197:       INFO = 0
198:       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
199:          INFO = -1
200:       ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
201:          INFO = -2
202:       ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
203:          INFO = -3
204:       ELSE IF( N.LT.0 ) THEN
205:          INFO = -4
206:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
207:          INFO = -6
208:       ELSE
209:          IF( VALEIG ) THEN
210:             IF( N.GT.0 .AND. VU.LE.VL )
211:      $         INFO = -8
212:          ELSE IF( INDEIG ) THEN
213:             IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN
214:                INFO = -9
215:             ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
216:                INFO = -10
217:             END IF
218:          END IF
219:       END IF
220:       IF( INFO.EQ.0 ) THEN
221:          IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
222:             INFO = -15
223:          END IF
224:       END IF
225: *
226:       IF( INFO.EQ.0 ) THEN
227:          IF( N.LE.1 ) THEN
228:             LWKMIN = 1
229:             WORK( 1 ) = LWKMIN
230:          ELSE
231:             LWKMIN = 2*N
232:             NB = ILAENV( 1, 'CHETRD', UPLO, N, -1, -1, -1 )
233:             NB = MAX( NB, ILAENV( 1, 'CUNMTR', UPLO, N, -1, -1, -1 ) )
234:             LWKOPT = MAX( 1, ( NB + 1 )*N )
235:             WORK( 1 ) = LWKOPT
236:          END IF
237: *
238:          IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY )
239:      $      INFO = -17
240:       END IF
241: *
242:       IF( INFO.NE.0 ) THEN
243:          CALL XERBLA( 'CHEEVX', -INFO )
244:          RETURN
245:       ELSE IF( LQUERY ) THEN
246:          RETURN
247:       END IF
248: *
249: *     Quick return if possible
250: *
251:       M = 0
252:       IF( N.EQ.0 ) THEN
253:          RETURN
254:       END IF
255: *
256:       IF( N.EQ.1 ) THEN
257:          IF( ALLEIG .OR. INDEIG ) THEN
258:             M = 1
259:             W( 1 ) = A( 1, 1 )
260:          ELSE IF( VALEIG ) THEN
261:             IF( VL.LT.REAL( A( 1, 1 ) ) .AND. VU.GE.REAL( A( 1, 1 ) ) )
262:      $           THEN
263:                M = 1
264:                W( 1 ) = A( 1, 1 )
265:             END IF
266:          END IF
267:          IF( WANTZ )
268:      $      Z( 1, 1 ) = CONE
269:          RETURN
270:       END IF
271: *
272: *     Get machine constants.
273: *
274:       SAFMIN = SLAMCH( 'Safe minimum' )
275:       EPS = SLAMCH( 'Precision' )
276:       SMLNUM = SAFMIN / EPS
277:       BIGNUM = ONE / SMLNUM
278:       RMIN = SQRT( SMLNUM )
279:       RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
280: *
281: *     Scale matrix to allowable range, if necessary.
282: *
283:       ISCALE = 0
284:       ABSTLL = ABSTOL
285:       IF( VALEIG ) THEN
286:          VLL = VL
287:          VUU = VU
288:       END IF
289:       ANRM = CLANHE( 'M', UPLO, N, A, LDA, RWORK )
290:       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
291:          ISCALE = 1
292:          SIGMA = RMIN / ANRM
293:       ELSE IF( ANRM.GT.RMAX ) THEN
294:          ISCALE = 1
295:          SIGMA = RMAX / ANRM
296:       END IF
297:       IF( ISCALE.EQ.1 ) THEN
298:          IF( LOWER ) THEN
299:             DO 10 J = 1, N
300:                CALL CSSCAL( N-J+1, SIGMA, A( J, J ), 1 )
301:    10       CONTINUE
302:          ELSE
303:             DO 20 J = 1, N
304:                CALL CSSCAL( J, SIGMA, A( 1, J ), 1 )
305:    20       CONTINUE
306:          END IF
307:          IF( ABSTOL.GT.0 )
308:      $      ABSTLL = ABSTOL*SIGMA
309:          IF( VALEIG ) THEN
310:             VLL = VL*SIGMA
311:             VUU = VU*SIGMA
312:          END IF
313:       END IF
314: *
315: *     Call CHETRD to reduce Hermitian matrix to tridiagonal form.
316: *
317:       INDD = 1
318:       INDE = INDD + N
319:       INDRWK = INDE + N
320:       INDTAU = 1
321:       INDWRK = INDTAU + N
322:       LLWORK = LWORK - INDWRK + 1
323:       CALL CHETRD( UPLO, N, A, LDA, RWORK( INDD ), RWORK( INDE ),
324:      $             WORK( INDTAU ), WORK( INDWRK ), LLWORK, IINFO )
325: *
326: *     If all eigenvalues are desired and ABSTOL is less than or equal to
327: *     zero, then call SSTERF or CUNGTR and CSTEQR.  If this fails for
328: *     some eigenvalue, then try SSTEBZ.
329: *
330:       TEST = .FALSE.
331:       IF( INDEIG ) THEN
332:          IF( IL.EQ.1 .AND. IU.EQ.N ) THEN
333:             TEST = .TRUE.
334:          END IF
335:       END IF
336:       IF( ( ALLEIG .OR. TEST ) .AND. ( ABSTOL.LE.ZERO ) ) THEN
337:          CALL SCOPY( N, RWORK( INDD ), 1, W, 1 )
338:          INDEE = INDRWK + 2*N
339:          IF( .NOT.WANTZ ) THEN
340:             CALL SCOPY( N-1, RWORK( INDE ), 1, RWORK( INDEE ), 1 )
341:             CALL SSTERF( N, W, RWORK( INDEE ), INFO )
342:          ELSE
343:             CALL CLACPY( 'A', N, N, A, LDA, Z, LDZ )
344:             CALL CUNGTR( UPLO, N, Z, LDZ, WORK( INDTAU ),
345:      $                   WORK( INDWRK ), LLWORK, IINFO )
346:             CALL SCOPY( N-1, RWORK( INDE ), 1, RWORK( INDEE ), 1 )
347:             CALL CSTEQR( JOBZ, N, W, RWORK( INDEE ), Z, LDZ,
348:      $                   RWORK( INDRWK ), INFO )
349:             IF( INFO.EQ.0 ) THEN
350:                DO 30 I = 1, N
351:                   IFAIL( I ) = 0
352:    30          CONTINUE
353:             END IF
354:          END IF
355:          IF( INFO.EQ.0 ) THEN
356:             M = N
357:             GO TO 40
358:          END IF
359:          INFO = 0
360:       END IF
361: *
362: *     Otherwise, call SSTEBZ and, if eigenvectors are desired, CSTEIN.
363: *
364:       IF( WANTZ ) THEN
365:          ORDER = 'B'
366:       ELSE
367:          ORDER = 'E'
368:       END IF
369:       INDIBL = 1
370:       INDISP = INDIBL + N
371:       INDIWK = INDISP + N
372:       CALL SSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL,
373:      $             RWORK( INDD ), RWORK( INDE ), M, NSPLIT, W,
374:      $             IWORK( INDIBL ), IWORK( INDISP ), RWORK( INDRWK ),
375:      $             IWORK( INDIWK ), INFO )
376: *
377:       IF( WANTZ ) THEN
378:          CALL CSTEIN( N, RWORK( INDD ), RWORK( INDE ), M, W,
379:      $                IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ,
380:      $                RWORK( INDRWK ), IWORK( INDIWK ), IFAIL, INFO )
381: *
382: *        Apply unitary matrix used in reduction to tridiagonal
383: *        form to eigenvectors returned by CSTEIN.
384: *
385:          CALL CUNMTR( 'L', UPLO, 'N', N, M, A, LDA, WORK( INDTAU ), Z,
386:      $                LDZ, WORK( INDWRK ), LLWORK, IINFO )
387:       END IF
388: *
389: *     If matrix was scaled, then rescale eigenvalues appropriately.
390: *
391:    40 CONTINUE
392:       IF( ISCALE.EQ.1 ) THEN
393:          IF( INFO.EQ.0 ) THEN
394:             IMAX = M
395:          ELSE
396:             IMAX = INFO - 1
397:          END IF
398:          CALL SSCAL( IMAX, ONE / SIGMA, W, 1 )
399:       END IF
400: *
401: *     If eigenvalues are not in order, then sort them, along with
402: *     eigenvectors.
403: *
404:       IF( WANTZ ) THEN
405:          DO 60 J = 1, M - 1
406:             I = 0
407:             TMP1 = W( J )
408:             DO 50 JJ = J + 1, M
409:                IF( W( JJ ).LT.TMP1 ) THEN
410:                   I = JJ
411:                   TMP1 = W( JJ )
412:                END IF
413:    50       CONTINUE
414: *
415:             IF( I.NE.0 ) THEN
416:                ITMP1 = IWORK( INDIBL+I-1 )
417:                W( I ) = W( J )
418:                IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 )
419:                W( J ) = TMP1
420:                IWORK( INDIBL+J-1 ) = ITMP1
421:                CALL CSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
422:                IF( INFO.NE.0 ) THEN
423:                   ITMP1 = IFAIL( I )
424:                   IFAIL( I ) = IFAIL( J )
425:                   IFAIL( J ) = ITMP1
426:                END IF
427:             END IF
428:    60    CONTINUE
429:       END IF
430: *
431: *     Set WORK(1) to optimal complex workspace size.
432: *
433:       WORK( 1 ) = LWKOPT
434: *
435:       RETURN
436: *
437: *     End of CHEEVX
438: *
439:       END
440: