1: *> \brief \b DTREVC3
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DTREVC3 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrevc3.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrevc3.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrevc3.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DTREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL,
22: * VR, LDVR, MM, M, WORK, LWORK, INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER HOWMNY, SIDE
26: * INTEGER INFO, LDT, LDVL, LDVR, LWORK, M, MM, N
27: * ..
28: * .. Array Arguments ..
29: * LOGICAL SELECT( * )
30: * DOUBLE PRECISION T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
31: * $ WORK( * )
32: * ..
33: *
34: *
35: *> \par Purpose:
36: * =============
37: *>
38: *> \verbatim
39: *>
40: *> DTREVC3 computes some or all of the right and/or left eigenvectors of
41: *> a real upper quasi-triangular matrix T.
42: *> Matrices of this type are produced by the Schur factorization of
43: *> a real general matrix: A = Q*T*Q**T, as computed by DHSEQR.
44: *>
45: *> The right eigenvector x and the left eigenvector y of T corresponding
46: *> to an eigenvalue w are defined by:
47: *>
48: *> T*x = w*x, (y**T)*T = w*(y**T)
49: *>
50: *> where y**T denotes the transpose of the vector y.
51: *> The eigenvalues are not input to this routine, but are read directly
52: *> from the diagonal blocks of T.
53: *>
54: *> This routine returns the matrices X and/or Y of right and left
55: *> eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
56: *> input matrix. If Q is the orthogonal factor that reduces a matrix
57: *> A to Schur form T, then Q*X and Q*Y are the matrices of right and
58: *> left eigenvectors of A.
59: *>
60: *> This uses a Level 3 BLAS version of the back transformation.
61: *> \endverbatim
62: *
63: * Arguments:
64: * ==========
65: *
66: *> \param[in] SIDE
67: *> \verbatim
68: *> SIDE is CHARACTER*1
69: *> = 'R': compute right eigenvectors only;
70: *> = 'L': compute left eigenvectors only;
71: *> = 'B': compute both right and left eigenvectors.
72: *> \endverbatim
73: *>
74: *> \param[in] HOWMNY
75: *> \verbatim
76: *> HOWMNY is CHARACTER*1
77: *> = 'A': compute all right and/or left eigenvectors;
78: *> = 'B': compute all right and/or left eigenvectors,
79: *> backtransformed by the matrices in VR and/or VL;
80: *> = 'S': compute selected right and/or left eigenvectors,
81: *> as indicated by the logical array SELECT.
82: *> \endverbatim
83: *>
84: *> \param[in,out] SELECT
85: *> \verbatim
86: *> SELECT is LOGICAL array, dimension (N)
87: *> If HOWMNY = 'S', SELECT specifies the eigenvectors to be
88: *> computed.
89: *> If w(j) is a real eigenvalue, the corresponding real
90: *> eigenvector is computed if SELECT(j) is .TRUE..
91: *> If w(j) and w(j+1) are the real and imaginary parts of a
92: *> complex eigenvalue, the corresponding complex eigenvector is
93: *> computed if either SELECT(j) or SELECT(j+1) is .TRUE., and
94: *> on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to
95: *> .FALSE..
96: *> Not referenced if HOWMNY = 'A' or 'B'.
97: *> \endverbatim
98: *>
99: *> \param[in] N
100: *> \verbatim
101: *> N is INTEGER
102: *> The order of the matrix T. N >= 0.
103: *> \endverbatim
104: *>
105: *> \param[in] T
106: *> \verbatim
107: *> T is DOUBLE PRECISION array, dimension (LDT,N)
108: *> The upper quasi-triangular matrix T in Schur canonical form.
109: *> \endverbatim
110: *>
111: *> \param[in] LDT
112: *> \verbatim
113: *> LDT is INTEGER
114: *> The leading dimension of the array T. LDT >= max(1,N).
115: *> \endverbatim
116: *>
117: *> \param[in,out] VL
118: *> \verbatim
119: *> VL is DOUBLE PRECISION array, dimension (LDVL,MM)
120: *> On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
121: *> contain an N-by-N matrix Q (usually the orthogonal matrix Q
122: *> of Schur vectors returned by DHSEQR).
123: *> On exit, if SIDE = 'L' or 'B', VL contains:
124: *> if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
125: *> if HOWMNY = 'B', the matrix Q*Y;
126: *> if HOWMNY = 'S', the left eigenvectors of T specified by
127: *> SELECT, stored consecutively in the columns
128: *> of VL, in the same order as their
129: *> eigenvalues.
130: *> A complex eigenvector corresponding to a complex eigenvalue
131: *> is stored in two consecutive columns, the first holding the
132: *> real part, and the second the imaginary part.
133: *> Not referenced if SIDE = 'R'.
134: *> \endverbatim
135: *>
136: *> \param[in] LDVL
137: *> \verbatim
138: *> LDVL is INTEGER
139: *> The leading dimension of the array VL.
140: *> LDVL >= 1, and if SIDE = 'L' or 'B', LDVL >= N.
141: *> \endverbatim
142: *>
143: *> \param[in,out] VR
144: *> \verbatim
145: *> VR is DOUBLE PRECISION array, dimension (LDVR,MM)
146: *> On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
147: *> contain an N-by-N matrix Q (usually the orthogonal matrix Q
148: *> of Schur vectors returned by DHSEQR).
149: *> On exit, if SIDE = 'R' or 'B', VR contains:
150: *> if HOWMNY = 'A', the matrix X of right eigenvectors of T;
151: *> if HOWMNY = 'B', the matrix Q*X;
152: *> if HOWMNY = 'S', the right eigenvectors of T specified by
153: *> SELECT, stored consecutively in the columns
154: *> of VR, in the same order as their
155: *> eigenvalues.
156: *> A complex eigenvector corresponding to a complex eigenvalue
157: *> is stored in two consecutive columns, the first holding the
158: *> real part and the second the imaginary part.
159: *> Not referenced if SIDE = 'L'.
160: *> \endverbatim
161: *>
162: *> \param[in] LDVR
163: *> \verbatim
164: *> LDVR is INTEGER
165: *> The leading dimension of the array VR.
166: *> LDVR >= 1, and if SIDE = 'R' or 'B', LDVR >= N.
167: *> \endverbatim
168: *>
169: *> \param[in] MM
170: *> \verbatim
171: *> MM is INTEGER
172: *> The number of columns in the arrays VL and/or VR. MM >= M.
173: *> \endverbatim
174: *>
175: *> \param[out] M
176: *> \verbatim
177: *> M is INTEGER
178: *> The number of columns in the arrays VL and/or VR actually
179: *> used to store the eigenvectors.
180: *> If HOWMNY = 'A' or 'B', M is set to N.
181: *> Each selected real eigenvector occupies one column and each
182: *> selected complex eigenvector occupies two columns.
183: *> \endverbatim
184: *>
185: *> \param[out] WORK
186: *> \verbatim
187: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
188: *> \endverbatim
189: *>
190: *> \param[in] LWORK
191: *> \verbatim
192: *> LWORK is INTEGER
193: *> The dimension of array WORK. LWORK >= max(1,3*N).
194: *> For optimum performance, LWORK >= N + 2*N*NB, where NB is
195: *> the optimal blocksize.
196: *>
197: *> If LWORK = -1, then a workspace query is assumed; the routine
198: *> only calculates the optimal size of the WORK array, returns
199: *> this value as the first entry of the WORK array, and no error
200: *> message related to LWORK is issued by XERBLA.
201: *> \endverbatim
202: *>
203: *> \param[out] INFO
204: *> \verbatim
205: *> INFO is INTEGER
206: *> = 0: successful exit
207: *> < 0: if INFO = -i, the i-th argument had an illegal value
208: *> \endverbatim
209: *
210: * Authors:
211: * ========
212: *
213: *> \author Univ. of Tennessee
214: *> \author Univ. of California Berkeley
215: *> \author Univ. of Colorado Denver
216: *> \author NAG Ltd.
217: *
218: *> \ingroup doubleOTHERcomputational
219: *
220: *> \par Further Details:
221: * =====================
222: *>
223: *> \verbatim
224: *>
225: *> The algorithm used in this program is basically backward (forward)
226: *> substitution, with scaling to make the the code robust against
227: *> possible overflow.
228: *>
229: *> Each eigenvector is normalized so that the element of largest
230: *> magnitude has magnitude 1; here the magnitude of a complex number
231: *> (x,y) is taken to be |x| + |y|.
232: *> \endverbatim
233: *>
234: * =====================================================================
235: SUBROUTINE DTREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL,
236: $ VR, LDVR, MM, M, WORK, LWORK, INFO )
237: IMPLICIT NONE
238: *
239: * -- LAPACK computational routine --
240: * -- LAPACK is a software package provided by Univ. of Tennessee, --
241: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
242: *
243: * .. Scalar Arguments ..
244: CHARACTER HOWMNY, SIDE
245: INTEGER INFO, LDT, LDVL, LDVR, LWORK, M, MM, N
246: * ..
247: * .. Array Arguments ..
248: LOGICAL SELECT( * )
249: DOUBLE PRECISION T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
250: $ WORK( * )
251: * ..
252: *
253: * =====================================================================
254: *
255: * .. Parameters ..
256: DOUBLE PRECISION ZERO, ONE
257: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
258: INTEGER NBMIN, NBMAX
259: PARAMETER ( NBMIN = 8, NBMAX = 128 )
260: * ..
261: * .. Local Scalars ..
262: LOGICAL ALLV, BOTHV, LEFTV, LQUERY, OVER, PAIR,
263: $ RIGHTV, SOMEV
264: INTEGER I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI,
265: $ IV, MAXWRK, NB, KI2
266: DOUBLE PRECISION BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE,
267: $ SMIN, SMLNUM, ULP, UNFL, VCRIT, VMAX, WI, WR,
268: $ XNORM
269: * ..
270: * .. External Functions ..
271: LOGICAL LSAME
272: INTEGER IDAMAX, ILAENV
273: DOUBLE PRECISION DDOT, DLAMCH
274: EXTERNAL LSAME, IDAMAX, ILAENV, DDOT, DLAMCH
275: * ..
276: * .. External Subroutines ..
277: EXTERNAL DAXPY, DCOPY, DGEMV, DLALN2, DSCAL, XERBLA,
278: $ DGEMM, DLASET, DLABAD, DLACPY
279: * ..
280: * .. Intrinsic Functions ..
281: INTRINSIC ABS, MAX, SQRT
282: * ..
283: * .. Local Arrays ..
284: DOUBLE PRECISION X( 2, 2 )
285: INTEGER ISCOMPLEX( NBMAX )
286: * ..
287: * .. Executable Statements ..
288: *
289: * Decode and test the input parameters
290: *
291: BOTHV = LSAME( SIDE, 'B' )
292: RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
293: LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
294: *
295: ALLV = LSAME( HOWMNY, 'A' )
296: OVER = LSAME( HOWMNY, 'B' )
297: SOMEV = LSAME( HOWMNY, 'S' )
298: *
299: INFO = 0
300: NB = ILAENV( 1, 'DTREVC', SIDE // HOWMNY, N, -1, -1, -1 )
301: MAXWRK = N + 2*N*NB
302: WORK(1) = MAXWRK
303: LQUERY = ( LWORK.EQ.-1 )
304: IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
305: INFO = -1
306: ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
307: INFO = -2
308: ELSE IF( N.LT.0 ) THEN
309: INFO = -4
310: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
311: INFO = -6
312: ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
313: INFO = -8
314: ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
315: INFO = -10
316: ELSE IF( LWORK.LT.MAX( 1, 3*N ) .AND. .NOT.LQUERY ) THEN
317: INFO = -14
318: ELSE
319: *
320: * Set M to the number of columns required to store the selected
321: * eigenvectors, standardize the array SELECT if necessary, and
322: * test MM.
323: *
324: IF( SOMEV ) THEN
325: M = 0
326: PAIR = .FALSE.
327: DO 10 J = 1, N
328: IF( PAIR ) THEN
329: PAIR = .FALSE.
330: SELECT( J ) = .FALSE.
331: ELSE
332: IF( J.LT.N ) THEN
333: IF( T( J+1, J ).EQ.ZERO ) THEN
334: IF( SELECT( J ) )
335: $ M = M + 1
336: ELSE
337: PAIR = .TRUE.
338: IF( SELECT( J ) .OR. SELECT( J+1 ) ) THEN
339: SELECT( J ) = .TRUE.
340: M = M + 2
341: END IF
342: END IF
343: ELSE
344: IF( SELECT( N ) )
345: $ M = M + 1
346: END IF
347: END IF
348: 10 CONTINUE
349: ELSE
350: M = N
351: END IF
352: *
353: IF( MM.LT.M ) THEN
354: INFO = -11
355: END IF
356: END IF
357: IF( INFO.NE.0 ) THEN
358: CALL XERBLA( 'DTREVC3', -INFO )
359: RETURN
360: ELSE IF( LQUERY ) THEN
361: RETURN
362: END IF
363: *
364: * Quick return if possible.
365: *
366: IF( N.EQ.0 )
367: $ RETURN
368: *
369: * Use blocked version of back-transformation if sufficient workspace.
370: * Zero-out the workspace to avoid potential NaN propagation.
371: *
372: IF( OVER .AND. LWORK .GE. N + 2*N*NBMIN ) THEN
373: NB = (LWORK - N) / (2*N)
374: NB = MIN( NB, NBMAX )
375: CALL DLASET( 'F', N, 1+2*NB, ZERO, ZERO, WORK, N )
376: ELSE
377: NB = 1
378: END IF
379: *
380: * Set the constants to control overflow.
381: *
382: UNFL = DLAMCH( 'Safe minimum' )
383: OVFL = ONE / UNFL
384: CALL DLABAD( UNFL, OVFL )
385: ULP = DLAMCH( 'Precision' )
386: SMLNUM = UNFL*( N / ULP )
387: BIGNUM = ( ONE-ULP ) / SMLNUM
388: *
389: * Compute 1-norm of each column of strictly upper triangular
390: * part of T to control overflow in triangular solver.
391: *
392: WORK( 1 ) = ZERO
393: DO 30 J = 2, N
394: WORK( J ) = ZERO
395: DO 20 I = 1, J - 1
396: WORK( J ) = WORK( J ) + ABS( T( I, J ) )
397: 20 CONTINUE
398: 30 CONTINUE
399: *
400: * Index IP is used to specify the real or complex eigenvalue:
401: * IP = 0, real eigenvalue,
402: * 1, first of conjugate complex pair: (wr,wi)
403: * -1, second of conjugate complex pair: (wr,wi)
404: * ISCOMPLEX array stores IP for each column in current block.
405: *
406: IF( RIGHTV ) THEN
407: *
408: * ============================================================
409: * Compute right eigenvectors.
410: *
411: * IV is index of column in current block.
412: * For complex right vector, uses IV-1 for real part and IV for complex part.
413: * Non-blocked version always uses IV=2;
414: * blocked version starts with IV=NB, goes down to 1 or 2.
415: * (Note the "0-th" column is used for 1-norms computed above.)
416: IV = 2
417: IF( NB.GT.2 ) THEN
418: IV = NB
419: END IF
420:
421: IP = 0
422: IS = M
423: DO 140 KI = N, 1, -1
424: IF( IP.EQ.-1 ) THEN
425: * previous iteration (ki+1) was second of conjugate pair,
426: * so this ki is first of conjugate pair; skip to end of loop
427: IP = 1
428: GO TO 140
429: ELSE IF( KI.EQ.1 ) THEN
430: * last column, so this ki must be real eigenvalue
431: IP = 0
432: ELSE IF( T( KI, KI-1 ).EQ.ZERO ) THEN
433: * zero on sub-diagonal, so this ki is real eigenvalue
434: IP = 0
435: ELSE
436: * non-zero on sub-diagonal, so this ki is second of conjugate pair
437: IP = -1
438: END IF
439:
440: IF( SOMEV ) THEN
441: IF( IP.EQ.0 ) THEN
442: IF( .NOT.SELECT( KI ) )
443: $ GO TO 140
444: ELSE
445: IF( .NOT.SELECT( KI-1 ) )
446: $ GO TO 140
447: END IF
448: END IF
449: *
450: * Compute the KI-th eigenvalue (WR,WI).
451: *
452: WR = T( KI, KI )
453: WI = ZERO
454: IF( IP.NE.0 )
455: $ WI = SQRT( ABS( T( KI, KI-1 ) ) )*
456: $ SQRT( ABS( T( KI-1, KI ) ) )
457: SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
458: *
459: IF( IP.EQ.0 ) THEN
460: *
461: * --------------------------------------------------------
462: * Real right eigenvector
463: *
464: WORK( KI + IV*N ) = ONE
465: *
466: * Form right-hand side.
467: *
468: DO 50 K = 1, KI - 1
469: WORK( K + IV*N ) = -T( K, KI )
470: 50 CONTINUE
471: *
472: * Solve upper quasi-triangular system:
473: * [ T(1:KI-1,1:KI-1) - WR ]*X = SCALE*WORK.
474: *
475: JNXT = KI - 1
476: DO 60 J = KI - 1, 1, -1
477: IF( J.GT.JNXT )
478: $ GO TO 60
479: J1 = J
480: J2 = J
481: JNXT = J - 1
482: IF( J.GT.1 ) THEN
483: IF( T( J, J-1 ).NE.ZERO ) THEN
484: J1 = J - 1
485: JNXT = J - 2
486: END IF
487: END IF
488: *
489: IF( J1.EQ.J2 ) THEN
490: *
491: * 1-by-1 diagonal block
492: *
493: CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
494: $ LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
495: $ ZERO, X, 2, SCALE, XNORM, IERR )
496: *
497: * Scale X(1,1) to avoid overflow when updating
498: * the right-hand side.
499: *
500: IF( XNORM.GT.ONE ) THEN
501: IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
502: X( 1, 1 ) = X( 1, 1 ) / XNORM
503: SCALE = SCALE / XNORM
504: END IF
505: END IF
506: *
507: * Scale if necessary
508: *
509: IF( SCALE.NE.ONE )
510: $ CALL DSCAL( KI, SCALE, WORK( 1+IV*N ), 1 )
511: WORK( J+IV*N ) = X( 1, 1 )
512: *
513: * Update right-hand side
514: *
515: CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
516: $ WORK( 1+IV*N ), 1 )
517: *
518: ELSE
519: *
520: * 2-by-2 diagonal block
521: *
522: CALL DLALN2( .FALSE., 2, 1, SMIN, ONE,
523: $ T( J-1, J-1 ), LDT, ONE, ONE,
524: $ WORK( J-1+IV*N ), N, WR, ZERO, X, 2,
525: $ SCALE, XNORM, IERR )
526: *
527: * Scale X(1,1) and X(2,1) to avoid overflow when
528: * updating the right-hand side.
529: *
530: IF( XNORM.GT.ONE ) THEN
531: BETA = MAX( WORK( J-1 ), WORK( J ) )
532: IF( BETA.GT.BIGNUM / XNORM ) THEN
533: X( 1, 1 ) = X( 1, 1 ) / XNORM
534: X( 2, 1 ) = X( 2, 1 ) / XNORM
535: SCALE = SCALE / XNORM
536: END IF
537: END IF
538: *
539: * Scale if necessary
540: *
541: IF( SCALE.NE.ONE )
542: $ CALL DSCAL( KI, SCALE, WORK( 1+IV*N ), 1 )
543: WORK( J-1+IV*N ) = X( 1, 1 )
544: WORK( J +IV*N ) = X( 2, 1 )
545: *
546: * Update right-hand side
547: *
548: CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
549: $ WORK( 1+IV*N ), 1 )
550: CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
551: $ WORK( 1+IV*N ), 1 )
552: END IF
553: 60 CONTINUE
554: *
555: * Copy the vector x or Q*x to VR and normalize.
556: *
557: IF( .NOT.OVER ) THEN
558: * ------------------------------
559: * no back-transform: copy x to VR and normalize.
560: CALL DCOPY( KI, WORK( 1 + IV*N ), 1, VR( 1, IS ), 1 )
561: *
562: II = IDAMAX( KI, VR( 1, IS ), 1 )
563: REMAX = ONE / ABS( VR( II, IS ) )
564: CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
565: *
566: DO 70 K = KI + 1, N
567: VR( K, IS ) = ZERO
568: 70 CONTINUE
569: *
570: ELSE IF( NB.EQ.1 ) THEN
571: * ------------------------------
572: * version 1: back-transform each vector with GEMV, Q*x.
573: IF( KI.GT.1 )
574: $ CALL DGEMV( 'N', N, KI-1, ONE, VR, LDVR,
575: $ WORK( 1 + IV*N ), 1, WORK( KI + IV*N ),
576: $ VR( 1, KI ), 1 )
577: *
578: II = IDAMAX( N, VR( 1, KI ), 1 )
579: REMAX = ONE / ABS( VR( II, KI ) )
580: CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
581: *
582: ELSE
583: * ------------------------------
584: * version 2: back-transform block of vectors with GEMM
585: * zero out below vector
586: DO K = KI + 1, N
587: WORK( K + IV*N ) = ZERO
588: END DO
589: ISCOMPLEX( IV ) = IP
590: * back-transform and normalization is done below
591: END IF
592: ELSE
593: *
594: * --------------------------------------------------------
595: * Complex right eigenvector.
596: *
597: * Initial solve
598: * [ ( T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I*WI) ]*X = 0.
599: * [ ( T(KI, KI-1) T(KI, KI) ) ]
600: *
601: IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN
602: WORK( KI-1 + (IV-1)*N ) = ONE
603: WORK( KI + (IV )*N ) = WI / T( KI-1, KI )
604: ELSE
605: WORK( KI-1 + (IV-1)*N ) = -WI / T( KI, KI-1 )
606: WORK( KI + (IV )*N ) = ONE
607: END IF
608: WORK( KI + (IV-1)*N ) = ZERO
609: WORK( KI-1 + (IV )*N ) = ZERO
610: *
611: * Form right-hand side.
612: *
613: DO 80 K = 1, KI - 2
614: WORK( K+(IV-1)*N ) = -WORK( KI-1+(IV-1)*N )*T(K,KI-1)
615: WORK( K+(IV )*N ) = -WORK( KI +(IV )*N )*T(K,KI )
616: 80 CONTINUE
617: *
618: * Solve upper quasi-triangular system:
619: * [ T(1:KI-2,1:KI-2) - (WR+i*WI) ]*X = SCALE*(WORK+i*WORK2)
620: *
621: JNXT = KI - 2
622: DO 90 J = KI - 2, 1, -1
623: IF( J.GT.JNXT )
624: $ GO TO 90
625: J1 = J
626: J2 = J
627: JNXT = J - 1
628: IF( J.GT.1 ) THEN
629: IF( T( J, J-1 ).NE.ZERO ) THEN
630: J1 = J - 1
631: JNXT = J - 2
632: END IF
633: END IF
634: *
635: IF( J1.EQ.J2 ) THEN
636: *
637: * 1-by-1 diagonal block
638: *
639: CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
640: $ LDT, ONE, ONE, WORK( J+(IV-1)*N ), N,
641: $ WR, WI, X, 2, SCALE, XNORM, IERR )
642: *
643: * Scale X(1,1) and X(1,2) to avoid overflow when
644: * updating the right-hand side.
645: *
646: IF( XNORM.GT.ONE ) THEN
647: IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
648: X( 1, 1 ) = X( 1, 1 ) / XNORM
649: X( 1, 2 ) = X( 1, 2 ) / XNORM
650: SCALE = SCALE / XNORM
651: END IF
652: END IF
653: *
654: * Scale if necessary
655: *
656: IF( SCALE.NE.ONE ) THEN
657: CALL DSCAL( KI, SCALE, WORK( 1+(IV-1)*N ), 1 )
658: CALL DSCAL( KI, SCALE, WORK( 1+(IV )*N ), 1 )
659: END IF
660: WORK( J+(IV-1)*N ) = X( 1, 1 )
661: WORK( J+(IV )*N ) = X( 1, 2 )
662: *
663: * Update the right-hand side
664: *
665: CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
666: $ WORK( 1+(IV-1)*N ), 1 )
667: CALL DAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1,
668: $ WORK( 1+(IV )*N ), 1 )
669: *
670: ELSE
671: *
672: * 2-by-2 diagonal block
673: *
674: CALL DLALN2( .FALSE., 2, 2, SMIN, ONE,
675: $ T( J-1, J-1 ), LDT, ONE, ONE,
676: $ WORK( J-1+(IV-1)*N ), N, WR, WI, X, 2,
677: $ SCALE, XNORM, IERR )
678: *
679: * Scale X to avoid overflow when updating
680: * the right-hand side.
681: *
682: IF( XNORM.GT.ONE ) THEN
683: BETA = MAX( WORK( J-1 ), WORK( J ) )
684: IF( BETA.GT.BIGNUM / XNORM ) THEN
685: REC = ONE / XNORM
686: X( 1, 1 ) = X( 1, 1 )*REC
687: X( 1, 2 ) = X( 1, 2 )*REC
688: X( 2, 1 ) = X( 2, 1 )*REC
689: X( 2, 2 ) = X( 2, 2 )*REC
690: SCALE = SCALE*REC
691: END IF
692: END IF
693: *
694: * Scale if necessary
695: *
696: IF( SCALE.NE.ONE ) THEN
697: CALL DSCAL( KI, SCALE, WORK( 1+(IV-1)*N ), 1 )
698: CALL DSCAL( KI, SCALE, WORK( 1+(IV )*N ), 1 )
699: END IF
700: WORK( J-1+(IV-1)*N ) = X( 1, 1 )
701: WORK( J +(IV-1)*N ) = X( 2, 1 )
702: WORK( J-1+(IV )*N ) = X( 1, 2 )
703: WORK( J +(IV )*N ) = X( 2, 2 )
704: *
705: * Update the right-hand side
706: *
707: CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
708: $ WORK( 1+(IV-1)*N ), 1 )
709: CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
710: $ WORK( 1+(IV-1)*N ), 1 )
711: CALL DAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1,
712: $ WORK( 1+(IV )*N ), 1 )
713: CALL DAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1,
714: $ WORK( 1+(IV )*N ), 1 )
715: END IF
716: 90 CONTINUE
717: *
718: * Copy the vector x or Q*x to VR and normalize.
719: *
720: IF( .NOT.OVER ) THEN
721: * ------------------------------
722: * no back-transform: copy x to VR and normalize.
723: CALL DCOPY( KI, WORK( 1+(IV-1)*N ), 1, VR(1,IS-1), 1 )
724: CALL DCOPY( KI, WORK( 1+(IV )*N ), 1, VR(1,IS ), 1 )
725: *
726: EMAX = ZERO
727: DO 100 K = 1, KI
728: EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+
729: $ ABS( VR( K, IS ) ) )
730: 100 CONTINUE
731: REMAX = ONE / EMAX
732: CALL DSCAL( KI, REMAX, VR( 1, IS-1 ), 1 )
733: CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
734: *
735: DO 110 K = KI + 1, N
736: VR( K, IS-1 ) = ZERO
737: VR( K, IS ) = ZERO
738: 110 CONTINUE
739: *
740: ELSE IF( NB.EQ.1 ) THEN
741: * ------------------------------
742: * version 1: back-transform each vector with GEMV, Q*x.
743: IF( KI.GT.2 ) THEN
744: CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
745: $ WORK( 1 + (IV-1)*N ), 1,
746: $ WORK( KI-1 + (IV-1)*N ), VR(1,KI-1), 1)
747: CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
748: $ WORK( 1 + (IV)*N ), 1,
749: $ WORK( KI + (IV)*N ), VR( 1, KI ), 1 )
750: ELSE
751: CALL DSCAL( N, WORK(KI-1+(IV-1)*N), VR(1,KI-1), 1)
752: CALL DSCAL( N, WORK(KI +(IV )*N), VR(1,KI ), 1)
753: END IF
754: *
755: EMAX = ZERO
756: DO 120 K = 1, N
757: EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+
758: $ ABS( VR( K, KI ) ) )
759: 120 CONTINUE
760: REMAX = ONE / EMAX
761: CALL DSCAL( N, REMAX, VR( 1, KI-1 ), 1 )
762: CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
763: *
764: ELSE
765: * ------------------------------
766: * version 2: back-transform block of vectors with GEMM
767: * zero out below vector
768: DO K = KI + 1, N
769: WORK( K + (IV-1)*N ) = ZERO
770: WORK( K + (IV )*N ) = ZERO
771: END DO
772: ISCOMPLEX( IV-1 ) = -IP
773: ISCOMPLEX( IV ) = IP
774: IV = IV - 1
775: * back-transform and normalization is done below
776: END IF
777: END IF
778:
779: IF( NB.GT.1 ) THEN
780: * --------------------------------------------------------
781: * Blocked version of back-transform
782: * For complex case, KI2 includes both vectors (KI-1 and KI)
783: IF( IP.EQ.0 ) THEN
784: KI2 = KI
785: ELSE
786: KI2 = KI - 1
787: END IF
788:
789: * Columns IV:NB of work are valid vectors.
790: * When the number of vectors stored reaches NB-1 or NB,
791: * or if this was last vector, do the GEMM
792: IF( (IV.LE.2) .OR. (KI2.EQ.1) ) THEN
793: CALL DGEMM( 'N', 'N', N, NB-IV+1, KI2+NB-IV, ONE,
794: $ VR, LDVR,
795: $ WORK( 1 + (IV)*N ), N,
796: $ ZERO,
797: $ WORK( 1 + (NB+IV)*N ), N )
798: * normalize vectors
799: DO K = IV, NB
800: IF( ISCOMPLEX(K).EQ.0 ) THEN
801: * real eigenvector
802: II = IDAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
803: REMAX = ONE / ABS( WORK( II + (NB+K)*N ) )
804: ELSE IF( ISCOMPLEX(K).EQ.1 ) THEN
805: * first eigenvector of conjugate pair
806: EMAX = ZERO
807: DO II = 1, N
808: EMAX = MAX( EMAX,
809: $ ABS( WORK( II + (NB+K )*N ) )+
810: $ ABS( WORK( II + (NB+K+1)*N ) ) )
811: END DO
812: REMAX = ONE / EMAX
813: * else if ISCOMPLEX(K).EQ.-1
814: * second eigenvector of conjugate pair
815: * reuse same REMAX as previous K
816: END IF
817: CALL DSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
818: END DO
819: CALL DLACPY( 'F', N, NB-IV+1,
820: $ WORK( 1 + (NB+IV)*N ), N,
821: $ VR( 1, KI2 ), LDVR )
822: IV = NB
823: ELSE
824: IV = IV - 1
825: END IF
826: END IF ! blocked back-transform
827: *
828: IS = IS - 1
829: IF( IP.NE.0 )
830: $ IS = IS - 1
831: 140 CONTINUE
832: END IF
833:
834: IF( LEFTV ) THEN
835: *
836: * ============================================================
837: * Compute left eigenvectors.
838: *
839: * IV is index of column in current block.
840: * For complex left vector, uses IV for real part and IV+1 for complex part.
841: * Non-blocked version always uses IV=1;
842: * blocked version starts with IV=1, goes up to NB-1 or NB.
843: * (Note the "0-th" column is used for 1-norms computed above.)
844: IV = 1
845: IP = 0
846: IS = 1
847: DO 260 KI = 1, N
848: IF( IP.EQ.1 ) THEN
849: * previous iteration (ki-1) was first of conjugate pair,
850: * so this ki is second of conjugate pair; skip to end of loop
851: IP = -1
852: GO TO 260
853: ELSE IF( KI.EQ.N ) THEN
854: * last column, so this ki must be real eigenvalue
855: IP = 0
856: ELSE IF( T( KI+1, KI ).EQ.ZERO ) THEN
857: * zero on sub-diagonal, so this ki is real eigenvalue
858: IP = 0
859: ELSE
860: * non-zero on sub-diagonal, so this ki is first of conjugate pair
861: IP = 1
862: END IF
863: *
864: IF( SOMEV ) THEN
865: IF( .NOT.SELECT( KI ) )
866: $ GO TO 260
867: END IF
868: *
869: * Compute the KI-th eigenvalue (WR,WI).
870: *
871: WR = T( KI, KI )
872: WI = ZERO
873: IF( IP.NE.0 )
874: $ WI = SQRT( ABS( T( KI, KI+1 ) ) )*
875: $ SQRT( ABS( T( KI+1, KI ) ) )
876: SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
877: *
878: IF( IP.EQ.0 ) THEN
879: *
880: * --------------------------------------------------------
881: * Real left eigenvector
882: *
883: WORK( KI + IV*N ) = ONE
884: *
885: * Form right-hand side.
886: *
887: DO 160 K = KI + 1, N
888: WORK( K + IV*N ) = -T( KI, K )
889: 160 CONTINUE
890: *
891: * Solve transposed quasi-triangular system:
892: * [ T(KI+1:N,KI+1:N) - WR ]**T * X = SCALE*WORK
893: *
894: VMAX = ONE
895: VCRIT = BIGNUM
896: *
897: JNXT = KI + 1
898: DO 170 J = KI + 1, N
899: IF( J.LT.JNXT )
900: $ GO TO 170
901: J1 = J
902: J2 = J
903: JNXT = J + 1
904: IF( J.LT.N ) THEN
905: IF( T( J+1, J ).NE.ZERO ) THEN
906: J2 = J + 1
907: JNXT = J + 2
908: END IF
909: END IF
910: *
911: IF( J1.EQ.J2 ) THEN
912: *
913: * 1-by-1 diagonal block
914: *
915: * Scale if necessary to avoid overflow when forming
916: * the right-hand side.
917: *
918: IF( WORK( J ).GT.VCRIT ) THEN
919: REC = ONE / VMAX
920: CALL DSCAL( N-KI+1, REC, WORK( KI+IV*N ), 1 )
921: VMAX = ONE
922: VCRIT = BIGNUM
923: END IF
924: *
925: WORK( J+IV*N ) = WORK( J+IV*N ) -
926: $ DDOT( J-KI-1, T( KI+1, J ), 1,
927: $ WORK( KI+1+IV*N ), 1 )
928: *
929: * Solve [ T(J,J) - WR ]**T * X = WORK
930: *
931: CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
932: $ LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
933: $ ZERO, X, 2, SCALE, XNORM, IERR )
934: *
935: * Scale if necessary
936: *
937: IF( SCALE.NE.ONE )
938: $ CALL DSCAL( N-KI+1, SCALE, WORK( KI+IV*N ), 1 )
939: WORK( J+IV*N ) = X( 1, 1 )
940: VMAX = MAX( ABS( WORK( J+IV*N ) ), VMAX )
941: VCRIT = BIGNUM / VMAX
942: *
943: ELSE
944: *
945: * 2-by-2 diagonal block
946: *
947: * Scale if necessary to avoid overflow when forming
948: * the right-hand side.
949: *
950: BETA = MAX( WORK( J ), WORK( J+1 ) )
951: IF( BETA.GT.VCRIT ) THEN
952: REC = ONE / VMAX
953: CALL DSCAL( N-KI+1, REC, WORK( KI+IV*N ), 1 )
954: VMAX = ONE
955: VCRIT = BIGNUM
956: END IF
957: *
958: WORK( J+IV*N ) = WORK( J+IV*N ) -
959: $ DDOT( J-KI-1, T( KI+1, J ), 1,
960: $ WORK( KI+1+IV*N ), 1 )
961: *
962: WORK( J+1+IV*N ) = WORK( J+1+IV*N ) -
963: $ DDOT( J-KI-1, T( KI+1, J+1 ), 1,
964: $ WORK( KI+1+IV*N ), 1 )
965: *
966: * Solve
967: * [ T(J,J)-WR T(J,J+1) ]**T * X = SCALE*( WORK1 )
968: * [ T(J+1,J) T(J+1,J+1)-WR ] ( WORK2 )
969: *
970: CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ),
971: $ LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
972: $ ZERO, X, 2, SCALE, XNORM, IERR )
973: *
974: * Scale if necessary
975: *
976: IF( SCALE.NE.ONE )
977: $ CALL DSCAL( N-KI+1, SCALE, WORK( KI+IV*N ), 1 )
978: WORK( J +IV*N ) = X( 1, 1 )
979: WORK( J+1+IV*N ) = X( 2, 1 )
980: *
981: VMAX = MAX( ABS( WORK( J +IV*N ) ),
982: $ ABS( WORK( J+1+IV*N ) ), VMAX )
983: VCRIT = BIGNUM / VMAX
984: *
985: END IF
986: 170 CONTINUE
987: *
988: * Copy the vector x or Q*x to VL and normalize.
989: *
990: IF( .NOT.OVER ) THEN
991: * ------------------------------
992: * no back-transform: copy x to VL and normalize.
993: CALL DCOPY( N-KI+1, WORK( KI + IV*N ), 1,
994: $ VL( KI, IS ), 1 )
995: *
996: II = IDAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
997: REMAX = ONE / ABS( VL( II, IS ) )
998: CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
999: *
1000: DO 180 K = 1, KI - 1
1001: VL( K, IS ) = ZERO
1002: 180 CONTINUE
1003: *
1004: ELSE IF( NB.EQ.1 ) THEN
1005: * ------------------------------
1006: * version 1: back-transform each vector with GEMV, Q*x.
1007: IF( KI.LT.N )
1008: $ CALL DGEMV( 'N', N, N-KI, ONE,
1009: $ VL( 1, KI+1 ), LDVL,
1010: $ WORK( KI+1 + IV*N ), 1,
1011: $ WORK( KI + IV*N ), VL( 1, KI ), 1 )
1012: *
1013: II = IDAMAX( N, VL( 1, KI ), 1 )
1014: REMAX = ONE / ABS( VL( II, KI ) )
1015: CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
1016: *
1017: ELSE
1018: * ------------------------------
1019: * version 2: back-transform block of vectors with GEMM
1020: * zero out above vector
1021: * could go from KI-NV+1 to KI-1
1022: DO K = 1, KI - 1
1023: WORK( K + IV*N ) = ZERO
1024: END DO
1025: ISCOMPLEX( IV ) = IP
1026: * back-transform and normalization is done below
1027: END IF
1028: ELSE
1029: *
1030: * --------------------------------------------------------
1031: * Complex left eigenvector.
1032: *
1033: * Initial solve:
1034: * [ ( T(KI,KI) T(KI,KI+1) )**T - (WR - I* WI) ]*X = 0.
1035: * [ ( T(KI+1,KI) T(KI+1,KI+1) ) ]
1036: *
1037: IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN
1038: WORK( KI + (IV )*N ) = WI / T( KI, KI+1 )
1039: WORK( KI+1 + (IV+1)*N ) = ONE
1040: ELSE
1041: WORK( KI + (IV )*N ) = ONE
1042: WORK( KI+1 + (IV+1)*N ) = -WI / T( KI+1, KI )
1043: END IF
1044: WORK( KI+1 + (IV )*N ) = ZERO
1045: WORK( KI + (IV+1)*N ) = ZERO
1046: *
1047: * Form right-hand side.
1048: *
1049: DO 190 K = KI + 2, N
1050: WORK( K+(IV )*N ) = -WORK( KI +(IV )*N )*T(KI, K)
1051: WORK( K+(IV+1)*N ) = -WORK( KI+1+(IV+1)*N )*T(KI+1,K)
1052: 190 CONTINUE
1053: *
1054: * Solve transposed quasi-triangular system:
1055: * [ T(KI+2:N,KI+2:N)**T - (WR-i*WI) ]*X = WORK1+i*WORK2
1056: *
1057: VMAX = ONE
1058: VCRIT = BIGNUM
1059: *
1060: JNXT = KI + 2
1061: DO 200 J = KI + 2, N
1062: IF( J.LT.JNXT )
1063: $ GO TO 200
1064: J1 = J
1065: J2 = J
1066: JNXT = J + 1
1067: IF( J.LT.N ) THEN
1068: IF( T( J+1, J ).NE.ZERO ) THEN
1069: J2 = J + 1
1070: JNXT = J + 2
1071: END IF
1072: END IF
1073: *
1074: IF( J1.EQ.J2 ) THEN
1075: *
1076: * 1-by-1 diagonal block
1077: *
1078: * Scale if necessary to avoid overflow when
1079: * forming the right-hand side elements.
1080: *
1081: IF( WORK( J ).GT.VCRIT ) THEN
1082: REC = ONE / VMAX
1083: CALL DSCAL( N-KI+1, REC, WORK(KI+(IV )*N), 1 )
1084: CALL DSCAL( N-KI+1, REC, WORK(KI+(IV+1)*N), 1 )
1085: VMAX = ONE
1086: VCRIT = BIGNUM
1087: END IF
1088: *
1089: WORK( J+(IV )*N ) = WORK( J+(IV)*N ) -
1090: $ DDOT( J-KI-2, T( KI+2, J ), 1,
1091: $ WORK( KI+2+(IV)*N ), 1 )
1092: WORK( J+(IV+1)*N ) = WORK( J+(IV+1)*N ) -
1093: $ DDOT( J-KI-2, T( KI+2, J ), 1,
1094: $ WORK( KI+2+(IV+1)*N ), 1 )
1095: *
1096: * Solve [ T(J,J)-(WR-i*WI) ]*(X11+i*X12)= WK+I*WK2
1097: *
1098: CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
1099: $ LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
1100: $ -WI, X, 2, SCALE, XNORM, IERR )
1101: *
1102: * Scale if necessary
1103: *
1104: IF( SCALE.NE.ONE ) THEN
1105: CALL DSCAL( N-KI+1, SCALE, WORK(KI+(IV )*N), 1)
1106: CALL DSCAL( N-KI+1, SCALE, WORK(KI+(IV+1)*N), 1)
1107: END IF
1108: WORK( J+(IV )*N ) = X( 1, 1 )
1109: WORK( J+(IV+1)*N ) = X( 1, 2 )
1110: VMAX = MAX( ABS( WORK( J+(IV )*N ) ),
1111: $ ABS( WORK( J+(IV+1)*N ) ), VMAX )
1112: VCRIT = BIGNUM / VMAX
1113: *
1114: ELSE
1115: *
1116: * 2-by-2 diagonal block
1117: *
1118: * Scale if necessary to avoid overflow when forming
1119: * the right-hand side elements.
1120: *
1121: BETA = MAX( WORK( J ), WORK( J+1 ) )
1122: IF( BETA.GT.VCRIT ) THEN
1123: REC = ONE / VMAX
1124: CALL DSCAL( N-KI+1, REC, WORK(KI+(IV )*N), 1 )
1125: CALL DSCAL( N-KI+1, REC, WORK(KI+(IV+1)*N), 1 )
1126: VMAX = ONE
1127: VCRIT = BIGNUM
1128: END IF
1129: *
1130: WORK( J +(IV )*N ) = WORK( J+(IV)*N ) -
1131: $ DDOT( J-KI-2, T( KI+2, J ), 1,
1132: $ WORK( KI+2+(IV)*N ), 1 )
1133: *
1134: WORK( J +(IV+1)*N ) = WORK( J+(IV+1)*N ) -
1135: $ DDOT( J-KI-2, T( KI+2, J ), 1,
1136: $ WORK( KI+2+(IV+1)*N ), 1 )
1137: *
1138: WORK( J+1+(IV )*N ) = WORK( J+1+(IV)*N ) -
1139: $ DDOT( J-KI-2, T( KI+2, J+1 ), 1,
1140: $ WORK( KI+2+(IV)*N ), 1 )
1141: *
1142: WORK( J+1+(IV+1)*N ) = WORK( J+1+(IV+1)*N ) -
1143: $ DDOT( J-KI-2, T( KI+2, J+1 ), 1,
1144: $ WORK( KI+2+(IV+1)*N ), 1 )
1145: *
1146: * Solve 2-by-2 complex linear equation
1147: * [ (T(j,j) T(j,j+1) )**T - (wr-i*wi)*I ]*X = SCALE*B
1148: * [ (T(j+1,j) T(j+1,j+1)) ]
1149: *
1150: CALL DLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ),
1151: $ LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
1152: $ -WI, X, 2, SCALE, XNORM, IERR )
1153: *
1154: * Scale if necessary
1155: *
1156: IF( SCALE.NE.ONE ) THEN
1157: CALL DSCAL( N-KI+1, SCALE, WORK(KI+(IV )*N), 1)
1158: CALL DSCAL( N-KI+1, SCALE, WORK(KI+(IV+1)*N), 1)
1159: END IF
1160: WORK( J +(IV )*N ) = X( 1, 1 )
1161: WORK( J +(IV+1)*N ) = X( 1, 2 )
1162: WORK( J+1+(IV )*N ) = X( 2, 1 )
1163: WORK( J+1+(IV+1)*N ) = X( 2, 2 )
1164: VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ),
1165: $ ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ),
1166: $ VMAX )
1167: VCRIT = BIGNUM / VMAX
1168: *
1169: END IF
1170: 200 CONTINUE
1171: *
1172: * Copy the vector x or Q*x to VL and normalize.
1173: *
1174: IF( .NOT.OVER ) THEN
1175: * ------------------------------
1176: * no back-transform: copy x to VL and normalize.
1177: CALL DCOPY( N-KI+1, WORK( KI + (IV )*N ), 1,
1178: $ VL( KI, IS ), 1 )
1179: CALL DCOPY( N-KI+1, WORK( KI + (IV+1)*N ), 1,
1180: $ VL( KI, IS+1 ), 1 )
1181: *
1182: EMAX = ZERO
1183: DO 220 K = KI, N
1184: EMAX = MAX( EMAX, ABS( VL( K, IS ) )+
1185: $ ABS( VL( K, IS+1 ) ) )
1186: 220 CONTINUE
1187: REMAX = ONE / EMAX
1188: CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
1189: CALL DSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 )
1190: *
1191: DO 230 K = 1, KI - 1
1192: VL( K, IS ) = ZERO
1193: VL( K, IS+1 ) = ZERO
1194: 230 CONTINUE
1195: *
1196: ELSE IF( NB.EQ.1 ) THEN
1197: * ------------------------------
1198: * version 1: back-transform each vector with GEMV, Q*x.
1199: IF( KI.LT.N-1 ) THEN
1200: CALL DGEMV( 'N', N, N-KI-1, ONE,
1201: $ VL( 1, KI+2 ), LDVL,
1202: $ WORK( KI+2 + (IV)*N ), 1,
1203: $ WORK( KI + (IV)*N ),
1204: $ VL( 1, KI ), 1 )
1205: CALL DGEMV( 'N', N, N-KI-1, ONE,
1206: $ VL( 1, KI+2 ), LDVL,
1207: $ WORK( KI+2 + (IV+1)*N ), 1,
1208: $ WORK( KI+1 + (IV+1)*N ),
1209: $ VL( 1, KI+1 ), 1 )
1210: ELSE
1211: CALL DSCAL( N, WORK(KI+ (IV )*N), VL(1, KI ), 1)
1212: CALL DSCAL( N, WORK(KI+1+(IV+1)*N), VL(1, KI+1), 1)
1213: END IF
1214: *
1215: EMAX = ZERO
1216: DO 240 K = 1, N
1217: EMAX = MAX( EMAX, ABS( VL( K, KI ) )+
1218: $ ABS( VL( K, KI+1 ) ) )
1219: 240 CONTINUE
1220: REMAX = ONE / EMAX
1221: CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
1222: CALL DSCAL( N, REMAX, VL( 1, KI+1 ), 1 )
1223: *
1224: ELSE
1225: * ------------------------------
1226: * version 2: back-transform block of vectors with GEMM
1227: * zero out above vector
1228: * could go from KI-NV+1 to KI-1
1229: DO K = 1, KI - 1
1230: WORK( K + (IV )*N ) = ZERO
1231: WORK( K + (IV+1)*N ) = ZERO
1232: END DO
1233: ISCOMPLEX( IV ) = IP
1234: ISCOMPLEX( IV+1 ) = -IP
1235: IV = IV + 1
1236: * back-transform and normalization is done below
1237: END IF
1238: END IF
1239:
1240: IF( NB.GT.1 ) THEN
1241: * --------------------------------------------------------
1242: * Blocked version of back-transform
1243: * For complex case, KI2 includes both vectors (KI and KI+1)
1244: IF( IP.EQ.0 ) THEN
1245: KI2 = KI
1246: ELSE
1247: KI2 = KI + 1
1248: END IF
1249:
1250: * Columns 1:IV of work are valid vectors.
1251: * When the number of vectors stored reaches NB-1 or NB,
1252: * or if this was last vector, do the GEMM
1253: IF( (IV.GE.NB-1) .OR. (KI2.EQ.N) ) THEN
1254: CALL DGEMM( 'N', 'N', N, IV, N-KI2+IV, ONE,
1255: $ VL( 1, KI2-IV+1 ), LDVL,
1256: $ WORK( KI2-IV+1 + (1)*N ), N,
1257: $ ZERO,
1258: $ WORK( 1 + (NB+1)*N ), N )
1259: * normalize vectors
1260: DO K = 1, IV
1261: IF( ISCOMPLEX(K).EQ.0) THEN
1262: * real eigenvector
1263: II = IDAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
1264: REMAX = ONE / ABS( WORK( II + (NB+K)*N ) )
1265: ELSE IF( ISCOMPLEX(K).EQ.1) THEN
1266: * first eigenvector of conjugate pair
1267: EMAX = ZERO
1268: DO II = 1, N
1269: EMAX = MAX( EMAX,
1270: $ ABS( WORK( II + (NB+K )*N ) )+
1271: $ ABS( WORK( II + (NB+K+1)*N ) ) )
1272: END DO
1273: REMAX = ONE / EMAX
1274: * else if ISCOMPLEX(K).EQ.-1
1275: * second eigenvector of conjugate pair
1276: * reuse same REMAX as previous K
1277: END IF
1278: CALL DSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
1279: END DO
1280: CALL DLACPY( 'F', N, IV,
1281: $ WORK( 1 + (NB+1)*N ), N,
1282: $ VL( 1, KI2-IV+1 ), LDVL )
1283: IV = 1
1284: ELSE
1285: IV = IV + 1
1286: END IF
1287: END IF ! blocked back-transform
1288: *
1289: IS = IS + 1
1290: IF( IP.NE.0 )
1291: $ IS = IS + 1
1292: 260 CONTINUE
1293: END IF
1294: *
1295: RETURN
1296: *
1297: * End of DTREVC3
1298: *
1299: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtrevc3.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack