zgeevx.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 557 行 · 第 1/3 页
HTML
557 行
CALL <a name="DLASCL.354"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, CSCALE, ANRM, 1, 1, DUM, 1, IERR )
ABNRM = DUM( 1 )
END IF
<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"> (CWorkspace: need 2*N, prefer N+N*NB)
</span><span class="comment">*</span><span class="comment"> (RWorkspace: none)
</span><span class="comment">*</span><span class="comment">
</span> ITAU = 1
IWRK = ITAU + N
CALL <a name="ZGEHRD.364"></a><a href="zgehrd.f.html#ZGEHRD.1">ZGEHRD</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="ZLACPY.373"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</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 unitary matrix in VL
</span><span class="comment">*</span><span class="comment"> (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
</span><span class="comment">*</span><span class="comment"> (RWorkspace: none)
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="ZUNGHR.379"></a><a href="zunghr.f.html#ZUNGHR.1">ZUNGHR</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
</span><span class="comment">*</span><span class="comment"> (CWorkspace: need 1, prefer HSWORK (see comments) )
</span><span class="comment">*</span><span class="comment"> (RWorkspace: none)
</span><span class="comment">*</span><span class="comment">
</span> IWRK = ITAU
CALL <a name="ZHSEQR.387"></a><a href="zhseqr.f.html#ZHSEQR.1">ZHSEQR</a>( <span class="string">'S'</span>, <span class="string">'V'</span>, N, ILO, IHI, A, LDA, W, 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="ZLACPY.396"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</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="ZLACPY.405"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</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 unitary matrix in VR
</span><span class="comment">*</span><span class="comment"> (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
</span><span class="comment">*</span><span class="comment"> (RWorkspace: none)
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="ZUNGHR.411"></a><a href="zunghr.f.html#ZUNGHR.1">ZUNGHR</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"> (CWorkspace: need 1, prefer HSWORK (see comments) )
</span><span class="comment">*</span><span class="comment"> (RWorkspace: none)
</span><span class="comment">*</span><span class="comment">
</span> IWRK = ITAU
CALL <a name="ZHSEQR.419"></a><a href="zhseqr.f.html#ZHSEQR.1">ZHSEQR</a>( <span class="string">'S'</span>, <span class="string">'V'</span>, N, ILO, IHI, A, LDA, W, 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"> If condition numbers desired, compute Schur form
</span><span class="comment">*</span><span class="comment">
</span> IF( WNTSNN ) THEN
JOB = <span class="string">'E'</span>
ELSE
JOB = <span class="string">'S'</span>
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> (CWorkspace: need 1, prefer HSWORK (see comments) )
</span><span class="comment">*</span><span class="comment"> (RWorkspace: none)
</span><span class="comment">*</span><span class="comment">
</span> IWRK = ITAU
CALL <a name="ZHSEQR.437"></a><a href="zhseqr.f.html#ZHSEQR.1">ZHSEQR</a>( JOB, <span class="string">'N'</span>, N, ILO, IHI, A, LDA, W, 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="ZHSEQR.441"></a><a href="zhseqr.f.html#ZHSEQR.1">ZHSEQR</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"> (CWorkspace: need 2*N)
</span><span class="comment">*</span><span class="comment"> (RWorkspace: need N)
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="ZTREVC.452"></a><a href="ztrevc.f.html#ZTREVC.1">ZTREVC</a>( SIDE, <span class="string">'B'</span>, SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
$ N, NOUT, WORK( IWRK ), RWORK, IERR )
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Compute condition numbers if desired
</span><span class="comment">*</span><span class="comment"> (CWorkspace: need N*N+2*N unless SENSE = 'E')
</span><span class="comment">*</span><span class="comment"> (RWorkspace: need 2*N unless SENSE = 'E')
</span><span class="comment">*</span><span class="comment">
</span> IF( .NOT.WNTSNN ) THEN
CALL <a name="ZTRSNA.461"></a><a href="ztrsna.f.html#ZTRSNA.1">ZTRSNA</a>( SENSE, <span class="string">'A'</span>, SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
$ RCONDE, RCONDV, N, NOUT, WORK( IWRK ), N, RWORK,
$ ICOND )
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">
</span> CALL <a name="ZGEBAK.470"></a><a href="zgebak.f.html#ZGEBAK.1">ZGEBAK</a>( BALANC, <span class="string">'L'</span>, N, ILO, IHI, SCALE, 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
SCL = ONE / DZNRM2( N, VL( 1, I ), 1 )
CALL ZDSCAL( N, SCL, VL( 1, I ), 1 )
DO 10 K = 1, N
RWORK( K ) = DBLE( VL( K, I ) )**2 +
$ DIMAG( VL( K, I ) )**2
10 CONTINUE
K = IDAMAX( N, RWORK, 1 )
TMP = DCONJG( VL( K, I ) ) / SQRT( RWORK( K ) )
CALL ZSCAL( N, TMP, VL( 1, I ), 1 )
VL( K, I ) = DCMPLX( DBLE( VL( K, I ) ), ZERO )
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">
</span> CALL <a name="ZGEBAK.493"></a><a href="zgebak.f.html#ZGEBAK.1">ZGEBAK</a>( BALANC, <span class="string">'R'</span>, N, ILO, IHI, SCALE, 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
SCL = ONE / DZNRM2( N, VR( 1, I ), 1 )
CALL ZDSCAL( N, SCL, VR( 1, I ), 1 )
DO 30 K = 1, N
RWORK( K ) = DBLE( VR( K, I ) )**2 +
$ DIMAG( VR( K, I ) )**2
30 CONTINUE
K = IDAMAX( N, RWORK, 1 )
TMP = DCONJG( VR( K, I ) ) / SQRT( RWORK( K ) )
CALL ZSCAL( N, TMP, VR( 1, I ), 1 )
VR( K, I ) = DCMPLX( DBLE( VR( K, I ) ), ZERO )
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="ZLASCL.516"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>( <span class="string">'G'</span>, 0, 0, CSCALE, ANRM, N-INFO, 1, W( INFO+1 ),
$ MAX( N-INFO, 1 ), IERR )
IF( INFO.EQ.0 ) THEN
IF( ( WNTSNV .OR. WNTSNB ) .AND. ICOND.EQ.0 )
$ CALL <a name="DLASCL.520"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, CSCALE, ANRM, N, 1, RCONDV, N,
$ IERR )
ELSE
CALL <a name="ZLASCL.523"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>( <span class="string">'G'</span>, 0, 0, CSCALE, ANRM, ILO-1, 1, W, 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="ZGEEVX.530"></a><a href="zgeevx.f.html#ZGEEVX.1">ZGEEVX</a>
</span><span class="comment">*</span><span class="comment">
</span> END
</pre>
</body>
</html>
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?