dtgsen.f.html

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

         DSUM = ONE
         CALL <a name="DLASSQ.544"></a><a href="dlassq.f.html#DLASSQ.1">DLASSQ</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="DLASSQ.553"></a><a href="dlassq.f.html#DLASSQ.1">DLASSQ</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
<span class="comment">*</span><span class="comment">
</span>      IF( WANTD ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute estimates of 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="DTGSYL.574"></a><a href="dtgsyl.f.html#DTGSYL.1">DTGSYL</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( 2*N1*N2+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="DTGSYL.581"></a><a href="dtgsyl.f.html#DTGSYL.1">DTGSYL</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( 2*N1*N2+1 ),
     $                   LWORK-2*N1*N2, IWORK, IERR )
         ELSE
<span class="comment">*</span><span class="comment">
</span><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="DLACN2.589"></a><a href="dlacn2.f.html#DLACN2.1">DLACN2</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="DLACN2.603"></a><a href="dlacn2.f.html#DLACN2.1">DLACN2</a>( MN2, WORK( MN2+1 ), WORK, IWORK, 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="DTGSYL.610"></a><a href="dtgsyl.f.html#DTGSYL.1">DTGSYL</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( 2*N1*N2+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="DTGSYL.619"></a><a href="dtgsyl.f.html#DTGSYL.1">DTGSYL</a>( <span class="string">'T'</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( 2*N1*N2+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="DLACN2.632"></a><a href="dlacn2.f.html#DLACN2.1">DLACN2</a>( MN2, WORK( MN2+1 ), WORK, IWORK, 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="DTGSYL.639"></a><a href="dtgsyl.f.html#DTGSYL.1">DTGSYL</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( 2*N1*N2+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="DTGSYL.648"></a><a href="dtgsyl.f.html#DTGSYL.1">DTGSYL</a>( <span class="string">'T'</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( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK,
     $                         IERR )
               END IF
               GO TO 50
            END IF
            DIF( 2 ) = DSCALE / DIF( 2 )
<span class="comment">*</span><span class="comment">
</span>         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>   60 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute generalized eigenvalues of reordered pair (A, B) and
</span><span class="comment">*</span><span class="comment">     normalize the generalized Schur form.
</span><span class="comment">*</span><span class="comment">
</span>      PAIR = .FALSE.
      DO 80 K = 1, N
         IF( PAIR ) THEN
            PAIR = .FALSE.
         ELSE
<span class="comment">*</span><span class="comment">
</span>            IF( K.LT.N ) THEN
               IF( A( K+1, K ).NE.ZERO ) THEN
                  PAIR = .TRUE.
               END IF
            END IF
<span class="comment">*</span><span class="comment">
</span>            IF( PAIR ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">             Compute the eigenvalue(s) at position K.
</span><span class="comment">*</span><span class="comment">
</span>               WORK( 1 ) = A( K, K )
               WORK( 2 ) = A( K+1, K )
               WORK( 3 ) = A( K, K+1 )
               WORK( 4 ) = A( K+1, K+1 )
               WORK( 5 ) = B( K, K )
               WORK( 6 ) = B( K+1, K )
               WORK( 7 ) = B( K, K+1 )
               WORK( 8 ) = B( K+1, K+1 )
               CALL <a name="DLAG2.690"></a><a href="dlag2.f.html#DLAG2.1">DLAG2</a>( WORK, 2, WORK( 5 ), 2, SMLNUM*EPS, BETA( K ),
     $                     BETA( K+1 ), ALPHAR( K ), ALPHAR( K+1 ),
     $                     ALPHAI( K ) )
               ALPHAI( K+1 ) = -ALPHAI( K )
<span class="comment">*</span><span class="comment">
</span>            ELSE
<span class="comment">*</span><span class="comment">
</span>               IF( SIGN( ONE, B( K, K ) ).LT.ZERO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 If B(K,K) is negative, make it positive
</span><span class="comment">*</span><span class="comment">
</span>                  DO 70 I = 1, N
                     A( K, I ) = -A( K, I )
                     B( K, I ) = -B( K, I )
                     Q( I, K ) = -Q( I, K )
   70             CONTINUE
               END IF
<span class="comment">*</span><span class="comment">
</span>               ALPHAR( K ) = A( K, K )
               ALPHAI( K ) = ZERO
               BETA( K ) = B( K, K )
<span class="comment">*</span><span class="comment">
</span>            END IF
         END IF
   80 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="DTGSEN.721"></a><a href="dtgsen.f.html#DTGSEN.1">DTGSEN</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>