zgges.f.html

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

</span><span class="comment">*</span><span class="comment">     Scale B if max element outside range [SMLNUM,BIGNUM]
</span><span class="comment">*</span><span class="comment">
</span>      BNRM = <a name="ZLANGE.320"></a><a href="zlange.f.html#ZLANGE.1">ZLANGE</a>( <span class="string">'M'</span>, N, N, B, LDB, RWORK )
      ILBSCL = .FALSE.
      IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
         BNRMTO = SMLNUM
         ILBSCL = .TRUE.
      ELSE IF( BNRM.GT.BIGNUM ) THEN
         BNRMTO = BIGNUM
         ILBSCL = .TRUE.
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( ILBSCL )
     $   CALL <a name="ZLASCL.331"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>( <span class="string">'G'</span>, 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Permute the matrix to make it more nearly triangular
</span><span class="comment">*</span><span class="comment">     (Real Workspace: need 6*N)
</span><span class="comment">*</span><span class="comment">
</span>      ILEFT = 1
      IRIGHT = N + 1
      IRWRK = IRIGHT + N
      CALL <a name="ZGGBAL.339"></a><a href="zggbal.f.html#ZGGBAL.1">ZGGBAL</a>( <span class="string">'P'</span>, N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
     $             RWORK( IRIGHT ), RWORK( IRWRK ), IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Reduce B to triangular form (QR decomposition of B)
</span><span class="comment">*</span><span class="comment">     (Complex Workspace: need N, prefer N*NB)
</span><span class="comment">*</span><span class="comment">
</span>      IROWS = IHI + 1 - ILO
      ICOLS = N + 1 - ILO
      ITAU = 1
      IWRK = ITAU + IROWS
      CALL <a name="ZGEQRF.349"></a><a href="zgeqrf.f.html#ZGEQRF.1">ZGEQRF</a>( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
     $             WORK( IWRK ), LWORK+1-IWRK, IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Apply the orthogonal transformation to matrix A
</span><span class="comment">*</span><span class="comment">     (Complex Workspace: need N, prefer N*NB)
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="ZUNMQR.355"></a><a href="zunmqr.f.html#ZUNMQR.1">ZUNMQR</a>( <span class="string">'L'</span>, <span class="string">'C'</span>, IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
     $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
     $             LWORK+1-IWRK, IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Initialize VSL
</span><span class="comment">*</span><span class="comment">     (Complex Workspace: need N, prefer N*NB)
</span><span class="comment">*</span><span class="comment">
</span>      IF( ILVSL ) THEN
         CALL <a name="ZLASET.363"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>( <span class="string">'Full'</span>, N, N, CZERO, CONE, VSL, LDVSL )
         IF( IROWS.GT.1 ) THEN
            CALL <a name="ZLACPY.365"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>( <span class="string">'L'</span>, IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
     $                   VSL( ILO+1, ILO ), LDVSL )
         END IF
         CALL <a name="ZUNGQR.368"></a><a href="zungqr.f.html#ZUNGQR.1">ZUNGQR</a>( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
     $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Initialize VSR
</span><span class="comment">*</span><span class="comment">
</span>      IF( ILVSR )
     $   CALL <a name="ZLASET.375"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>( <span class="string">'Full'</span>, N, N, CZERO, CONE, VSR, LDVSR )
<span class="comment">*</span><span class="comment">
</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>      CALL <a name="ZGGHRD.380"></a><a href="zgghrd.f.html#ZGGHRD.1">ZGGHRD</a>( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
     $             LDVSL, VSR, LDVSR, IERR )
<span class="comment">*</span><span class="comment">
</span>      SDIM = 0
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Perform QZ algorithm, computing Schur vectors if desired
</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
      CALL <a name="ZHGEQZ.390"></a><a href="zhgeqz.f.html#ZHGEQZ.1">ZHGEQZ</a>( <span class="string">'S'</span>, JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
     $             ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK( IWRK ),
     $             LWORK+1-IWRK, RWORK( IRWRK ), 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 30
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Sort eigenvalues ALPHA/BETA if desired
</span><span class="comment">*</span><span class="comment">     (Workspace: none needed)
</span><span class="comment">*</span><span class="comment">
</span>      IF( WANTST ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Undo scaling on eigenvalues before selecting
</span><span class="comment">*</span><span class="comment">
</span>         IF( ILASCL )
     $      CALL <a name="ZLASCL.412"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRM, ANRMTO, N, 1, ALPHA, N, IERR )
         IF( ILBSCL )
     $      CALL <a name="ZLASCL.414"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>( <span class="string">'G'</span>, 0, 0, BNRM, BNRMTO, N, 1, BETA, N, IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Select eigenvalues
</span><span class="comment">*</span><span class="comment">
</span>         DO 10 I = 1, N
            BWORK( I ) = SELCTG( ALPHA( I ), BETA( I ) )
   10    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         CALL <a name="ZTGSEN.422"></a><a href="ztgsen.f.html#ZTGSEN.1">ZTGSEN</a>( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB, ALPHA,
     $                BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PVSL, PVSR,
     $                DIF, WORK( IWRK ), LWORK-IWRK+1, IDUM, 1, IERR )
         IF( IERR.EQ.1 )
     $      INFO = N + 3
<span class="comment">*</span><span class="comment">
</span>      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Apply back-permutation to VSL and VSR
</span><span class="comment">*</span><span class="comment">     (Workspace: none needed)
</span><span class="comment">*</span><span class="comment">
</span>      IF( ILVSL )
     $   CALL <a name="ZGGBAK.434"></a><a href="zggbak.f.html#ZGGBAK.1">ZGGBAK</a>( <span class="string">'P'</span>, <span class="string">'L'</span>, N, ILO, IHI, RWORK( ILEFT ),
     $                RWORK( IRIGHT ), N, VSL, LDVSL, IERR )
      IF( ILVSR )
     $   CALL <a name="ZGGBAK.437"></a><a href="zggbak.f.html#ZGGBAK.1">ZGGBAK</a>( <span class="string">'P'</span>, <span class="string">'R'</span>, N, ILO, IHI, RWORK( ILEFT ),
     $                RWORK( IRIGHT ), N, VSR, LDVSR, IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Undo scaling
</span><span class="comment">*</span><span class="comment">
</span>      IF( ILASCL ) THEN
         CALL <a name="ZLASCL.443"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>( <span class="string">'U'</span>, 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
         CALL <a name="ZLASCL.444"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( ILBSCL ) THEN
         CALL <a name="ZLASCL.448"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>( <span class="string">'U'</span>, 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
         CALL <a name="ZLASCL.449"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</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>      IF( WANTST ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Check if reordering is correct
</span><span class="comment">*</span><span class="comment">
</span>         LASTSL = .TRUE.
         SDIM = 0
         DO 20 I = 1, N
            CURSL = SELCTG( ALPHA( I ), BETA( I ) )
            IF( CURSL )
     $         SDIM = SDIM + 1
            IF( CURSL .AND. .NOT.LASTSL )
     $         INFO = N + 2
            LASTSL = CURSL
   20    CONTINUE
<span class="comment">*</span><span class="comment">
</span>      END IF
<span class="comment">*</span><span class="comment">
</span>   30 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      WORK( 1 ) = LWKOPT
<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="ZGGES.475"></a><a href="zgges.f.html#ZGGES.1">ZGGES</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>