dtgevc.f
来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 1,148 行 · 第 1/3 页
F
1,148 行
TEMP = ZERO
TEMP2 = ZERO
IF( S( J, J-1 ).EQ.ZERO ) THEN
IEND = J - 1
ELSE
IEND = J - 2
END IF
DO 30 I = 1, IEND
TEMP = TEMP + ABS( S( I, J ) )
TEMP2 = TEMP2 + ABS( P( I, J ) )
30 CONTINUE
WORK( J ) = TEMP
WORK( N+J ) = TEMP2
DO 40 I = IEND + 1, MIN( J+1, N )
TEMP = TEMP + ABS( S( I, J ) )
TEMP2 = TEMP2 + ABS( P( I, J ) )
40 CONTINUE
ANORM = MAX( ANORM, TEMP )
BNORM = MAX( BNORM, TEMP2 )
50 CONTINUE
*
ASCALE = ONE / MAX( ANORM, SAFMIN )
BSCALE = ONE / MAX( BNORM, SAFMIN )
*
* Left eigenvectors
*
IF( COMPL ) THEN
IEIG = 0
*
* Main loop over eigenvalues
*
ILCPLX = .FALSE.
DO 220 JE = 1, N
*
* Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
* (b) this would be the second of a complex pair.
* Check for complex eigenvalue, so as to be sure of which
* entry(-ies) of SELECT to look at.
*
IF( ILCPLX ) THEN
ILCPLX = .FALSE.
GO TO 220
END IF
NW = 1
IF( JE.LT.N ) THEN
IF( S( JE+1, JE ).NE.ZERO ) THEN
ILCPLX = .TRUE.
NW = 2
END IF
END IF
IF( ILALL ) THEN
ILCOMP = .TRUE.
ELSE IF( ILCPLX ) THEN
ILCOMP = SELECT( JE ) .OR. SELECT( JE+1 )
ELSE
ILCOMP = SELECT( JE )
END IF
IF( .NOT.ILCOMP )
$ GO TO 220
*
* Decide if (a) singular pencil, (b) real eigenvalue, or
* (c) complex eigenvalue.
*
IF( .NOT.ILCPLX ) THEN
IF( ABS( S( JE, JE ) ).LE.SAFMIN .AND.
$ ABS( P( JE, JE ) ).LE.SAFMIN ) THEN
*
* Singular matrix pencil -- return unit eigenvector
*
IEIG = IEIG + 1
DO 60 JR = 1, N
VL( JR, IEIG ) = ZERO
60 CONTINUE
VL( IEIG, IEIG ) = ONE
GO TO 220
END IF
END IF
*
* Clear vector
*
DO 70 JR = 1, NW*N
WORK( 2*N+JR ) = ZERO
70 CONTINUE
* T
* Compute coefficients in ( a A - b B ) y = 0
* a is ACOEF
* b is BCOEFR + i*BCOEFI
*
IF( .NOT.ILCPLX ) THEN
*
* Real eigenvalue
*
TEMP = ONE / MAX( ABS( S( JE, JE ) )*ASCALE,
$ ABS( P( JE, JE ) )*BSCALE, SAFMIN )
SALFAR = ( TEMP*S( JE, JE ) )*ASCALE
SBETA = ( TEMP*P( JE, JE ) )*BSCALE
ACOEF = SBETA*ASCALE
BCOEFR = SALFAR*BSCALE
BCOEFI = ZERO
*
* Scale to avoid underflow
*
SCALE = ONE
LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEF ).LT.SMALL
LSB = ABS( SALFAR ).GE.SAFMIN .AND. ABS( BCOEFR ).LT.
$ SMALL
IF( LSA )
$ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
IF( LSB )
$ SCALE = MAX( SCALE, ( SMALL / ABS( SALFAR ) )*
$ MIN( BNORM, BIG ) )
IF( LSA .OR. LSB ) THEN
SCALE = MIN( SCALE, ONE /
$ ( SAFMIN*MAX( ONE, ABS( ACOEF ),
$ ABS( BCOEFR ) ) ) )
IF( LSA ) THEN
ACOEF = ASCALE*( SCALE*SBETA )
ELSE
ACOEF = SCALE*ACOEF
END IF
IF( LSB ) THEN
BCOEFR = BSCALE*( SCALE*SALFAR )
ELSE
BCOEFR = SCALE*BCOEFR
END IF
END IF
ACOEFA = ABS( ACOEF )
BCOEFA = ABS( BCOEFR )
*
* First component is 1
*
WORK( 2*N+JE ) = ONE
XMAX = ONE
ELSE
*
* Complex eigenvalue
*
CALL DLAG2( S( JE, JE ), LDS, P( JE, JE ), LDP,
$ SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2,
$ BCOEFI )
BCOEFI = -BCOEFI
IF( BCOEFI.EQ.ZERO ) THEN
INFO = JE
RETURN
END IF
*
* Scale to avoid over/underflow
*
ACOEFA = ABS( ACOEF )
BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
SCALE = ONE
IF( ACOEFA*ULP.LT.SAFMIN .AND. ACOEFA.GE.SAFMIN )
$ SCALE = ( SAFMIN / ULP ) / ACOEFA
IF( BCOEFA*ULP.LT.SAFMIN .AND. BCOEFA.GE.SAFMIN )
$ SCALE = MAX( SCALE, ( SAFMIN / ULP ) / BCOEFA )
IF( SAFMIN*ACOEFA.GT.ASCALE )
$ SCALE = ASCALE / ( SAFMIN*ACOEFA )
IF( SAFMIN*BCOEFA.GT.BSCALE )
$ SCALE = MIN( SCALE, BSCALE / ( SAFMIN*BCOEFA ) )
IF( SCALE.NE.ONE ) THEN
ACOEF = SCALE*ACOEF
ACOEFA = ABS( ACOEF )
BCOEFR = SCALE*BCOEFR
BCOEFI = SCALE*BCOEFI
BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
END IF
*
* Compute first two components of eigenvector
*
TEMP = ACOEF*S( JE+1, JE )
TEMP2R = ACOEF*S( JE, JE ) - BCOEFR*P( JE, JE )
TEMP2I = -BCOEFI*P( JE, JE )
IF( ABS( TEMP ).GT.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN
WORK( 2*N+JE ) = ONE
WORK( 3*N+JE ) = ZERO
WORK( 2*N+JE+1 ) = -TEMP2R / TEMP
WORK( 3*N+JE+1 ) = -TEMP2I / TEMP
ELSE
WORK( 2*N+JE+1 ) = ONE
WORK( 3*N+JE+1 ) = ZERO
TEMP = ACOEF*S( JE, JE+1 )
WORK( 2*N+JE ) = ( BCOEFR*P( JE+1, JE+1 )-ACOEF*
$ S( JE+1, JE+1 ) ) / TEMP
WORK( 3*N+JE ) = BCOEFI*P( JE+1, JE+1 ) / TEMP
END IF
XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ),
$ ABS( WORK( 2*N+JE+1 ) )+ABS( WORK( 3*N+JE+1 ) ) )
END IF
*
DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
*
* T
* Triangular solve of (a A - b B) y = 0
*
* T
* (rowwise in (a A - b B) , or columnwise in (a A - b B) )
*
IL2BY2 = .FALSE.
*
DO 160 J = JE + NW, N
IF( IL2BY2 ) THEN
IL2BY2 = .FALSE.
GO TO 160
END IF
*
NA = 1
BDIAG( 1 ) = P( J, J )
IF( J.LT.N ) THEN
IF( S( J+1, J ).NE.ZERO ) THEN
IL2BY2 = .TRUE.
BDIAG( 2 ) = P( J+1, J+1 )
NA = 2
END IF
END IF
*
* Check whether scaling is necessary for dot products
*
XSCALE = ONE / MAX( ONE, XMAX )
TEMP = MAX( WORK( J ), WORK( N+J ),
$ ACOEFA*WORK( J )+BCOEFA*WORK( N+J ) )
IF( IL2BY2 )
$ TEMP = MAX( TEMP, WORK( J+1 ), WORK( N+J+1 ),
$ ACOEFA*WORK( J+1 )+BCOEFA*WORK( N+J+1 ) )
IF( TEMP.GT.BIGNUM*XSCALE ) THEN
DO 90 JW = 0, NW - 1
DO 80 JR = JE, J - 1
WORK( ( JW+2 )*N+JR ) = XSCALE*
$ WORK( ( JW+2 )*N+JR )
80 CONTINUE
90 CONTINUE
XMAX = XMAX*XSCALE
END IF
*
* Compute dot products
*
* j-1
* SUM = sum conjg( a*S(k,j) - b*P(k,j) )*x(k)
* k=je
*
* To reduce the op count, this is done as
*
* _ j-1 _ j-1
* a*conjg( sum S(k,j)*x(k) ) - b*conjg( sum P(k,j)*x(k) )
* k=je k=je
*
* which may cause underflow problems if A or B are close
* to underflow. (E.g., less than SMALL.)
*
*
* A series of compiler directives to defeat vectorization
* for the next loop
*
*$PL$ CMCHAR=' '
CDIR$ NEXTSCALAR
C$DIR SCALAR
CDIR$ NEXT SCALAR
CVD$L NOVECTOR
CDEC$ NOVECTOR
CVD$ NOVECTOR
*VDIR NOVECTOR
*VOCL LOOP,SCALAR
CIBM PREFER SCALAR
*$PL$ CMCHAR='*'
*
DO 120 JW = 1, NW
*
*$PL$ CMCHAR=' '
CDIR$ NEXTSCALAR
C$DIR SCALAR
CDIR$ NEXT SCALAR
CVD$L NOVECTOR
CDEC$ NOVECTOR
CVD$ NOVECTOR
*VDIR NOVECTOR
*VOCL LOOP,SCALAR
CIBM PREFER SCALAR
*$PL$ CMCHAR='*'
*
DO 110 JA = 1, NA
SUMS( JA, JW ) = ZERO
SUMP( JA, JW ) = ZERO
*
DO 100 JR = JE, J - 1
SUMS( JA, JW ) = SUMS( JA, JW ) +
$ S( JR, J+JA-1 )*
$ WORK( ( JW+1 )*N+JR )
SUMP( JA, JW ) = SUMP( JA, JW ) +
$ P( JR, J+JA-1 )*
$ WORK( ( JW+1 )*N+JR )
100 CONTINUE
110 CONTINUE
120 CONTINUE
*
*$PL$ CMCHAR=' '
CDIR$ NEXTSCALAR
C$DIR SCALAR
CDIR$ NEXT SCALAR
CVD$L NOVECTOR
CDEC$ NOVECTOR
CVD$ NOVECTOR
*VDIR NOVECTOR
*VOCL LOOP,SCALAR
CIBM PREFER SCALAR
*$PL$ CMCHAR='*'
*
DO 130 JA = 1, NA
IF( ILCPLX ) THEN
SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) +
$ BCOEFR*SUMP( JA, 1 ) -
$ BCOEFI*SUMP( JA, 2 )
SUM( JA, 2 ) = -ACOEF*SUMS( JA, 2 ) +
$ BCOEFR*SUMP( JA, 2 ) +
$ BCOEFI*SUMP( JA, 1 )
ELSE
SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) +
$ BCOEFR*SUMP( JA, 1 )
END IF
130 CONTINUE
*
* T
* Solve ( a A - b B ) y = SUM(,)
* with scaling and perturbation of the denominator
*
CALL DLALN2( .TRUE., NA, NW, DMIN, ACOEF, S( J, J ), LDS,
$ BDIAG( 1 ), BDIAG( 2 ), SUM, 2, BCOEFR,
$ BCOEFI, WORK( 2*N+J ), N, SCALE, TEMP,
$ IINFO )
IF( SCALE.LT.ONE ) THEN
DO 150 JW = 0, NW - 1
DO 140 JR = JE, J - 1
WORK( ( JW+2 )*N+JR ) = SCALE*
$ WORK( ( JW+2 )*N+JR )
140 CONTINUE
150 CONTINUE
XMAX = SCALE*XMAX
END IF
XMAX = MAX( XMAX, TEMP )
160 CONTINUE
*
* Copy eigenvector to VL, back transforming if
* HOWMNY='B'.
*
IEIG = IEIG + 1
IF( ILBACK ) THEN
DO 170 JW = 0, NW - 1
CALL DGEMV( 'N', N, N+1-JE, ONE, VL( 1, JE ), LDVL,
$ WORK( ( JW+2 )*N+JE ), 1, ZERO,
$ WORK( ( JW+4 )*N+1 ), 1 )
170 CONTINUE
CALL DLACPY( ' ', N, NW, WORK( 4*N+1 ), N, VL( 1, JE ),
$ LDVL )
IBEG = 1
ELSE
CALL DLACPY( ' ', N, NW, WORK( 2*N+1 ), N, VL( 1, IEIG ),
$ LDVL )
IBEG = JE
END IF
*
* Scale eigenvector
*
XMAX = ZERO
IF( ILCPLX ) THEN
DO 180 J = IBEG, N
XMAX = MAX( XMAX, ABS( VL( J, IEIG ) )+
$ ABS( VL( J, IEIG+1 ) ) )
180 CONTINUE
ELSE
DO 190 J = IBEG, N
XMAX = MAX( XMAX, ABS( VL( J, IEIG ) ) )
190 CONTINUE
END IF
*
IF( XMAX.GT.SAFMIN ) THEN
XSCALE = ONE / XMAX
*
DO 210 JW = 0, NW - 1
DO 200 JR = IBEG, N
VL( JR, IEIG+JW ) = XSCALE*VL( JR, IEIG+JW )
200 CONTINUE
210 CONTINUE
END IF
IEIG = IEIG + NW - 1
*
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?