zggevx.f.html

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

</span><span class="comment">*</span><span class="comment">     Reduce to generalized Hessenberg form
</span><span class="comment">*</span><span class="comment">     (Workspace: none needed)
</span><span class="comment">*</span><span class="comment">
</span>      IF( ILV .OR. .NOT.WANTSN ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Eigenvectors requested -- work on whole matrix.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="ZGGHRD.499"></a><a href="zgghrd.f.html#ZGGHRD.1">ZGGHRD</a>( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
     $                LDVL, VR, LDVR, IERR )
      ELSE
         CALL <a name="ZGGHRD.502"></a><a href="zgghrd.f.html#ZGGHRD.1">ZGGHRD</a>( <span class="string">'N'</span>, <span class="string">'N'</span>, IROWS, 1, IROWS, A( ILO, ILO ), LDA,
     $                B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Perform QZ algorithm (Compute eigenvalues, and optionally, the
</span><span class="comment">*</span><span class="comment">     Schur forms and Schur vectors)
</span><span class="comment">*</span><span class="comment">     (Complex Workspace: need N)
</span><span class="comment">*</span><span class="comment">     (Real Workspace: need N)
</span><span class="comment">*</span><span class="comment">
</span>      IWRK = ITAU
      IF( ILV .OR. .NOT.WANTSN ) THEN
         CHTEMP = <span class="string">'S'</span>
      ELSE
         CHTEMP = <span class="string">'E'</span>
      END IF
<span class="comment">*</span><span class="comment">
</span>      CALL <a name="ZHGEQZ.518"></a><a href="zhgeqz.f.html#ZHGEQZ.1">ZHGEQZ</a>( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
     $             ALPHA, BETA, VL, LDVL, VR, LDVR, WORK( IWRK ),
     $             LWORK+1-IWRK, RWORK, IERR )
      IF( IERR.NE.0 ) THEN
         IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
            INFO = IERR
         ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
            INFO = IERR - N
         ELSE
            INFO = N + 1
         END IF
         GO TO 90
      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">     <a name="ZTGEVC.533"></a><a href="ztgevc.f.html#ZTGEVC.1">ZTGEVC</a>: (Complex Workspace: need 2*N )
</span><span class="comment">*</span><span class="comment">             (Real Workspace:    need 2*N )
</span><span class="comment">*</span><span class="comment">     <a name="ZTGSNA.535"></a><a href="ztgsna.f.html#ZTGSNA.1">ZTGSNA</a>: (Complex Workspace: need 2*N*N if SENSE='V' or 'B')
</span><span class="comment">*</span><span class="comment">             (Integer Workspace: need N+2 )
</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="ZTGEVC.550"></a><a href="ztgevc.f.html#ZTGEVC.1">ZTGEVC</a>( CHTEMP, <span class="string">'B'</span>, LDUMMA, N, A, LDA, B, LDB, VL,
     $                   LDVL, VR, LDVR, N, IN, WORK( IWRK ), RWORK,
     $                   IERR )
            IF( IERR.NE.0 ) THEN
               INFO = N + 2
               GO TO 90
            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.561"></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.562"></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
</span><span class="comment">*</span><span class="comment">           re-calculate eigenvectors and estimate the condition numbers
</span><span class="comment">*</span><span class="comment">           one at a time.
</span><span class="comment">*</span><span class="comment">
</span>            DO 20 I = 1, N
<span class="comment">*</span><span class="comment">
</span>               DO 10 J = 1, N
                  BWORK( J ) = .FALSE.
   10          CONTINUE
               BWORK( I ) = .TRUE.
<span class="comment">*</span><span class="comment">
</span>               IWRK = N + 1
               IWRK1 = IWRK + N
<span class="comment">*</span><span class="comment">
</span>               IF( WANTSE .OR. WANTSB ) THEN
                  CALL <a name="ZTGEVC.581"></a><a href="ztgevc.f.html#ZTGEVC.1">ZTGEVC</a>( <span class="string">'B'</span>, <span class="string">'S'</span>, BWORK, N, A, LDA, B, LDB,
     $                         WORK( 1 ), N, WORK( IWRK ), N, 1, M,
     $                         WORK( IWRK1 ), RWORK, IERR )
                  IF( IERR.NE.0 ) THEN
                     INFO = N + 2
                     GO TO 90
                  END IF
               END IF
<span class="comment">*</span><span class="comment">
</span>               CALL <a name="ZTGSNA.590"></a><a href="ztgsna.f.html#ZTGSNA.1">ZTGSNA</a>( SENSE, <span class="string">'S'</span>, BWORK, N, A, LDA, B, LDB,
     $                      WORK( 1 ), N, WORK( IWRK ), N, RCONDE( I ),
     $                      RCONDV( I ), 1, 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="ZGGBAK.603"></a><a href="zggbak.f.html#ZGGBAK.1">ZGGBAK</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 50 JC = 1, N
            TEMP = ZERO
            DO 30 JR = 1, N
               TEMP = MAX( TEMP, ABS1( VL( JR, JC ) ) )
   30       CONTINUE
            IF( TEMP.LT.SMLNUM )
     $         GO TO 50
            TEMP = ONE / TEMP
            DO 40 JR = 1, N
               VL( JR, JC ) = VL( JR, JC )*TEMP
   40       CONTINUE
   50    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( ILVR ) THEN
         CALL <a name="ZGGBAK.621"></a><a href="zggbak.f.html#ZGGBAK.1">ZGGBAK</a>( BALANC, <span class="string">'R'</span>, N, ILO, IHI, LSCALE, RSCALE, N, VR,
     $                LDVR, IERR )
         DO 80 JC = 1, N
            TEMP = ZERO
            DO 60 JR = 1, N
               TEMP = MAX( TEMP, ABS1( VR( JR, JC ) ) )
   60       CONTINUE
            IF( TEMP.LT.SMLNUM )
     $         GO TO 80
            TEMP = ONE / TEMP
            DO 70 JR = 1, N
               VR( JR, JC ) = VR( JR, JC )*TEMP
   70       CONTINUE
   80    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 )
     $   CALL <a name="ZLASCL.640"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
<span class="comment">*</span><span class="comment">
</span>      IF( ILBSCL )
     $   CALL <a name="ZLASCL.643"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>( <span class="string">'G'</span>, 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
<span class="comment">*</span><span class="comment">
</span>   90 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="ZGGEVX.650"></a><a href="zggevx.f.html#ZGGEVX.1">ZGGEVX</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>