Annotation of rpl/lapack/lapack/zheevr.f, revision 1.8
1.1 bertrand 1: SUBROUTINE ZHEEVR( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU,
2: $ ABSTOL, M, W, Z, LDZ, ISUPPZ, WORK, LWORK,
3: $ RWORK, LRWORK, IWORK, LIWORK, INFO )
4: *
1.5 bertrand 5: * -- LAPACK driver routine (version 3.2.2) --
1.1 bertrand 6: * -- LAPACK is a software package provided by Univ. of Tennessee, --
7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.5 bertrand 8: * June 2010
1.1 bertrand 9: *
10: * .. Scalar Arguments ..
11: CHARACTER JOBZ, RANGE, UPLO
12: INTEGER IL, INFO, IU, LDA, LDZ, LIWORK, LRWORK, LWORK,
13: $ M, N
14: DOUBLE PRECISION ABSTOL, VL, VU
15: * ..
16: * .. Array Arguments ..
17: INTEGER ISUPPZ( * ), IWORK( * )
18: DOUBLE PRECISION RWORK( * ), W( * )
19: COMPLEX*16 A( LDA, * ), WORK( * ), Z( LDZ, * )
20: * ..
21: *
22: * Purpose
23: * =======
24: *
25: * ZHEEVR computes selected eigenvalues and, optionally, eigenvectors
26: * of a complex Hermitian matrix A. Eigenvalues and eigenvectors can
27: * be selected by specifying either a range of values or a range of
28: * indices for the desired eigenvalues.
29: *
30: * ZHEEVR first reduces the matrix A to tridiagonal form T with a call
31: * to ZHETRD. Then, whenever possible, ZHEEVR calls ZSTEMR to compute
32: * eigenspectrum using Relatively Robust Representations. ZSTEMR
33: * computes eigenvalues by the dqds algorithm, while orthogonal
34: * eigenvectors are computed from various "good" L D L^T representations
35: * (also known as Relatively Robust Representations). Gram-Schmidt
36: * orthogonalization is avoided as far as possible. More specifically,
37: * the various steps of the algorithm are as follows.
38: *
39: * For each unreduced block (submatrix) of T,
40: * (a) Compute T - sigma I = L D L^T, so that L and D
41: * define all the wanted eigenvalues to high relative accuracy.
42: * This means that small relative changes in the entries of D and L
43: * cause only small relative changes in the eigenvalues and
44: * eigenvectors. The standard (unfactored) representation of the
45: * tridiagonal matrix T does not have this property in general.
46: * (b) Compute the eigenvalues to suitable accuracy.
47: * If the eigenvectors are desired, the algorithm attains full
48: * accuracy of the computed eigenvalues only right before
49: * the corresponding vectors have to be computed, see steps c) and d).
50: * (c) For each cluster of close eigenvalues, select a new
51: * shift close to the cluster, find a new factorization, and refine
52: * the shifted eigenvalues to suitable accuracy.
53: * (d) For each eigenvalue with a large enough relative separation compute
54: * the corresponding eigenvector by forming a rank revealing twisted
55: * factorization. Go back to (c) for any clusters that remain.
56: *
57: * The desired accuracy of the output can be specified by the input
58: * parameter ABSTOL.
59: *
60: * For more details, see DSTEMR's documentation and:
61: * - Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representations
62: * to compute orthogonal eigenvectors of symmetric tridiagonal matrices,"
63: * Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004.
64: * - Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and
65: * Relative Gaps," SIAM Journal on Matrix Analysis and Applications, Vol. 25,
66: * 2004. Also LAPACK Working Note 154.
67: * - Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric
68: * tridiagonal eigenvalue/eigenvector problem",
69: * Computer Science Division Technical Report No. UCB/CSD-97-971,
70: * UC Berkeley, May 1997.
71: *
72: *
73: * Note 1 : ZHEEVR calls ZSTEMR when the full spectrum is requested
74: * on machines which conform to the ieee-754 floating point standard.
75: * ZHEEVR calls DSTEBZ and ZSTEIN on non-ieee machines and
76: * when partial spectrum requests are made.
77: *
78: * Normal execution of ZSTEMR may create NaNs and infinities and
79: * hence may abort due to a floating point exception in environments
80: * which do not handle NaNs and infinities in the ieee standard default
81: * manner.
82: *
83: * Arguments
84: * =========
85: *
86: * JOBZ (input) CHARACTER*1
87: * = 'N': Compute eigenvalues only;
88: * = 'V': Compute eigenvalues and eigenvectors.
89: *
90: * RANGE (input) CHARACTER*1
91: * = 'A': all eigenvalues will be found.
92: * = 'V': all eigenvalues in the half-open interval (VL,VU]
93: * will be found.
94: * = 'I': the IL-th through IU-th eigenvalues will be found.
95: ********** For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and
96: ********** ZSTEIN are called
97: *
98: * UPLO (input) CHARACTER*1
99: * = 'U': Upper triangle of A is stored;
100: * = 'L': Lower triangle of A is stored.
101: *
102: * N (input) INTEGER
103: * The order of the matrix A. N >= 0.
104: *
105: * A (input/output) COMPLEX*16 array, dimension (LDA, N)
106: * On entry, the Hermitian matrix A. If UPLO = 'U', the
107: * leading N-by-N upper triangular part of A contains the
108: * upper triangular part of the matrix A. If UPLO = 'L',
109: * the leading N-by-N lower triangular part of A contains
110: * the lower triangular part of the matrix A.
111: * On exit, the lower triangle (if UPLO='L') or the upper
112: * triangle (if UPLO='U') of A, including the diagonal, is
113: * destroyed.
114: *
115: * LDA (input) INTEGER
116: * The leading dimension of the array A. LDA >= max(1,N).
117: *
118: * VL (input) DOUBLE PRECISION
119: * VU (input) DOUBLE PRECISION
120: * If RANGE='V', the lower and upper bounds of the interval to
121: * be searched for eigenvalues. VL < VU.
122: * Not referenced if RANGE = 'A' or 'I'.
123: *
124: * IL (input) INTEGER
125: * IU (input) INTEGER
126: * If RANGE='I', the indices (in ascending order) of the
127: * smallest and largest eigenvalues to be returned.
128: * 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
129: * Not referenced if RANGE = 'A' or 'V'.
130: *
131: * ABSTOL (input) DOUBLE PRECISION
132: * The absolute error tolerance for the eigenvalues.
133: * An approximate eigenvalue is accepted as converged
134: * when it is determined to lie in an interval [a,b]
135: * of width less than or equal to
136: *
137: * ABSTOL + EPS * max( |a|,|b| ) ,
138: *
139: * where EPS is the machine precision. If ABSTOL is less than
140: * or equal to zero, then EPS*|T| will be used in its place,
141: * where |T| is the 1-norm of the tridiagonal matrix obtained
142: * by reducing A to tridiagonal form.
143: *
144: * See "Computing Small Singular Values of Bidiagonal Matrices
145: * with Guaranteed High Relative Accuracy," by Demmel and
146: * Kahan, LAPACK Working Note #3.
147: *
148: * If high relative accuracy is important, set ABSTOL to
149: * DLAMCH( 'Safe minimum' ). Doing so will guarantee that
150: * eigenvalues are computed to high relative accuracy when
151: * possible in future releases. The current code does not
152: * make any guarantees about high relative accuracy, but
153: * furutre releases will. See J. Barlow and J. Demmel,
154: * "Computing Accurate Eigensystems of Scaled Diagonally
155: * Dominant Matrices", LAPACK Working Note #7, for a discussion
156: * of which matrices define their eigenvalues to high relative
157: * accuracy.
158: *
159: * M (output) INTEGER
160: * The total number of eigenvalues found. 0 <= M <= N.
161: * If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
162: *
163: * W (output) DOUBLE PRECISION array, dimension (N)
164: * The first M elements contain the selected eigenvalues in
165: * ascending order.
166: *
167: * Z (output) COMPLEX*16 array, dimension (LDZ, max(1,M))
168: * If JOBZ = 'V', then if INFO = 0, the first M columns of Z
169: * contain the orthonormal eigenvectors of the matrix A
170: * corresponding to the selected eigenvalues, with the i-th
171: * column of Z holding the eigenvector associated with W(i).
172: * If JOBZ = 'N', then Z is not referenced.
173: * Note: the user must ensure that at least max(1,M) columns are
174: * supplied in the array Z; if RANGE = 'V', the exact value of M
175: * is not known in advance and an upper bound must be used.
176: *
177: * LDZ (input) INTEGER
178: * The leading dimension of the array Z. LDZ >= 1, and if
179: * JOBZ = 'V', LDZ >= max(1,N).
180: *
181: * ISUPPZ (output) INTEGER array, dimension ( 2*max(1,M) )
182: * The support of the eigenvectors in Z, i.e., the indices
183: * indicating the nonzero elements in Z. The i-th eigenvector
184: * is nonzero only in elements ISUPPZ( 2*i-1 ) through
185: * ISUPPZ( 2*i ).
186: ********** Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
187: *
188: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
189: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
190: *
191: * LWORK (input) INTEGER
192: * The length of the array WORK. LWORK >= max(1,2*N).
193: * For optimal efficiency, LWORK >= (NB+1)*N,
194: * where NB is the max of the blocksize for ZHETRD and for
195: * ZUNMTR as returned by ILAENV.
196: *
197: * If LWORK = -1, then a workspace query is assumed; the routine
198: * only calculates the optimal sizes of the WORK, RWORK and
199: * IWORK arrays, returns these values as the first entries of
200: * the WORK, RWORK and IWORK arrays, and no error message
201: * related to LWORK or LRWORK or LIWORK is issued by XERBLA.
202: *
203: * RWORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
204: * On exit, if INFO = 0, RWORK(1) returns the optimal
205: * (and minimal) LRWORK.
206: *
207: * LRWORK (input) INTEGER
208: * The length of the array RWORK. LRWORK >= max(1,24*N).
209: *
210: * If LRWORK = -1, then a workspace query is assumed; the
211: * routine only calculates the optimal sizes of the WORK, RWORK
212: * and IWORK arrays, returns these values as the first entries
213: * of the WORK, RWORK and IWORK arrays, and no error message
214: * related to LWORK or LRWORK or LIWORK is issued by XERBLA.
215: *
216: * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
217: * On exit, if INFO = 0, IWORK(1) returns the optimal
218: * (and minimal) LIWORK.
219: *
220: * LIWORK (input) INTEGER
221: * The dimension of the array IWORK. LIWORK >= max(1,10*N).
222: *
223: * If LIWORK = -1, then a workspace query is assumed; the
224: * routine only calculates the optimal sizes of the WORK, RWORK
225: * and IWORK arrays, returns these values as the first entries
226: * of the WORK, RWORK and IWORK arrays, and no error message
227: * related to LWORK or LRWORK or LIWORK is issued by XERBLA.
228: *
229: * INFO (output) INTEGER
230: * = 0: successful exit
231: * < 0: if INFO = -i, the i-th argument had an illegal value
232: * > 0: Internal error
233: *
234: * Further Details
235: * ===============
236: *
237: * Based on contributions by
238: * Inderjit Dhillon, IBM Almaden, USA
239: * Osni Marques, LBNL/NERSC, USA
240: * Ken Stanley, Computer Science Division, University of
241: * California at Berkeley, USA
242: * Jason Riedy, Computer Science Division, University of
243: * California at Berkeley, USA
244: *
245: * =====================================================================
246: *
247: * .. Parameters ..
248: DOUBLE PRECISION ZERO, ONE, TWO
249: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
250: * ..
251: * .. Local Scalars ..
252: LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
253: $ WANTZ, TRYRAC
254: CHARACTER ORDER
255: INTEGER I, IEEEOK, IINFO, IMAX, INDIBL, INDIFL, INDISP,
256: $ INDIWO, INDRD, INDRDD, INDRE, INDREE, INDRWK,
257: $ INDTAU, INDWK, INDWKN, ISCALE, ITMP1, J, JJ,
258: $ LIWMIN, LLWORK, LLRWORK, LLWRKN, LRWMIN,
259: $ LWKOPT, LWMIN, NB, NSPLIT
260: DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
261: $ SIGMA, SMLNUM, TMP1, VLL, VUU
262: * ..
263: * .. External Functions ..
264: LOGICAL LSAME
265: INTEGER ILAENV
266: DOUBLE PRECISION DLAMCH, ZLANSY
267: EXTERNAL LSAME, ILAENV, DLAMCH, ZLANSY
268: * ..
269: * .. External Subroutines ..
270: EXTERNAL DCOPY, DSCAL, DSTEBZ, DSTERF, XERBLA, ZDSCAL,
271: $ ZHETRD, ZSTEMR, ZSTEIN, ZSWAP, ZUNMTR
272: * ..
273: * .. Intrinsic Functions ..
274: INTRINSIC DBLE, MAX, MIN, SQRT
275: * ..
276: * .. Executable Statements ..
277: *
278: * Test the input parameters.
279: *
280: IEEEOK = ILAENV( 10, 'ZHEEVR', 'N', 1, 2, 3, 4 )
281: *
282: LOWER = LSAME( UPLO, 'L' )
283: WANTZ = LSAME( JOBZ, 'V' )
284: ALLEIG = LSAME( RANGE, 'A' )
285: VALEIG = LSAME( RANGE, 'V' )
286: INDEIG = LSAME( RANGE, 'I' )
287: *
288: LQUERY = ( ( LWORK.EQ.-1 ) .OR. ( LRWORK.EQ.-1 ) .OR.
289: $ ( LIWORK.EQ.-1 ) )
290: *
291: LRWMIN = MAX( 1, 24*N )
292: LIWMIN = MAX( 1, 10*N )
293: LWMIN = MAX( 1, 2*N )
294: *
295: INFO = 0
296: IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
297: INFO = -1
298: ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
299: INFO = -2
300: ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
301: INFO = -3
302: ELSE IF( N.LT.0 ) THEN
303: INFO = -4
304: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
305: INFO = -6
306: ELSE
307: IF( VALEIG ) THEN
308: IF( N.GT.0 .AND. VU.LE.VL )
309: $ INFO = -8
310: ELSE IF( INDEIG ) THEN
311: IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN
312: INFO = -9
313: ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
314: INFO = -10
315: END IF
316: END IF
317: END IF
318: IF( INFO.EQ.0 ) THEN
319: IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
320: INFO = -15
321: END IF
322: END IF
323: *
324: IF( INFO.EQ.0 ) THEN
325: NB = ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 )
326: NB = MAX( NB, ILAENV( 1, 'ZUNMTR', UPLO, N, -1, -1, -1 ) )
327: LWKOPT = MAX( ( NB+1 )*N, LWMIN )
328: WORK( 1 ) = LWKOPT
329: RWORK( 1 ) = LRWMIN
330: IWORK( 1 ) = LIWMIN
331: *
332: IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
333: INFO = -18
334: ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
335: INFO = -20
336: ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
337: INFO = -22
338: END IF
339: END IF
340: *
341: IF( INFO.NE.0 ) THEN
342: CALL XERBLA( 'ZHEEVR', -INFO )
343: RETURN
344: ELSE IF( LQUERY ) THEN
345: RETURN
346: END IF
347: *
348: * Quick return if possible
349: *
350: M = 0
351: IF( N.EQ.0 ) THEN
352: WORK( 1 ) = 1
353: RETURN
354: END IF
355: *
356: IF( N.EQ.1 ) THEN
357: WORK( 1 ) = 2
358: IF( ALLEIG .OR. INDEIG ) THEN
359: M = 1
360: W( 1 ) = DBLE( A( 1, 1 ) )
361: ELSE
362: IF( VL.LT.DBLE( A( 1, 1 ) ) .AND. VU.GE.DBLE( A( 1, 1 ) ) )
363: $ THEN
364: M = 1
365: W( 1 ) = DBLE( A( 1, 1 ) )
366: END IF
367: END IF
1.5 bertrand 368: IF( WANTZ ) THEN
369: Z( 1, 1 ) = ONE
370: ISUPPZ( 1 ) = 1
371: ISUPPZ( 2 ) = 1
372: END IF
1.1 bertrand 373: RETURN
374: END IF
375: *
376: * Get machine constants.
377: *
378: SAFMIN = DLAMCH( 'Safe minimum' )
379: EPS = DLAMCH( 'Precision' )
380: SMLNUM = SAFMIN / EPS
381: BIGNUM = ONE / SMLNUM
382: RMIN = SQRT( SMLNUM )
383: RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
384: *
385: * Scale matrix to allowable range, if necessary.
386: *
387: ISCALE = 0
388: ABSTLL = ABSTOL
389: IF (VALEIG) THEN
390: VLL = VL
391: VUU = VU
392: END IF
393: ANRM = ZLANSY( 'M', UPLO, N, A, LDA, RWORK )
394: IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
395: ISCALE = 1
396: SIGMA = RMIN / ANRM
397: ELSE IF( ANRM.GT.RMAX ) THEN
398: ISCALE = 1
399: SIGMA = RMAX / ANRM
400: END IF
401: IF( ISCALE.EQ.1 ) THEN
402: IF( LOWER ) THEN
403: DO 10 J = 1, N
404: CALL ZDSCAL( N-J+1, SIGMA, A( J, J ), 1 )
405: 10 CONTINUE
406: ELSE
407: DO 20 J = 1, N
408: CALL ZDSCAL( J, SIGMA, A( 1, J ), 1 )
409: 20 CONTINUE
410: END IF
411: IF( ABSTOL.GT.0 )
412: $ ABSTLL = ABSTOL*SIGMA
413: IF( VALEIG ) THEN
414: VLL = VL*SIGMA
415: VUU = VU*SIGMA
416: END IF
417: END IF
418:
419: * Initialize indices into workspaces. Note: The IWORK indices are
420: * used only if DSTERF or ZSTEMR fail.
421:
422: * WORK(INDTAU:INDTAU+N-1) stores the complex scalar factors of the
423: * elementary reflectors used in ZHETRD.
424: INDTAU = 1
425: * INDWK is the starting offset of the remaining complex workspace,
426: * and LLWORK is the remaining complex workspace size.
427: INDWK = INDTAU + N
428: LLWORK = LWORK - INDWK + 1
429:
430: * RWORK(INDRD:INDRD+N-1) stores the real tridiagonal's diagonal
431: * entries.
432: INDRD = 1
433: * RWORK(INDRE:INDRE+N-1) stores the off-diagonal entries of the
434: * tridiagonal matrix from ZHETRD.
435: INDRE = INDRD + N
436: * RWORK(INDRDD:INDRDD+N-1) is a copy of the diagonal entries over
437: * -written by ZSTEMR (the DSTERF path copies the diagonal to W).
438: INDRDD = INDRE + N
439: * RWORK(INDREE:INDREE+N-1) is a copy of the off-diagonal entries over
440: * -written while computing the eigenvalues in DSTERF and ZSTEMR.
441: INDREE = INDRDD + N
442: * INDRWK is the starting offset of the left-over real workspace, and
443: * LLRWORK is the remaining workspace size.
444: INDRWK = INDREE + N
445: LLRWORK = LRWORK - INDRWK + 1
446:
447: * IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in DSTEBZ and
448: * stores the block indices of each of the M<=N eigenvalues.
449: INDIBL = 1
450: * IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in DSTEBZ and
451: * stores the starting and finishing indices of each block.
452: INDISP = INDIBL + N
453: * IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors
454: * that corresponding to eigenvectors that fail to converge in
455: * DSTEIN. This information is discarded; if any fail, the driver
456: * returns INFO > 0.
457: INDIFL = INDISP + N
458: * INDIWO is the offset of the remaining integer workspace.
459: INDIWO = INDISP + N
460:
461: *
462: * Call ZHETRD to reduce Hermitian matrix to tridiagonal form.
463: *
464: CALL ZHETRD( UPLO, N, A, LDA, RWORK( INDRD ), RWORK( INDRE ),
465: $ WORK( INDTAU ), WORK( INDWK ), LLWORK, IINFO )
466: *
467: * If all eigenvalues are desired
468: * then call DSTERF or ZSTEMR and ZUNMTR.
469: *
470: TEST = .FALSE.
471: IF( INDEIG ) THEN
472: IF( IL.EQ.1 .AND. IU.EQ.N ) THEN
473: TEST = .TRUE.
474: END IF
475: END IF
476: IF( ( ALLEIG.OR.TEST ) .AND. ( IEEEOK.EQ.1 ) ) THEN
477: IF( .NOT.WANTZ ) THEN
478: CALL DCOPY( N, RWORK( INDRD ), 1, W, 1 )
479: CALL DCOPY( N-1, RWORK( INDRE ), 1, RWORK( INDREE ), 1 )
480: CALL DSTERF( N, W, RWORK( INDREE ), INFO )
481: ELSE
482: CALL DCOPY( N-1, RWORK( INDRE ), 1, RWORK( INDREE ), 1 )
483: CALL DCOPY( N, RWORK( INDRD ), 1, RWORK( INDRDD ), 1 )
484: *
485: IF (ABSTOL .LE. TWO*N*EPS) THEN
486: TRYRAC = .TRUE.
487: ELSE
488: TRYRAC = .FALSE.
489: END IF
490: CALL ZSTEMR( JOBZ, 'A', N, RWORK( INDRDD ),
491: $ RWORK( INDREE ), VL, VU, IL, IU, M, W,
492: $ Z, LDZ, N, ISUPPZ, TRYRAC,
493: $ RWORK( INDRWK ), LLRWORK,
494: $ IWORK, LIWORK, INFO )
495: *
496: * Apply unitary matrix used in reduction to tridiagonal
497: * form to eigenvectors returned by ZSTEIN.
498: *
499: IF( WANTZ .AND. INFO.EQ.0 ) THEN
500: INDWKN = INDWK
501: LLWRKN = LWORK - INDWKN + 1
502: CALL ZUNMTR( 'L', UPLO, 'N', N, M, A, LDA,
503: $ WORK( INDTAU ), Z, LDZ, WORK( INDWKN ),
504: $ LLWRKN, IINFO )
505: END IF
506: END IF
507: *
508: *
509: IF( INFO.EQ.0 ) THEN
510: M = N
511: GO TO 30
512: END IF
513: INFO = 0
514: END IF
515: *
516: * Otherwise, call DSTEBZ and, if eigenvectors are desired, ZSTEIN.
517: * Also call DSTEBZ and ZSTEIN if ZSTEMR fails.
518: *
519: IF( WANTZ ) THEN
520: ORDER = 'B'
521: ELSE
522: ORDER = 'E'
523: END IF
524:
525: CALL DSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL,
526: $ RWORK( INDRD ), RWORK( INDRE ), M, NSPLIT, W,
527: $ IWORK( INDIBL ), IWORK( INDISP ), RWORK( INDRWK ),
528: $ IWORK( INDIWO ), INFO )
529: *
530: IF( WANTZ ) THEN
531: CALL ZSTEIN( N, RWORK( INDRD ), RWORK( INDRE ), M, W,
532: $ IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ,
533: $ RWORK( INDRWK ), IWORK( INDIWO ), IWORK( INDIFL ),
534: $ INFO )
535: *
536: * Apply unitary matrix used in reduction to tridiagonal
537: * form to eigenvectors returned by ZSTEIN.
538: *
539: INDWKN = INDWK
540: LLWRKN = LWORK - INDWKN + 1
541: CALL ZUNMTR( 'L', UPLO, 'N', N, M, A, LDA, WORK( INDTAU ), Z,
542: $ LDZ, WORK( INDWKN ), LLWRKN, IINFO )
543: END IF
544: *
545: * If matrix was scaled, then rescale eigenvalues appropriately.
546: *
547: 30 CONTINUE
548: IF( ISCALE.EQ.1 ) THEN
549: IF( INFO.EQ.0 ) THEN
550: IMAX = M
551: ELSE
552: IMAX = INFO - 1
553: END IF
554: CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
555: END IF
556: *
557: * If eigenvalues are not in order, then sort them, along with
558: * eigenvectors.
559: *
560: IF( WANTZ ) THEN
561: DO 50 J = 1, M - 1
562: I = 0
563: TMP1 = W( J )
564: DO 40 JJ = J + 1, M
565: IF( W( JJ ).LT.TMP1 ) THEN
566: I = JJ
567: TMP1 = W( JJ )
568: END IF
569: 40 CONTINUE
570: *
571: IF( I.NE.0 ) THEN
572: ITMP1 = IWORK( INDIBL+I-1 )
573: W( I ) = W( J )
574: IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 )
575: W( J ) = TMP1
576: IWORK( INDIBL+J-1 ) = ITMP1
577: CALL ZSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
578: END IF
579: 50 CONTINUE
580: END IF
581: *
582: * Set WORK(1) to optimal workspace size.
583: *
584: WORK( 1 ) = LWKOPT
585: RWORK( 1 ) = LRWMIN
586: IWORK( 1 ) = LIWMIN
587: *
588: RETURN
589: *
590: * End of ZHEEVR
591: *
592: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zheevr.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack