ztgsen.f.html

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

         I = N1 + 1
         CALL <a name="ZLACPY.493"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>( <span class="string">'Full'</span>, N1, N2, A( 1, I ), LDA, WORK, N1 )
         CALL <a name="ZLACPY.494"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>( <span class="string">'Full'</span>, N1, N2, B( 1, I ), LDB, WORK( N1*N2+1 ),
     $                N1 )
         IJB = 0
         CALL <a name="ZTGSYL.497"></a><a href="ztgsyl.f.html#ZTGSYL.1">ZTGSYL</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="ZLASSQ.507"></a><a href="zlassq.f.html#ZLASSQ.1">ZLASSQ</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="ZLASSQ.516"></a><a href="zlassq.f.html#ZLASSQ.1">ZLASSQ</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="ZTGSYL.536"></a><a href="ztgsyl.f.html#ZTGSYL.1">ZTGSYL</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="ZTGSYL.543"></a><a href="ztgsyl.f.html#ZTGSYL.1">ZTGSYL</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="ZLACN2.550"></a><a href="zlacn2.f.html#ZLACN2.1">ZLACN2</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="ZLACN2.564"></a><a href="zlacn2.f.html#ZLACN2.1">ZLACN2</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="ZTGSYL.571"></a><a href="ztgsyl.f.html#ZTGSYL.1">ZTGSYL</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="ZTGSYL.580"></a><a href="ztgsyl.f.html#ZTGSYL.1">ZTGSYL</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="ZLACN2.593"></a><a href="zlacn2.f.html#ZLACN2.1">ZLACN2</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="ZTGSYL.600"></a><a href="ztgsyl.f.html#ZTGSYL.1">ZTGSYL</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="ZTGSYL.609"></a><a href="ztgsyl.f.html#ZTGSYL.1">ZTGSYL</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 = DCONJG( B( K, K ) / DSCALE )
            TEMP2 = B( K, K ) / DSCALE
            B( K, K ) = DSCALE
            CALL ZSCAL( N-K, TEMP1, B( K, K+1 ), LDB )
            CALL ZSCAL( N-K+1, TEMP1, A( K, K ), LDA )
            IF( WANTQ )
     $         CALL ZSCAL( N, TEMP2, Q( 1, K ), 1 )
         ELSE
            B( K, K ) = DCMPLX( 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="ZTGSEN.651"></a><a href="ztgsen.f.html#ZTGSEN.1">ZTGSEN</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>