dggev.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 514 行 · 第 1/3 页
HTML
514 行
CALL <a name="DLASET.326"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'Full'</span>, N, N, ZERO, ONE, VL, LDVL )
IF( IROWS.GT.1 ) THEN
CALL <a name="DLACPY.328"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'L'</span>, IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
$ VL( ILO+1, ILO ), LDVL )
END IF
CALL <a name="DORGQR.331"></a><a href="dorgqr.f.html#DORGQR.1">DORGQR</a>( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
$ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Initialize VR
</span><span class="comment">*</span><span class="comment">
</span> IF( ILVR )
$ CALL <a name="DLASET.338"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'Full'</span>, N, N, ZERO, ONE, VR, LDVR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Reduce to generalized Hessenberg form
</span><span class="comment">*</span><span class="comment"> (Workspace: none needed)
</span><span class="comment">*</span><span class="comment">
</span> IF( ILV ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Eigenvectors requested -- work on whole matrix.
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="DGGHRD.347"></a><a href="dgghrd.f.html#DGGHRD.1">DGGHRD</a>( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
$ LDVL, VR, LDVR, IERR )
ELSE
CALL <a name="DGGHRD.350"></a><a href="dgghrd.f.html#DGGHRD.1">DGGHRD</a>( <span class="string">'N'</span>, <span class="string">'N'</span>, IROWS, 1, IROWS, A( ILO, ILO ), LDA,
$ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR )
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Perform QZ algorithm (Compute eigenvalues, and optionally, the
</span><span class="comment">*</span><span class="comment"> Schur forms and Schur vectors)
</span><span class="comment">*</span><span class="comment"> (Workspace: need N)
</span><span class="comment">*</span><span class="comment">
</span> IWRK = ITAU
IF( ILV ) THEN
CHTEMP = <span class="string">'S'</span>
ELSE
CHTEMP = <span class="string">'E'</span>
END IF
CALL <a name="DHGEQZ.364"></a><a href="dhgeqz.f.html#DHGEQZ.1">DHGEQZ</a>( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
$ ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR,
$ WORK( IWRK ), LWORK+1-IWRK, IERR )
IF( IERR.NE.0 ) THEN
IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
INFO = IERR
ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
INFO = IERR - N
ELSE
INFO = N + 1
END IF
GO TO 110
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Compute Eigenvectors
</span><span class="comment">*</span><span class="comment"> (Workspace: need 6*N)
</span><span class="comment">*</span><span class="comment">
</span> 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
CALL <a name="DTGEVC.391"></a><a href="dtgevc.f.html#DTGEVC.1">DTGEVC</a>( CHTEMP, <span class="string">'B'</span>, LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
$ VR, LDVR, N, IN, WORK( IWRK ), IERR )
IF( IERR.NE.0 ) THEN
INFO = N + 2
GO TO 110
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="DGGBAK.402"></a><a href="dggbak.f.html#DGGBAK.1">DGGBAK</a>( <span class="string">'P'</span>, <span class="string">'L'</span>, N, ILO, IHI, WORK( ILEFT ),
$ WORK( IRIGHT ), N, VL, LDVL, IERR )
DO 50 JC = 1, N
IF( ALPHAI( JC ).LT.ZERO )
$ GO TO 50
TEMP = ZERO
IF( ALPHAI( JC ).EQ.ZERO ) THEN
DO 10 JR = 1, N
TEMP = MAX( TEMP, ABS( VL( JR, JC ) ) )
10 CONTINUE
ELSE
DO 20 JR = 1, N
TEMP = MAX( TEMP, ABS( VL( JR, JC ) )+
$ ABS( VL( JR, JC+1 ) ) )
20 CONTINUE
END IF
IF( TEMP.LT.SMLNUM )
$ GO TO 50
TEMP = ONE / TEMP
IF( ALPHAI( JC ).EQ.ZERO ) THEN
DO 30 JR = 1, N
VL( JR, JC ) = VL( JR, JC )*TEMP
30 CONTINUE
ELSE
DO 40 JR = 1, N
VL( JR, JC ) = VL( JR, JC )*TEMP
VL( JR, JC+1 ) = VL( JR, JC+1 )*TEMP
40 CONTINUE
END IF
50 CONTINUE
END IF
IF( ILVR ) THEN
CALL <a name="DGGBAK.434"></a><a href="dggbak.f.html#DGGBAK.1">DGGBAK</a>( <span class="string">'P'</span>, <span class="string">'R'</span>, N, ILO, IHI, WORK( ILEFT ),
$ WORK( IRIGHT ), N, VR, LDVR, IERR )
DO 100 JC = 1, N
IF( ALPHAI( JC ).LT.ZERO )
$ GO TO 100
TEMP = ZERO
IF( ALPHAI( JC ).EQ.ZERO ) THEN
DO 60 JR = 1, N
TEMP = MAX( TEMP, ABS( VR( JR, JC ) ) )
60 CONTINUE
ELSE
DO 70 JR = 1, N
TEMP = MAX( TEMP, ABS( VR( JR, JC ) )+
$ ABS( VR( JR, JC+1 ) ) )
70 CONTINUE
END IF
IF( TEMP.LT.SMLNUM )
$ GO TO 100
TEMP = ONE / TEMP
IF( ALPHAI( JC ).EQ.ZERO ) THEN
DO 80 JR = 1, N
VR( JR, JC ) = VR( JR, JC )*TEMP
80 CONTINUE
ELSE
DO 90 JR = 1, N
VR( JR, JC ) = VR( JR, JC )*TEMP
VR( JR, JC+1 ) = VR( JR, JC+1 )*TEMP
90 CONTINUE
END IF
100 CONTINUE
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> End of eigenvector calculation
</span><span class="comment">*</span><span class="comment">
</span> 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="DLASCL.473"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
CALL <a name="DLASCL.474"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</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="DLASCL.478"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</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> 110 CONTINUE
<span class="comment">*</span><span class="comment">
</span> 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="DGGEV.487"></a><a href="dggev.f.html#DGGEV.1">DGGEV</a>
</span><span class="comment">*</span><span class="comment">
</span> END
</pre>
</body>
</html>
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?