OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
dtrevc3.f
Go to the documentation of this file.
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
dlabad
subroutine dlabad(small, large)
DLABAD
Definition dlabad.f:74
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
dlaset
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
dlaln2
subroutine dlaln2(ltrans, na, nw, smin, ca, a, lda, d1, d2, b, ldb, wr, wi, x, ldx, scale, xnorm, info)
DLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form.
Definition dlaln2.f:218
dtrevc3
subroutine dtrevc3(side, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, mm, m, work, lwork, info)
DTREVC3
Definition dtrevc3.f:237
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dgemv
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:156
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:187
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21