1: *> \brief <b> DSBEVX_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices</b>
2: *
3: * @precisions fortran d -> s
4: *
5: * =========== DOCUMENTATION ===========
6: *
7: * Online html documentation available at
8: * http://www.netlib.org/lapack/explore-html/
9: *
10: *> \htmlonly
11: *> Download DSBEVX_2STAGE + dependencies
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsbevx_2stage.f">
13: *> [TGZ]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsbevx_2stage.f">
15: *> [ZIP]</a>
16: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsbevx_2stage.f">
17: *> [TXT]</a>
18: *> \endhtmlonly
19: *
20: * Definition:
21: * ===========
22: *
23: * SUBROUTINE DSBEVX_2STAGE( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q,
24: * LDQ, VL, VU, IL, IU, ABSTOL, M, W, Z,
25: * LDZ, WORK, LWORK, IWORK, IFAIL, INFO )
26: *
27: * IMPLICIT NONE
28: *
29: * .. Scalar Arguments ..
30: * CHARACTER JOBZ, RANGE, UPLO
31: * INTEGER IL, INFO, IU, KD, LDAB, LDQ, LDZ, M, N, LWORK
32: * DOUBLE PRECISION ABSTOL, VL, VU
33: * ..
34: * .. Array Arguments ..
35: * INTEGER IFAIL( * ), IWORK( * )
36: * DOUBLE PRECISION AB( LDAB, * ), Q( LDQ, * ), W( * ), WORK( * ),
37: * $ Z( LDZ, * )
38: * ..
39: *
40: *
41: *> \par Purpose:
42: * =============
43: *>
44: *> \verbatim
45: *>
46: *> DSBEVX_2STAGE computes selected eigenvalues and, optionally, eigenvectors
47: *> of a real symmetric band matrix A using the 2stage technique for
48: *> the reduction to tridiagonal. Eigenvalues and eigenvectors can
49: *> be selected by specifying either a range of values or a range of
50: *> indices for the desired eigenvalues.
51: *> \endverbatim
52: *
53: * Arguments:
54: * ==========
55: *
56: *> \param[in] JOBZ
57: *> \verbatim
58: *> JOBZ is CHARACTER*1
59: *> = 'N': Compute eigenvalues only;
60: *> = 'V': Compute eigenvalues and eigenvectors.
61: *> Not available in this release.
62: *> \endverbatim
63: *>
64: *> \param[in] RANGE
65: *> \verbatim
66: *> RANGE is CHARACTER*1
67: *> = 'A': all eigenvalues will be found;
68: *> = 'V': all eigenvalues in the half-open interval (VL,VU]
69: *> will be found;
70: *> = 'I': the IL-th through IU-th eigenvalues will be found.
71: *> \endverbatim
72: *>
73: *> \param[in] UPLO
74: *> \verbatim
75: *> UPLO is CHARACTER*1
76: *> = 'U': Upper triangle of A is stored;
77: *> = 'L': Lower triangle of A is stored.
78: *> \endverbatim
79: *>
80: *> \param[in] N
81: *> \verbatim
82: *> N is INTEGER
83: *> The order of the matrix A. N >= 0.
84: *> \endverbatim
85: *>
86: *> \param[in] KD
87: *> \verbatim
88: *> KD is INTEGER
89: *> The number of superdiagonals of the matrix A if UPLO = 'U',
90: *> or the number of subdiagonals if UPLO = 'L'. KD >= 0.
91: *> \endverbatim
92: *>
93: *> \param[in,out] AB
94: *> \verbatim
95: *> AB is DOUBLE PRECISION array, dimension (LDAB, N)
96: *> On entry, the upper or lower triangle of the symmetric band
97: *> matrix A, stored in the first KD+1 rows of the array. The
98: *> j-th column of A is stored in the j-th column of the array AB
99: *> as follows:
100: *> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
101: *> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
102: *>
103: *> On exit, AB is overwritten by values generated during the
104: *> reduction to tridiagonal form. If UPLO = 'U', the first
105: *> superdiagonal and the diagonal of the tridiagonal matrix T
106: *> are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
107: *> the diagonal and first subdiagonal of T are returned in the
108: *> first two rows of AB.
109: *> \endverbatim
110: *>
111: *> \param[in] LDAB
112: *> \verbatim
113: *> LDAB is INTEGER
114: *> The leading dimension of the array AB. LDAB >= KD + 1.
115: *> \endverbatim
116: *>
117: *> \param[out] Q
118: *> \verbatim
119: *> Q is DOUBLE PRECISION array, dimension (LDQ, N)
120: *> If JOBZ = 'V', the N-by-N orthogonal matrix used in the
121: *> reduction to tridiagonal form.
122: *> If JOBZ = 'N', the array Q is not referenced.
123: *> \endverbatim
124: *>
125: *> \param[in] LDQ
126: *> \verbatim
127: *> LDQ is INTEGER
128: *> The leading dimension of the array Q. If JOBZ = 'V', then
129: *> LDQ >= max(1,N).
130: *> \endverbatim
131: *>
132: *> \param[in] VL
133: *> \verbatim
134: *> VL is DOUBLE PRECISION
135: *> If RANGE='V', the lower bound of the interval to
136: *> be searched for eigenvalues. VL < VU.
137: *> Not referenced if RANGE = 'A' or 'I'.
138: *> \endverbatim
139: *>
140: *> \param[in] VU
141: *> \verbatim
142: *> VU is DOUBLE PRECISION
143: *> If RANGE='V', the upper bound of the interval to
144: *> be searched for eigenvalues. VL < VU.
145: *> Not referenced if RANGE = 'A' or 'I'.
146: *> \endverbatim
147: *>
148: *> \param[in] IL
149: *> \verbatim
150: *> IL is INTEGER
151: *> If RANGE='I', the index of the
152: *> smallest eigenvalue to be returned.
153: *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
154: *> Not referenced if RANGE = 'A' or 'V'.
155: *> \endverbatim
156: *>
157: *> \param[in] IU
158: *> \verbatim
159: *> IU is INTEGER
160: *> If RANGE='I', the index of the
161: *> largest eigenvalue to be returned.
162: *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
163: *> Not referenced if RANGE = 'A' or 'V'.
164: *> \endverbatim
165: *>
166: *> \param[in] ABSTOL
167: *> \verbatim
168: *> ABSTOL is DOUBLE PRECISION
169: *> The absolute error tolerance for the eigenvalues.
170: *> An approximate eigenvalue is accepted as converged
171: *> when it is determined to lie in an interval [a,b]
172: *> of width less than or equal to
173: *>
174: *> ABSTOL + EPS * max( |a|,|b| ) ,
175: *>
176: *> where EPS is the machine precision. If ABSTOL is less than
177: *> or equal to zero, then EPS*|T| will be used in its place,
178: *> where |T| is the 1-norm of the tridiagonal matrix obtained
179: *> by reducing AB to tridiagonal form.
180: *>
181: *> Eigenvalues will be computed most accurately when ABSTOL is
182: *> set to twice the underflow threshold 2*DLAMCH('S'), not zero.
183: *> If this routine returns with INFO>0, indicating that some
184: *> eigenvectors did not converge, try setting ABSTOL to
185: *> 2*DLAMCH('S').
186: *>
187: *> See "Computing Small Singular Values of Bidiagonal Matrices
188: *> with Guaranteed High Relative Accuracy," by Demmel and
189: *> Kahan, LAPACK Working Note #3.
190: *> \endverbatim
191: *>
192: *> \param[out] M
193: *> \verbatim
194: *> M is INTEGER
195: *> The total number of eigenvalues found. 0 <= M <= N.
196: *> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
197: *> \endverbatim
198: *>
199: *> \param[out] W
200: *> \verbatim
201: *> W is DOUBLE PRECISION array, dimension (N)
202: *> The first M elements contain the selected eigenvalues in
203: *> ascending order.
204: *> \endverbatim
205: *>
206: *> \param[out] Z
207: *> \verbatim
208: *> Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M))
209: *> If JOBZ = 'V', then if INFO = 0, the first M columns of Z
210: *> contain the orthonormal eigenvectors of the matrix A
211: *> corresponding to the selected eigenvalues, with the i-th
212: *> column of Z holding the eigenvector associated with W(i).
213: *> If an eigenvector fails to converge, then that column of Z
214: *> contains the latest approximation to the eigenvector, and the
215: *> index of the eigenvector is returned in IFAIL.
216: *> If JOBZ = 'N', then Z is not referenced.
217: *> Note: the user must ensure that at least max(1,M) columns are
218: *> supplied in the array Z; if RANGE = 'V', the exact value of M
219: *> is not known in advance and an upper bound must be used.
220: *> \endverbatim
221: *>
222: *> \param[in] LDZ
223: *> \verbatim
224: *> LDZ is INTEGER
225: *> The leading dimension of the array Z. LDZ >= 1, and if
226: *> JOBZ = 'V', LDZ >= max(1,N).
227: *> \endverbatim
228: *>
229: *> \param[out] WORK
230: *> \verbatim
231: *> WORK is DOUBLE PRECISION array, dimension (LWORK)
232: *> \endverbatim
233: *>
234: *> \param[in] LWORK
235: *> \verbatim
236: *> LWORK is INTEGER
237: *> The length of the array WORK. LWORK >= 1, when N <= 1;
238: *> otherwise
239: *> If JOBZ = 'N' and N > 1, LWORK must be queried.
240: *> LWORK = MAX(1, 7*N, dimension) where
241: *> dimension = (2KD+1)*N + KD*NTHREADS + 2*N
242: *> where KD is the size of the band.
243: *> NTHREADS is the number of threads used when
244: *> openMP compilation is enabled, otherwise =1.
245: *> If JOBZ = 'V' and N > 1, LWORK must be queried. Not yet available
246: *>
247: *> If LWORK = -1, then a workspace query is assumed; the routine
248: *> only calculates the optimal size of the WORK array, returns
249: *> this value as the first entry of the WORK array, and no error
250: *> message related to LWORK is issued by XERBLA.
251: *> \endverbatim
252: *>
253: *> \param[out] IWORK
254: *> \verbatim
255: *> IWORK is INTEGER array, dimension (5*N)
256: *> \endverbatim
257: *>
258: *> \param[out] IFAIL
259: *> \verbatim
260: *> IFAIL is INTEGER array, dimension (N)
261: *> If JOBZ = 'V', then if INFO = 0, the first M elements of
262: *> IFAIL are zero. If INFO > 0, then IFAIL contains the
263: *> indices of the eigenvectors that failed to converge.
264: *> If JOBZ = 'N', then IFAIL is not referenced.
265: *> \endverbatim
266: *>
267: *> \param[out] INFO
268: *> \verbatim
269: *> INFO is INTEGER
270: *> = 0: successful exit.
271: *> < 0: if INFO = -i, the i-th argument had an illegal value.
272: *> > 0: if INFO = i, then i eigenvectors failed to converge.
273: *> Their indices are stored in array IFAIL.
274: *> \endverbatim
275: *
276: * Authors:
277: * ========
278: *
279: *> \author Univ. of Tennessee
280: *> \author Univ. of California Berkeley
281: *> \author Univ. of Colorado Denver
282: *> \author NAG Ltd.
283: *
284: *> \date June 2016
285: *
286: *> \ingroup doubleOTHEReigen
287: *
288: *> \par Further Details:
289: * =====================
290: *>
291: *> \verbatim
292: *>
293: *> All details about the 2stage techniques are available in:
294: *>
295: *> Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
296: *> Parallel reduction to condensed forms for symmetric eigenvalue problems
297: *> using aggregated fine-grained and memory-aware kernels. In Proceedings
298: *> of 2011 International Conference for High Performance Computing,
299: *> Networking, Storage and Analysis (SC '11), New York, NY, USA,
300: *> Article 8 , 11 pages.
301: *> http://doi.acm.org/10.1145/2063384.2063394
302: *>
303: *> A. Haidar, J. Kurzak, P. Luszczek, 2013.
304: *> An improved parallel singular value algorithm and its implementation
305: *> for multicore hardware, In Proceedings of 2013 International Conference
306: *> for High Performance Computing, Networking, Storage and Analysis (SC '13).
307: *> Denver, Colorado, USA, 2013.
308: *> Article 90, 12 pages.
309: *> http://doi.acm.org/10.1145/2503210.2503292
310: *>
311: *> A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
312: *> A novel hybrid CPU-GPU generalized eigensolver for electronic structure
313: *> calculations based on fine-grained memory aware tasks.
314: *> International Journal of High Performance Computing Applications.
315: *> Volume 28 Issue 2, Pages 196-209, May 2014.
316: *> http://hpc.sagepub.com/content/28/2/196
317: *>
318: *> \endverbatim
319: *
320: * =====================================================================
321: SUBROUTINE DSBEVX_2STAGE( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q,
322: $ LDQ, VL, VU, IL, IU, ABSTOL, M, W, Z,
323: $ LDZ, WORK, LWORK, IWORK, IFAIL, INFO )
324: *
325: IMPLICIT NONE
326: *
327: * -- LAPACK driver routine (version 3.8.0) --
328: * -- LAPACK is a software package provided by Univ. of Tennessee, --
329: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
330: * June 2016
331: *
332: * .. Scalar Arguments ..
333: CHARACTER JOBZ, RANGE, UPLO
334: INTEGER IL, INFO, IU, KD, LDAB, LDQ, LDZ, M, N, LWORK
335: DOUBLE PRECISION ABSTOL, VL, VU
336: * ..
337: * .. Array Arguments ..
338: INTEGER IFAIL( * ), IWORK( * )
339: DOUBLE PRECISION AB( LDAB, * ), Q( LDQ, * ), W( * ), WORK( * ),
340: $ Z( LDZ, * )
341: * ..
342: *
343: * =====================================================================
344: *
345: * .. Parameters ..
346: DOUBLE PRECISION ZERO, ONE
347: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
348: * ..
349: * .. Local Scalars ..
350: LOGICAL ALLEIG, INDEIG, LOWER, TEST, VALEIG, WANTZ,
351: $ LQUERY
352: CHARACTER ORDER
353: INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
354: $ INDISP, INDIWO, INDWRK, ISCALE, ITMP1, J, JJ,
355: $ LLWORK, LWMIN, LHTRD, LWTRD, IB, INDHOUS,
356: $ NSPLIT
357: DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
358: $ SIGMA, SMLNUM, TMP1, VLL, VUU
359: * ..
360: * .. External Functions ..
361: LOGICAL LSAME
362: INTEGER ILAENV2STAGE
363: DOUBLE PRECISION DLAMCH, DLANSB
364: EXTERNAL LSAME, DLAMCH, DLANSB, ILAENV2STAGE
365: * ..
366: * .. External Subroutines ..
367: EXTERNAL DCOPY, DGEMV, DLACPY, DLASCL, DSCAL,
368: $ DSTEBZ, DSTEIN, DSTEQR, DSTERF, DSWAP, XERBLA,
369: $ DSYTRD_SB2ST
370: * ..
371: * .. Intrinsic Functions ..
372: INTRINSIC MAX, MIN, SQRT
373: * ..
374: * .. Executable Statements ..
375: *
376: * Test the input parameters.
377: *
378: WANTZ = LSAME( JOBZ, 'V' )
379: ALLEIG = LSAME( RANGE, 'A' )
380: VALEIG = LSAME( RANGE, 'V' )
381: INDEIG = LSAME( RANGE, 'I' )
382: LOWER = LSAME( UPLO, 'L' )
383: LQUERY = ( LWORK.EQ.-1 )
384: *
385: INFO = 0
386: IF( .NOT.( LSAME( JOBZ, 'N' ) ) ) THEN
387: INFO = -1
388: ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
389: INFO = -2
390: ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
391: INFO = -3
392: ELSE IF( N.LT.0 ) THEN
393: INFO = -4
394: ELSE IF( KD.LT.0 ) THEN
395: INFO = -5
396: ELSE IF( LDAB.LT.KD+1 ) THEN
397: INFO = -7
398: ELSE IF( WANTZ .AND. LDQ.LT.MAX( 1, N ) ) THEN
399: INFO = -9
400: ELSE
401: IF( VALEIG ) THEN
402: IF( N.GT.0 .AND. VU.LE.VL )
403: $ INFO = -11
404: ELSE IF( INDEIG ) THEN
405: IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN
406: INFO = -12
407: ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
408: INFO = -13
409: END IF
410: END IF
411: END IF
412: IF( INFO.EQ.0 ) THEN
413: IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) )
414: $ INFO = -18
415: END IF
416: *
417: IF( INFO.EQ.0 ) THEN
418: IF( N.LE.1 ) THEN
419: LWMIN = 1
420: WORK( 1 ) = LWMIN
421: ELSE
422: IB = ILAENV2STAGE( 2, 'DSYTRD_SB2ST', JOBZ,
423: $ N, KD, -1, -1 )
424: LHTRD = ILAENV2STAGE( 3, 'DSYTRD_SB2ST', JOBZ,
425: $ N, KD, IB, -1 )
426: LWTRD = ILAENV2STAGE( 4, 'DSYTRD_SB2ST', JOBZ,
427: $ N, KD, IB, -1 )
428: LWMIN = 2*N + LHTRD + LWTRD
429: WORK( 1 ) = LWMIN
430: ENDIF
431: *
432: IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY )
433: $ INFO = -20
434: END IF
435: *
436: IF( INFO.NE.0 ) THEN
437: CALL XERBLA( 'DSBEVX_2STAGE ', -INFO )
438: RETURN
439: ELSE IF( LQUERY ) THEN
440: RETURN
441: END IF
442: *
443: * Quick return if possible
444: *
445: M = 0
446: IF( N.EQ.0 )
447: $ RETURN
448: *
449: IF( N.EQ.1 ) THEN
450: M = 1
451: IF( LOWER ) THEN
452: TMP1 = AB( 1, 1 )
453: ELSE
454: TMP1 = AB( KD+1, 1 )
455: END IF
456: IF( VALEIG ) THEN
457: IF( .NOT.( VL.LT.TMP1 .AND. VU.GE.TMP1 ) )
458: $ M = 0
459: END IF
460: IF( M.EQ.1 ) THEN
461: W( 1 ) = TMP1
462: IF( WANTZ )
463: $ Z( 1, 1 ) = ONE
464: END IF
465: RETURN
466: END IF
467: *
468: * Get machine constants.
469: *
470: SAFMIN = DLAMCH( 'Safe minimum' )
471: EPS = DLAMCH( 'Precision' )
472: SMLNUM = SAFMIN / EPS
473: BIGNUM = ONE / SMLNUM
474: RMIN = SQRT( SMLNUM )
475: RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
476: *
477: * Scale matrix to allowable range, if necessary.
478: *
479: ISCALE = 0
480: ABSTLL = ABSTOL
481: IF( VALEIG ) THEN
482: VLL = VL
483: VUU = VU
484: ELSE
485: VLL = ZERO
486: VUU = ZERO
487: END IF
488: ANRM = DLANSB( 'M', UPLO, N, KD, AB, LDAB, WORK )
489: IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
490: ISCALE = 1
491: SIGMA = RMIN / ANRM
492: ELSE IF( ANRM.GT.RMAX ) THEN
493: ISCALE = 1
494: SIGMA = RMAX / ANRM
495: END IF
496: IF( ISCALE.EQ.1 ) THEN
497: IF( LOWER ) THEN
498: CALL DLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
499: ELSE
500: CALL DLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
501: END IF
502: IF( ABSTOL.GT.0 )
503: $ ABSTLL = ABSTOL*SIGMA
504: IF( VALEIG ) THEN
505: VLL = VL*SIGMA
506: VUU = VU*SIGMA
507: END IF
508: END IF
509: *
510: * Call DSYTRD_SB2ST to reduce symmetric band matrix to tridiagonal form.
511: *
512: INDD = 1
513: INDE = INDD + N
514: INDHOUS = INDE + N
515: INDWRK = INDHOUS + LHTRD
516: LLWORK = LWORK - INDWRK + 1
517: *
518: CALL DSYTRD_SB2ST( "N", JOBZ, UPLO, N, KD, AB, LDAB, WORK( INDD ),
519: $ WORK( INDE ), WORK( INDHOUS ), LHTRD,
520: $ WORK( INDWRK ), LLWORK, IINFO )
521: *
522: * If all eigenvalues are desired and ABSTOL is less than or equal
523: * to zero, then call DSTERF or SSTEQR. If this fails for some
524: * eigenvalue, then try DSTEBZ.
525: *
526: TEST = .FALSE.
527: IF (INDEIG) THEN
528: IF (IL.EQ.1 .AND. IU.EQ.N) THEN
529: TEST = .TRUE.
530: END IF
531: END IF
532: IF ((ALLEIG .OR. TEST) .AND. (ABSTOL.LE.ZERO)) THEN
533: CALL DCOPY( N, WORK( INDD ), 1, W, 1 )
534: INDEE = INDWRK + 2*N
535: IF( .NOT.WANTZ ) THEN
536: CALL DCOPY( N-1, WORK( INDE ), 1, WORK( INDEE ), 1 )
537: CALL DSTERF( N, W, WORK( INDEE ), INFO )
538: ELSE
539: CALL DLACPY( 'A', N, N, Q, LDQ, Z, LDZ )
540: CALL DCOPY( N-1, WORK( INDE ), 1, WORK( INDEE ), 1 )
541: CALL DSTEQR( JOBZ, N, W, WORK( INDEE ), Z, LDZ,
542: $ WORK( INDWRK ), INFO )
543: IF( INFO.EQ.0 ) THEN
544: DO 10 I = 1, N
545: IFAIL( I ) = 0
546: 10 CONTINUE
547: END IF
548: END IF
549: IF( INFO.EQ.0 ) THEN
550: M = N
551: GO TO 30
552: END IF
553: INFO = 0
554: END IF
555: *
556: * Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN.
557: *
558: IF( WANTZ ) THEN
559: ORDER = 'B'
560: ELSE
561: ORDER = 'E'
562: END IF
563: INDIBL = 1
564: INDISP = INDIBL + N
565: INDIWO = INDISP + N
566: CALL DSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL,
567: $ WORK( INDD ), WORK( INDE ), M, NSPLIT, W,
568: $ IWORK( INDIBL ), IWORK( INDISP ), WORK( INDWRK ),
569: $ IWORK( INDIWO ), INFO )
570: *
571: IF( WANTZ ) THEN
572: CALL DSTEIN( N, WORK( INDD ), WORK( INDE ), M, W,
573: $ IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ,
574: $ WORK( INDWRK ), IWORK( INDIWO ), IFAIL, INFO )
575: *
576: * Apply orthogonal matrix used in reduction to tridiagonal
577: * form to eigenvectors returned by DSTEIN.
578: *
579: DO 20 J = 1, M
580: CALL DCOPY( N, Z( 1, J ), 1, WORK( 1 ), 1 )
581: CALL DGEMV( 'N', N, N, ONE, Q, LDQ, WORK, 1, ZERO,
582: $ Z( 1, J ), 1 )
583: 20 CONTINUE
584: END IF
585: *
586: * If matrix was scaled, then rescale eigenvalues appropriately.
587: *
588: 30 CONTINUE
589: IF( ISCALE.EQ.1 ) THEN
590: IF( INFO.EQ.0 ) THEN
591: IMAX = M
592: ELSE
593: IMAX = INFO - 1
594: END IF
595: CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
596: END IF
597: *
598: * If eigenvalues are not in order, then sort them, along with
599: * eigenvectors.
600: *
601: IF( WANTZ ) THEN
602: DO 50 J = 1, M - 1
603: I = 0
604: TMP1 = W( J )
605: DO 40 JJ = J + 1, M
606: IF( W( JJ ).LT.TMP1 ) THEN
607: I = JJ
608: TMP1 = W( JJ )
609: END IF
610: 40 CONTINUE
611: *
612: IF( I.NE.0 ) THEN
613: ITMP1 = IWORK( INDIBL+I-1 )
614: W( I ) = W( J )
615: IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 )
616: W( J ) = TMP1
617: IWORK( INDIBL+J-1 ) = ITMP1
618: CALL DSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
619: IF( INFO.NE.0 ) THEN
620: ITMP1 = IFAIL( I )
621: IFAIL( I ) = IFAIL( J )
622: IFAIL( J ) = ITMP1
623: END IF
624: END IF
625: 50 CONTINUE
626: END IF
627: *
628: * Set WORK(1) to optimal workspace size.
629: *
630: WORK( 1 ) = LWMIN
631: *
632: RETURN
633: *
634: * End of DSBEVX_2STAGE
635: *
636: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsbevx_2stage.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack