sggevx.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 741 行 · 第 1/4 页
HTML
741 行
END IF
GO TO 130
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Compute Eigenvectors and estimate condition numbers if desired
</span><span class="comment">*</span><span class="comment"> (Workspace: <a name="STGEVC.545"></a><a href="stgevc.f.html#STGEVC.1">STGEVC</a>: need 6*N
</span><span class="comment">*</span><span class="comment"> <a name="STGSNA.546"></a><a href="stgsna.f.html#STGSNA.1">STGSNA</a>: need 2*N*(N+2)+16 if SENSE = 'V' or 'B',
</span><span class="comment">*</span><span class="comment"> need N otherwise )
</span><span class="comment">*</span><span class="comment">
</span> IF( ILV .OR. .NOT.WANTSN ) THEN
IF( ILV ) THEN
IF( ILVL ) THEN
IF( ILVR ) THEN
CHTEMP = <span class="string">'B'</span>
ELSE
CHTEMP = <span class="string">'L'</span>
END IF
ELSE
CHTEMP = <span class="string">'R'</span>
END IF
<span class="comment">*</span><span class="comment">
</span> CALL <a name="STGEVC.561"></a><a href="stgevc.f.html#STGEVC.1">STGEVC</a>( CHTEMP, <span class="string">'B'</span>, LDUMMA, N, A, LDA, B, LDB, VL,
$ LDVL, VR, LDVR, N, IN, WORK, IERR )
IF( IERR.NE.0 ) THEN
INFO = N + 2
GO TO 130
END IF
END IF
<span class="comment">*</span><span class="comment">
</span> IF( .NOT.WANTSN ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> compute eigenvectors (<a name="STGEVC.571"></a><a href="stgevc.f.html#STGEVC.1">STGEVC</a>) and estimate condition
</span><span class="comment">*</span><span class="comment"> numbers (<a name="STGSNA.572"></a><a href="stgsna.f.html#STGSNA.1">STGSNA</a>). Note that the definition of the condition
</span><span class="comment">*</span><span class="comment"> number is not invariant under transformation (u,v) to
</span><span class="comment">*</span><span class="comment"> (Q*u, Z*v), where (u,v) are eigenvectors of the generalized
</span><span class="comment">*</span><span class="comment"> Schur form (S,T), Q and Z are orthogonal matrices. In order
</span><span class="comment">*</span><span class="comment"> to avoid using extra 2*N*N workspace, we have to recalculate
</span><span class="comment">*</span><span class="comment"> eigenvectors and estimate one condition numbers at a time.
</span><span class="comment">*</span><span class="comment">
</span> PAIR = .FALSE.
DO 20 I = 1, N
<span class="comment">*</span><span class="comment">
</span> IF( PAIR ) THEN
PAIR = .FALSE.
GO TO 20
END IF
MM = 1
IF( I.LT.N ) THEN
IF( A( I+1, I ).NE.ZERO ) THEN
PAIR = .TRUE.
MM = 2
END IF
END IF
<span class="comment">*</span><span class="comment">
</span> DO 10 J = 1, N
BWORK( J ) = .FALSE.
10 CONTINUE
IF( MM.EQ.1 ) THEN
BWORK( I ) = .TRUE.
ELSE IF( MM.EQ.2 ) THEN
BWORK( I ) = .TRUE.
BWORK( I+1 ) = .TRUE.
END IF
<span class="comment">*</span><span class="comment">
</span> IWRK = MM*N + 1
IWRK1 = IWRK + MM*N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Compute a pair of left and right eigenvectors.
</span><span class="comment">*</span><span class="comment"> (compute workspace: need up to 4*N + 6*N)
</span><span class="comment">*</span><span class="comment">
</span> IF( WANTSE .OR. WANTSB ) THEN
CALL <a name="STGEVC.611"></a><a href="stgevc.f.html#STGEVC.1">STGEVC</a>( <span class="string">'B'</span>, <span class="string">'S'</span>, BWORK, N, A, LDA, B, LDB,
$ WORK( 1 ), N, WORK( IWRK ), N, MM, M,
$ WORK( IWRK1 ), IERR )
IF( IERR.NE.0 ) THEN
INFO = N + 2
GO TO 130
END IF
END IF
<span class="comment">*</span><span class="comment">
</span> CALL <a name="STGSNA.620"></a><a href="stgsna.f.html#STGSNA.1">STGSNA</a>( SENSE, <span class="string">'S'</span>, BWORK, N, A, LDA, B, LDB,
$ WORK( 1 ), N, WORK( IWRK ), N, RCONDE( I ),
$ RCONDV( I ), MM, M, WORK( IWRK1 ),
$ LWORK-IWRK1+1, IWORK, IERR )
<span class="comment">*</span><span class="comment">
</span> 20 CONTINUE
END IF
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Undo balancing on VL and VR and normalization
</span><span class="comment">*</span><span class="comment"> (Workspace: none needed)
</span><span class="comment">*</span><span class="comment">
</span> IF( ILVL ) THEN
CALL <a name="SGGBAK.633"></a><a href="sggbak.f.html#SGGBAK.1">SGGBAK</a>( BALANC, <span class="string">'L'</span>, N, ILO, IHI, LSCALE, RSCALE, N, VL,
$ LDVL, IERR )
<span class="comment">*</span><span class="comment">
</span> DO 70 JC = 1, N
IF( ALPHAI( JC ).LT.ZERO )
$ GO TO 70
TEMP = ZERO
IF( ALPHAI( JC ).EQ.ZERO ) THEN
DO 30 JR = 1, N
TEMP = MAX( TEMP, ABS( VL( JR, JC ) ) )
30 CONTINUE
ELSE
DO 40 JR = 1, N
TEMP = MAX( TEMP, ABS( VL( JR, JC ) )+
$ ABS( VL( JR, JC+1 ) ) )
40 CONTINUE
END IF
IF( TEMP.LT.SMLNUM )
$ GO TO 70
TEMP = ONE / TEMP
IF( ALPHAI( JC ).EQ.ZERO ) THEN
DO 50 JR = 1, N
VL( JR, JC ) = VL( JR, JC )*TEMP
50 CONTINUE
ELSE
DO 60 JR = 1, N
VL( JR, JC ) = VL( JR, JC )*TEMP
VL( JR, JC+1 ) = VL( JR, JC+1 )*TEMP
60 CONTINUE
END IF
70 CONTINUE
END IF
IF( ILVR ) THEN
CALL <a name="SGGBAK.666"></a><a href="sggbak.f.html#SGGBAK.1">SGGBAK</a>( BALANC, <span class="string">'R'</span>, N, ILO, IHI, LSCALE, RSCALE, N, VR,
$ LDVR, IERR )
DO 120 JC = 1, N
IF( ALPHAI( JC ).LT.ZERO )
$ GO TO 120
TEMP = ZERO
IF( ALPHAI( JC ).EQ.ZERO ) THEN
DO 80 JR = 1, N
TEMP = MAX( TEMP, ABS( VR( JR, JC ) ) )
80 CONTINUE
ELSE
DO 90 JR = 1, N
TEMP = MAX( TEMP, ABS( VR( JR, JC ) )+
$ ABS( VR( JR, JC+1 ) ) )
90 CONTINUE
END IF
IF( TEMP.LT.SMLNUM )
$ GO TO 120
TEMP = ONE / TEMP
IF( ALPHAI( JC ).EQ.ZERO ) THEN
DO 100 JR = 1, N
VR( JR, JC ) = VR( JR, JC )*TEMP
100 CONTINUE
ELSE
DO 110 JR = 1, N
VR( JR, JC ) = VR( JR, JC )*TEMP
VR( JR, JC+1 ) = VR( JR, JC+1 )*TEMP
110 CONTINUE
END IF
120 CONTINUE
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Undo scaling if necessary
</span><span class="comment">*</span><span class="comment">
</span> IF( ILASCL ) THEN
CALL <a name="SLASCL.701"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
CALL <a name="SLASCL.702"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
END IF
<span class="comment">*</span><span class="comment">
</span> IF( ILBSCL ) THEN
CALL <a name="SLASCL.706"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
END IF
<span class="comment">*</span><span class="comment">
</span> 130 CONTINUE
WORK( 1 ) = MAXWRK
<span class="comment">*</span><span class="comment">
</span> RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> End of <a name="SGGEVX.714"></a><a href="sggevx.f.html#SGGEVX.1">SGGEVX</a>
</span><span class="comment">*</span><span class="comment">
</span> END
</pre>
</body>
</html>
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?