OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
strevc3.f
Go to the documentation of this file.
1*> \brief \b STREVC3
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download STREVC3 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/strevc3.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/strevc3.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/strevc3.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE STREVC3( 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* REAL T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
31* $ WORK( * )
32* ..
33*
34*
35*> \par Purpose:
36* =============
37*>
38*> \verbatim
39*>
40*> STREVC3 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 SHSEQR.
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 REAL 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 REAL 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 SHSEQR).
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 REAL 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 SHSEQR).
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 REAL 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 realOTHERcomputational
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 strevc3( 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 REAL T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
250 $ work( * )
251* ..
252*
253* =====================================================================
254*
255* .. Parameters ..
256 REAL ZERO, ONE
257 parameter( zero = 0.0e+0, one = 1.0e+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 REAL 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 ISAMAX, ILAENV
273 REAL SDOT, SLAMCH
274 EXTERNAL lsame, isamax, ilaenv, sdot, slamch
275* ..
276* .. External Subroutines ..
277 EXTERNAL saxpy, scopy, sgemv, slaln2, sscal, xerbla,
278 $ slacpy, sgemm, slabad, slaset
279* ..
280* .. Intrinsic Functions ..
281 INTRINSIC abs, max, sqrt
282* ..
283* .. Local Arrays ..
284 REAL 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, 'strevc', SIDE // HOWMNY, N, -1, -1, -1 )
301 MAXWRK = N + 2*N*NB
302 WORK(1) = MAXWRK
303.EQ. LQUERY = ( LWORK-1 )
304.NOT..AND..NOT. IF( RIGHTV LEFTV ) THEN
305 INFO = -1
306.NOT..AND..NOT..AND..NOT. ELSE IF( ALLV OVER SOMEV ) THEN
307 INFO = -2
308.LT. ELSE IF( N0 ) THEN
309 INFO = -4
310.LT. ELSE IF( LDTMAX( 1, N ) ) THEN
311 INFO = -6
312.LT..OR..AND..LT. ELSE IF( LDVL1 ( LEFTV LDVLN ) ) THEN
313 INFO = -8
314.LT..OR..AND..LT. ELSE IF( LDVR1 ( RIGHTV LDVRN ) ) THEN
315 INFO = -10
316.LT..AND..NOT. ELSE IF( LWORKMAX( 1, 3*N ) 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.LT. IF( JN ) THEN
333.EQ. IF( T( J+1, J )ZERO ) THEN
334 IF( SELECT( J ) )
335 $ M = M + 1
336 ELSE
337 PAIR = .TRUE.
338.OR. IF( SELECT( J ) 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.LT. IF( MMM ) THEN
354 INFO = -11
355 END IF
356 END IF
357.NE. IF( INFO0 ) THEN
358 CALL XERBLA( 'strevc3', -INFO )
359 RETURN
360 ELSE IF( LQUERY ) THEN
361 RETURN
362 END IF
363*
364* Quick return if possible.
365*
366.EQ. IF( N0 )
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.AND..GE. IF( OVER LWORK N + 2*N*NBMIN ) THEN
373 NB = (LWORK - N) / (2*N)
374 NB = MIN( NB, NBMAX )
375 CALL SLASET( '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 = SLAMCH( 'safe minimum' )
383 OVFL = ONE / UNFL
384 CALL SLABAD( UNFL, OVFL )
385 ULP = SLAMCH( '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.GT. IF( NB2 ) THEN
418 IV = NB
419 END IF
420
421 IP = 0
422 IS = M
423 DO 140 KI = N, 1, -1
424.EQ. IF( IP-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.EQ. ELSE IF( KI1 ) THEN
430* last column, so this ki must be real eigenvalue
431 IP = 0
432.EQ. ELSE IF( T( KI, KI-1 )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.EQ. IF( IP0 ) THEN
442.NOT. IF( SELECT( KI ) )
443 $ GO TO 140
444 ELSE
445.NOT. IF( 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.NE. IF( IP0 )
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.EQ. IF( IP0 ) 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.GT. IF( JJNXT )
478 $ GO TO 60
479 J1 = J
480 J2 = J
481 JNXT = J - 1
482.GT. IF( J1 ) THEN
483.NE. IF( T( J, J-1 )ZERO ) THEN
484 J1 = J - 1
485 JNXT = J - 2
486 END IF
487 END IF
488*
489.EQ. IF( J1J2 ) THEN
490*
491* 1-by-1 diagonal block
492*
493 CALL SLALN2( .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.GT. IF( XNORMONE ) THEN
501.GT. IF( WORK( J )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.NE. IF( SCALEONE )
510 $ CALL SSCAL( KI, SCALE, WORK( 1+IV*N ), 1 )
511 WORK( J+IV*N ) = X( 1, 1 )
512*
513* Update right-hand side
514*
515 CALL SAXPY( 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 SLALN2( .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.GT. IF( XNORMONE ) THEN
531 BETA = MAX( WORK( J-1 ), WORK( J ) )
532.GT. IF( BETABIGNUM / 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.NE. IF( SCALEONE )
542 $ CALL SSCAL( 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 SAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
549 $ WORK( 1+IV*N ), 1 )
550 CALL SAXPY( 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.NOT. IF( OVER ) THEN
558* ------------------------------
559* no back-transform: copy x to VR and normalize.
560 CALL SCOPY( KI, WORK( 1 + IV*N ), 1, VR( 1, IS ), 1 )
561*
562 II = ISAMAX( KI, VR( 1, IS ), 1 )
563 REMAX = ONE / ABS( VR( II, IS ) )
564 CALL SSCAL( KI, REMAX, VR( 1, IS ), 1 )
565*
566 DO 70 K = KI + 1, N
567 VR( K, IS ) = ZERO
568 70 CONTINUE
569*
570.EQ. ELSE IF( NB1 ) THEN
571* ------------------------------
572* version 1: back-transform each vector with GEMV, Q*x.
573.GT. IF( KI1 )
574 $ CALL SGEMV( '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 = ISAMAX( N, VR( 1, KI ), 1 )
579 REMAX = ONE / ABS( VR( II, KI ) )
580 CALL SSCAL( 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.GE. IF( ABS( T( KI-1, KI ) )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.GT. IF( JJNXT )
624 $ GO TO 90
625 J1 = J
626 J2 = J
627 JNXT = J - 1
628.GT. IF( J1 ) THEN
629.NE. IF( T( J, J-1 )ZERO ) THEN
630 J1 = J - 1
631 JNXT = J - 2
632 END IF
633 END IF
634*
635.EQ. IF( J1J2 ) THEN
636*
637* 1-by-1 diagonal block
638*
639 CALL SLALN2( .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.GT. IF( XNORMONE ) THEN
647.GT. IF( WORK( J )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.NE. IF( SCALEONE ) THEN
657 CALL SSCAL( KI, SCALE, WORK( 1+(IV-1)*N ), 1 )
658 CALL SSCAL( 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 SAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
666 $ WORK( 1+(IV-1)*N ), 1 )
667 CALL SAXPY( 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 SLALN2( .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.GT. IF( XNORMONE ) THEN
683 BETA = MAX( WORK( J-1 ), WORK( J ) )
684.GT. IF( BETABIGNUM / 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.NE. IF( SCALEONE ) THEN
697 CALL SSCAL( KI, SCALE, WORK( 1+(IV-1)*N ), 1 )
698 CALL SSCAL( 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 SAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
708 $ WORK( 1+(IV-1)*N ), 1 )
709 CALL SAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
710 $ WORK( 1+(IV-1)*N ), 1 )
711 CALL SAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1,
712 $ WORK( 1+(IV )*N ), 1 )
713 CALL SAXPY( 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.NOT. IF( OVER ) THEN
721* ------------------------------
722* no back-transform: copy x to VR and normalize.
723 CALL SCOPY( KI, WORK( 1+(IV-1)*N ), 1, VR(1,IS-1), 1 )
724 CALL SCOPY( 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 SSCAL( KI, REMAX, VR( 1, IS-1 ), 1 )
733 CALL SSCAL( 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.EQ. ELSE IF( NB1 ) THEN
741* ------------------------------
742* version 1: back-transform each vector with GEMV, Q*x.
743.GT. IF( KI2 ) THEN
744 CALL SGEMV( '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 SGEMV( '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 SSCAL( N, WORK(KI-1+(IV-1)*N), VR(1,KI-1), 1)
752 CALL SSCAL( 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 SSCAL( N, REMAX, VR( 1, KI-1 ), 1 )
762 CALL SSCAL( 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.GT. IF( NB1 ) THEN
780* --------------------------------------------------------
781* Blocked version of back-transform
782* For complex case, KI2 includes both vectors (KI-1 and KI)
783.EQ. IF( IP0 ) 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.LE..OR..EQ. IF( (IV2) (KI21) ) THEN
793 CALL SGEMM( '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.EQ. IF( ISCOMPLEX(K)0 ) THEN
801* real eigenvector
802 II = ISAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
803 REMAX = ONE / ABS( WORK( II + (NB+K)*N ) )
804.EQ. ELSE IF( ISCOMPLEX(K)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 SSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
818 END DO
819 CALL SLACPY( '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.NE. IF( IP0 )
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.EQ. IF( IP1 ) 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.EQ. ELSE IF( KIN ) THEN
854* last column, so this ki must be real eigenvalue
855 IP = 0
856.EQ. ELSE IF( T( KI+1, KI )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.NOT. IF( 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.NE. IF( IP0 )
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.EQ. IF( IP0 ) 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.LT. IF( JJNXT )
900 $ GO TO 170
901 J1 = J
902 J2 = J
903 JNXT = J + 1
904.LT. IF( JN ) THEN
905.NE. IF( T( J+1, J )ZERO ) THEN
906 J2 = J + 1
907 JNXT = J + 2
908 END IF
909 END IF
910*
911.EQ. IF( J1J2 ) 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.GT. IF( WORK( J )VCRIT ) THEN
919 REC = ONE / VMAX
920 CALL SSCAL( 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 $ SDOT( 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 SLALN2( .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.NE. IF( SCALEONE )
938 $ CALL SSCAL( 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.GT. IF( BETAVCRIT ) THEN
952 REC = ONE / VMAX
953 CALL SSCAL( 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 $ SDOT( 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 $ SDOT( 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 SLALN2( .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.NE. IF( SCALEONE )
977 $ CALL SSCAL( 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.NOT. IF( OVER ) THEN
991* ------------------------------
992* no back-transform: copy x to VL and normalize.
993 CALL SCOPY( N-KI+1, WORK( KI + IV*N ), 1,
994 $ VL( KI, IS ), 1 )
995*
996 II = ISAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
997 REMAX = ONE / ABS( VL( II, IS ) )
998 CALL SSCAL( 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.EQ. ELSE IF( NB1 ) THEN
1005* ------------------------------
1006* version 1: back-transform each vector with GEMV, Q*x.
1007.LT. IF( KIN )
1008 $ CALL SGEMV( '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 = ISAMAX( N, VL( 1, KI ), 1 )
1014 REMAX = ONE / ABS( VL( II, KI ) )
1015 CALL SSCAL( 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.GE. IF( ABS( T( KI, KI+1 ) )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.LT. IF( JJNXT )
1063 $ GO TO 200
1064 J1 = J
1065 J2 = J
1066 JNXT = J + 1
1067.LT. IF( JN ) THEN
1068.NE. IF( T( J+1, J )ZERO ) THEN
1069 J2 = J + 1
1070 JNXT = J + 2
1071 END IF
1072 END IF
1073*
1074.EQ. IF( J1J2 ) 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.GT. IF( WORK( J )VCRIT ) THEN
1082 REC = ONE / VMAX
1083 CALL SSCAL( N-KI+1, REC, WORK(KI+(IV )*N), 1 )
1084 CALL SSCAL( 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 $ SDOT( 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 $ SDOT( 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 SLALN2( .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.NE. IF( SCALEONE ) THEN
1105 CALL SSCAL( N-KI+1, SCALE, WORK(KI+(IV )*N), 1)
1106 CALL SSCAL( 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.GT. IF( BETAVCRIT ) THEN
1123 REC = ONE / VMAX
1124 CALL SSCAL( N-KI+1, REC, WORK(KI+(IV )*N), 1 )
1125 CALL SSCAL( 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 $ SDOT( 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 $ SDOT( 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 $ SDOT( 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 $ SDOT( 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 SLALN2( .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.NE. IF( SCALEONE ) THEN
1157 CALL SSCAL( N-KI+1, SCALE, WORK(KI+(IV )*N), 1)
1158 CALL SSCAL( 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.NOT. IF( OVER ) THEN
1175* ------------------------------
1176* no back-transform: copy x to VL and normalize.
1177 CALL SCOPY( N-KI+1, WORK( KI + (IV )*N ), 1,
1178 $ VL( KI, IS ), 1 )
1179 CALL SCOPY( 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 SSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
1189 CALL SSCAL( 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.EQ. ELSE IF( NB1 ) THEN
1197* ------------------------------
1198* version 1: back-transform each vector with GEMV, Q*x.
1199.LT. IF( KIN-1 ) THEN
1200 CALL SGEMV( '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 SGEMV( '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 SSCAL( N, WORK(KI+ (IV )*N), VL(1, KI ), 1)
1212 CALL SSCAL( 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 SSCAL( N, REMAX, VL( 1, KI ), 1 )
1222 CALL SSCAL( 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.GT. IF( NB1 ) THEN
1241* --------------------------------------------------------
1242* Blocked version of back-transform
1243* For complex case, KI2 includes both vectors (KI and KI+1)
1244.EQ. IF( IP0 ) 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.GE..OR..EQ. IF( (IVNB-1) (KI2N) ) THEN
1254 CALL SGEMM( '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.EQ. IF( ISCOMPLEX(K)0) THEN
1262* real eigenvector
1263 II = ISAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
1264 REMAX = ONE / ABS( WORK( II + (NB+K)*N ) )
1265.EQ. ELSE IF( ISCOMPLEX(K)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 SSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
1279 END DO
1280 CALL SLACPY( '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.NE. IF( IP0 )
1291 $ IS = IS + 1
1292 260 CONTINUE
1293 END IF
1294*
1295 RETURN
1296*
1297* End of STREVC3
1298*
1299 END
slabad
subroutine slabad(small, large)
SLABAD
Definition slabad.f:74
slaset
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition slaset.f:110
slacpy
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:103
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
slaln2
subroutine slaln2(ltrans, na, nw, smin, ca, a, lda, d1, d2, b, ldb, wr, wi, x, ldx, scale, xnorm, info)
SLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form.
Definition slaln2.f:218
strevc3
subroutine strevc3(side, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, mm, m, work, lwork, info)
STREVC3
Definition strevc3.f:237
strevc
subroutine strevc(side, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, mm, m, work, info)
STREVC
Definition strevc.f:222
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
saxpy
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
Definition saxpy.f:89
sgemv
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
Definition sgemv.f:156
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:187
max
#define max(a, b)
Definition macros.h:21