Annotation of rpl/lapack/lapack/dstemr.f, revision 1.3
1.1 bertrand 1: SUBROUTINE DSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
2: $ M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK,
3: $ IWORK, LIWORK, INFO )
4: IMPLICIT NONE
5: *
6: * -- LAPACK computational routine (version 3.2.1) --
7: * -- LAPACK is a software package provided by Univ. of Tennessee, --
8: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
9: * -- April 2009 --
10: *
11: * .. Scalar Arguments ..
12: CHARACTER JOBZ, RANGE
13: LOGICAL TRYRAC
14: INTEGER IL, INFO, IU, LDZ, NZC, LIWORK, LWORK, M, N
15: DOUBLE PRECISION VL, VU
16: * ..
17: * .. Array Arguments ..
18: INTEGER ISUPPZ( * ), IWORK( * )
19: DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
20: DOUBLE PRECISION Z( LDZ, * )
21: * ..
22: *
23: * Purpose
24: * =======
25: *
26: * DSTEMR computes selected eigenvalues and, optionally, eigenvectors
27: * of a real symmetric tridiagonal matrix T. Any such unreduced matrix has
28: * a well defined set of pairwise different real eigenvalues, the corresponding
29: * real eigenvectors are pairwise orthogonal.
30: *
31: * The spectrum may be computed either completely or partially by specifying
32: * either an interval (VL,VU] or a range of indices IL:IU for the desired
33: * eigenvalues.
34: *
35: * Depending on the number of desired eigenvalues, these are computed either
36: * by bisection or the dqds algorithm. Numerically orthogonal eigenvectors are
37: * computed by the use of various suitable L D L^T factorizations near clusters
38: * of close eigenvalues (referred to as RRRs, Relatively Robust
39: * Representations). An informal sketch of the algorithm follows.
40: *
41: * For each unreduced block (submatrix) of T,
42: * (a) Compute T - sigma I = L D L^T, so that L and D
43: * define all the wanted eigenvalues to high relative accuracy.
44: * This means that small relative changes in the entries of D and L
45: * cause only small relative changes in the eigenvalues and
46: * eigenvectors. The standard (unfactored) representation of the
47: * tridiagonal matrix T does not have this property in general.
48: * (b) Compute the eigenvalues to suitable accuracy.
49: * If the eigenvectors are desired, the algorithm attains full
50: * accuracy of the computed eigenvalues only right before
51: * the corresponding vectors have to be computed, see steps c) and d).
52: * (c) For each cluster of close eigenvalues, select a new
53: * shift close to the cluster, find a new factorization, and refine
54: * the shifted eigenvalues to suitable accuracy.
55: * (d) For each eigenvalue with a large enough relative separation compute
56: * the corresponding eigenvector by forming a rank revealing twisted
57: * factorization. Go back to (c) for any clusters that remain.
58: *
59: * For more details, see:
60: * - Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representations
61: * to compute orthogonal eigenvectors of symmetric tridiagonal matrices,"
62: * Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004.
63: * - Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and
64: * Relative Gaps," SIAM Journal on Matrix Analysis and Applications, Vol. 25,
65: * 2004. Also LAPACK Working Note 154.
66: * - Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric
67: * tridiagonal eigenvalue/eigenvector problem",
68: * Computer Science Division Technical Report No. UCB/CSD-97-971,
69: * UC Berkeley, May 1997.
70: *
71: * Further Details
72: * 1.DSTEMR works only on machines which follow IEEE-754
73: * floating-point standard in their handling of infinities and NaNs.
74: * This permits the use of efficient inner loops avoiding a check for
75: * zero divisors.
76: *
77: * Arguments
78: * =========
79: *
80: * JOBZ (input) CHARACTER*1
81: * = 'N': Compute eigenvalues only;
82: * = 'V': Compute eigenvalues and eigenvectors.
83: *
84: * RANGE (input) CHARACTER*1
85: * = 'A': all eigenvalues will be found.
86: * = 'V': all eigenvalues in the half-open interval (VL,VU]
87: * will be found.
88: * = 'I': the IL-th through IU-th eigenvalues will be found.
89: *
90: * N (input) INTEGER
91: * The order of the matrix. N >= 0.
92: *
93: * D (input/output) DOUBLE PRECISION array, dimension (N)
94: * On entry, the N diagonal elements of the tridiagonal matrix
95: * T. On exit, D is overwritten.
96: *
97: * E (input/output) DOUBLE PRECISION array, dimension (N)
98: * On entry, the (N-1) subdiagonal elements of the tridiagonal
99: * matrix T in elements 1 to N-1 of E. E(N) need not be set on
100: * input, but is used internally as workspace.
101: * On exit, E is overwritten.
102: *
103: * VL (input) DOUBLE PRECISION
104: * VU (input) DOUBLE PRECISION
105: * If RANGE='V', the lower and upper bounds of the interval to
106: * be searched for eigenvalues. VL < VU.
107: * Not referenced if RANGE = 'A' or 'I'.
108: *
109: * IL (input) INTEGER
110: * IU (input) INTEGER
111: * If RANGE='I', the indices (in ascending order) of the
112: * smallest and largest eigenvalues to be returned.
113: * 1 <= IL <= IU <= N, if N > 0.
114: * Not referenced if RANGE = 'A' or 'V'.
115: *
116: * M (output) INTEGER
117: * The total number of eigenvalues found. 0 <= M <= N.
118: * If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
119: *
120: * W (output) DOUBLE PRECISION array, dimension (N)
121: * The first M elements contain the selected eigenvalues in
122: * ascending order.
123: *
124: * Z (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
125: * If JOBZ = 'V', and if INFO = 0, then the first M columns of Z
126: * contain the orthonormal eigenvectors of the matrix T
127: * corresponding to the selected eigenvalues, with the i-th
128: * column of Z holding the eigenvector associated with W(i).
129: * If JOBZ = 'N', then Z is not referenced.
130: * Note: the user must ensure that at least max(1,M) columns are
131: * supplied in the array Z; if RANGE = 'V', the exact value of M
132: * is not known in advance and can be computed with a workspace
133: * query by setting NZC = -1, see below.
134: *
135: * LDZ (input) INTEGER
136: * The leading dimension of the array Z. LDZ >= 1, and if
137: * JOBZ = 'V', then LDZ >= max(1,N).
138: *
139: * NZC (input) INTEGER
140: * The number of eigenvectors to be held in the array Z.
141: * If RANGE = 'A', then NZC >= max(1,N).
142: * If RANGE = 'V', then NZC >= the number of eigenvalues in (VL,VU].
143: * If RANGE = 'I', then NZC >= IU-IL+1.
144: * If NZC = -1, then a workspace query is assumed; the
145: * routine calculates the number of columns of the array Z that
146: * are needed to hold the eigenvectors.
147: * This value is returned as the first entry of the Z array, and
148: * no error message related to NZC is issued by XERBLA.
149: *
150: * ISUPPZ (output) INTEGER ARRAY, dimension ( 2*max(1,M) )
151: * The support of the eigenvectors in Z, i.e., the indices
152: * indicating the nonzero elements in Z. The i-th computed eigenvector
153: * is nonzero only in elements ISUPPZ( 2*i-1 ) through
154: * ISUPPZ( 2*i ). This is relevant in the case when the matrix
155: * is split. ISUPPZ is only accessed when JOBZ is 'V' and N > 0.
156: *
157: * TRYRAC (input/output) LOGICAL
158: * If TRYRAC.EQ..TRUE., indicates that the code should check whether
159: * the tridiagonal matrix defines its eigenvalues to high relative
160: * accuracy. If so, the code uses relative-accuracy preserving
161: * algorithms that might be (a bit) slower depending on the matrix.
162: * If the matrix does not define its eigenvalues to high relative
163: * accuracy, the code can uses possibly faster algorithms.
164: * If TRYRAC.EQ..FALSE., the code is not required to guarantee
165: * relatively accurate eigenvalues and can use the fastest possible
166: * techniques.
167: * On exit, a .TRUE. TRYRAC will be set to .FALSE. if the matrix
168: * does not define its eigenvalues to high relative accuracy.
169: *
170: * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
171: * On exit, if INFO = 0, WORK(1) returns the optimal
172: * (and minimal) LWORK.
173: *
174: * LWORK (input) INTEGER
175: * The dimension of the array WORK. LWORK >= max(1,18*N)
176: * if JOBZ = 'V', and LWORK >= max(1,12*N) if JOBZ = 'N'.
177: * If LWORK = -1, then a workspace query is assumed; the routine
178: * only calculates the optimal size of the WORK array, returns
179: * this value as the first entry of the WORK array, and no error
180: * message related to LWORK is issued by XERBLA.
181: *
182: * IWORK (workspace/output) INTEGER array, dimension (LIWORK)
183: * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
184: *
185: * LIWORK (input) INTEGER
186: * The dimension of the array IWORK. LIWORK >= max(1,10*N)
187: * if the eigenvectors are desired, and LIWORK >= max(1,8*N)
188: * if only the eigenvalues are to be computed.
189: * If LIWORK = -1, then a workspace query is assumed; the
190: * routine only calculates the optimal size of the IWORK array,
191: * returns this value as the first entry of the IWORK array, and
192: * no error message related to LIWORK is issued by XERBLA.
193: *
194: * INFO (output) INTEGER
195: * On exit, INFO
196: * = 0: successful exit
197: * < 0: if INFO = -i, the i-th argument had an illegal value
198: * > 0: if INFO = 1X, internal error in DLARRE,
199: * if INFO = 2X, internal error in DLARRV.
200: * Here, the digit X = ABS( IINFO ) < 10, where IINFO is
201: * the nonzero error code returned by DLARRE or
202: * DLARRV, respectively.
203: *
204: *
205: * Further Details
206: * ===============
207: *
208: * Based on contributions by
209: * Beresford Parlett, University of California, Berkeley, USA
210: * Jim Demmel, University of California, Berkeley, USA
211: * Inderjit Dhillon, University of Texas, Austin, USA
212: * Osni Marques, LBNL/NERSC, USA
213: * Christof Voemel, University of California, Berkeley, USA
214: *
215: * =====================================================================
216: *
217: * .. Parameters ..
218: DOUBLE PRECISION ZERO, ONE, FOUR, MINRGP
219: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0,
220: $ FOUR = 4.0D0,
221: $ MINRGP = 1.0D-3 )
222: * ..
223: * .. Local Scalars ..
224: LOGICAL ALLEIG, INDEIG, LQUERY, VALEIG, WANTZ, ZQUERY
225: INTEGER I, IBEGIN, IEND, IFIRST, IIL, IINDBL, IINDW,
226: $ IINDWK, IINFO, IINSPL, IIU, ILAST, IN, INDD,
227: $ INDE2, INDERR, INDGP, INDGRS, INDWRK, ITMP,
228: $ ITMP2, J, JBLK, JJ, LIWMIN, LWMIN, NSPLIT,
229: $ NZCMIN, OFFSET, WBEGIN, WEND
230: DOUBLE PRECISION BIGNUM, CS, EPS, PIVMIN, R1, R2, RMAX, RMIN,
231: $ RTOL1, RTOL2, SAFMIN, SCALE, SMLNUM, SN,
232: $ THRESH, TMP, TNRM, WL, WU
233: * ..
234: * ..
235: * .. External Functions ..
236: LOGICAL LSAME
237: DOUBLE PRECISION DLAMCH, DLANST
238: EXTERNAL LSAME, DLAMCH, DLANST
239: * ..
240: * .. External Subroutines ..
241: EXTERNAL DCOPY, DLAE2, DLAEV2, DLARRC, DLARRE, DLARRJ,
242: $ DLARRR, DLARRV, DLASRT, DSCAL, DSWAP, XERBLA
243: * ..
244: * .. Intrinsic Functions ..
245: INTRINSIC MAX, MIN, SQRT
246:
247:
248: * ..
249: * .. Executable Statements ..
250: *
251: * Test the input parameters.
252: *
253: WANTZ = LSAME( JOBZ, 'V' )
254: ALLEIG = LSAME( RANGE, 'A' )
255: VALEIG = LSAME( RANGE, 'V' )
256: INDEIG = LSAME( RANGE, 'I' )
257: *
258: LQUERY = ( ( LWORK.EQ.-1 ).OR.( LIWORK.EQ.-1 ) )
259: ZQUERY = ( NZC.EQ.-1 )
260:
261: * DSTEMR needs WORK of size 6*N, IWORK of size 3*N.
262: * In addition, DLARRE needs WORK of size 6*N, IWORK of size 5*N.
263: * Furthermore, DLARRV needs WORK of size 12*N, IWORK of size 7*N.
264: IF( WANTZ ) THEN
265: LWMIN = 18*N
266: LIWMIN = 10*N
267: ELSE
268: * need less workspace if only the eigenvalues are wanted
269: LWMIN = 12*N
270: LIWMIN = 8*N
271: ENDIF
272:
273: WL = ZERO
274: WU = ZERO
275: IIL = 0
276: IIU = 0
277:
278: IF( VALEIG ) THEN
279: * We do not reference VL, VU in the cases RANGE = 'I','A'
280: * The interval (WL, WU] contains all the wanted eigenvalues.
281: * It is either given by the user or computed in DLARRE.
282: WL = VL
283: WU = VU
284: ELSEIF( INDEIG ) THEN
285: * We do not reference IL, IU in the cases RANGE = 'V','A'
286: IIL = IL
287: IIU = IU
288: ENDIF
289: *
290: INFO = 0
291: IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
292: INFO = -1
293: ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
294: INFO = -2
295: ELSE IF( N.LT.0 ) THEN
296: INFO = -3
297: ELSE IF( VALEIG .AND. N.GT.0 .AND. WU.LE.WL ) THEN
298: INFO = -7
299: ELSE IF( INDEIG .AND. ( IIL.LT.1 .OR. IIL.GT.N ) ) THEN
300: INFO = -8
301: ELSE IF( INDEIG .AND. ( IIU.LT.IIL .OR. IIU.GT.N ) ) THEN
302: INFO = -9
303: ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
304: INFO = -13
305: ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
306: INFO = -17
307: ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
308: INFO = -19
309: END IF
310: *
311: * Get machine constants.
312: *
313: SAFMIN = DLAMCH( 'Safe minimum' )
314: EPS = DLAMCH( 'Precision' )
315: SMLNUM = SAFMIN / EPS
316: BIGNUM = ONE / SMLNUM
317: RMIN = SQRT( SMLNUM )
318: RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
319: *
320: IF( INFO.EQ.0 ) THEN
321: WORK( 1 ) = LWMIN
322: IWORK( 1 ) = LIWMIN
323: *
324: IF( WANTZ .AND. ALLEIG ) THEN
325: NZCMIN = N
326: ELSE IF( WANTZ .AND. VALEIG ) THEN
327: CALL DLARRC( 'T', N, VL, VU, D, E, SAFMIN,
328: $ NZCMIN, ITMP, ITMP2, INFO )
329: ELSE IF( WANTZ .AND. INDEIG ) THEN
330: NZCMIN = IIU-IIL+1
331: ELSE
332: * WANTZ .EQ. FALSE.
333: NZCMIN = 0
334: ENDIF
335: IF( ZQUERY .AND. INFO.EQ.0 ) THEN
336: Z( 1,1 ) = NZCMIN
337: ELSE IF( NZC.LT.NZCMIN .AND. .NOT.ZQUERY ) THEN
338: INFO = -14
339: END IF
340: END IF
341:
342: IF( INFO.NE.0 ) THEN
343: *
344: CALL XERBLA( 'DSTEMR', -INFO )
345: *
346: RETURN
347: ELSE IF( LQUERY .OR. ZQUERY ) THEN
348: RETURN
349: END IF
350: *
351: * Handle N = 0, 1, and 2 cases immediately
352: *
353: M = 0
354: IF( N.EQ.0 )
355: $ RETURN
356: *
357: IF( N.EQ.1 ) THEN
358: IF( ALLEIG .OR. INDEIG ) THEN
359: M = 1
360: W( 1 ) = D( 1 )
361: ELSE
362: IF( WL.LT.D( 1 ) .AND. WU.GE.D( 1 ) ) THEN
363: M = 1
364: W( 1 ) = D( 1 )
365: END IF
366: END IF
367: IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN
368: Z( 1, 1 ) = ONE
369: ISUPPZ(1) = 1
370: ISUPPZ(2) = 1
371: END IF
372: RETURN
373: END IF
374: *
375: IF( N.EQ.2 ) THEN
376: IF( .NOT.WANTZ ) THEN
377: CALL DLAE2( D(1), E(1), D(2), R1, R2 )
378: ELSE IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN
379: CALL DLAEV2( D(1), E(1), D(2), R1, R2, CS, SN )
380: END IF
381: IF( ALLEIG.OR.
382: $ (VALEIG.AND.(R2.GT.WL).AND.
383: $ (R2.LE.WU)).OR.
384: $ (INDEIG.AND.(IIL.EQ.1)) ) THEN
385: M = M+1
386: W( M ) = R2
387: IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN
388: Z( 1, M ) = -SN
389: Z( 2, M ) = CS
390: * Note: At most one of SN and CS can be zero.
391: IF (SN.NE.ZERO) THEN
392: IF (CS.NE.ZERO) THEN
393: ISUPPZ(2*M-1) = 1
394: ISUPPZ(2*M-1) = 2
395: ELSE
396: ISUPPZ(2*M-1) = 1
397: ISUPPZ(2*M-1) = 1
398: END IF
399: ELSE
400: ISUPPZ(2*M-1) = 2
401: ISUPPZ(2*M) = 2
402: END IF
403: ENDIF
404: ENDIF
405: IF( ALLEIG.OR.
406: $ (VALEIG.AND.(R1.GT.WL).AND.
407: $ (R1.LE.WU)).OR.
408: $ (INDEIG.AND.(IIU.EQ.2)) ) THEN
409: M = M+1
410: W( M ) = R1
411: IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN
412: Z( 1, M ) = CS
413: Z( 2, M ) = SN
414: * Note: At most one of SN and CS can be zero.
415: IF (SN.NE.ZERO) THEN
416: IF (CS.NE.ZERO) THEN
417: ISUPPZ(2*M-1) = 1
418: ISUPPZ(2*M-1) = 2
419: ELSE
420: ISUPPZ(2*M-1) = 1
421: ISUPPZ(2*M-1) = 1
422: END IF
423: ELSE
424: ISUPPZ(2*M-1) = 2
425: ISUPPZ(2*M) = 2
426: END IF
427: ENDIF
428: ENDIF
429: RETURN
430: END IF
431:
432: * Continue with general N
433:
434: INDGRS = 1
435: INDERR = 2*N + 1
436: INDGP = 3*N + 1
437: INDD = 4*N + 1
438: INDE2 = 5*N + 1
439: INDWRK = 6*N + 1
440: *
441: IINSPL = 1
442: IINDBL = N + 1
443: IINDW = 2*N + 1
444: IINDWK = 3*N + 1
445: *
446: * Scale matrix to allowable range, if necessary.
447: * The allowable range is related to the PIVMIN parameter; see the
448: * comments in DLARRD. The preference for scaling small values
449: * up is heuristic; we expect users' matrices not to be close to the
450: * RMAX threshold.
451: *
452: SCALE = ONE
453: TNRM = DLANST( 'M', N, D, E )
454: IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
455: SCALE = RMIN / TNRM
456: ELSE IF( TNRM.GT.RMAX ) THEN
457: SCALE = RMAX / TNRM
458: END IF
459: IF( SCALE.NE.ONE ) THEN
460: CALL DSCAL( N, SCALE, D, 1 )
461: CALL DSCAL( N-1, SCALE, E, 1 )
462: TNRM = TNRM*SCALE
463: IF( VALEIG ) THEN
464: * If eigenvalues in interval have to be found,
465: * scale (WL, WU] accordingly
466: WL = WL*SCALE
467: WU = WU*SCALE
468: ENDIF
469: END IF
470: *
471: * Compute the desired eigenvalues of the tridiagonal after splitting
472: * into smaller subblocks if the corresponding off-diagonal elements
473: * are small
474: * THRESH is the splitting parameter for DLARRE
475: * A negative THRESH forces the old splitting criterion based on the
476: * size of the off-diagonal. A positive THRESH switches to splitting
477: * which preserves relative accuracy.
478: *
479: IF( TRYRAC ) THEN
480: * Test whether the matrix warrants the more expensive relative approach.
481: CALL DLARRR( N, D, E, IINFO )
482: ELSE
483: * The user does not care about relative accurately eigenvalues
484: IINFO = -1
485: ENDIF
486: * Set the splitting criterion
487: IF (IINFO.EQ.0) THEN
488: THRESH = EPS
489: ELSE
490: THRESH = -EPS
491: * relative accuracy is desired but T does not guarantee it
492: TRYRAC = .FALSE.
493: ENDIF
494: *
495: IF( TRYRAC ) THEN
496: * Copy original diagonal, needed to guarantee relative accuracy
497: CALL DCOPY(N,D,1,WORK(INDD),1)
498: ENDIF
499: * Store the squares of the offdiagonal values of T
500: DO 5 J = 1, N-1
501: WORK( INDE2+J-1 ) = E(J)**2
502: 5 CONTINUE
503:
504: * Set the tolerance parameters for bisection
505: IF( .NOT.WANTZ ) THEN
506: * DLARRE computes the eigenvalues to full precision.
507: RTOL1 = FOUR * EPS
508: RTOL2 = FOUR * EPS
509: ELSE
510: * DLARRE computes the eigenvalues to less than full precision.
511: * DLARRV will refine the eigenvalue approximations, and we can
512: * need less accurate initial bisection in DLARRE.
513: * Note: these settings do only affect the subset case and DLARRE
514: RTOL1 = SQRT(EPS)
515: RTOL2 = MAX( SQRT(EPS)*5.0D-3, FOUR * EPS )
516: ENDIF
517: CALL DLARRE( RANGE, N, WL, WU, IIL, IIU, D, E,
518: $ WORK(INDE2), RTOL1, RTOL2, THRESH, NSPLIT,
519: $ IWORK( IINSPL ), M, W, WORK( INDERR ),
520: $ WORK( INDGP ), IWORK( IINDBL ),
521: $ IWORK( IINDW ), WORK( INDGRS ), PIVMIN,
522: $ WORK( INDWRK ), IWORK( IINDWK ), IINFO )
523: IF( IINFO.NE.0 ) THEN
524: INFO = 10 + ABS( IINFO )
525: RETURN
526: END IF
527: * Note that if RANGE .NE. 'V', DLARRE computes bounds on the desired
528: * part of the spectrum. All desired eigenvalues are contained in
529: * (WL,WU]
530:
531:
532: IF( WANTZ ) THEN
533: *
534: * Compute the desired eigenvectors corresponding to the computed
535: * eigenvalues
536: *
537: CALL DLARRV( N, WL, WU, D, E,
538: $ PIVMIN, IWORK( IINSPL ), M,
539: $ 1, M, MINRGP, RTOL1, RTOL2,
540: $ W, WORK( INDERR ), WORK( INDGP ), IWORK( IINDBL ),
541: $ IWORK( IINDW ), WORK( INDGRS ), Z, LDZ,
542: $ ISUPPZ, WORK( INDWRK ), IWORK( IINDWK ), IINFO )
543: IF( IINFO.NE.0 ) THEN
544: INFO = 20 + ABS( IINFO )
545: RETURN
546: END IF
547: ELSE
548: * DLARRE computes eigenvalues of the (shifted) root representation
549: * DLARRV returns the eigenvalues of the unshifted matrix.
550: * However, if the eigenvectors are not desired by the user, we need
551: * to apply the corresponding shifts from DLARRE to obtain the
552: * eigenvalues of the original matrix.
553: DO 20 J = 1, M
554: ITMP = IWORK( IINDBL+J-1 )
555: W( J ) = W( J ) + E( IWORK( IINSPL+ITMP-1 ) )
556: 20 CONTINUE
557: END IF
558: *
559:
560: IF ( TRYRAC ) THEN
561: * Refine computed eigenvalues so that they are relatively accurate
562: * with respect to the original matrix T.
563: IBEGIN = 1
564: WBEGIN = 1
565: DO 39 JBLK = 1, IWORK( IINDBL+M-1 )
566: IEND = IWORK( IINSPL+JBLK-1 )
567: IN = IEND - IBEGIN + 1
568: WEND = WBEGIN - 1
569: * check if any eigenvalues have to be refined in this block
570: 36 CONTINUE
571: IF( WEND.LT.M ) THEN
572: IF( IWORK( IINDBL+WEND ).EQ.JBLK ) THEN
573: WEND = WEND + 1
574: GO TO 36
575: END IF
576: END IF
577: IF( WEND.LT.WBEGIN ) THEN
578: IBEGIN = IEND + 1
579: GO TO 39
580: END IF
581:
582: OFFSET = IWORK(IINDW+WBEGIN-1)-1
583: IFIRST = IWORK(IINDW+WBEGIN-1)
584: ILAST = IWORK(IINDW+WEND-1)
585: RTOL2 = FOUR * EPS
586: CALL DLARRJ( IN,
587: $ WORK(INDD+IBEGIN-1), WORK(INDE2+IBEGIN-1),
588: $ IFIRST, ILAST, RTOL2, OFFSET, W(WBEGIN),
589: $ WORK( INDERR+WBEGIN-1 ),
590: $ WORK( INDWRK ), IWORK( IINDWK ), PIVMIN,
591: $ TNRM, IINFO )
592: IBEGIN = IEND + 1
593: WBEGIN = WEND + 1
594: 39 CONTINUE
595: ENDIF
596: *
597: * If matrix was scaled, then rescale eigenvalues appropriately.
598: *
599: IF( SCALE.NE.ONE ) THEN
600: CALL DSCAL( M, ONE / SCALE, W, 1 )
601: END IF
602: *
603: * If eigenvalues are not in increasing order, then sort them,
604: * possibly along with eigenvectors.
605: *
606: IF( NSPLIT.GT.1 ) THEN
607: IF( .NOT. WANTZ ) THEN
608: CALL DLASRT( 'I', M, W, IINFO )
609: IF( IINFO.NE.0 ) THEN
610: INFO = 3
611: RETURN
612: END IF
613: ELSE
614: DO 60 J = 1, M - 1
615: I = 0
616: TMP = W( J )
617: DO 50 JJ = J + 1, M
618: IF( W( JJ ).LT.TMP ) THEN
619: I = JJ
620: TMP = W( JJ )
621: END IF
622: 50 CONTINUE
623: IF( I.NE.0 ) THEN
624: W( I ) = W( J )
625: W( J ) = TMP
626: IF( WANTZ ) THEN
627: CALL DSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
628: ITMP = ISUPPZ( 2*I-1 )
629: ISUPPZ( 2*I-1 ) = ISUPPZ( 2*J-1 )
630: ISUPPZ( 2*J-1 ) = ITMP
631: ITMP = ISUPPZ( 2*I )
632: ISUPPZ( 2*I ) = ISUPPZ( 2*J )
633: ISUPPZ( 2*J ) = ITMP
634: END IF
635: END IF
636: 60 CONTINUE
637: END IF
638: ENDIF
639: *
640: *
641: WORK( 1 ) = LWMIN
642: IWORK( 1 ) = LIWMIN
643: RETURN
644: *
645: * End of DSTEMR
646: *
647: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dstemr.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack