cggesx.f.html

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

</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="CLASET.438"></a><a href="claset.f.html#CLASET.1">CLASET</a>( <span class="string">'Full'</span>, N, N, CZERO, CONE, VSL, LDVSL )
         IF( IROWS.GT.1 ) THEN
            CALL <a name="CLACPY.440"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</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="CUNGQR.443"></a><a href="cungqr.f.html#CUNGQR.1">CUNGQR</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="CLASET.450"></a><a href="claset.f.html#CLASET.1">CLASET</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="CGGHRD.455"></a><a href="cgghrd.f.html#CGGHRD.1">CGGHRD</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="CHGEQZ.465"></a><a href="chgeqz.f.html#CHGEQZ.1">CHGEQZ</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 40
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Sort eigenvalues ALPHA/BETA and compute the reciprocal of
</span><span class="comment">*</span><span class="comment">     condition number(s)
</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 SELCTGing
</span><span class="comment">*</span><span class="comment">
</span>         IF( ILASCL )
     $      CALL <a name="CLASCL.487"></a><a href="clascl.f.html#CLASCL.1">CLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
         IF( ILBSCL )
     $      CALL <a name="CLASCL.489"></a><a href="clascl.f.html#CLASCL.1">CLASCL</a>( <span class="string">'G'</span>, 0, 0, BNRMTO, BNRM, 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><span class="comment">*</span><span class="comment">        Reorder eigenvalues, transform Generalized Schur vectors, and
</span><span class="comment">*</span><span class="comment">        compute reciprocal condition numbers
</span><span class="comment">*</span><span class="comment">        (Complex Workspace: If IJOB &gt;= 1, need MAX(1, 2*SDIM*(N-SDIM))
</span><span class="comment">*</span><span class="comment">                            otherwise, need 1 )
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CTGSEN.502"></a><a href="ctgsen.f.html#CTGSEN.1">CTGSEN</a>( IJOB, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB,
     $                ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PL, PR,
     $                DIF, WORK( IWRK ), LWORK-IWRK+1, IWORK, LIWORK,
     $                IERR )
<span class="comment">*</span><span class="comment">
</span>         IF( IJOB.GE.1 )
     $      MAXWRK = MAX( MAXWRK, 2*SDIM*( N-SDIM ) )
         IF( IERR.EQ.-21 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">            not enough complex workspace
</span><span class="comment">*</span><span class="comment">
</span>            INFO = -21
         ELSE
            IF( IJOB.EQ.1 .OR. IJOB.EQ.4 ) THEN
               RCONDE( 1 ) = PL
               RCONDE( 2 ) = PR
            END IF
            IF( IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN
               RCONDV( 1 ) = DIF( 1 )
               RCONDV( 2 ) = DIF( 2 )
            END IF
            IF( IERR.EQ.1 )
     $         INFO = N + 3
         END IF
<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 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="CGGBAK.533"></a><a href="cggbak.f.html#CGGBAK.1">CGGBAK</a>( <span class="string">'P'</span>, <span class="string">'L'</span>, N, ILO, IHI, RWORK( ILEFT ),
     $                RWORK( IRIGHT ), N, VSL, LDVSL, IERR )
<span class="comment">*</span><span class="comment">
</span>      IF( ILVSR )
     $   CALL <a name="CGGBAK.537"></a><a href="cggbak.f.html#CGGBAK.1">CGGBAK</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="CLASCL.543"></a><a href="clascl.f.html#CLASCL.1">CLASCL</a>( <span class="string">'U'</span>, 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
         CALL <a name="CLASCL.544"></a><a href="clascl.f.html#CLASCL.1">CLASCL</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="CLASCL.548"></a><a href="clascl.f.html#CLASCL.1">CLASCL</a>( <span class="string">'U'</span>, 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
         CALL <a name="CLASCL.549"></a><a href="clascl.f.html#CLASCL.1">CLASCL</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 30 I = 1, N
            CURSL = SELCTG( ALPHA( I ), BETA( I ) )
            IF( CURSL )
     $         SDIM = SDIM + 1
            IF( CURSL .AND. .NOT.LASTSL )
     $         INFO = N + 2
            LASTSL = CURSL
   30    CONTINUE
<span class="comment">*</span><span class="comment">
</span>      END IF
<span class="comment">*</span><span class="comment">
</span>   40 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      WORK( 1 ) = MAXWRK
      IWORK( 1 ) = LIWMIN
<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="CGGESX.576"></a><a href="cggesx.f.html#CGGESX.1">CGGESX</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>