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