stgex2.f.html

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

         CALL SROT( 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 SROT( N, Z( 1, J1 ), 1, Z( 1, J1+1 ), 1, IR( 1, 1 ),
     $                 IR( 2, 1 ) )
         IF( WANTQ )
     $      CALL SROT( 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="SLACPY.317"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'Full'</span>, N1, N2, T( 1, N1+1 ), LDST, LI, LDST )
         CALL <a name="SLACPY.318"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'Full'</span>, N1, N2, S( 1, N1+1 ), LDST,
     $                IR( N2+1, N1+1 ), LDST )
         CALL <a name="STGSY2.320"></a><a href="stgsy2.f.html#STGSY2.1">STGSY2</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 SSCAL( N1, -ONE, LI( 1, I ), 1 )
            LI( N1+I, I ) = SCALE
   10    CONTINUE
         CALL <a name="SGEQR2.337"></a><a href="sgeqr2.f.html#SGEQR2.1">SGEQR2</a>( M, N2, LI, LDST, TAUL, WORK, LINFO )
         IF( LINFO.NE.0 )
     $      GO TO 70
         CALL <a name="SORG2R.340"></a><a href="sorg2r.f.html#SORG2R.1">SORG2R</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="SGERQ2.353"></a><a href="sgerq2.f.html#SGERQ2.1">SGERQ2</a>( N1, M, IR( N2+1, 1 ), LDST, TAUR, WORK, LINFO )
         IF( LINFO.NE.0 )
     $      GO TO 70
         CALL <a name="SORGR2.356"></a><a href="sorgr2.f.html#SORGR2.1">SORGR2</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 SGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, M, M, M, ONE, LI, LDST, S, LDST, ZERO,
     $               WORK, M )
         CALL SGEMM( <span class="string">'N'</span>, <span class="string">'T'</span>, M, M, M, ONE, WORK, M, IR, LDST, ZERO, S,
     $               LDST )
         CALL SGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, M, M, M, ONE, LI, LDST, T, LDST, ZERO,
     $               WORK, M )
         CALL SGEMM( <span class="string">'N'</span>, <span class="string">'T'</span>, M, M, M, ONE, WORK, M, IR, LDST, ZERO, T,
     $               LDST )
         CALL <a name="SLACPY.370"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'F'</span>, M, M, S, LDST, SCPY, LDST )
         CALL <a name="SLACPY.371"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'F'</span>, M, M, T, LDST, TCPY, LDST )
         CALL <a name="SLACPY.372"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'F'</span>, M, M, IR, LDST, IRCOP, LDST )
         CALL <a name="SLACPY.373"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</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="SGERQ2.378"></a><a href="sgerq2.f.html#SGERQ2.1">SGERQ2</a>( M, M, T, LDST, TAUR, WORK, LINFO )
         IF( LINFO.NE.0 )
     $      GO TO 70
         CALL <a name="SORMR2.381"></a><a href="sormr2.f.html#SORMR2.1">SORMR2</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="SORMR2.385"></a><a href="sormr2.f.html#SORMR2.1">SORMR2</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="SLASSQ.395"></a><a href="slassq.f.html#SLASSQ.1">SLASSQ</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="SGEQR2.402"></a><a href="sgeqr2.f.html#SGEQR2.1">SGEQR2</a>( M, M, TCPY, LDST, TAUL, WORK, LINFO )
         IF( LINFO.NE.0 )
     $      GO TO 70
         CALL <a name="SORM2R.405"></a><a href="sorm2r.f.html#SORM2R.1">SORM2R</a>( <span class="string">'L'</span>, <span class="string">'T'</span>, M, M, M, TCPY, LDST, TAUL, SCPY, LDST,
     $                WORK, INFO )
         CALL <a name="SORM2R.407"></a><a href="sorm2r.f.html#SORM2R.1">SORM2R</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="SLASSQ.417"></a><a href="slassq.f.html#SLASSQ.1">SLASSQ</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="SLACPY.426"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'F'</span>, M, M, SCPY, LDST, S, LDST )
            CALL <a name="SLACPY.427"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'F'</span>, M, M, TCPY, LDST, T, LDST )
            CALL <a name="SLACPY.428"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'F'</span>, M, M, IRCOP, LDST, IR, LDST )
            CALL <a name="SLACPY.429"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</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="SLASET.436"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'Lower'</span>, M-1, M-1, ZERO, ZERO, T(2,1), LDST )
<span class="comment">*</span><span class="comment">