dgeev.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 448 行 · 第 1/3 页
HTML
448 行
</span><span class="comment">*</span><span class="comment"> (Workspace: need N+1, prefer N+HSWORK (see comments) )
</span><span class="comment">*</span><span class="comment">
</span> IWRK = ITAU
CALL <a name="DHSEQR.285"></a><a href="dhseqr.f.html#DHSEQR.1">DHSEQR</a>( <span class="string">'S'</span>, <span class="string">'V'</span>, N, ILO, IHI, A, LDA, WR, WI, VL, LDVL,
$ WORK( IWRK ), LWORK-IWRK+1, INFO )
<span class="comment">*</span><span class="comment">
</span> IF( WANTVR ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Want left and right eigenvectors
</span><span class="comment">*</span><span class="comment"> Copy Schur vectors to VR
</span><span class="comment">*</span><span class="comment">
</span> SIDE = <span class="string">'B'</span>
CALL <a name="DLACPY.294"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'F'</span>, N, N, VL, LDVL, VR, LDVR )
END IF
<span class="comment">*</span><span class="comment">
</span> ELSE IF( WANTVR ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Want right eigenvectors
</span><span class="comment">*</span><span class="comment"> Copy Householder vectors to VR
</span><span class="comment">*</span><span class="comment">
</span> SIDE = <span class="string">'R'</span>
CALL <a name="DLACPY.303"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'L'</span>, N, N, A, LDA, VR, LDVR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Generate orthogonal matrix in VR
</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="DORGHR.308"></a><a href="dorghr.f.html#DORGHR.1">DORGHR</a>( N, ILO, IHI, VR, LDVR, 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 VR
</span><span class="comment">*</span><span class="comment"> (Workspace: need N+1, prefer N+HSWORK (see comments) )
</span><span class="comment">*</span><span class="comment">
</span> IWRK = ITAU
CALL <a name="DHSEQR.315"></a><a href="dhseqr.f.html#DHSEQR.1">DHSEQR</a>( <span class="string">'S'</span>, <span class="string">'V'</span>, N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
$ WORK( IWRK ), LWORK-IWRK+1, INFO )
<span class="comment">*</span><span class="comment">
</span> ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Compute eigenvalues only
</span><span class="comment">*</span><span class="comment"> (Workspace: need N+1, prefer N+HSWORK (see comments) )
</span><span class="comment">*</span><span class="comment">
</span> IWRK = ITAU
CALL <a name="DHSEQR.324"></a><a href="dhseqr.f.html#DHSEQR.1">DHSEQR</a>( <span class="string">'E'</span>, <span class="string">'N'</span>, N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
$ WORK( IWRK ), LWORK-IWRK+1, INFO )
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> If INFO > 0 from <a name="DHSEQR.328"></a><a href="dhseqr.f.html#DHSEQR.1">DHSEQR</a>, then quit
</span><span class="comment">*</span><span class="comment">
</span> IF( INFO.GT.0 )
$ GO TO 50
<span class="comment">*</span><span class="comment">
</span> IF( WANTVL .OR. WANTVR ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Compute left and/or right eigenvectors
</span><span class="comment">*</span><span class="comment"> (Workspace: need 4*N)
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="DTREVC.338"></a><a href="dtrevc.f.html#DTREVC.1">DTREVC</a>( SIDE, <span class="string">'B'</span>, SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
$ N, NOUT, WORK( IWRK ), IERR )
END IF
<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"> Undo balancing of left eigenvectors
</span><span class="comment">*</span><span class="comment"> (Workspace: need N)
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="DGEBAK.347"></a><a href="dgebak.f.html#DGEBAK.1">DGEBAK</a>( <span class="string">'B'</span>, <span class="string">'L'</span>, N, ILO, IHI, WORK( IBAL ), N, VL, LDVL,
$ IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Normalize left eigenvectors and make largest component real
</span><span class="comment">*</span><span class="comment">
</span> DO 20 I = 1, N
IF( WI( I ).EQ.ZERO ) THEN
SCL = ONE / DNRM2( N, VL( 1, I ), 1 )
CALL DSCAL( N, SCL, VL( 1, I ), 1 )
ELSE IF( WI( I ).GT.ZERO ) THEN
SCL = ONE / <a name="DLAPY2.357"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>( DNRM2( N, VL( 1, I ), 1 ),
$ DNRM2( N, VL( 1, I+1 ), 1 ) )
CALL DSCAL( N, SCL, VL( 1, I ), 1 )
CALL DSCAL( N, SCL, VL( 1, I+1 ), 1 )
DO 10 K = 1, N
WORK( IWRK+K-1 ) = VL( K, I )**2 + VL( K, I+1 )**2
10 CONTINUE
K = IDAMAX( N, WORK( IWRK ), 1 )
CALL <a name="DLARTG.365"></a><a href="dlartg.f.html#DLARTG.1">DLARTG</a>( VL( K, I ), VL( K, I+1 ), CS, SN, R )
CALL DROT( N, VL( 1, I ), 1, VL( 1, I+1 ), 1, CS, SN )
VL( K, I+1 ) = ZERO
END IF
20 CONTINUE
END IF
<span class="comment">*</span><span class="comment">
</span> IF( WANTVR ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Undo balancing of right eigenvectors
</span><span class="comment">*</span><span class="comment"> (Workspace: need N)
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="DGEBAK.377"></a><a href="dgebak.f.html#DGEBAK.1">DGEBAK</a>( <span class="string">'B'</span>, <span class="string">'R'</span>, N, ILO, IHI, WORK( IBAL ), N, VR, LDVR,
$ IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Normalize right eigenvectors and make largest component real
</span><span class="comment">*</span><span class="comment">
</span> DO 40 I = 1, N
IF( WI( I ).EQ.ZERO ) THEN
SCL = ONE / DNRM2( N, VR( 1, I ), 1 )
CALL DSCAL( N, SCL, VR( 1, I ), 1 )
ELSE IF( WI( I ).GT.ZERO ) THEN
SCL = ONE / <a name="DLAPY2.387"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>( DNRM2( N, VR( 1, I ), 1 ),
$ DNRM2( N, VR( 1, I+1 ), 1 ) )
CALL DSCAL( N, SCL, VR( 1, I ), 1 )
CALL DSCAL( N, SCL, VR( 1, I+1 ), 1 )
DO 30 K = 1, N
WORK( IWRK+K-1 ) = VR( K, I )**2 + VR( K, I+1 )**2
30 CONTINUE
K = IDAMAX( N, WORK( IWRK ), 1 )
CALL <a name="DLARTG.395"></a><a href="dlartg.f.html#DLARTG.1">DLARTG</a>( VR( K, I ), VR( K, I+1 ), CS, SN, R )
CALL DROT( N, VR( 1, I ), 1, VR( 1, I+1 ), 1, CS, SN )
VR( K, I+1 ) = ZERO
END IF
40 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> 50 CONTINUE
IF( SCALEA ) THEN
CALL <a name="DLASCL.406"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, CSCALE, ANRM, N-INFO, 1, WR( INFO+1 ),
$ MAX( N-INFO, 1 ), IERR )
CALL <a name="DLASCL.408"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, CSCALE, ANRM, N-INFO, 1, WI( INFO+1 ),
$ MAX( N-INFO, 1 ), IERR )
IF( INFO.GT.0 ) THEN
CALL <a name="DLASCL.411"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, CSCALE, ANRM, ILO-1, 1, WR, N,
$ IERR )
CALL <a name="DLASCL.413"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N,
$ IERR )
END IF
END IF
<span class="comment">*</span><span class="comment">
</span> WORK( 1 ) = MAXWRK
RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> End of <a name="DGEEV.421"></a><a href="dgeev.f.html#DGEEV.1">DGEEV</a>
</span><span class="comment">*</span><span class="comment">
</span> END
</pre>
</body>
</html>
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?