dtgevc.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 1,038 行 · 第 1/5 页
HTML
1,038 行
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> (3) v(j) := s / C(j,j)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Step 2 is sometimes called the "dot product" step, since it is an
</span><span class="comment">*</span><span class="comment"> inner product between the j-th row and the portion of the eigenvector
</span><span class="comment">*</span><span class="comment"> that has been computed so far.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> The "columnwise" method consists basically in doing the sums
</span><span class="comment">*</span><span class="comment"> for all the rows in parallel. As each v(j) is computed, the
</span><span class="comment">*</span><span class="comment"> contribution of v(j) times the j-th column of C is added to the
</span><span class="comment">*</span><span class="comment"> partial sums. Since FORTRAN arrays are stored columnwise, this has
</span><span class="comment">*</span><span class="comment"> the advantage that at each step, the elements of C that are accessed
</span><span class="comment">*</span><span class="comment"> are adjacent to one another, whereas with the rowwise method, the
</span><span class="comment">*</span><span class="comment"> elements accessed at a step are spaced LDS (and LDP) words apart.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> When finding left eigenvectors, the matrix in question is the
</span><span class="comment">*</span><span class="comment"> transpose of the one in storage, so the rowwise method then
</span><span class="comment">*</span><span class="comment"> actually accesses columns of A and B at each step, and so is the
</span><span class="comment">*</span><span class="comment"> preferred method.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> .. Parameters ..
</span> DOUBLE PRECISION ZERO, ONE, SAFETY
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0,
$ SAFETY = 1.0D+2 )
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Local Scalars ..
</span> LOGICAL COMPL, COMPR, IL2BY2, ILABAD, ILALL, ILBACK,
$ ILBBAD, ILCOMP, ILCPLX, LSA, LSB
INTEGER I, IBEG, IEIG, IEND, IHWMNY, IINFO, IM, ISIDE,
$ J, JA, JC, JE, JR, JW, NA, NW
DOUBLE PRECISION ACOEF, ACOEFA, ANORM, ASCALE, BCOEFA, BCOEFI,
$ BCOEFR, BIG, BIGNUM, BNORM, BSCALE, CIM2A,
$ CIM2B, CIMAGA, CIMAGB, CRE2A, CRE2B, CREALA,
$ CREALB, DMIN, SAFMIN, SALFAR, SBETA, SCALE,
$ SMALL, TEMP, TEMP2, TEMP2I, TEMP2R, ULP, XMAX,
$ XSCALE
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Local Arrays ..
</span> DOUBLE PRECISION BDIAG( 2 ), SUM( 2, 2 ), SUMS( 2, 2 ),
$ SUMP( 2, 2 )
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. External Functions ..
</span> LOGICAL <a name="LSAME.234"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
DOUBLE PRECISION <a name="DLAMCH.235"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>
EXTERNAL <a name="LSAME.236"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="DLAMCH.236"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. External Subroutines ..
</span> EXTERNAL DGEMV, <a name="DLABAD.239"></a><a href="dlabad.f.html#DLABAD.1">DLABAD</a>, <a name="DLACPY.239"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>, <a name="DLAG2.239"></a><a href="dlag2.f.html#DLAG2.1">DLAG2</a>, <a name="DLALN2.239"></a><a href="dlaln2.f.html#DLALN2.1">DLALN2</a>, <a name="XERBLA.239"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Intrinsic Functions ..
</span> INTRINSIC ABS, MAX, MIN
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Decode and Test the input parameters
</span><span class="comment">*</span><span class="comment">
</span> IF( <a name="LSAME.248"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( HOWMNY, <span class="string">'A'</span> ) ) THEN
IHWMNY = 1
ILALL = .TRUE.
ILBACK = .FALSE.
ELSE IF( <a name="LSAME.252"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( HOWMNY, <span class="string">'S'</span> ) ) THEN
IHWMNY = 2
ILALL = .FALSE.
ILBACK = .FALSE.
ELSE IF( <a name="LSAME.256"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( HOWMNY, <span class="string">'B'</span> ) ) THEN
IHWMNY = 3
ILALL = .TRUE.
ILBACK = .TRUE.
ELSE
IHWMNY = -1
ILALL = .TRUE.
END IF
<span class="comment">*</span><span class="comment">
</span> IF( <a name="LSAME.265"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( SIDE, <span class="string">'R'</span> ) ) THEN
ISIDE = 1
COMPL = .FALSE.
COMPR = .TRUE.
ELSE IF( <a name="LSAME.269"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( SIDE, <span class="string">'L'</span> ) ) THEN
ISIDE = 2
COMPL = .TRUE.
COMPR = .FALSE.
ELSE IF( <a name="LSAME.273"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( SIDE, <span class="string">'B'</span> ) ) THEN
ISIDE = 3
COMPL = .TRUE.
COMPR = .TRUE.
ELSE
ISIDE = -1
END IF
<span class="comment">*</span><span class="comment">
</span> INFO = 0
IF( ISIDE.LT.0 ) THEN
INFO = -1
ELSE IF( IHWMNY.LT.0 ) THEN
INFO = -2
ELSE IF( N.LT.0 ) THEN
INFO = -4
ELSE IF( LDS.LT.MAX( 1, N ) ) THEN
INFO = -6
ELSE IF( LDP.LT.MAX( 1, N ) ) THEN
INFO = -8
END IF
IF( INFO.NE.0 ) THEN
CALL <a name="XERBLA.294"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DTGEVC.294"></a><a href="dtgevc.f.html#DTGEVC.1">DTGEVC</a>'</span>, -INFO )
RETURN
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Count the number of eigenvectors to be computed
</span><span class="comment">*</span><span class="comment">
</span> IF( .NOT.ILALL ) THEN
IM = 0
ILCPLX = .FALSE.
DO 10 J = 1, N
IF( ILCPLX ) THEN
ILCPLX = .FALSE.
GO TO 10
END IF
IF( J.LT.N ) THEN
IF( S( J+1, J ).NE.ZERO )
$ ILCPLX = .TRUE.
END IF
IF( ILCPLX ) THEN
IF( SELECT( J ) .OR. SELECT( J+1 ) )
$ IM = IM + 2
ELSE
IF( SELECT( J ) )
$ IM = IM + 1
END IF
10 CONTINUE
ELSE
IM = N
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Check 2-by-2 diagonal blocks of A, B
</span><span class="comment">*</span><span class="comment">
</span> ILABAD = .FALSE.
ILBBAD = .FALSE.
DO 20 J = 1, N - 1
IF( S( J+1, J ).NE.ZERO ) THEN
IF( P( J, J ).EQ.ZERO .OR. P( J+1, J+1 ).EQ.ZERO .OR.
$ P( J, J+1 ).NE.ZERO )ILBBAD = .TRUE.
IF( J.LT.N-1 ) THEN
IF( S( J+2, J+1 ).NE.ZERO )
$ ILABAD = .TRUE.
END IF
END IF
20 CONTINUE
<span class="comment">*</span><span class="comment">
</span> IF( ILABAD ) THEN
INFO = -5
ELSE IF( ILBBAD ) THEN
INFO = -7
ELSE IF( COMPL .AND. LDVL.LT.N .OR. LDVL.LT.1 ) THEN
INFO = -10
ELSE IF( COMPR .AND. LDVR.LT.N .OR. LDVR.LT.1 ) THEN
INFO = -12
ELSE IF( MM.LT.IM ) THEN
INFO = -13
END IF
IF( INFO.NE.0 ) THEN
CALL <a name="XERBLA.351"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DTGEVC.351"></a><a href="dtgevc.f.html#DTGEVC.1">DTGEVC</a>'</span>, -INFO )
RETURN
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span> M = IM
IF( N.EQ.0 )
$ RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Machine Constants
</span><span class="comment">*</span><span class="comment">
</span> SAFMIN = <a name="DLAMCH.363"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Safe minimum'</span> )
BIG = ONE / SAFMIN
CALL <a name="DLABAD.365"></a><a href="dlabad.f.html#DLABAD.1">DLABAD</a>( SAFMIN, BIG )
ULP = <a name="DLAMCH.366"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Epsilon'</span> )*<a name="DLAMCH.366"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Base'</span> )
SMALL = SAFMIN*N / ULP
BIG = ONE / SMALL
BIGNUM = ONE / ( SAFMIN*N )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Compute the 1-norm of each column of the strictly upper triangular
</span><span class="comment">*</span><span class="comment"> part (i.e., excluding all elements belonging to the diagonal
</span><span class="comment">*</span><span class="comment"> blocks) of A and B to check for possible overflow in the
</span><span class="comment">*</span><span class="comment"> triangular solver.
</span><span class="comment">*</span><span class="comment">
</span> ANORM = ABS( S( 1, 1 ) )
IF( N.GT.1 )
$ ANORM = ANORM + ABS( S( 2, 1 ) )
BNORM = ABS( P( 1, 1 ) )
WORK( 1 ) = ZERO
WORK( N+1 ) = ZERO
<span class="comment">*</span><span class="comment">
</span> DO 50 J = 2, N
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 )
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?