strevc.f
来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 982 行 · 第 1/3 页
F
982 行
J1 = J
J2 = J
JNXT = J - 1
IF( J.GT.1 ) THEN
IF( T( J, J-1 ).NE.ZERO ) THEN
J1 = J - 1
JNXT = J - 2
END IF
END IF
*
IF( J1.EQ.J2 ) THEN
*
* 1-by-1 diagonal block
*
CALL SLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
$ LDT, ONE, ONE, WORK( J+N ), N, WR,
$ ZERO, X, 2, SCALE, XNORM, IERR )
*
* Scale X(1,1) to avoid overflow when updating
* the right-hand side.
*
IF( XNORM.GT.ONE ) THEN
IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
X( 1, 1 ) = X( 1, 1 ) / XNORM
SCALE = SCALE / XNORM
END IF
END IF
*
* Scale if necessary
*
IF( SCALE.NE.ONE )
$ CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 )
WORK( J+N ) = X( 1, 1 )
*
* Update right-hand side
*
CALL SAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
$ WORK( 1+N ), 1 )
*
ELSE
*
* 2-by-2 diagonal block
*
CALL SLALN2( .FALSE., 2, 1, SMIN, ONE,
$ T( J-1, J-1 ), LDT, ONE, ONE,
$ WORK( J-1+N ), N, WR, ZERO, X, 2,
$ SCALE, XNORM, IERR )
*
* Scale X(1,1) and X(2,1) to avoid overflow when
* updating the right-hand side.
*
IF( XNORM.GT.ONE ) THEN
BETA = MAX( WORK( J-1 ), WORK( J ) )
IF( BETA.GT.BIGNUM / XNORM ) THEN
X( 1, 1 ) = X( 1, 1 ) / XNORM
X( 2, 1 ) = X( 2, 1 ) / XNORM
SCALE = SCALE / XNORM
END IF
END IF
*
* Scale if necessary
*
IF( SCALE.NE.ONE )
$ CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 )
WORK( J-1+N ) = X( 1, 1 )
WORK( J+N ) = X( 2, 1 )
*
* Update right-hand side
*
CALL SAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
$ WORK( 1+N ), 1 )
CALL SAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
$ WORK( 1+N ), 1 )
END IF
60 CONTINUE
*
* Copy the vector x or Q*x to VR and normalize.
*
IF( .NOT.OVER ) THEN
CALL SCOPY( KI, WORK( 1+N ), 1, VR( 1, IS ), 1 )
*
II = ISAMAX( KI, VR( 1, IS ), 1 )
REMAX = ONE / ABS( VR( II, IS ) )
CALL SSCAL( KI, REMAX, VR( 1, IS ), 1 )
*
DO 70 K = KI + 1, N
VR( K, IS ) = ZERO
70 CONTINUE
ELSE
IF( KI.GT.1 )
$ CALL SGEMV( 'N', N, KI-1, ONE, VR, LDVR,
$ WORK( 1+N ), 1, WORK( KI+N ),
$ VR( 1, KI ), 1 )
*
II = ISAMAX( N, VR( 1, KI ), 1 )
REMAX = ONE / ABS( VR( II, KI ) )
CALL SSCAL( N, REMAX, VR( 1, KI ), 1 )
END IF
*
ELSE
*
* Complex right eigenvector.
*
* Initial solve
* [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0.
* [ (T(KI,KI-1) T(KI,KI) ) ]
*
IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN
WORK( KI-1+N ) = ONE
WORK( KI+N2 ) = WI / T( KI-1, KI )
ELSE
WORK( KI-1+N ) = -WI / T( KI, KI-1 )
WORK( KI+N2 ) = ONE
END IF
WORK( KI+N ) = ZERO
WORK( KI-1+N2 ) = ZERO
*
* Form right-hand side
*
DO 80 K = 1, KI - 2
WORK( K+N ) = -WORK( KI-1+N )*T( K, KI-1 )
WORK( K+N2 ) = -WORK( KI+N2 )*T( K, KI )
80 CONTINUE
*
* Solve upper quasi-triangular system:
* (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2)
*
JNXT = KI - 2
DO 90 J = KI - 2, 1, -1
IF( J.GT.JNXT )
$ GO TO 90
J1 = J
J2 = J
JNXT = J - 1
IF( J.GT.1 ) THEN
IF( T( J, J-1 ).NE.ZERO ) THEN
J1 = J - 1
JNXT = J - 2
END IF
END IF
*
IF( J1.EQ.J2 ) THEN
*
* 1-by-1 diagonal block
*
CALL SLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
$ LDT, ONE, ONE, WORK( J+N ), N, WR, WI,
$ X, 2, SCALE, XNORM, IERR )
*
* Scale X(1,1) and X(1,2) to avoid overflow when
* updating the right-hand side.
*
IF( XNORM.GT.ONE ) THEN
IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
X( 1, 1 ) = X( 1, 1 ) / XNORM
X( 1, 2 ) = X( 1, 2 ) / XNORM
SCALE = SCALE / XNORM
END IF
END IF
*
* Scale if necessary
*
IF( SCALE.NE.ONE ) THEN
CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 )
CALL SSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
END IF
WORK( J+N ) = X( 1, 1 )
WORK( J+N2 ) = X( 1, 2 )
*
* Update the right-hand side
*
CALL SAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
$ WORK( 1+N ), 1 )
CALL SAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1,
$ WORK( 1+N2 ), 1 )
*
ELSE
*
* 2-by-2 diagonal block
*
CALL SLALN2( .FALSE., 2, 2, SMIN, ONE,
$ T( J-1, J-1 ), LDT, ONE, ONE,
$ WORK( J-1+N ), N, WR, WI, X, 2, SCALE,
$ XNORM, IERR )
*
* Scale X to avoid overflow when updating
* the right-hand side.
*
IF( XNORM.GT.ONE ) THEN
BETA = MAX( WORK( J-1 ), WORK( J ) )
IF( BETA.GT.BIGNUM / XNORM ) THEN
REC = ONE / XNORM
X( 1, 1 ) = X( 1, 1 )*REC
X( 1, 2 ) = X( 1, 2 )*REC
X( 2, 1 ) = X( 2, 1 )*REC
X( 2, 2 ) = X( 2, 2 )*REC
SCALE = SCALE*REC
END IF
END IF
*
* Scale if necessary
*
IF( SCALE.NE.ONE ) THEN
CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 )
CALL SSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
END IF
WORK( J-1+N ) = X( 1, 1 )
WORK( J+N ) = X( 2, 1 )
WORK( J-1+N2 ) = X( 1, 2 )
WORK( J+N2 ) = X( 2, 2 )
*
* Update the right-hand side
*
CALL SAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
$ WORK( 1+N ), 1 )
CALL SAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
$ WORK( 1+N ), 1 )
CALL SAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1,
$ WORK( 1+N2 ), 1 )
CALL SAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1,
$ WORK( 1+N2 ), 1 )
END IF
90 CONTINUE
*
* Copy the vector x or Q*x to VR and normalize.
*
IF( .NOT.OVER ) THEN
CALL SCOPY( KI, WORK( 1+N ), 1, VR( 1, IS-1 ), 1 )
CALL SCOPY( KI, WORK( 1+N2 ), 1, VR( 1, IS ), 1 )
*
EMAX = ZERO
DO 100 K = 1, KI
EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+
$ ABS( VR( K, IS ) ) )
100 CONTINUE
*
REMAX = ONE / EMAX
CALL SSCAL( KI, REMAX, VR( 1, IS-1 ), 1 )
CALL SSCAL( KI, REMAX, VR( 1, IS ), 1 )
*
DO 110 K = KI + 1, N
VR( K, IS-1 ) = ZERO
VR( K, IS ) = ZERO
110 CONTINUE
*
ELSE
*
IF( KI.GT.2 ) THEN
CALL SGEMV( 'N', N, KI-2, ONE, VR, LDVR,
$ WORK( 1+N ), 1, WORK( KI-1+N ),
$ VR( 1, KI-1 ), 1 )
CALL SGEMV( 'N', N, KI-2, ONE, VR, LDVR,
$ WORK( 1+N2 ), 1, WORK( KI+N2 ),
$ VR( 1, KI ), 1 )
ELSE
CALL SSCAL( N, WORK( KI-1+N ), VR( 1, KI-1 ), 1 )
CALL SSCAL( N, WORK( KI+N2 ), VR( 1, KI ), 1 )
END IF
*
EMAX = ZERO
DO 120 K = 1, N
EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+
$ ABS( VR( K, KI ) ) )
120 CONTINUE
REMAX = ONE / EMAX
CALL SSCAL( N, REMAX, VR( 1, KI-1 ), 1 )
CALL SSCAL( N, REMAX, VR( 1, KI ), 1 )
END IF
END IF
*
IS = IS - 1
IF( IP.NE.0 )
$ IS = IS - 1
130 CONTINUE
IF( IP.EQ.1 )
$ IP = 0
IF( IP.EQ.-1 )
$ IP = 1
140 CONTINUE
END IF
*
IF( LEFTV ) THEN
*
* Compute left eigenvectors.
*
IP = 0
IS = 1
DO 260 KI = 1, N
*
IF( IP.EQ.-1 )
$ GO TO 250
IF( KI.EQ.N )
$ GO TO 150
IF( T( KI+1, KI ).EQ.ZERO )
$ GO TO 150
IP = 1
*
150 CONTINUE
IF( SOMEV ) THEN
IF( .NOT.SELECT( KI ) )
$ GO TO 250
END IF
*
* Compute the KI-th eigenvalue (WR,WI).
*
WR = T( KI, KI )
WI = ZERO
IF( IP.NE.0 )
$ WI = SQRT( ABS( T( KI, KI+1 ) ) )*
$ SQRT( ABS( T( KI+1, KI ) ) )
SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
*
IF( IP.EQ.0 ) THEN
*
* Real left eigenvector.
*
WORK( KI+N ) = ONE
*
* Form right-hand side
*
DO 160 K = KI + 1, N
WORK( K+N ) = -T( KI, K )
160 CONTINUE
*
* Solve the quasi-triangular system:
* (T(KI+1:N,KI+1:N) - WR)'*X = SCALE*WORK
*
VMAX = ONE
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?