dtgex2.f.html

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

         CALL DROT( N-J1+1, B( J1, J1 ), LDB, B( J1+1, J1 ), LDB,
     $              LI( 1, 1 ), LI( 2, 1 ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Set  N1-by-N2 (2,1) - blocks to ZERO.
</span><span class="comment">*</span><span class="comment">
</span>         A( J1+1, J1 ) = ZERO
         B( J1+1, J1 ) = ZERO
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Accumulate transformations into Q and Z if requested.
</span><span class="comment">*</span><span class="comment">
</span>         IF( WANTZ )
     $      CALL DROT( N, Z( 1, J1 ), 1, Z( 1, J1+1 ), 1, IR( 1, 1 ),
     $                 IR( 2, 1 ) )
         IF( WANTQ )
     $      CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J1+1 ), 1, LI( 1, 1 ),
     $                 LI( 2, 1 ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Exit with INFO = 0 if swap was successfully performed.
</span><span class="comment">*</span><span class="comment">
</span>         RETURN
<span class="comment">*</span><span class="comment">
</span>      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        CASE 2: Swap 1-by-1 and 2-by-2 blocks, or 2-by-2
</span><span class="comment">*</span><span class="comment">                and 2-by-2 blocks.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Solve the generalized Sylvester equation
</span><span class="comment">*</span><span class="comment">                 S11 * R - L * S22 = SCALE * S12
</span><span class="comment">*</span><span class="comment">                 T11 * R - L * T22 = SCALE * T12
</span><span class="comment">*</span><span class="comment">        for R and L. Solutions in LI and IR.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="DLACPY.317"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'Full'</span>, N1, N2, T( 1, N1+1 ), LDST, LI, LDST )
         CALL <a name="DLACPY.318"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'Full'</span>, N1, N2, S( 1, N1+1 ), LDST,
     $                IR( N2+1, N1+1 ), LDST )
         CALL <a name="DTGSY2.320"></a><a href="dtgsy2.f.html#DTGSY2.1">DTGSY2</a>( <span class="string">'N'</span>, 0, N1, N2, S, LDST, S( N1+1, N1+1 ), LDST,
     $                IR( N2+1, N1+1 ), LDST, T, LDST, T( N1+1, N1+1 ),
     $                LDST, LI, LDST, SCALE, DSUM, DSCALE, IWORK, IDUM,
     $                LINFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute orthogonal matrix QL:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    QL' * LI = [ TL ]
</span><span class="comment">*</span><span class="comment">                               [ 0  ]
</span><span class="comment">*</span><span class="comment">        where
</span><span class="comment">*</span><span class="comment">                    LI =  [      -L              ]
</span><span class="comment">*</span><span class="comment">                          [ SCALE * identity(N2) ]
</span><span class="comment">*</span><span class="comment">
</span>         DO 10 I = 1, N2
            CALL DSCAL( N1, -ONE, LI( 1, I ), 1 )
            LI( N1+I, I ) = SCALE
   10    CONTINUE
         CALL <a name="DGEQR2.337"></a><a href="dgeqr2.f.html#DGEQR2.1">DGEQR2</a>( M, N2, LI, LDST, TAUL, WORK, LINFO )
         IF( LINFO.NE.0 )
     $      GO TO 70
         CALL <a name="DORG2R.340"></a><a href="dorg2r.f.html#DORG2R.1">DORG2R</a>( M, M, N2, LI, LDST, TAUL, WORK, LINFO )
         IF( LINFO.NE.0 )
     $      GO TO 70
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute orthogonal matrix RQ:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    IR * RQ' =   [ 0  TR],
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">         where IR = [ SCALE * identity(N1), R ]
</span><span class="comment">*</span><span class="comment">
</span>         DO 20 I = 1, N1
            IR( N2+I, I ) = SCALE
   20    CONTINUE
         CALL <a name="DGERQ2.353"></a><a href="dgerq2.f.html#DGERQ2.1">DGERQ2</a>( N1, M, IR( N2+1, 1 ), LDST, TAUR, WORK, LINFO )
         IF( LINFO.NE.0 )
     $      GO TO 70
         CALL <a name="DORGR2.356"></a><a href="dorgr2.f.html#DORGR2.1">DORGR2</a>( M, M, N1, IR, LDST, TAUR, WORK, LINFO )
         IF( LINFO.NE.0 )
     $      GO TO 70
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Perform the swapping tentatively:
</span><span class="comment">*</span><span class="comment">
</span>         CALL DGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, M, M, M, ONE, LI, LDST, S, LDST, ZERO,
     $               WORK, M )
         CALL DGEMM( <span class="string">'N'</span>, <span class="string">'T'</span>, M, M, M, ONE, WORK, M, IR, LDST, ZERO, S,
     $               LDST )
         CALL DGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, M, M, M, ONE, LI, LDST, T, LDST, ZERO,
     $               WORK, M )
         CALL DGEMM( <span class="string">'N'</span>, <span class="string">'T'</span>, M, M, M, ONE, WORK, M, IR, LDST, ZERO, T,
     $               LDST )
         CALL <a name="DLACPY.370"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'F'</span>, M, M, S, LDST, SCPY, LDST )
         CALL <a name="DLACPY.371"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'F'</span>, M, M, T, LDST, TCPY, LDST )
         CALL <a name="DLACPY.372"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'F'</span>, M, M, IR, LDST, IRCOP, LDST )
         CALL <a name="DLACPY.373"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'F'</span>, M, M, LI, LDST, LICOP, LDST )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Triangularize the B-part by an RQ factorization.
