ctgsen.f.html

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

         I = N1 + 1
         CALL <a name="CLACPY.490"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'Full'</span>, N1, N2, A( 1, I ), LDA, WORK, N1 )
         CALL <a name="CLACPY.491"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'Full'</span>, N1, N2, B( 1, I ), LDB, WORK( N1*N2+1 ),
     $                N1 )
         IJB = 0
         CALL <a name="CTGSYL.494"></a><a href="ctgsyl.f.html#CTGSYL.1">CTGSYL</a>( <span class="string">'N'</span>, IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK,
     $                N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ), N1,
     $                DSCALE, DIF( 1 ), WORK( N1*N2*2+1 ),
     $                LWORK-2*N1*N2, IWORK, IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Estimate the reciprocal of norms of &quot;projections&quot; onto
</span><span class="comment">*</span><span class="comment">        left and right eigenspaces
</span><span class="comment">*</span><span class="comment">
</span>         RDSCAL = ZERO
         DSUM = ONE
         CALL <a name="CLASSQ.504"></a><a href="classq.f.html#CLASSQ.1">CLASSQ</a>( N1*N2, WORK, 1, RDSCAL, DSUM )
         PL = RDSCAL*SQRT( DSUM )
         IF( PL.EQ.ZERO ) THEN
            PL = ONE
         ELSE
            PL = DSCALE / ( SQRT( DSCALE*DSCALE / PL+PL )*SQRT( PL ) )
         END IF
         RDSCAL = ZERO
         DSUM = ONE
         CALL <a name="CLASSQ.513"></a><a href="classq.f.html#CLASSQ.1">CLASSQ</a>( N1*N2, WORK( N1*N2+1 ), 1, RDSCAL, DSUM )
         PR = RDSCAL*SQRT( DSUM )
         IF( PR.EQ.ZERO ) THEN
            PR = ONE
         ELSE
            PR = DSCALE / ( SQRT( DSCALE*DSCALE / PR+PR )*SQRT( PR ) )
         END IF
      END IF
      IF( WANTD ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute estimates Difu and Difl.
</span><span class="comment">*</span><span class="comment">
</span>         IF( WANTD1 ) THEN
            N1 = M
            N2 = N - M
            I = N1 + 1
            IJB = IDIFJB
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Frobenius norm-based Difu estimate.
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="CTGSYL.533"></a><a href="ctgsyl.f.html#CTGSYL.1">CTGSYL</a>( <span class="string">'N'</span>, IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK,
     $                   N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ),
     $                   N1, DSCALE, DIF( 1 ), WORK( N1*N2*2+1 ),
     $                   LWORK-2*N1*N2, IWORK, IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Frobenius norm-based Difl estimate.
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="CTGSYL.540"></a><a href="ctgsyl.f.html#CTGSYL.1">CTGSYL</a>( <span class="string">'N'</span>, IJB, N2, N1, A( I, I ), LDA, A, LDA, WORK,
     $                   N2, B( I, I ), LDB, B, LDB, WORK( N1*N2+1 ),
     $                   N2, DSCALE, DIF( 2 ), WORK( N1*N2*2+1 ),
     $                   LWORK-2*N1*N2, IWORK, IERR )
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute 1-norm-based estimates of Difu and Difl using
</span><span class="comment">*</span><span class="comment">           reversed communication with <a name="CLACN2.547"></a><a href="clacn2.f.html#CLACN2.1">CLACN2</a>. In each step a
</span><span class="comment">*</span><span class="comment">           generalized Sylvester equation or a transposed variant
</span><span class="comment">*</span><span class="comment">           is solved.
</span><span class="comment">*</span><span class="comment">
</span>            KASE = 0
            N1 = M
            N2 = N - M
            I = N1 + 1
            IJB = 0
            MN2 = 2*N1*N2
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           1-norm-based estimate of Difu.
</span><span class="comment">*</span><span class="comment">
</span>   40       CONTINUE
            CALL <a name="CLACN2.561"></a><a href="clacn2.f.html#CLACN2.1">CLACN2</a>( MN2, WORK( MN2+1 ), WORK, DIF( 1 ), KASE,
     $                   ISAVE )
            IF( KASE.NE.0 ) THEN
               IF( KASE.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Solve generalized Sylvester equation
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="CTGSYL.568"></a><a href="ctgsyl.f.html#CTGSYL.1">CTGSYL</a>( <span class="string">'N'</span>, IJB, N1, N2, A, LDA, A( I, I ), LDA,
     $                         WORK, N1, B, LDB, B( I, I ), LDB,
     $                         WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ),
     $                         WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK,
     $                         IERR )
               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Solve the transposed variant.
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="CTGSYL.577"></a><a href="ctgsyl.f.html#CTGSYL.1">CTGSYL</a>( <span class="string">'C'</span>, IJB, N1, N2, A, LDA, A( I, I ), LDA,
     $                         WORK, N1, B, LDB, B( I, I ), LDB,
     $                         WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ),
     $                         WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK,
     $                         IERR )
               END IF
               GO TO 40
            END IF
            DIF( 1 ) = DSCALE / DIF( 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           1-norm-based estimate of Difl.
</span><span class="comment">*</span><span class="comment">
</span>   50       CONTINUE
            CALL <a name="CLACN2.590"></a><a href="clacn2.f.html#CLACN2.1">CLACN2</a>( MN2, WORK( MN2+1 ), WORK, DIF( 2 ), KASE,
     $                   ISAVE )
            IF( KASE.NE.0 ) THEN
               IF( KASE.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Solve generalized Sylvester equation
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="CTGSYL.597"></a><a href="ctgsyl.f.html#CTGSYL.1">CTGSYL</a>( <span class="string">'N'</span>, IJB, N2, N1, A( I, I ), LDA, A, LDA,
     $                         WORK, N2, B( I, I ), LDB, B, LDB,
     $                         WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ),
     $                         WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK,
     $                         IERR )
               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Solve the transposed variant.
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="CTGSYL.606"></a><a href="ctgsyl.f.html#CTGSYL.1">CTGSYL</a>( <span class="string">'C'</span>, IJB, N2, N1, A( I, I ), LDA, A, LDA,
     $                         WORK, N2, B, LDB, B( I, I ), LDB,
     $                         WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ),
     $                         WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK,
     $                         IERR )
               END IF
               GO TO 50
            END IF
            DIF( 2 ) = DSCALE / DIF( 2 )
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If B(K,K) is complex, make it real and positive (normalization
</span><span class="comment">*</span><span class="comment">     of the generalized Schur form) and Store the generalized 
</span><span class="comment">*</span><span class="comment">     eigenvalues of reordered pair (A, B)
</span><span class="comment">*</span><span class="comment">
</span>      DO 60 K = 1, N
         DSCALE = ABS( B( K, K ) )
         IF( DSCALE.GT.SAFMIN ) THEN
            TEMP1 = CONJG( B( K, K ) / DSCALE )
            TEMP2 = B( K, K ) / DSCALE
            B( K, K ) = DSCALE
            CALL CSCAL( N-K, TEMP1, B( K, K+1 ), LDB )
            CALL CSCAL( N-K+1, TEMP1, A( K, K ), LDA )
            IF( WANTQ )
     $         CALL CSCAL( N, TEMP2, Q( 1, K ), 1 )
         ELSE
            B( K, K ) = CMPLX( ZERO, ZERO )
         END IF
<span class="comment">*</span><span class="comment">
</span>         ALPHA( K ) = A( K, K )
         BETA( K ) = B( K, K )
<span class="comment">*</span><span class="comment">
</span>   60 CONTINUE
<span class="comment">*</span><span class="comment">
</span>   70 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      WORK( 1 ) = LWMIN
      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="CTGSEN.648"></a><a href="ctgsen.f.html#CTGSEN.1">CTGSEN</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>