strevc.f
来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 982 行 · 第 1/3 页
F
982 行
VCRIT = BIGNUM
*
JNXT = KI + 1
DO 170 J = KI + 1, N
IF( J.LT.JNXT )
$ GO TO 170
J1 = J
J2 = J
JNXT = J + 1
IF( J.LT.N ) THEN
IF( T( J+1, J ).NE.ZERO ) THEN
J2 = J + 1
JNXT = J + 2
END IF
END IF
*
IF( J1.EQ.J2 ) THEN
*
* 1-by-1 diagonal block
*
* Scale if necessary to avoid overflow when forming
* the right-hand side.
*
IF( WORK( J ).GT.VCRIT ) THEN
REC = ONE / VMAX
CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
VMAX = ONE
VCRIT = BIGNUM
END IF
*
WORK( J+N ) = WORK( J+N ) -
$ SDOT( J-KI-1, T( KI+1, J ), 1,
$ WORK( KI+1+N ), 1 )
*
* Solve (T(J,J)-WR)'*X = WORK
*
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 if necessary
*
IF( SCALE.NE.ONE )
$ CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
WORK( J+N ) = X( 1, 1 )
VMAX = MAX( ABS( WORK( J+N ) ), VMAX )
VCRIT = BIGNUM / VMAX
*
ELSE
*
* 2-by-2 diagonal block
*
* Scale if necessary to avoid overflow when forming
* the right-hand side.
*
BETA = MAX( WORK( J ), WORK( J+1 ) )
IF( BETA.GT.VCRIT ) THEN
REC = ONE / VMAX
CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
VMAX = ONE
VCRIT = BIGNUM
END IF
*
WORK( J+N ) = WORK( J+N ) -
$ SDOT( J-KI-1, T( KI+1, J ), 1,
$ WORK( KI+1+N ), 1 )
*
WORK( J+1+N ) = WORK( J+1+N ) -
$ SDOT( J-KI-1, T( KI+1, J+1 ), 1,
$ WORK( KI+1+N ), 1 )
*
* Solve
* [T(J,J)-WR T(J,J+1) ]'* X = SCALE*( WORK1 )
* [T(J+1,J) T(J+1,J+1)-WR] ( WORK2 )
*
CALL SLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ),
$ LDT, ONE, ONE, WORK( J+N ), N, WR,
$ ZERO, X, 2, SCALE, XNORM, IERR )
*
* Scale if necessary
*
IF( SCALE.NE.ONE )
$ CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
WORK( J+N ) = X( 1, 1 )
WORK( J+1+N ) = X( 2, 1 )
*
VMAX = MAX( ABS( WORK( J+N ) ),
$ ABS( WORK( J+1+N ) ), VMAX )
VCRIT = BIGNUM / VMAX
*
END IF
170 CONTINUE
*
* Copy the vector x or Q*x to VL and normalize.
*
IF( .NOT.OVER ) THEN
CALL SCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
*
II = ISAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
REMAX = ONE / ABS( VL( II, IS ) )
CALL SSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
*
DO 180 K = 1, KI - 1
VL( K, IS ) = ZERO
180 CONTINUE
*
ELSE
*
IF( KI.LT.N )
$ CALL SGEMV( 'N', N, N-KI, ONE, VL( 1, KI+1 ), LDVL,
$ WORK( KI+1+N ), 1, WORK( KI+N ),
$ VL( 1, KI ), 1 )
*
II = ISAMAX( N, VL( 1, KI ), 1 )
REMAX = ONE / ABS( VL( II, KI ) )
CALL SSCAL( N, REMAX, VL( 1, KI ), 1 )
*
END IF
*
ELSE
*
* Complex left eigenvector.
*
* Initial solve:
* ((T(KI,KI) T(KI,KI+1) )' - (WR - I* WI))*X = 0.
* ((T(KI+1,KI) T(KI+1,KI+1)) )
*
IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN
WORK( KI+N ) = WI / T( KI, KI+1 )
WORK( KI+1+N2 ) = ONE
ELSE
WORK( KI+N ) = ONE
WORK( KI+1+N2 ) = -WI / T( KI+1, KI )
END IF
WORK( KI+1+N ) = ZERO
WORK( KI+N2 ) = ZERO
*
* Form right-hand side
*
DO 190 K = KI + 2, N
WORK( K+N ) = -WORK( KI+N )*T( KI, K )
WORK( K+N2 ) = -WORK( KI+1+N2 )*T( KI+1, K )
190 CONTINUE
*
* Solve complex quasi-triangular system:
* ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2
*
VMAX = ONE
VCRIT = BIGNUM
*
JNXT = KI + 2
DO 200 J = KI + 2, N
IF( J.LT.JNXT )
$ GO TO 200
J1 = J
J2 = J
JNXT = J + 1
IF( J.LT.N ) THEN
IF( T( J+1, J ).NE.ZERO ) THEN
J2 = J + 1
JNXT = J + 2
END IF
END IF
*
IF( J1.EQ.J2 ) THEN
*
* 1-by-1 diagonal block
*
* Scale if necessary to avoid overflow when
* forming the right-hand side elements.
*
IF( WORK( J ).GT.VCRIT ) THEN
REC = ONE / VMAX
CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
CALL SSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
VMAX = ONE
VCRIT = BIGNUM
END IF
*
WORK( J+N ) = WORK( J+N ) -
$ SDOT( J-KI-2, T( KI+2, J ), 1,
$ WORK( KI+2+N ), 1 )
WORK( J+N2 ) = WORK( J+N2 ) -
$ SDOT( J-KI-2, T( KI+2, J ), 1,
$ WORK( KI+2+N2 ), 1 )
*
* Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2
*
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 if necessary
*
IF( SCALE.NE.ONE ) THEN
CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
CALL SSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
END IF
WORK( J+N ) = X( 1, 1 )
WORK( J+N2 ) = X( 1, 2 )
VMAX = MAX( ABS( WORK( J+N ) ),
$ ABS( WORK( J+N2 ) ), VMAX )
VCRIT = BIGNUM / VMAX
*
ELSE
*
* 2-by-2 diagonal block
*
* Scale if necessary to avoid overflow when forming
* the right-hand side elements.
*
BETA = MAX( WORK( J ), WORK( J+1 ) )
IF( BETA.GT.VCRIT ) THEN
REC = ONE / VMAX
CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
CALL SSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
VMAX = ONE
VCRIT = BIGNUM
END IF
*
WORK( J+N ) = WORK( J+N ) -
$ SDOT( J-KI-2, T( KI+2, J ), 1,
$ WORK( KI+2+N ), 1 )
*
WORK( J+N2 ) = WORK( J+N2 ) -
$ SDOT( J-KI-2, T( KI+2, J ), 1,
$ WORK( KI+2+N2 ), 1 )
*
WORK( J+1+N ) = WORK( J+1+N ) -
$ SDOT( J-KI-2, T( KI+2, J+1 ), 1,
$ WORK( KI+2+N ), 1 )
*
WORK( J+1+N2 ) = WORK( J+1+N2 ) -
$ SDOT( J-KI-2, T( KI+2, J+1 ), 1,
$ WORK( KI+2+N2 ), 1 )
*
* Solve 2-by-2 complex linear equation
* ([T(j,j) T(j,j+1) ]'-(wr-i*wi)*I)*X = SCALE*B
* ([T(j+1,j) T(j+1,j+1)] )
*
CALL SLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ),
$ LDT, ONE, ONE, WORK( J+N ), N, WR,
$ -WI, X, 2, SCALE, XNORM, IERR )
*
* Scale if necessary
*
IF( SCALE.NE.ONE ) THEN
CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
CALL SSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
END IF
WORK( J+N ) = X( 1, 1 )
WORK( J+N2 ) = X( 1, 2 )
WORK( J+1+N ) = X( 2, 1 )
WORK( J+1+N2 ) = X( 2, 2 )
VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ),
$ ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ), VMAX )
VCRIT = BIGNUM / VMAX
*
END IF
200 CONTINUE
*
* Copy the vector x or Q*x to VL and normalize.
*
IF( .NOT.OVER ) THEN
CALL SCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
CALL SCOPY( N-KI+1, WORK( KI+N2 ), 1, VL( KI, IS+1 ),
$ 1 )
*
EMAX = ZERO
DO 220 K = KI, N
EMAX = MAX( EMAX, ABS( VL( K, IS ) )+
$ ABS( VL( K, IS+1 ) ) )
220 CONTINUE
REMAX = ONE / EMAX
CALL SSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
CALL SSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 )
*
DO 230 K = 1, KI - 1
VL( K, IS ) = ZERO
VL( K, IS+1 ) = ZERO
230 CONTINUE
ELSE
IF( KI.LT.N-1 ) THEN
CALL SGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
$ LDVL, WORK( KI+2+N ), 1, WORK( KI+N ),
$ VL( 1, KI ), 1 )
CALL SGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
$ LDVL, WORK( KI+2+N2 ), 1,
$ WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
ELSE
CALL SSCAL( N, WORK( KI+N ), VL( 1, KI ), 1 )
CALL SSCAL( N, WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
END IF
*
EMAX = ZERO
DO 240 K = 1, N
EMAX = MAX( EMAX, ABS( VL( K, KI ) )+
$ ABS( VL( K, KI+1 ) ) )
240 CONTINUE
REMAX = ONE / EMAX
CALL SSCAL( N, REMAX, VL( 1, KI ), 1 )
CALL SSCAL( N, REMAX, VL( 1, KI+1 ), 1 )
*
END IF
*
END IF
*
IS = IS + 1
IF( IP.NE.0 )
$ IS = IS + 1
250 CONTINUE
IF( IP.EQ.-1 )
$ IP = 0
IF( IP.EQ.1 )
$ IP = -1
*
260 CONTINUE
*
END IF
*
RETURN
*
* End of STREVC
*
END
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?