LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
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 trevc3
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, 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 = max( 1, 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 ulp = dlamch( 'Precision' )
385 smlnum = unfl*( n / ulp )
386 bignum = ( one-ulp ) / smlnum
387*
388* Compute 1-norm of each column of strictly upper triangular
389* part of T to control overflow in triangular solver.
390*
391 work( 1 ) = zero
392 DO 30 j = 2, n
393 work( j ) = zero
394 DO 20 i = 1, j - 1
395 work( j ) = work( j ) + abs( t( i, j ) )
396 20 CONTINUE
397 30 CONTINUE
398*
399* Index IP is used to specify the real or complex eigenvalue:
400* IP = 0, real eigenvalue,
401* 1, first of conjugate complex pair: (wr,wi)
402* -1, second of conjugate complex pair: (wr,wi)
403* ISCOMPLEX array stores IP for each column in current block.
404*
405 IF( rightv ) THEN
406*
407* ============================================================
408* Compute right eigenvectors.
409*
410* IV is index of column in current block.
411* For complex right vector, uses IV-1 for real part and IV for complex part.
412* Non-blocked version always uses IV=2;
413* blocked version starts with IV=NB, goes down to 1 or 2.
414* (Note the "0-th" column is used for 1-norms computed above.)
415 iv = 2
416 IF( nb.GT.2 ) THEN
417 iv = nb
418 END IF
419
420 ip = 0
421 is = m
422 DO 140 ki = n, 1, -1
423 IF( ip.EQ.-1 ) THEN
424* previous iteration (ki+1) was second of conjugate pair,
425* so this ki is first of conjugate pair; skip to end of loop
426 ip = 1
427 GO TO 140
428 ELSE IF( ki.EQ.1 ) THEN
429* last column, so this ki must be real eigenvalue
430 ip = 0
431 ELSE IF( t( ki, ki-1 ).EQ.zero ) THEN
432* zero on sub-diagonal, so this ki is real eigenvalue
433 ip = 0
434 ELSE
435* non-zero on sub-diagonal, so this ki is second of conjugate pair
436 ip = -1
437 END IF
438
439 IF( somev ) THEN
440 IF( ip.EQ.0 ) THEN
441 IF( .NOT.SELECT( ki ) )
442 $ GO TO 140
443 ELSE
444 IF( .NOT.SELECT( ki-1 ) )
445 $ GO TO 140
446 END IF
447 END IF
448*
449* Compute the KI-th eigenvalue (WR,WI).
450*
451 wr = t( ki, ki )
452 wi = zero
453 IF( ip.NE.0 )
454 $ wi = sqrt( abs( t( ki, ki-1 ) ) )*
455 $ sqrt( abs( t( ki-1, ki ) ) )
456 smin = max( ulp*( abs( wr )+abs( wi ) ), smlnum )
457*
458 IF( ip.EQ.0 ) THEN
459*
460* --------------------------------------------------------
461* Real right eigenvector
462*
463 work( ki + iv*n ) = one
464*
465* Form right-hand side.
466*
467 DO 50 k = 1, ki - 1
468 work( k + iv*n ) = -t( k, ki )
469 50 CONTINUE
470*
471* Solve upper quasi-triangular system:
472* [ T(1:KI-1,1:KI-1) - WR ]*X = SCALE*WORK.
473*
474 jnxt = ki - 1
475 DO 60 j = ki - 1, 1, -1
476 IF( j.GT.jnxt )
477 $ GO TO 60
478 j1 = j
479 j2 = j
480 jnxt = j - 1
481 IF( j.GT.1 ) THEN
482 IF( t( j, j-1 ).NE.zero ) THEN
483 j1 = j - 1
484 jnxt = j - 2
485 END IF
486 END IF
487*
488 IF( j1.EQ.j2 ) THEN
489*
490* 1-by-1 diagonal block
491*
492 CALL dlaln2( .false., 1, 1, smin, one, t( j, j ),
493 $ ldt, one, one, work( j+iv*n ), n, wr,
494 $ zero, x, 2, scale, xnorm, ierr )
495*
496* Scale X(1,1) to avoid overflow when updating
497* the right-hand side.
498*
499 IF( xnorm.GT.one ) THEN
500 IF( work( j ).GT.bignum / xnorm ) THEN
501 x( 1, 1 ) = x( 1, 1 ) / xnorm
502 scale = scale / xnorm
503 END IF
504 END IF
505*
506* Scale if necessary
507*
508 IF( scale.NE.one )
509 $ CALL dscal( ki, scale, work( 1+iv*n ), 1 )
510 work( j+iv*n ) = x( 1, 1 )
511*
512* Update right-hand side
513*
514 CALL daxpy( j-1, -x( 1, 1 ), t( 1, j ), 1,
515 $ work( 1+iv*n ), 1 )
516*
517 ELSE
518*
519* 2-by-2 diagonal block
520*
521 CALL dlaln2( .false., 2, 1, smin, one,
522 $ t( j-1, j-1 ), ldt, one, one,
523 $ work( j-1+iv*n ), n, wr, zero, x, 2,
524 $ scale, xnorm, ierr )
525*
526* Scale X(1,1) and X(2,1) to avoid overflow when
527* updating the right-hand side.
528*
529 IF( xnorm.GT.one ) THEN
530 beta = max( work( j-1 ), work( j ) )
531 IF( beta.GT.bignum / xnorm ) THEN
532 x( 1, 1 ) = x( 1, 1 ) / xnorm
533 x( 2, 1 ) = x( 2, 1 ) / xnorm
534 scale = scale / xnorm
535 END IF
536 END IF
537*
538* Scale if necessary
539*
540 IF( scale.NE.one )
541 $ CALL dscal( ki, scale, work( 1+iv*n ), 1 )
542 work( j-1+iv*n ) = x( 1, 1 )
543 work( j +iv*n ) = x( 2, 1 )
544*
545* Update right-hand side
546*
547 CALL daxpy( j-2, -x( 1, 1 ), t( 1, j-1 ), 1,
548 $ work( 1+iv*n ), 1 )
549 CALL daxpy( j-2, -x( 2, 1 ), t( 1, j ), 1,
550 $ work( 1+iv*n ), 1 )
551 END IF
552 60 CONTINUE
553*
554* Copy the vector x or Q*x to VR and normalize.
555*
556 IF( .NOT.over ) THEN
557* ------------------------------
558* no back-transform: copy x to VR and normalize.
559 CALL dcopy( ki, work( 1 + iv*n ), 1, vr( 1, is ), 1 )
560*
561 ii = idamax( ki, vr( 1, is ), 1 )
562 remax = one / abs( vr( ii, is ) )
563 CALL dscal( ki, remax, vr( 1, is ), 1 )
564*
565 DO 70 k = ki + 1, n
566 vr( k, is ) = zero
567 70 CONTINUE
568*
569 ELSE IF( nb.EQ.1 ) THEN
570* ------------------------------
571* version 1: back-transform each vector with GEMV, Q*x.
572 IF( ki.GT.1 )
573 $ CALL dgemv( 'N', n, ki-1, one, vr, ldvr,
574 $ work( 1 + iv*n ), 1, work( ki + iv*n ),
575 $ vr( 1, ki ), 1 )
576*
577 ii = idamax( n, vr( 1, ki ), 1 )
578 remax = one / abs( vr( ii, ki ) )
579 CALL dscal( n, remax, vr( 1, ki ), 1 )
580*
581 ELSE
582* ------------------------------
583* version 2: back-transform block of vectors with GEMM
584* zero out below vector
585 DO k = ki + 1, n
586 work( k + iv*n ) = zero
587 END DO
588 iscomplex( iv ) = ip
589* back-transform and normalization is done below
590 END IF
591 ELSE
592*
593* --------------------------------------------------------
594* Complex right eigenvector.
595*
596* Initial solve
597* [ ( T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I*WI) ]*X = 0.
598* [ ( T(KI, KI-1) T(KI, KI) ) ]
599*
600 IF( abs( t( ki-1, ki ) ).GE.abs( t( ki, ki-1 ) ) ) THEN
601 work( ki-1 + (iv-1)*n ) = one
602 work( ki + (iv )*n ) = wi / t( ki-1, ki )
603 ELSE
604 work( ki-1 + (iv-1)*n ) = -wi / t( ki, ki-1 )
605 work( ki + (iv )*n ) = one
606 END IF
607 work( ki + (iv-1)*n ) = zero
608 work( ki-1 + (iv )*n ) = zero
609*
610* Form right-hand side.
611*
612 DO 80 k = 1, ki - 2
613 work( k+(iv-1)*n ) = -work( ki-1+(iv-1)*n )*t(k,ki-1)
614 work( k+(iv )*n ) = -work( ki +(iv )*n )*t(k,ki )
615 80 CONTINUE
616*
617* Solve upper quasi-triangular system:
618* [ T(1:KI-2,1:KI-2) - (WR+i*WI) ]*X = SCALE*(WORK+i*WORK2)
619*
620 jnxt = ki - 2
621 DO 90 j = ki - 2, 1, -1
622 IF( j.GT.jnxt )
623 $ GO TO 90
624 j1 = j
625 j2 = j
626 jnxt = j - 1
627 IF( j.GT.1 ) THEN
628 IF( t( j, j-1 ).NE.zero ) THEN
629 j1 = j - 1
630 jnxt = j - 2
631 END IF
632 END IF
633*
634 IF( j1.EQ.j2 ) THEN
635*
636* 1-by-1 diagonal block
637*
638 CALL dlaln2( .false., 1, 2, smin, one, t( j, j ),
639 $ ldt, one, one, work( j+(iv-1)*n ), n,
640 $ wr, wi, x, 2, scale, xnorm, ierr )
641*
642* Scale X(1,1) and X(1,2) to avoid overflow when
643* updating the right-hand side.
644*
645 IF( xnorm.GT.one ) THEN
646 IF( work( j ).GT.bignum / xnorm ) THEN
647 x( 1, 1 ) = x( 1, 1 ) / xnorm
648 x( 1, 2 ) = x( 1, 2 ) / xnorm
649 scale = scale / xnorm
650 END IF
651 END IF
652*
653* Scale if necessary
654*
655 IF( scale.NE.one ) THEN
656 CALL dscal( ki, scale, work( 1+(iv-1)*n ), 1 )
657 CALL dscal( ki, scale, work( 1+(iv )*n ), 1 )
658 END IF
659 work( j+(iv-1)*n ) = x( 1, 1 )
660 work( j+(iv )*n ) = x( 1, 2 )
661*
662* Update the right-hand side
663*
664 CALL daxpy( j-1, -x( 1, 1 ), t( 1, j ), 1,
665 $ work( 1+(iv-1)*n ), 1 )
666 CALL daxpy( j-1, -x( 1, 2 ), t( 1, j ), 1,
667 $ work( 1+(iv )*n ), 1 )
668*
669 ELSE
670*
671* 2-by-2 diagonal block
672*
673 CALL dlaln2( .false., 2, 2, smin, one,
674 $ t( j-1, j-1 ), ldt, one, one,
675 $ work( j-1+(iv-1)*n ), n, wr, wi, x, 2,
676 $ scale, xnorm, ierr )
677*
678* Scale X to avoid overflow when updating
679* the right-hand side.
680*
681 IF( xnorm.GT.one ) THEN
682 beta = max( work( j-1 ), work( j ) )
683 IF( beta.GT.bignum / xnorm ) THEN
684 rec = one / xnorm
685 x( 1, 1 ) = x( 1, 1 )*rec
686 x( 1, 2 ) = x( 1, 2 )*rec
687 x( 2, 1 ) = x( 2, 1 )*rec
688 x( 2, 2 ) = x( 2, 2 )*rec
689 scale = scale*rec
690 END IF
691 END IF
692*
693* Scale if necessary
694*
695 IF( scale.NE.one ) THEN
696 CALL dscal( ki, scale, work( 1+(iv-1)*n ), 1 )
697 CALL dscal( ki, scale, work( 1+(iv )*n ), 1 )
698 END IF
699 work( j-1+(iv-1)*n ) = x( 1, 1 )
700 work( j +(iv-1)*n ) = x( 2, 1 )
701 work( j-1+(iv )*n ) = x( 1, 2 )
702 work( j +(iv )*n ) = x( 2, 2 )
703*
704* Update the right-hand side
705*
706 CALL daxpy( j-2, -x( 1, 1 ), t( 1, j-1 ), 1,
707 $ work( 1+(iv-1)*n ), 1 )
708 CALL daxpy( j-2, -x( 2, 1 ), t( 1, j ), 1,
709 $ work( 1+(iv-1)*n ), 1 )
710 CALL daxpy( j-2, -x( 1, 2 ), t( 1, j-1 ), 1,
711 $ work( 1+(iv )*n ), 1 )
712 CALL daxpy( j-2, -x( 2, 2 ), t( 1, j ), 1,
713 $ work( 1+(iv )*n ), 1 )
714 END IF
715 90 CONTINUE
716*
717* Copy the vector x or Q*x to VR and normalize.
718*
719 IF( .NOT.over ) THEN
720* ------------------------------
721* no back-transform: copy x to VR and normalize.
722 CALL dcopy( ki, work( 1+(iv-1)*n ), 1, vr(1,is-1), 1 )
723 CALL dcopy( ki, work( 1+(iv )*n ), 1, vr(1,is ), 1 )
724*
725 emax = zero
726 DO 100 k = 1, ki
727 emax = max( emax, abs( vr( k, is-1 ) )+
728 $ abs( vr( k, is ) ) )
729 100 CONTINUE
730 remax = one / emax
731 CALL dscal( ki, remax, vr( 1, is-1 ), 1 )
732 CALL dscal( ki, remax, vr( 1, is ), 1 )
733*
734 DO 110 k = ki + 1, n
735 vr( k, is-1 ) = zero
736 vr( k, is ) = zero
737 110 CONTINUE
738*
739 ELSE IF( nb.EQ.1 ) THEN
740* ------------------------------
741* version 1: back-transform each vector with GEMV, Q*x.
742 IF( ki.GT.2 ) THEN
743 CALL dgemv( 'N', n, ki-2, one, vr, ldvr,
744 $ work( 1 + (iv-1)*n ), 1,
745 $ work( ki-1 + (iv-1)*n ), vr(1,ki-1), 1)
746 CALL dgemv( 'N', n, ki-2, one, vr, ldvr,
747 $ work( 1 + (iv)*n ), 1,
748 $ work( ki + (iv)*n ), vr( 1, ki ), 1 )
749 ELSE
750 CALL dscal( n, work(ki-1+(iv-1)*n), vr(1,ki-1), 1)
751 CALL dscal( n, work(ki +(iv )*n), vr(1,ki ), 1)
752 END IF
753*
754 emax = zero
755 DO 120 k = 1, n
756 emax = max( emax, abs( vr( k, ki-1 ) )+
757 $ abs( vr( k, ki ) ) )
758 120 CONTINUE
759 remax = one / emax
760 CALL dscal( n, remax, vr( 1, ki-1 ), 1 )
761 CALL dscal( n, remax, vr( 1, ki ), 1 )
762*
763 ELSE
764* ------------------------------
765* version 2: back-transform block of vectors with GEMM
766* zero out below vector
767 DO k = ki + 1, n
768 work( k + (iv-1)*n ) = zero
769 work( k + (iv )*n ) = zero
770 END DO
771 iscomplex( iv-1 ) = -ip
772 iscomplex( iv ) = ip
773 iv = iv - 1
774* back-transform and normalization is done below
775 END IF
776 END IF
777
778 IF( nb.GT.1 ) THEN
779* --------------------------------------------------------
780* Blocked version of back-transform
781* For complex case, KI2 includes both vectors (KI-1 and KI)
782 IF( ip.EQ.0 ) THEN
783 ki2 = ki
784 ELSE
785 ki2 = ki - 1
786 END IF
787
788* Columns IV:NB of work are valid vectors.
789* When the number of vectors stored reaches NB-1 or NB,
790* or if this was last vector, do the GEMM
791 IF( (iv.LE.2) .OR. (ki2.EQ.1) ) THEN
792 CALL dgemm( 'N', 'N', n, nb-iv+1, ki2+nb-iv, one,
793 $ vr, ldvr,
794 $ work( 1 + (iv)*n ), n,
795 $ zero,
796 $ work( 1 + (nb+iv)*n ), n )
797* normalize vectors
798 DO k = iv, nb
799 IF( iscomplex(k).EQ.0 ) THEN
800* real eigenvector
801 ii = idamax( n, work( 1 + (nb+k)*n ), 1 )
802 remax = one / abs( work( ii + (nb+k)*n ) )
803 ELSE IF( iscomplex(k).EQ.1 ) THEN
804* first eigenvector of conjugate pair
805 emax = zero
806 DO ii = 1, n
807 emax = max( emax,
808 $ abs( work( ii + (nb+k )*n ) )+
809 $ abs( work( ii + (nb+k+1)*n ) ) )
810 END DO
811 remax = one / emax
812* else if ISCOMPLEX(K).EQ.-1
813* second eigenvector of conjugate pair
814* reuse same REMAX as previous K
815 END IF
816 CALL dscal( n, remax, work( 1 + (nb+k)*n ), 1 )
817 END DO
818 CALL dlacpy( 'F', n, nb-iv+1,
819 $ work( 1 + (nb+iv)*n ), n,
820 $ vr( 1, ki2 ), ldvr )
821 iv = nb
822 ELSE
823 iv = iv - 1
824 END IF
825 END IF ! blocked back-transform
826*
827 is = is - 1
828 IF( ip.NE.0 )
829 $ is = is - 1
830 140 CONTINUE
831 END IF
832
833 IF( leftv ) THEN
834*
835* ============================================================
836* Compute left eigenvectors.
837*
838* IV is index of column in current block.
839* For complex left vector, uses IV for real part and IV+1 for complex part.
840* Non-blocked version always uses IV=1;
841* blocked version starts with IV=1, goes up to NB-1 or NB.
842* (Note the "0-th" column is used for 1-norms computed above.)
843 iv = 1
844 ip = 0
845 is = 1
846 DO 260 ki = 1, n
847 IF( ip.EQ.1 ) THEN
848* previous iteration (ki-1) was first of conjugate pair,
849* so this ki is second of conjugate pair; skip to end of loop
850 ip = -1
851 GO TO 260
852 ELSE IF( ki.EQ.n ) THEN
853* last column, so this ki must be real eigenvalue
854 ip = 0
855 ELSE IF( t( ki+1, ki ).EQ.zero ) THEN
856* zero on sub-diagonal, so this ki is real eigenvalue
857 ip = 0
858 ELSE
859* non-zero on sub-diagonal, so this ki is first of conjugate pair
860 ip = 1
861 END IF
862*
863 IF( somev ) THEN
864 IF( .NOT.SELECT( ki ) )
865 $ GO TO 260
866 END IF
867*
868* Compute the KI-th eigenvalue (WR,WI).
869*
870 wr = t( ki, ki )
871 wi = zero
872 IF( ip.NE.0 )
873 $ wi = sqrt( abs( t( ki, ki+1 ) ) )*
874 $ sqrt( abs( t( ki+1, ki ) ) )
875 smin = max( ulp*( abs( wr )+abs( wi ) ), smlnum )
876*
877 IF( ip.EQ.0 ) THEN
878*
879* --------------------------------------------------------
880* Real left eigenvector
881*
882 work( ki + iv*n ) = one
883*
884* Form right-hand side.
885*
886 DO 160 k = ki + 1, n
887 work( k + iv*n ) = -t( ki, k )
888 160 CONTINUE
889*
890* Solve transposed quasi-triangular system:
891* [ T(KI+1:N,KI+1:N) - WR ]**T * X = SCALE*WORK
892*
893 vmax = one
894 vcrit = bignum
895*
896 jnxt = ki + 1
897 DO 170 j = ki + 1, n
898 IF( j.LT.jnxt )
899 $ GO TO 170
900 j1 = j
901 j2 = j
902 jnxt = j + 1
903 IF( j.LT.n ) THEN
904 IF( t( j+1, j ).NE.zero ) THEN
905 j2 = j + 1
906 jnxt = j + 2
907 END IF
908 END IF
909*
910 IF( j1.EQ.j2 ) THEN
911*
912* 1-by-1 diagonal block
913*
914* Scale if necessary to avoid overflow when forming
915* the right-hand side.
916*
917 IF( work( j ).GT.vcrit ) THEN
918 rec = one / vmax
919 CALL dscal( n-ki+1, rec, work( ki+iv*n ), 1 )
920 vmax = one
921 vcrit = bignum
922 END IF
923*
924 work( j+iv*n ) = work( j+iv*n ) -
925 $ ddot( j-ki-1, t( ki+1, j ), 1,
926 $ work( ki+1+iv*n ), 1 )
927*
928* Solve [ T(J,J) - WR ]**T * X = WORK
929*
930 CALL dlaln2( .false., 1, 1, smin, one, t( j, j ),
931 $ ldt, one, one, work( j+iv*n ), n, wr,
932 $ zero, x, 2, scale, xnorm, ierr )
933*
934* Scale if necessary
935*
936 IF( scale.NE.one )
937 $ CALL dscal( n-ki+1, scale, work( ki+iv*n ), 1 )
938 work( j+iv*n ) = x( 1, 1 )
939 vmax = max( abs( work( j+iv*n ) ), vmax )
940 vcrit = bignum / vmax
941*
942 ELSE
943*
944* 2-by-2 diagonal block
945*
946* Scale if necessary to avoid overflow when forming
947* the right-hand side.
948*
949 beta = max( work( j ), work( j+1 ) )
950 IF( beta.GT.vcrit ) THEN
951 rec = one / vmax
952 CALL dscal( n-ki+1, rec, work( ki+iv*n ), 1 )
953 vmax = one
954 vcrit = bignum
955 END IF
956*
957 work( j+iv*n ) = work( j+iv*n ) -
958 $ ddot( j-ki-1, t( ki+1, j ), 1,
959 $ work( ki+1+iv*n ), 1 )
960*
961 work( j+1+iv*n ) = work( j+1+iv*n ) -
962 $ ddot( j-ki-1, t( ki+1, j+1 ), 1,
963 $ work( ki+1+iv*n ), 1 )
964*
965* Solve
966* [ T(J,J)-WR T(J,J+1) ]**T * X = SCALE*( WORK1 )
967* [ T(J+1,J) T(J+1,J+1)-WR ] ( WORK2 )
968*
969 CALL dlaln2( .true., 2, 1, smin, one, t( j, j ),
970 $ ldt, one, one, work( j+iv*n ), n, wr,
971 $ zero, x, 2, scale, xnorm, ierr )
972*
973* Scale if necessary
974*
975 IF( scale.NE.one )
976 $ CALL dscal( n-ki+1, scale, work( ki+iv*n ), 1 )
977 work( j +iv*n ) = x( 1, 1 )
978 work( j+1+iv*n ) = x( 2, 1 )
979*
980 vmax = max( abs( work( j +iv*n ) ),
981 $ abs( work( j+1+iv*n ) ), vmax )
982 vcrit = bignum / vmax
983*
984 END IF
985 170 CONTINUE
986*
987* Copy the vector x or Q*x to VL and normalize.
988*
989 IF( .NOT.over ) THEN
990* ------------------------------
991* no back-transform: copy x to VL and normalize.
992 CALL dcopy( n-ki+1, work( ki + iv*n ), 1,
993 $ vl( ki, is ), 1 )
994*
995 ii = idamax( n-ki+1, vl( ki, is ), 1 ) + ki - 1
996 remax = one / abs( vl( ii, is ) )
997 CALL dscal( n-ki+1, remax, vl( ki, is ), 1 )
998*
999 DO 180 k = 1, ki - 1
1000 vl( k, is ) = zero
1001 180 CONTINUE
1002*
1003 ELSE IF( nb.EQ.1 ) THEN
1004* ------------------------------
1005* version 1: back-transform each vector with GEMV, Q*x.
1006 IF( ki.LT.n )
1007 $ CALL dgemv( 'N', n, n-ki, one,
1008 $ vl( 1, ki+1 ), ldvl,
1009 $ work( ki+1 + iv*n ), 1,
1010 $ work( ki + iv*n ), vl( 1, ki ), 1 )
1011*
1012 ii = idamax( n, vl( 1, ki ), 1 )
1013 remax = one / abs( vl( ii, ki ) )
1014 CALL dscal( n, remax, vl( 1, ki ), 1 )
1015*
1016 ELSE
1017* ------------------------------
1018* version 2: back-transform block of vectors with GEMM
1019* zero out above vector
1020* could go from KI-NV+1 to KI-1
1021 DO k = 1, ki - 1
1022 work( k + iv*n ) = zero
1023 END DO
1024 iscomplex( iv ) = ip
1025* back-transform and normalization is done below
1026 END IF
1027 ELSE
1028*
1029* --------------------------------------------------------
1030* Complex left eigenvector.
1031*
1032* Initial solve:
1033* [ ( T(KI,KI) T(KI,KI+1) )**T - (WR - I* WI) ]*X = 0.
1034* [ ( T(KI+1,KI) T(KI+1,KI+1) ) ]
1035*
1036 IF( abs( t( ki, ki+1 ) ).GE.abs( t( ki+1, ki ) ) ) THEN
1037 work( ki + (iv )*n ) = wi / t( ki, ki+1 )
1038 work( ki+1 + (iv+1)*n ) = one
1039 ELSE
1040 work( ki + (iv )*n ) = one
1041 work( ki+1 + (iv+1)*n ) = -wi / t( ki+1, ki )
1042 END IF
1043 work( ki+1 + (iv )*n ) = zero
1044 work( ki + (iv+1)*n ) = zero
1045*
1046* Form right-hand side.
1047*
1048 DO 190 k = ki + 2, n
1049 work( k+(iv )*n ) = -work( ki +(iv )*n )*t(ki, k)
1050 work( k+(iv+1)*n ) = -work( ki+1+(iv+1)*n )*t(ki+1,k)
1051 190 CONTINUE
1052*
1053* Solve transposed quasi-triangular system:
1054* [ T(KI+2:N,KI+2:N)**T - (WR-i*WI) ]*X = WORK1+i*WORK2
1055*
1056 vmax = one
1057 vcrit = bignum
1058*
1059 jnxt = ki + 2
1060 DO 200 j = ki + 2, n
1061 IF( j.LT.jnxt )
1062 $ GO TO 200
1063 j1 = j
1064 j2 = j
1065 jnxt = j + 1
1066 IF( j.LT.n ) THEN
1067 IF( t( j+1, j ).NE.zero ) THEN
1068 j2 = j + 1
1069 jnxt = j + 2
1070 END IF
1071 END IF
1072*
1073 IF( j1.EQ.j2 ) THEN
1074*
1075* 1-by-1 diagonal block
1076*
1077* Scale if necessary to avoid overflow when
1078* forming the right-hand side elements.
1079*
1080 IF( work( j ).GT.vcrit ) THEN
1081 rec = one / vmax
1082 CALL dscal( n-ki+1, rec, work(ki+(iv )*n), 1 )
1083 CALL dscal( n-ki+1, rec, work(ki+(iv+1)*n), 1 )
1084 vmax = one
1085 vcrit = bignum
1086 END IF
1087*
1088 work( j+(iv )*n ) = work( j+(iv)*n ) -
1089 $ ddot( j-ki-2, t( ki+2, j ), 1,
1090 $ work( ki+2+(iv)*n ), 1 )
1091 work( j+(iv+1)*n ) = work( j+(iv+1)*n ) -
1092 $ ddot( j-ki-2, t( ki+2, j ), 1,
1093 $ work( ki+2+(iv+1)*n ), 1 )
1094*
1095* Solve [ T(J,J)-(WR-i*WI) ]*(X11+i*X12)= WK+I*WK2
1096*
1097 CALL dlaln2( .false., 1, 2, smin, one, t( j, j ),
1098 $ ldt, one, one, work( j+iv*n ), n, wr,
1099 $ -wi, x, 2, scale, xnorm, ierr )
1100*
1101* Scale if necessary
1102*
1103 IF( scale.NE.one ) THEN
1104 CALL dscal( n-ki+1, scale, work(ki+(iv )*n), 1)
1105 CALL dscal( n-ki+1, scale, work(ki+(iv+1)*n), 1)
1106 END IF
1107 work( j+(iv )*n ) = x( 1, 1 )
1108 work( j+(iv+1)*n ) = x( 1, 2 )
1109 vmax = max( abs( work( j+(iv )*n ) ),
1110 $ abs( work( j+(iv+1)*n ) ), vmax )
1111 vcrit = bignum / vmax
1112*
1113 ELSE
1114*
1115* 2-by-2 diagonal block
1116*
1117* Scale if necessary to avoid overflow when forming
1118* the right-hand side elements.
1119*
1120 beta = max( work( j ), work( j+1 ) )
1121 IF( beta.GT.vcrit ) THEN
1122 rec = one / vmax
1123 CALL dscal( n-ki+1, rec, work(ki+(iv )*n), 1 )
1124 CALL dscal( n-ki+1, rec, work(ki+(iv+1)*n), 1 )
1125 vmax = one
1126 vcrit = bignum
1127 END IF
1128*
1129 work( j +(iv )*n ) = work( j+(iv)*n ) -
1130 $ ddot( j-ki-2, t( ki+2, j ), 1,
1131 $ work( ki+2+(iv)*n ), 1 )
1132*
1133 work( j +(iv+1)*n ) = work( j+(iv+1)*n ) -
1134 $ ddot( j-ki-2, t( ki+2, j ), 1,
1135 $ work( ki+2+(iv+1)*n ), 1 )
1136*
1137 work( j+1+(iv )*n ) = work( j+1+(iv)*n ) -
1138 $ ddot( j-ki-2, t( ki+2, j+1 ), 1,
1139 $ work( ki+2+(iv)*n ), 1 )
1140*
1141 work( j+1+(iv+1)*n ) = work( j+1+(iv+1)*n ) -
1142 $ ddot( j-ki-2, t( ki+2, j+1 ), 1,
1143 $ work( ki+2+(iv+1)*n ), 1 )
1144*
1145* Solve 2-by-2 complex linear equation
1146* [ (T(j,j) T(j,j+1) )**T - (wr-i*wi)*I ]*X = SCALE*B
1147* [ (T(j+1,j) T(j+1,j+1)) ]
1148*
1149 CALL dlaln2( .true., 2, 2, smin, one, t( j, j ),
1150 $ ldt, one, one, work( j+iv*n ), n, wr,
1151 $ -wi, x, 2, scale, xnorm, ierr )
1152*
1153* Scale if necessary
1154*
1155 IF( scale.NE.one ) THEN
1156 CALL dscal( n-ki+1, scale, work(ki+(iv )*n), 1)
1157 CALL dscal( n-ki+1, scale, work(ki+(iv+1)*n), 1)
1158 END IF
1159 work( j +(iv )*n ) = x( 1, 1 )
1160 work( j +(iv+1)*n ) = x( 1, 2 )
1161 work( j+1+(iv )*n ) = x( 2, 1 )
1162 work( j+1+(iv+1)*n ) = x( 2, 2 )
1163 vmax = max( abs( x( 1, 1 ) ), abs( x( 1, 2 ) ),
1164 $ abs( x( 2, 1 ) ), abs( x( 2, 2 ) ),
1165 $ vmax )
1166 vcrit = bignum / vmax
1167*
1168 END IF
1169 200 CONTINUE
1170*
1171* Copy the vector x or Q*x to VL and normalize.
1172*
1173 IF( .NOT.over ) THEN
1174* ------------------------------
1175* no back-transform: copy x to VL and normalize.
1176 CALL dcopy( n-ki+1, work( ki + (iv )*n ), 1,
1177 $ vl( ki, is ), 1 )
1178 CALL dcopy( n-ki+1, work( ki + (iv+1)*n ), 1,
1179 $ vl( ki, is+1 ), 1 )
1180*
1181 emax = zero
1182 DO 220 k = ki, n
1183 emax = max( emax, abs( vl( k, is ) )+
1184 $ abs( vl( k, is+1 ) ) )
1185 220 CONTINUE
1186 remax = one / emax
1187 CALL dscal( n-ki+1, remax, vl( ki, is ), 1 )
1188 CALL dscal( n-ki+1, remax, vl( ki, is+1 ), 1 )
1189*
1190 DO 230 k = 1, ki - 1
1191 vl( k, is ) = zero
1192 vl( k, is+1 ) = zero
1193 230 CONTINUE
1194*
1195 ELSE IF( nb.EQ.1 ) THEN
1196* ------------------------------
1197* version 1: back-transform each vector with GEMV, Q*x.
1198 IF( ki.LT.n-1 ) THEN
1199 CALL dgemv( 'N', n, n-ki-1, one,
1200 $ vl( 1, ki+2 ), ldvl,
1201 $ work( ki+2 + (iv)*n ), 1,
1202 $ work( ki + (iv)*n ),
1203 $ vl( 1, ki ), 1 )
1204 CALL dgemv( 'N', n, n-ki-1, one,
1205 $ vl( 1, ki+2 ), ldvl,
1206 $ work( ki+2 + (iv+1)*n ), 1,
1207 $ work( ki+1 + (iv+1)*n ),
1208 $ vl( 1, ki+1 ), 1 )
1209 ELSE
1210 CALL dscal( n, work(ki+ (iv )*n), vl(1, ki ), 1)
1211 CALL dscal( n, work(ki+1+(iv+1)*n), vl(1, ki+1), 1)
1212 END IF
1213*
1214 emax = zero
1215 DO 240 k = 1, n
1216 emax = max( emax, abs( vl( k, ki ) )+
1217 $ abs( vl( k, ki+1 ) ) )
1218 240 CONTINUE
1219 remax = one / emax
1220 CALL dscal( n, remax, vl( 1, ki ), 1 )
1221 CALL dscal( n, remax, vl( 1, ki+1 ), 1 )
1222*
1223 ELSE
1224* ------------------------------
1225* version 2: back-transform block of vectors with GEMM
1226* zero out above vector
1227* could go from KI-NV+1 to KI-1
1228 DO k = 1, ki - 1
1229 work( k + (iv )*n ) = zero
1230 work( k + (iv+1)*n ) = zero
1231 END DO
1232 iscomplex( iv ) = ip
1233 iscomplex( iv+1 ) = -ip
1234 iv = iv + 1
1235* back-transform and normalization is done below
1236 END IF
1237 END IF
1238
1239 IF( nb.GT.1 ) THEN
1240* --------------------------------------------------------
1241* Blocked version of back-transform
1242* For complex case, KI2 includes both vectors (KI and KI+1)
1243 IF( ip.EQ.0 ) THEN
1244 ki2 = ki
1245 ELSE
1246 ki2 = ki + 1
1247 END IF
1248
1249* Columns 1:IV of work are valid vectors.
1250* When the number of vectors stored reaches NB-1 or NB,
1251* or if this was last vector, do the GEMM
1252 IF( (iv.GE.nb-1) .OR. (ki2.EQ.n) ) THEN
1253 CALL dgemm( 'N', 'N', n, iv, n-ki2+iv, one,
1254 $ vl( 1, ki2-iv+1 ), ldvl,
1255 $ work( ki2-iv+1 + (1)*n ), n,
1256 $ zero,
1257 $ work( 1 + (nb+1)*n ), n )
1258* normalize vectors
1259 DO k = 1, iv
1260 IF( iscomplex(k).EQ.0) THEN
1261* real eigenvector
1262 ii = idamax( n, work( 1 + (nb+k)*n ), 1 )
1263 remax = one / abs( work( ii + (nb+k)*n ) )
1264 ELSE IF( iscomplex(k).EQ.1) THEN
1265* first eigenvector of conjugate pair
1266 emax = zero
1267 DO ii = 1, n
1268 emax = max( emax,
1269 $ abs( work( ii + (nb+k )*n ) )+
1270 $ abs( work( ii + (nb+k+1)*n ) ) )
1271 END DO
1272 remax = one / emax
1273* else if ISCOMPLEX(K).EQ.-1
1274* second eigenvector of conjugate pair
1275* reuse same REMAX as previous K
1276 END IF
1277 CALL dscal( n, remax, work( 1 + (nb+k)*n ), 1 )
1278 END DO
1279 CALL dlacpy( 'F', n, iv,
1280 $ work( 1 + (nb+1)*n ), n,
1281 $ vl( 1, ki2-iv+1 ), ldvl )
1282 iv = 1
1283 ELSE
1284 iv = iv + 1
1285 END IF
1286 END IF ! blocked back-transform
1287*
1288 is = is + 1
1289 IF( ip.NE.0 )
1290 $ is = is + 1
1291 260 CONTINUE
1292 END IF
1293*
1294 RETURN
1295*
1296* End of DTREVC3
1297*
1298 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
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
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
dgemv
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
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
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
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
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
dtrevc3
subroutine dtrevc3(side, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, mm, m, work, lwork, info)
DTREVC3
Definition dtrevc3.f:237