3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download STREVC + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/strevc.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/strevc.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/strevc.f">
21 * SUBROUTINE STREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
22 * LDVR, MM, M, WORK, INFO )
24 * .. Scalar Arguments ..
25 * CHARACTER HOWMNY, SIDE
26 * INTEGER INFO, LDT, LDVL, LDVR, M, MM, N
28 * .. Array Arguments ..
30 * REAL T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
40 *> STREVC 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.
45 *> The right eigenvector x and the left eigenvector y of T corresponding
46 *> to an eigenvalue w are defined by:
48 *> T*x = w*x, (y**T)*T = w*(y**T)
50 *> where y**T denotes the transpose of y.
51 *> The eigenvalues are not input to this routine, but are read directly
52 *> from the diagonal blocks of T.
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.
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.
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.
82 *> \param[in,out] SELECT
84 *> SELECT is LOGICAL array, dimension (N)
85 *> If HOWMNY = 'S', SELECT specifies the eigenvectors to be
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
94 *> Not referenced if HOWMNY = 'A' or 'B'.
100 *> The order of the matrix T. N >= 0.
105 *> T is REAL array, dimension (LDT,N)
106 *> The upper quasi-triangular matrix T in Schur canonical form.
112 *> The leading dimension of the array T. LDT >= max(1,N).
117 *> VL is REAL 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 SHSEQR).
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
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'.
137 *> The leading dimension of the array VL. LDVL >= 1, and if
138 *> SIDE = 'L' or 'B', LDVL >= N.
143 *> VR is REAL 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 SHSEQR).
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
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'.
163 *> The leading dimension of the array VR. LDVR >= 1, and if
164 *> SIDE = 'R' or 'B', LDVR >= N.
170 *> The number of columns in the arrays VL and/or VR. MM >= M.
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.
185 *> WORK is REAL array, dimension (3*N)
191 *> = 0: successful exit
192 *> < 0: if INFO = -i, the i-th argument had an illegal value
198 *> \author Univ. of Tennessee
199 *> \author Univ. of California Berkeley
200 *> \author Univ. of Colorado Denver
203 *> \date November 2011
205 *> \ingroup realOTHERcomputational
207 *> \par Further Details:
208 * =====================
212 *> The algorithm used in this program is basically backward (forward)
213 *> substitution, with scaling to make the the code robust against
214 *> possible overflow.
216 *> Each eigenvector is normalized so that the element of largest
217 *> magnitude has magnitude 1; here the magnitude of a complex number
218 *> (x,y) is taken to be |x| + |y|.
221 * =====================================================================
222 SUBROUTINE STREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
223 $ LDVR, MM, M, WORK, INFO )
225 * -- LAPACK computational routine (version 3.4.0) --
226 * -- LAPACK is a software package provided by Univ. of Tennessee, --
227 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
230 * .. Scalar Arguments ..
231 CHARACTER HOWMNY, SIDE
232 INTEGER INFO, LDT, LDVL, LDVR, M, MM, N
234 * .. Array Arguments ..
236 REAL T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
240 * =====================================================================
244 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
246 * .. Local Scalars ..
247 LOGICAL ALLV, BOTHV, LEFTV, OVER, PAIR, RIGHTV, SOMEV
248 INTEGER I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI, N2
249 REAL BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE,
250 $ SMIN, SMLNUM, ULP, UNFL, VCRIT, VMAX, WI, WR,
253 * .. External Functions ..
257 EXTERNAL LSAME, ISAMAX, SDOT, SLAMCH
259 * .. External Subroutines ..
260 EXTERNAL SAXPY, SCOPY, SGEMV, SLABAD, SLALN2, SSCAL,
263 * .. Intrinsic Functions ..
264 INTRINSIC ABS, MAX, SQRT
269 * .. Executable Statements ..
271 * Decode and test the input parameters
273 BOTHV = LSAME( SIDE, 'B' )
274 RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
275 LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
277 ALLV = LSAME( HOWMNY, 'A' )
278 OVER = LSAME( HOWMNY, 'B' )
279 SOMEV = LSAME( HOWMNY, 'S' )
282 IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
284 ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
286 ELSE IF( N.LT.0 ) THEN
288 ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
290 ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
292 ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
296 * Set M to the number of columns required to store the selected
297 * eigenvectors, standardize the array SELECT if necessary, and
306 SELECT( J ) = .FALSE.
309 IF( T( J+1, J ).EQ.ZERO ) THEN
314 IF( SELECT( J ) .OR. SELECT( J+1 ) ) THEN
334 CALL XERBLA( 'STREVC', -INFO )
338 * Quick return if possible.
343 * Set the constants to control overflow.
345 UNFL = SLAMCH( 'Safe minimum' )
347 CALL SLABAD( UNFL, OVFL )
348 ULP = SLAMCH( 'Precision' )
349 SMLNUM = UNFL*( N / ULP )
350 BIGNUM = ( ONE-ULP ) / SMLNUM
352 * Compute 1-norm of each column of strictly upper triangular
353 * part of T to control overflow in triangular solver.
359 WORK( J ) = WORK( J ) + ABS( T( I, J ) )
363 * Index IP is used to specify the real or complex eigenvalue:
364 * IP = 0, real eigenvalue,
365 * 1, first of conjugate complex pair: (wr,wi)
366 * -1, second of conjugate complex pair: (wr,wi)
372 * Compute right eigenvectors.
382 IF( T( KI, KI-1 ).EQ.ZERO )
389 IF( .NOT.SELECT( KI ) )
392 IF( .NOT.SELECT( KI-1 ) )
397 * Compute the KI-th eigenvalue (WR,WI).
402 $ WI = SQRT( ABS( T( KI, KI-1 ) ) )*
403 $ SQRT( ABS( T( KI-1, KI ) ) )
404 SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
408 * Real right eigenvector
412 * Form right-hand side
415 WORK( K+N ) = -T( K, KI )
418 * Solve the upper quasi-triangular system:
419 * (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK.
422 DO 60 J = KI - 1, 1, -1
429 IF( T( J, J-1 ).NE.ZERO ) THEN
437 * 1-by-1 diagonal block
439 CALL SLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
440 $ LDT, ONE, ONE, WORK( J+N ), N, WR,
441 $ ZERO, X, 2, SCALE, XNORM, IERR )
443 * Scale X(1,1) to avoid overflow when updating
444 * the right-hand side.
446 IF( XNORM.GT.ONE ) THEN
447 IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
448 X( 1, 1 ) = X( 1, 1 ) / XNORM
449 SCALE = SCALE / XNORM
456 $ CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 )
457 WORK( J+N ) = X( 1, 1 )
459 * Update right-hand side
461 CALL SAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
466 * 2-by-2 diagonal block
468 CALL SLALN2( .FALSE., 2, 1, SMIN, ONE,
469 $ T( J-1, J-1 ), LDT, ONE, ONE,
470 $ WORK( J-1+N ), N, WR, ZERO, X, 2,
471 $ SCALE, XNORM, IERR )
473 * Scale X(1,1) and X(2,1) to avoid overflow when
474 * updating the right-hand side.
476 IF( XNORM.GT.ONE ) THEN
477 BETA = MAX( WORK( J-1 ), WORK( J ) )
478 IF( BETA.GT.BIGNUM / XNORM ) THEN
479 X( 1, 1 ) = X( 1, 1 ) / XNORM
480 X( 2, 1 ) = X( 2, 1 ) / XNORM
481 SCALE = SCALE / XNORM
488 $ CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 )
489 WORK( J-1+N ) = X( 1, 1 )
490 WORK( J+N ) = X( 2, 1 )
492 * Update right-hand side
494 CALL SAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
496 CALL SAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
501 * Copy the vector x or Q*x to VR and normalize.
504 CALL SCOPY( KI, WORK( 1+N ), 1, VR( 1, IS ), 1 )
506 II = ISAMAX( KI, VR( 1, IS ), 1 )
507 REMAX = ONE / ABS( VR( II, IS ) )
508 CALL SSCAL( KI, REMAX, VR( 1, IS ), 1 )
515 $ CALL SGEMV( 'N', N, KI-1, ONE, VR, LDVR,
516 $ WORK( 1+N ), 1, WORK( KI+N ),
519 II = ISAMAX( N, VR( 1, KI ), 1 )
520 REMAX = ONE / ABS( VR( II, KI ) )
521 CALL SSCAL( N, REMAX, VR( 1, KI ), 1 )
526 * Complex right eigenvector.
529 * [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0.
530 * [ (T(KI,KI-1) T(KI,KI) ) ]
532 IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN
534 WORK( KI+N2 ) = WI / T( KI-1, KI )
536 WORK( KI-1+N ) = -WI / T( KI, KI-1 )
540 WORK( KI-1+N2 ) = ZERO
542 * Form right-hand side
545 WORK( K+N ) = -WORK( KI-1+N )*T( K, KI-1 )
546 WORK( K+N2 ) = -WORK( KI+N2 )*T( K, KI )
549 * Solve upper quasi-triangular system:
550 * (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2)
553 DO 90 J = KI - 2, 1, -1
560 IF( T( J, J-1 ).NE.ZERO ) THEN
568 * 1-by-1 diagonal block
570 CALL SLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
571 $ LDT, ONE, ONE, WORK( J+N ), N, WR, WI,
572 $ X, 2, SCALE, XNORM, IERR )
574 * Scale X(1,1) and X(1,2) to avoid overflow when
575 * updating the right-hand side.
577 IF( XNORM.GT.ONE ) THEN
578 IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
579 X( 1, 1 ) = X( 1, 1 ) / XNORM
580 X( 1, 2 ) = X( 1, 2 ) / XNORM
581 SCALE = SCALE / XNORM
587 IF( SCALE.NE.ONE ) THEN
588 CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 )
589 CALL SSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
591 WORK( J+N ) = X( 1, 1 )
592 WORK( J+N2 ) = X( 1, 2 )
594 * Update the right-hand side
596 CALL SAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
598 CALL SAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1,
603 * 2-by-2 diagonal block
605 CALL SLALN2( .FALSE., 2, 2, SMIN, ONE,
606 $ T( J-1, J-1 ), LDT, ONE, ONE,
607 $ WORK( J-1+N ), N, WR, WI, X, 2, SCALE,
610 * Scale X to avoid overflow when updating
611 * the right-hand side.
613 IF( XNORM.GT.ONE ) THEN
614 BETA = MAX( WORK( J-1 ), WORK( J ) )
615 IF( BETA.GT.BIGNUM / XNORM ) THEN
617 X( 1, 1 ) = X( 1, 1 )*REC
618 X( 1, 2 ) = X( 1, 2 )*REC
619 X( 2, 1 ) = X( 2, 1 )*REC
620 X( 2, 2 ) = X( 2, 2 )*REC
627 IF( SCALE.NE.ONE ) THEN
628 CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 )
629 CALL SSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
631 WORK( J-1+N ) = X( 1, 1 )
632 WORK( J+N ) = X( 2, 1 )
633 WORK( J-1+N2 ) = X( 1, 2 )
634 WORK( J+N2 ) = X( 2, 2 )
636 * Update the right-hand side
638 CALL SAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
640 CALL SAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
642 CALL SAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1,
644 CALL SAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1,
649 * Copy the vector x or Q*x to VR and normalize.
652 CALL SCOPY( KI, WORK( 1+N ), 1, VR( 1, IS-1 ), 1 )
653 CALL SCOPY( KI, WORK( 1+N2 ), 1, VR( 1, IS ), 1 )
657 EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+
658 $ ABS( VR( K, IS ) ) )
662 CALL SSCAL( KI, REMAX, VR( 1, IS-1 ), 1 )
663 CALL SSCAL( KI, REMAX, VR( 1, IS ), 1 )
673 CALL SGEMV( 'N', N, KI-2, ONE, VR, LDVR,
674 $ WORK( 1+N ), 1, WORK( KI-1+N ),
676 CALL SGEMV( 'N', N, KI-2, ONE, VR, LDVR,
677 $ WORK( 1+N2 ), 1, WORK( KI+N2 ),
680 CALL SSCAL( N, WORK( KI-1+N ), VR( 1, KI-1 ), 1 )
681 CALL SSCAL( N, WORK( KI+N2 ), VR( 1, KI ), 1 )
686 EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+
687 $ ABS( VR( K, KI ) ) )
690 CALL SSCAL( N, REMAX, VR( 1, KI-1 ), 1 )
691 CALL SSCAL( N, REMAX, VR( 1, KI ), 1 )
708 * Compute left eigenvectors.
718 IF( T( KI+1, KI ).EQ.ZERO )
724 IF( .NOT.SELECT( KI ) )
728 * Compute the KI-th eigenvalue (WR,WI).
733 $ WI = SQRT( ABS( T( KI, KI+1 ) ) )*
734 $ SQRT( ABS( T( KI+1, KI ) ) )
735 SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
739 * Real left eigenvector.
743 * Form right-hand side
746 WORK( K+N ) = -T( KI, K )
749 * Solve the quasi-triangular system:
750 * (T(KI+1:N,KI+1:N) - WR)**T*X = SCALE*WORK
763 IF( T( J+1, J ).NE.ZERO ) THEN
771 * 1-by-1 diagonal block
773 * Scale if necessary to avoid overflow when forming
774 * the right-hand side.
776 IF( WORK( J ).GT.VCRIT ) THEN
778 CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
783 WORK( J+N ) = WORK( J+N ) -
784 $ SDOT( J-KI-1, T( KI+1, J ), 1,
785 $ WORK( KI+1+N ), 1 )
787 * Solve (T(J,J)-WR)**T*X = WORK
789 CALL SLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
790 $ LDT, ONE, ONE, WORK( J+N ), N, WR,
791 $ ZERO, X, 2, SCALE, XNORM, IERR )
796 $ CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
797 WORK( J+N ) = X( 1, 1 )
798 VMAX = MAX( ABS( WORK( J+N ) ), VMAX )
799 VCRIT = BIGNUM / VMAX
803 * 2-by-2 diagonal block
805 * Scale if necessary to avoid overflow when forming
806 * the right-hand side.
808 BETA = MAX( WORK( J ), WORK( J+1 ) )
809 IF( BETA.GT.VCRIT ) THEN
811 CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
816 WORK( J+N ) = WORK( J+N ) -
817 $ SDOT( J-KI-1, T( KI+1, J ), 1,
818 $ WORK( KI+1+N ), 1 )
820 WORK( J+1+N ) = WORK( J+1+N ) -
821 $ SDOT( J-KI-1, T( KI+1, J+1 ), 1,
822 $ WORK( KI+1+N ), 1 )
825 * [T(J,J)-WR T(J,J+1) ]**T* X = SCALE*( WORK1 )
826 * [T(J+1,J) T(J+1,J+1)-WR] ( WORK2 )
828 CALL SLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ),
829 $ LDT, ONE, ONE, WORK( J+N ), N, WR,
830 $ ZERO, X, 2, SCALE, XNORM, IERR )
835 $ CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
836 WORK( J+N ) = X( 1, 1 )
837 WORK( J+1+N ) = X( 2, 1 )
839 VMAX = MAX( ABS( WORK( J+N ) ),
840 $ ABS( WORK( J+1+N ) ), VMAX )
841 VCRIT = BIGNUM / VMAX
846 * Copy the vector x or Q*x to VL and normalize.
849 CALL SCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
851 II = ISAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
852 REMAX = ONE / ABS( VL( II, IS ) )
853 CALL SSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
862 $ CALL SGEMV( 'N', N, N-KI, ONE, VL( 1, KI+1 ), LDVL,
863 $ WORK( KI+1+N ), 1, WORK( KI+N ),
866 II = ISAMAX( N, VL( 1, KI ), 1 )
867 REMAX = ONE / ABS( VL( II, KI ) )
868 CALL SSCAL( N, REMAX, VL( 1, KI ), 1 )
874 * Complex left eigenvector.
877 * ((T(KI,KI) T(KI,KI+1) )**T - (WR - I* WI))*X = 0.
878 * ((T(KI+1,KI) T(KI+1,KI+1)) )
880 IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN
881 WORK( KI+N ) = WI / T( KI, KI+1 )
882 WORK( KI+1+N2 ) = ONE
885 WORK( KI+1+N2 ) = -WI / T( KI+1, KI )
887 WORK( KI+1+N ) = ZERO
890 * Form right-hand side
893 WORK( K+N ) = -WORK( KI+N )*T( KI, K )
894 WORK( K+N2 ) = -WORK( KI+1+N2 )*T( KI+1, K )
897 * Solve complex quasi-triangular system:
898 * ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2
911 IF( T( J+1, J ).NE.ZERO ) THEN
919 * 1-by-1 diagonal block
921 * Scale if necessary to avoid overflow when
922 * forming the right-hand side elements.
924 IF( WORK( J ).GT.VCRIT ) THEN
926 CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
927 CALL SSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
932 WORK( J+N ) = WORK( J+N ) -
933 $ SDOT( J-KI-2, T( KI+2, J ), 1,
934 $ WORK( KI+2+N ), 1 )
935 WORK( J+N2 ) = WORK( J+N2 ) -
936 $ SDOT( J-KI-2, T( KI+2, J ), 1,
937 $ WORK( KI+2+N2 ), 1 )
939 * Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2
941 CALL SLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
942 $ LDT, ONE, ONE, WORK( J+N ), N, WR,
943 $ -WI, X, 2, SCALE, XNORM, IERR )
947 IF( SCALE.NE.ONE ) THEN
948 CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
949 CALL SSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
951 WORK( J+N ) = X( 1, 1 )
952 WORK( J+N2 ) = X( 1, 2 )
953 VMAX = MAX( ABS( WORK( J+N ) ),
954 $ ABS( WORK( J+N2 ) ), VMAX )
955 VCRIT = BIGNUM / VMAX
959 * 2-by-2 diagonal block
961 * Scale if necessary to avoid overflow when forming
962 * the right-hand side elements.
964 BETA = MAX( WORK( J ), WORK( J+1 ) )
965 IF( BETA.GT.VCRIT ) THEN
967 CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
968 CALL SSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
973 WORK( J+N ) = WORK( J+N ) -
974 $ SDOT( J-KI-2, T( KI+2, J ), 1,
975 $ WORK( KI+2+N ), 1 )
977 WORK( J+N2 ) = WORK( J+N2 ) -
978 $ SDOT( J-KI-2, T( KI+2, J ), 1,
979 $ WORK( KI+2+N2 ), 1 )
981 WORK( J+1+N ) = WORK( J+1+N ) -
982 $ SDOT( J-KI-2, T( KI+2, J+1 ), 1,
983 $ WORK( KI+2+N ), 1 )
985 WORK( J+1+N2 ) = WORK( J+1+N2 ) -
986 $ SDOT( J-KI-2, T( KI+2, J+1 ), 1,
987 $ WORK( KI+2+N2 ), 1 )
989 * Solve 2-by-2 complex linear equation
990 * ([T(j,j) T(j,j+1) ]**T-(wr-i*wi)*I)*X = SCALE*B
991 * ([T(j+1,j) T(j+1,j+1)] )
993 CALL SLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ),
994 $ LDT, ONE, ONE, WORK( J+N ), N, WR,
995 $ -WI, X, 2, SCALE, XNORM, IERR )
999 IF( SCALE.NE.ONE ) THEN
1000 CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
1001 CALL SSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
1003 WORK( J+N ) = X( 1, 1 )
1004 WORK( J+N2 ) = X( 1, 2 )
1005 WORK( J+1+N ) = X( 2, 1 )
1006 WORK( J+1+N2 ) = X( 2, 2 )
1007 VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ),
1008 $ ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ), VMAX )
1009 VCRIT = BIGNUM / VMAX
1014 * Copy the vector x or Q*x to VL and normalize.
1016 IF( .NOT.OVER ) THEN
1017 CALL SCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
1018 CALL SCOPY( N-KI+1, WORK( KI+N2 ), 1, VL( KI, IS+1 ),
1023 EMAX = MAX( EMAX, ABS( VL( K, IS ) )+
1024 $ ABS( VL( K, IS+1 ) ) )
1027 CALL SSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
1028 CALL SSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 )
1030 DO 230 K = 1, KI - 1
1032 VL( K, IS+1 ) = ZERO
1035 IF( KI.LT.N-1 ) THEN
1036 CALL SGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
1037 $ LDVL, WORK( KI+2+N ), 1, WORK( KI+N ),
1039 CALL SGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
1040 $ LDVL, WORK( KI+2+N2 ), 1,
1041 $ WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
1043 CALL SSCAL( N, WORK( KI+N ), VL( 1, KI ), 1 )
1044 CALL SSCAL( N, WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
1049 EMAX = MAX( EMAX, ABS( VL( K, KI ) )+
1050 $ ABS( VL( K, KI+1 ) ) )
1053 CALL SSCAL( N, REMAX, VL( 1, KI ), 1 )
1054 CALL SSCAL( N, REMAX, VL( 1, KI+1 ), 1 )