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