dgegs.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 463 行 · 第 1/3 页
HTML
463 行
</span> ILEFT = 1
IRIGHT = N + 1
IWORK = IRIGHT + N
CALL <a name="DGGBAL.295"></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 10
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Reduce B to triangular form, and initialize VSL and/or VSR
</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
ICOLS = N + 1 - ILO
ITAU = IWORK
IWORK = ITAU + IROWS
CALL <a name="DGEQRF.310"></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 10
END IF
<span class="comment">*</span><span class="comment">
</span> CALL <a name="DORMQR.319"></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 10
END IF
<span class="comment">*</span><span class="comment">
</span> IF( ILVSL ) THEN
CALL <a name="DLASET.330"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'Full'</span>, N, N, ZERO, ONE, VSL, LDVSL )
CALL <a name="DLACPY.331"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'L'</span>, IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
$ VSL( ILO+1, ILO ), LDVSL )
CALL <a name="DORGQR.333"></a><a href="dorgqr.f.html#DORGQR.1">DORGQR</a>( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
$ 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 10
END IF
END IF
<span class="comment">*</span><span class="comment">
</span> IF( ILVSR )
$ CALL <a name="DLASET.345"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'Full'</span>, N, N, ZERO, ONE, VSR, LDVSR )
<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> CALL <a name="DGGHRD.349"></a><a href="dgghrd.f.html#DGGHRD.1">DGGHRD</a>( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
$ LDVSL, VSR, LDVSR, IINFO )
IF( IINFO.NE.0 ) THEN
INFO = N + 5
GO TO 10
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Perform QZ algorithm, computing Schur vectors if desired
</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
CALL <a name="DHGEQZ.361"></a><a href="dhgeqz.f.html#DHGEQZ.1">DHGEQZ</a>( <span class="string">'S'</span>, JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
$ ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
$ 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 10
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Apply permutation to VSL and VSR
</span><span class="comment">*</span><span class="comment">
</span> IF( ILVSL ) THEN
CALL <a name="DGGBAK.380"></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, VSL, LDVSL, IINFO )
IF( IINFO.NE.0 ) THEN
INFO = N + 7
GO TO 10
END IF
END IF
IF( ILVSR ) THEN
CALL <a name="DGGBAK.388"></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, VSR, LDVSR, IINFO )
IF( IINFO.NE.0 ) THEN
INFO = N + 8
GO TO 10
END IF
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Undo scaling
</span><span class="comment">*</span><span class="comment">
</span> IF( ILASCL ) THEN
CALL <a name="DLASCL.399"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'H'</span>, -1, -1, ANRMTO, ANRM, N, N, A, LDA, IINFO )
IF( IINFO.NE.0 ) THEN
INFO = N + 9
RETURN
END IF
CALL <a name="DLASCL.404"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, -1, -1, ANRMTO, ANRM, N, 1, ALPHAR, N,
$ IINFO )
IF( IINFO.NE.0 ) THEN
INFO = N + 9
RETURN
END IF
CALL <a name="DLASCL.410"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, -1, -1, ANRMTO, ANRM, N, 1, ALPHAI, N,
$ IINFO )
IF( IINFO.NE.0 ) THEN
INFO = N + 9
RETURN
END IF
END IF
<span class="comment">*</span><span class="comment">
</span> IF( ILBSCL ) THEN
CALL <a name="DLASCL.419"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'U'</span>, -1, -1, BNRMTO, BNRM, N, N, B, LDB, IINFO )
IF( IINFO.NE.0 ) THEN
INFO = N + 9
RETURN
END IF
CALL <a name="DLASCL.424"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, -1, -1, BNRMTO, BNRM, N, 1, BETA, N, IINFO )
IF( IINFO.NE.0 ) THEN
INFO = N + 9
RETURN
END IF
END IF
<span class="comment">*</span><span class="comment">
</span> 10 CONTINUE
WORK( 1 ) = LWKOPT
<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="DGEGS.436"></a><a href="dgegs.f.html#DGEGS.1">DGEGS</a>
</span><span class="comment">*</span><span class="comment">
</span> END
</pre>
</body>
</html>
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?