dgegv.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 690 行 · 第 1/4 页
HTML
690 行
SAFMIN = <a name="DLAMCH.328"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'S'</span> )
SAFMIN = SAFMIN + SAFMIN
SAFMAX = ONE / SAFMIN
ONEPLS = ONE + ( 4*EPS )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Scale A
</span><span class="comment">*</span><span class="comment">
</span> ANRM = <a name="DLANGE.335"></a><a href="dlange.f.html#DLANGE.1">DLANGE</a>( <span class="string">'M'</span>, N, N, A, LDA, WORK )
ANRM1 = ANRM
ANRM2 = ONE
IF( ANRM.LT.ONE ) THEN
IF( SAFMAX*ANRM.LT.ONE ) THEN
ANRM1 = SAFMIN
ANRM2 = SAFMAX*ANRM
END IF
END IF
<span class="comment">*</span><span class="comment">
</span> IF( ANRM.GT.ZERO ) THEN
CALL <a name="DLASCL.346"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, -1, -1, ANRM, ONE, N, N, A, LDA, IINFO )
IF( IINFO.NE.0 ) THEN
INFO = N + 10
RETURN
END IF
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Scale B
</span><span class="comment">*</span><span class="comment">
</span> BNRM = <a name="DLANGE.355"></a><a href="dlange.f.html#DLANGE.1">DLANGE</a>( <span class="string">'M'</span>, N, N, B, LDB, WORK )
BNRM1 = BNRM
BNRM2 = ONE
IF( BNRM.LT.ONE ) THEN
IF( SAFMAX*BNRM.LT.ONE ) THEN
BNRM1 = SAFMIN
BNRM2 = SAFMAX*BNRM
END IF
END IF
<span class="comment">*</span><span class="comment">
</span> IF( BNRM.GT.ZERO ) THEN
CALL <a name="DLASCL.366"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, -1, -1, BNRM, ONE, N, N, B, LDB, IINFO )
IF( IINFO.NE.0 ) THEN
INFO = N + 10
RETURN
END IF
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Permute the matrix to make it more nearly triangular
</span><span class="comment">*</span><span class="comment"> Workspace layout: (8*N words -- "work" requires 6*N words)
</span><span class="comment">*</span><span class="comment"> left_permutation, right_permutation, work...
</span><span class="comment">*</span><span class="comment">
</span> ILEFT = 1
IRIGHT = N + 1
IWORK = IRIGHT + N
CALL <a name="DGGBAL.380"></a><a href="dggbal.f.html#DGGBAL.1">DGGBAL</a>( <span class="string">'P'</span>, N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
$ WORK( IRIGHT ), WORK( IWORK ), IINFO )
IF( IINFO.NE.0 ) THEN
INFO = N + 1
GO TO 120
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Reduce B to triangular form, and initialize VL and/or VR
</span><span class="comment">*</span><span class="comment"> Workspace layout: ("work..." must have at least N words)
</span><span class="comment">*</span><span class="comment"> left_permutation, right_permutation, tau, work...
</span><span class="comment">*</span><span class="comment">
</span> IROWS = IHI + 1 - ILO
IF( ILV ) THEN
ICOLS = N + 1 - ILO
ELSE
ICOLS = IROWS
END IF
ITAU = IWORK
IWORK = ITAU + IROWS
CALL <a name="DGEQRF.399"></a><a href="dgeqrf.f.html#DGEQRF.1">DGEQRF</a>( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
$ WORK( IWORK ), LWORK+1-IWORK, IINFO )
IF( IINFO.GE.0 )
$ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
IF( IINFO.NE.0 ) THEN
INFO = N + 2
GO TO 120
END IF
<span class="comment">*</span><span class="comment">
</span> CALL <a name="DORMQR.408"></a><a href="dormqr.f.html#DORMQR.1">DORMQR</a>( <span class="string">'L'</span>, <span class="string">'T'</span>, IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
$ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWORK ),
$ LWORK+1-IWORK, IINFO )
IF( IINFO.GE.0 )
$ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
IF( IINFO.NE.0 ) THEN
INFO = N + 3
GO TO 120
END IF
<span class="comment">*</span><span class="comment">
</span> IF( ILVL ) THEN
CALL <a name="DLASET.419"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'Full'</span>, N, N, ZERO, ONE, VL, LDVL )
CALL <a name="DLACPY.420"></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 )
CALL <a name="DORGQR.422"></a><a href="dorgqr.f.html#DORGQR.1">DORGQR</a>( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
$ WORK( ITAU ), WORK( IWORK ), LWORK+1-IWORK,
$ IINFO )
IF( IINFO.GE.0 )
$ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
IF( IINFO.NE.0 ) THEN
INFO = N + 4
GO TO 120
END IF
END IF
<span class="comment">*</span><span class="comment">
</span> IF( ILVR )
$ CALL <a name="DLASET.434"></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">
</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.442"></a><a href="dgghrd.f.html#DGGHRD.1">DGGHRD</a>( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
$ LDVL, VR, LDVR, IINFO )
ELSE
CALL <a name="DGGHRD.445"></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, IINFO )
END IF
IF( IINFO.NE.0 ) THEN
INFO = N + 5
GO TO 120
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Perform QZ algorithm
</span><span class="comment">*</span><span class="comment"> Workspace layout: ("work..." must have at least 1 word)
</span><span class="comment">*</span><span class="comment"> left_permutation, right_permutation, work...
</span><span class="comment">*</span><span class="comment">
</span> IWORK = ITAU
IF( ILV ) THEN
CHTEMP = <span class="string">'S'</span>
ELSE
CHTEMP = <span class="string">'E'</span>
END IF
CALL <a name="DHGEQZ.463"></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( IWORK ), LWORK+1-IWORK, IINFO )
IF( IINFO.GE.0 )
$ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
IF( IINFO.NE.0 ) THEN
IF( IINFO.GT.0 .AND. IINFO.LE.N ) THEN
INFO = IINFO
ELSE IF( IINFO.GT.N .AND. IINFO.LE.2*N ) THEN
INFO = IINFO - N
ELSE
INFO = N + 6
END IF
GO TO 120
END IF
<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"> Compute Eigenvectors (<a name="DTGEVC.481"></a><a href="dtgevc.f.html#DTGEVC.1">DTGEVC</a> requires 6*N words of workspace)
</span><span class="comment">*</span><span class="comment">
</span> 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="DTGEVC.493"></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( IWORK ), IINFO )
IF( IINFO.NE.0 ) THEN
INFO = N + 7
GO TO 120
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Undo balancing on VL and VR, rescale
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?