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