sgeev.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 448 行 · 第 1/3 页
HTML
448 行
</span> EXTERNAL <a name="SGEBAK.132"></a><a href="sgebak.f.html#SGEBAK.1">SGEBAK</a>, <a name="SGEBAL.132"></a><a href="sgebal.f.html#SGEBAL.1">SGEBAL</a>, <a name="SGEHRD.132"></a><a href="sgehrd.f.html#SGEHRD.1">SGEHRD</a>, <a name="SHSEQR.132"></a><a href="shseqr.f.html#SHSEQR.1">SHSEQR</a>, <a name="SLABAD.132"></a><a href="slabad.f.html#SLABAD.1">SLABAD</a>, <a name="SLACPY.132"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>,
$ <a name="SLARTG.133"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>, <a name="SLASCL.133"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>, <a name="SORGHR.133"></a><a href="sorghr.f.html#SORGHR.1">SORGHR</a>, SROT, SSCAL, <a name="STREVC.133"></a><a href="strevc.f.html#STREVC.1">STREVC</a>,
$ <a name="XERBLA.134"></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"> .. External Functions ..
</span> LOGICAL <a name="LSAME.137"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
INTEGER <a name="ILAENV.138"></a><a href="hfy-index.html#ILAENV">ILAENV</a>, ISAMAX
REAL <a name="SLAMCH.139"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, <a name="SLANGE.139"></a><a href="slange.f.html#SLANGE.1">SLANGE</a>, <a name="SLAPY2.139"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>, SNRM2
EXTERNAL <a name="LSAME.140"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="ILAENV.140"></a><a href="hfy-index.html#ILAENV">ILAENV</a>, ISAMAX, <a name="SLAMCH.140"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, <a name="SLANGE.140"></a><a href="slange.f.html#SLANGE.1">SLANGE</a>, <a name="SLAPY2.140"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>,
$ SNRM2
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Intrinsic Functions ..
</span> INTRINSIC MAX, SQRT
<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"> Test the input arguments
</span><span class="comment">*</span><span class="comment">
</span> INFO = 0
LQUERY = ( LWORK.EQ.-1 )
WANTVL = <a name="LSAME.152"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBVL, <span class="string">'V'</span> )
WANTVR = <a name="LSAME.153"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBVR, <span class="string">'V'</span> )
IF( ( .NOT.WANTVL ) .AND. ( .NOT.<a name="LSAME.154"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBVL, <span class="string">'N'</span> ) ) ) THEN
INFO = -1
ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.<a name="LSAME.156"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBVR, <span class="string">'N'</span> ) ) ) THEN
INFO = -2
ELSE IF( N.LT.0 ) THEN
INFO = -3
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
INFO = -5
ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
INFO = -9
ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
INFO = -11
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Compute workspace
</span><span class="comment">*</span><span class="comment"> (Note: Comments in the code beginning "Workspace:" describe the
</span><span class="comment">*</span><span class="comment"> minimal amount of workspace needed at that point in the code,
</span><span class="comment">*</span><span class="comment"> as well as the preferred amount for good performance.
</span><span class="comment">*</span><span class="comment"> NB refers to the optimal block size for the immediately
</span><span class="comment">*</span><span class="comment"> following subroutine, as returned by <a name="ILAENV.173"></a><a href="hfy-index.html#ILAENV">ILAENV</a>.
</span><span class="comment">*</span><span class="comment"> HSWORK refers to the workspace preferred by <a name="SHSEQR.174"></a><a href="shseqr.f.html#SHSEQR.1">SHSEQR</a>, as
</span><span class="comment">*</span><span class="comment"> calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
</span><span class="comment">*</span><span class="comment"> the worst case.)
</span><span class="comment">*</span><span class="comment">
</span> IF( INFO.EQ.0 ) THEN
IF( N.EQ.0 ) THEN
MINWRK = 1
MAXWRK = 1
ELSE
MAXWRK = 2*N + N*<a name="ILAENV.183"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="SGEHRD.183"></a><a href="sgehrd.f.html#SGEHRD.1">SGEHRD</a>'</span>, <span class="string">' '</span>, N, 1, N, 0 )
IF( WANTVL ) THEN
MINWRK = 4*N
MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*<a name="ILAENV.186"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1,
$ <span class="string">'<a name="SORGHR.187"></a><a href="sorghr.f.html#SORGHR.1">SORGHR</a>'</span>, <span class="string">' '</span>, N, 1, N, -1 ) )
CALL <a name="SHSEQR.188"></a><a href="shseqr.f.html#SHSEQR.1">SHSEQR</a>( <span class="string">'S'</span>, <span class="string">'V'</span>, N, 1, N, A, LDA, WR, WI, VL, LDVL,
$ WORK, -1, INFO )
HSWORK = WORK( 1 )
MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
MAXWRK = MAX( MAXWRK, 4*N )
ELSE IF( WANTVR ) THEN
MINWRK = 4*N
MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*<a name="ILAENV.195"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1,
$ <span class="string">'<a name="SORGHR.196"></a><a href="sorghr.f.html#SORGHR.1">SORGHR</a>'</span>, <span class="string">' '</span>, N, 1, N, -1 ) )
CALL <a name="SHSEQR.197"></a><a href="shseqr.f.html#SHSEQR.1">SHSEQR</a>( <span class="string">'S'</span>, <span class="string">'V'</span>, N, 1, N, A, LDA, WR, WI, VR, LDVR,
$ WORK, -1, INFO )
HSWORK = WORK( 1 )
MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
MAXWRK = MAX( MAXWRK, 4*N )
ELSE
MINWRK = 3*N
CALL <a name="SHSEQR.204"></a><a href="shseqr.f.html#SHSEQR.1">SHSEQR</a>( <span class="string">'E'</span>, <span class="string">'N'</span>, N, 1, N, A, LDA, WR, WI, VR, LDVR,
$ WORK, -1, INFO )
HSWORK = WORK( 1 )
MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
END IF
MAXWRK = MAX( MAXWRK, MINWRK )
END IF
WORK( 1 ) = MAXWRK
<span class="comment">*</span><span class="comment">
</span> IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
INFO = -13
END IF
END IF
<span class="comment">*</span><span class="comment">
</span> IF( INFO.NE.0 ) THEN
CALL <a name="XERBLA.219"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="SGEEV.219"></a><a href="sgeev.f.html#SGEEV.1">SGEEV</a> '</span>, -INFO )
RETURN
ELSE IF( LQUERY ) THEN
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> IF( N.EQ.0 )
$ RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Get machine constants
</span><span class="comment">*</span><span class="comment">
</span> EPS = <a name="SLAMCH.232"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'P'</span> )
SMLNUM = <a name="SLAMCH.233"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'S'</span> )
BIGNUM = ONE / SMLNUM
CALL <a name="SLABAD.235"></a><a href="slabad.f.html#SLABAD.1">SLABAD</a>( SMLNUM, BIGNUM )
SMLNUM = SQRT( SMLNUM ) / EPS
BIGNUM = ONE / SMLNUM
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Scale A if max element outside range [SMLNUM,BIGNUM]
</span><span class="comment">*</span><span class="comment">
</span> ANRM = <a name="SLANGE.241"></a><a href="slange.f.html#SLANGE.1">SLANGE</a>( <span class="string">'M'</span>, N, N, A, LDA, DUM )
SCALEA = .FALSE.
IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
SCALEA = .TRUE.
CSCALE = SMLNUM
ELSE IF( ANRM.GT.BIGNUM ) THEN
SCALEA = .TRUE.
CSCALE = BIGNUM
END IF
IF( SCALEA )
$ CALL <a name="SLASCL.251"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Balance the matrix
</span><span class="comment">*</span><span class="comment"> (Workspace: need N)
</span><span class="comment">*</span><span class="comment">
</span> IBAL = 1
CALL <a name="SGEBAL.257"></a><a href="sgebal.f.html#SGEBAL.1">SGEBAL</a>( <span class="string">'B'</span>, N, A, LDA, ILO, IHI, WORK( IBAL ), IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Reduce to upper Hessenberg form
</span><span class="comment">*</span><span class="comment"> (Workspace: need 3*N, prefer 2*N+N*NB)
</span><span class="comment">*</span><span class="comment">
</span> ITAU = IBAL + N
IWRK = ITAU + N
CALL <a name="SGEHRD.264"></a><a href="sgehrd.f.html#SGEHRD.1">SGEHRD</a>( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
$ LWORK-IWRK+1, IERR )
<span class="comment">*</span><span class="comment">
</span> IF( WANTVL ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Want left eigenvectors
</span><span class="comment">*</span><span class="comment"> Copy Householder vectors to VL
</span><span class="comment">*</span><span class="comment">
</span> SIDE = <span class="string">'L'</span>
CALL <a name="SLACPY.273"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'L'</span>, N, N, A, LDA, VL, LDVL )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Generate orthogonal matrix in VL
</span><span class="comment">*</span><span class="comment"> (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="SORGHR.278"></a><a href="sorghr.f.html#SORGHR.1">SORGHR</a>( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
$ LWORK-IWRK+1, IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Perform QR iteration, accumulating Schur vectors in VL
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?