dggevx.f.html

来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 743 行 · 第 1/4 页

HTML
743
字号
         ELSE
            INFO = N + 1
         END IF
         GO TO 130
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute Eigenvectors and estimate condition numbers if desired
</span><span class="comment">*</span><span class="comment">     (Workspace: <a name="DTGEVC.547"></a><a href="dtgevc.f.html#DTGEVC.1">DTGEVC</a>: need 6*N
</span><span class="comment">*</span><span class="comment">                 <a name="DTGSNA.548"></a><a href="dtgsna.f.html#DTGSNA.1">DTGSNA</a>: need 2*N*(N+2)+16 if SENSE = 'V' or 'B',
</span><span class="comment">*</span><span class="comment">                         need N otherwise )
</span><span class="comment">*</span><span class="comment">
</span>      IF( ILV .OR. .NOT.WANTSN ) THEN
         IF( ILV ) THEN
            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.563"></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, IERR )
            IF( IERR.NE.0 ) THEN
               INFO = N + 2
               GO TO 130
            END IF
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( .NOT.WANTSN ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           compute eigenvectors (<a name="DTGEVC.573"></a><a href="dtgevc.f.html#DTGEVC.1">DTGEVC</a>) and estimate condition
</span><span class="comment">*</span><span class="comment">           numbers (<a name="DTGSNA.574"></a><a href="dtgsna.f.html#DTGSNA.1">DTGSNA</a>). Note that the definition of the condition
</span><span class="comment">*</span><span class="comment">           number is not invariant under transformation (u,v) to
</span><span class="comment">*</span><span class="comment">           (Q*u, Z*v), where (u,v) are eigenvectors of the generalized
</span><span class="comment">*</span><span class="comment">           Schur form (S,T), Q and Z are orthogonal matrices. In order
</span><span class="comment">*</span><span class="comment">           to avoid using extra 2*N*N workspace, we have to recalculate
</span><span class="comment">*</span><span class="comment">           eigenvectors and estimate one condition numbers at a time.
</span><span class="comment">*</span><span class="comment">
</span>            PAIR = .FALSE.
            DO 20 I = 1, N
<span class="comment">*</span><span class="comment">
</span>               IF( PAIR ) THEN
                  PAIR = .FALSE.
                  GO TO 20
               END IF
               MM = 1
               IF( I.LT.N ) THEN
                  IF( A( I+1, I ).NE.ZERO ) THEN
                     PAIR = .TRUE.
                     MM = 2
                  END IF
               END IF
<span class="comment">*</span><span class="comment">
</span>               DO 10 J = 1, N
                  BWORK( J ) = .FALSE.
   10          CONTINUE
               IF( MM.EQ.1 ) THEN
                  BWORK( I ) = .TRUE.
               ELSE IF( MM.EQ.2 ) THEN
                  BWORK( I ) = .TRUE.
                  BWORK( I+1 ) = .TRUE.
               END IF
<span class="comment">*</span><span class="comment">
</span>               IWRK = MM*N + 1
               IWRK1 = IWRK + MM*N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Compute a pair of left and right eigenvectors.
</span><span class="comment">*</span><span class="comment">              (compute workspace: need up to 4*N + 6*N)
</span><span class="comment">*</span><span class="comment">
</span>               IF( WANTSE .OR. WANTSB ) THEN
                  CALL <a name="DTGEVC.613"></a><a href="dtgevc.f.html#DTGEVC.1">DTGEVC</a>( <span class="string">'B'</span>, <span class="string">'S'</span>, BWORK, N, A, LDA, B, LDB,
     $                         WORK( 1 ), N, WORK( IWRK ), N, MM, M,
     $                         WORK( IWRK1 ), IERR )
                  IF( IERR.NE.0 ) THEN
                     INFO = N + 2
                     GO TO 130
                  END IF
               END IF
<span class="comment">*</span><span class="comment">
</span>               CALL <a name="DTGSNA.622"></a><a href="dtgsna.f.html#DTGSNA.1">DTGSNA</a>( SENSE, <span class="string">'S'</span>, BWORK, N, A, LDA, B, LDB,
     $                      WORK( 1 ), N, WORK( IWRK ), N, RCONDE( I ),
     $                      RCONDV( I ), MM, M, WORK( IWRK1 ),
     $                      LWORK-IWRK1+1, IWORK, IERR )
<span class="comment">*</span><span class="comment">
</span>   20       CONTINUE
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Undo balancing on VL and VR and normalization
</span><span class="comment">*</span><span class="comment">     (Workspace: none needed)
</span><span class="comment">*</span><span class="comment">
</span>      IF( ILVL ) THEN
         CALL <a name="DGGBAK.635"></a><a href="dggbak.f.html#DGGBAK.1">DGGBAK</a>( BALANC, <span class="string">'L'</span>, N, ILO, IHI, LSCALE, RSCALE, N, VL,
     $                LDVL, IERR )
<span class="comment">*</span><span class="comment">
</span>         DO 70 JC = 1, N
            IF( ALPHAI( JC ).LT.ZERO )
     $         GO TO 70
            TEMP = ZERO
            IF( ALPHAI( JC ).EQ.ZERO ) THEN
               DO 30 JR = 1, N
                  TEMP = MAX( TEMP, ABS( VL( JR, JC ) ) )
   30          CONTINUE
            ELSE
               DO 40 JR = 1, N
                  TEMP = MAX( TEMP, ABS( VL( JR, JC ) )+
     $                   ABS( VL( JR, JC+1 ) ) )
   40          CONTINUE
            END IF
            IF( TEMP.LT.SMLNUM )
     $         GO TO 70
            TEMP = ONE / TEMP
            IF( ALPHAI( JC ).EQ.ZERO ) THEN
               DO 50 JR = 1, N
                  VL( JR, JC ) = VL( JR, JC )*TEMP
   50          CONTINUE
            ELSE
               DO 60 JR = 1, N
                  VL( JR, JC ) = VL( JR, JC )*TEMP
                  VL( JR, JC+1 ) = VL( JR, JC+1 )*TEMP
   60          CONTINUE
            END IF
   70    CONTINUE
      END IF
      IF( ILVR ) THEN
         CALL <a name="DGGBAK.668"></a><a href="dggbak.f.html#DGGBAK.1">DGGBAK</a>( BALANC, <span class="string">'R'</span>, N, ILO, IHI, LSCALE, RSCALE, N, VR,
     $                LDVR, IERR )
         DO 120 JC = 1, N
            IF( ALPHAI( JC ).LT.ZERO )
     $         GO TO 120
            TEMP = ZERO
            IF( ALPHAI( JC ).EQ.ZERO ) THEN
               DO 80 JR = 1, N
                  TEMP = MAX( TEMP, ABS( VR( JR, JC ) ) )
   80          CONTINUE
            ELSE
               DO 90 JR = 1, N
                  TEMP = MAX( TEMP, ABS( VR( JR, JC ) )+
     $                   ABS( VR( JR, JC+1 ) ) )
   90          CONTINUE
            END IF
            IF( TEMP.LT.SMLNUM )
     $         GO TO 120
            TEMP = ONE / TEMP
            IF( ALPHAI( JC ).EQ.ZERO ) THEN
               DO 100 JR = 1, N
                  VR( JR, JC ) = VR( JR, JC )*TEMP
  100          CONTINUE
            ELSE
               DO 110 JR = 1, N
                  VR( JR, JC ) = VR( JR, JC )*TEMP
                  VR( JR, JC+1 ) = VR( JR, JC+1 )*TEMP
  110          CONTINUE
            END IF
  120    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>      IF( ILASCL ) THEN
         CALL <a name="DLASCL.703"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
         CALL <a name="DLASCL.704"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( ILBSCL ) THEN
         CALL <a name="DLASCL.708"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
      END IF
<span class="comment">*</span><span class="comment">
</span>  130 CONTINUE
      WORK( 1 ) = MAXWRK
<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="DGGEVX.716"></a><a href="dggevx.f.html#DGGEVX.1">DGGEVX</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>

⌨️ 快捷键说明

复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?