</span><span class="comment">*</span><span class="comment">        Apply transformation (from left) to A-part, giving S.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="DGERQ2.378"></a><a href="dgerq2.f.html#DGERQ2.1">DGERQ2</a>( M, M, T, LDST, TAUR, WORK, LINFO )
         IF( LINFO.NE.0 )
     $      GO TO 70
         CALL <a name="DORMR2.381"></a><a href="dormr2.f.html#DORMR2.1">DORMR2</a>( <span class="string">'R'</span>, <span class="string">'T'</span>, M, M, M, T, LDST, TAUR, S, LDST, WORK,
     $                LINFO )
         IF( LINFO.NE.0 )
     $      GO TO 70
         CALL <a name="DORMR2.385"></a><a href="dormr2.f.html#DORMR2.1">DORMR2</a>( <span class="string">'L'</span>, <span class="string">'N'</span>, M, M, M, T, LDST, TAUR, IR, LDST, WORK,
     $                LINFO )
         IF( LINFO.NE.0 )
     $      GO TO 70
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute F-norm(S21) in BRQA21. (T21 is 0.)
</span><span class="comment">*</span><span class="comment">
</span>         DSCALE = ZERO
         DSUM = ONE
         DO 30 I = 1, N2
            CALL <a name="DLASSQ.395"></a><a href="dlassq.f.html#DLASSQ.1">DLASSQ</a>( N1, S( N2+1, I ), 1, DSCALE, DSUM )
   30    CONTINUE
         BRQA21 = DSCALE*SQRT( DSUM )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Triangularize the B-part by a QR factorization.
</span><span class="comment">*</span><span class="comment">        Apply transformation (from right) to A-part, giving S.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="DGEQR2.402"></a><a href="dgeqr2.f.html#DGEQR2.1">DGEQR2</a>( M, M, TCPY, LDST, TAUL, WORK, LINFO )
         IF( LINFO.NE.0 )
     $      GO TO 70
         CALL <a name="DORM2R.405"></a><a href="dorm2r.f.html#DORM2R.1">DORM2R</a>( <span class="string">'L'</span>, <span class="string">'T'</span>, M, M, M, TCPY, LDST, TAUL, SCPY, LDST,
     $                WORK, INFO )
         CALL <a name="DORM2R.407"></a><a href="dorm2r.f.html#DORM2R.1">DORM2R</a>( <span class="string">'R'</span>, <span class="string">'N'</span>, M, M, M, TCPY, LDST, TAUL, LICOP, LDST,
     $                WORK, INFO )
         IF( LINFO.NE.0 )
     $      GO TO 70
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute F-norm(S21) in BQRA21. (T21 is 0.)
</span><span class="comment">*</span><span class="comment">
</span>         DSCALE = ZERO
         DSUM = ONE
         DO 40 I = 1, N2
            CALL <a name="DLASSQ.417"></a><a href="dlassq.f.html#DLASSQ.1">DLASSQ</a>( N1, SCPY( N2+1, I ), 1, DSCALE, DSUM )
   40    CONTINUE
         BQRA21 = DSCALE*SQRT( DSUM )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Decide which method to use.
</span><span class="comment">*</span><span class="comment">          Weak stability test:
</span><span class="comment">*</span><span class="comment">             F-norm(S21) &lt;= O(EPS * F-norm((S, T)))
</span><span class="comment">*</span><span class="comment">
</span>         IF( BQRA21.LE.BRQA21 .AND. BQRA21.LE.THRESH ) THEN
            CALL <a name="DLACPY.426"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'F'</span>, M, M, SCPY, LDST, S, LDST )
            CALL <a name="DLACPY.427"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'F'</span>, M, M, TCPY, LDST, T, LDST )
            CALL <a name="DLACPY.428"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'F'</span>, M, M, IRCOP, LDST, IR, LDST )
            CALL <a name="DLACPY.429"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'F'</span>, M, M, LICOP, LDST, LI, LDST )
         ELSE IF( BRQA21.GE.THRESH ) THEN
            GO TO 70
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Set lower triangle of B-part to zero
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="DLASET.436"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'Lower'</span>, M-1, M-1, ZERO, ZERO, T(2,1), LDST )
<span class="comment">*</span><span class="comment">