stgex2.f.html

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

</span>         IF( WANDS ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Strong stability test:
</span><span class="comment">*</span><span class="comment">              F-norm((A-QL*S*QR', B-QL*T*QR')) &lt;= O(EPS*F-norm((A,B)))
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="SLACPY.443"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'Full'</span>, M, M, A( J1, J1 ), LDA, WORK( M*M+1 ),
     $                   M )
            CALL SGEMM( <span class="string">'N'</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">'N'</span>, M, M, M, -ONE, WORK, M, IR, LDST, ONE,
     $                  WORK( M*M+1 ), M )
            DSCALE = ZERO
            DSUM = ONE
            CALL <a name="SLASSQ.451"></a><a href="slassq.f.html#SLASSQ.1">SLASSQ</a>( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )
<span class="comment">*</span><span class="comment">
</span>            CALL <a name="SLACPY.453"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'Full'</span>, M, M, B( J1, J1 ), LDB, WORK( M*M+1 ),
     $                   M )
            CALL SGEMM( <span class="string">'N'</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">'N'</span>, M, M, M, -ONE, WORK, M, IR, LDST, ONE,
     $                  WORK( M*M+1 ), M )
            CALL <a name="SLASSQ.459"></a><a href="slassq.f.html#SLASSQ.1">SLASSQ</a>( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )
            SS = DSCALE*SQRT( DSUM )
            STRONG = ( SS.LE.THRESH )
            IF( .NOT.STRONG )
     $         GO TO 70
<span class="comment">*</span><span class="comment">
</span>         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        If the swap is accepted (&quot;weakly&quot; and &quot;strongly&quot;), apply the
</span><span class="comment">*</span><span class="comment">        transformations and set N1-by-N2 (2,1)-block to zero.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SLASET.470"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'Full'</span>, N1, N2, ZERO, ZERO, S(N2+1,1), LDST )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        copy back M-by-M diagonal block starting at index J1 of (A, B)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SLACPY.474"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'F'</span>, M, M, S, LDST, A( J1, J1 ), LDA )
         CALL <a name="SLACPY.475"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'F'</span>, M, M, T, LDST, B( J1, J1 ), LDB )
         CALL <a name="SLASET.476"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'Full'</span>, LDST, LDST, ZERO, ZERO, T, LDST )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Standardize existing 2-by-2 blocks.
</span><span class="comment">*</span><span class="comment">
</span>         DO 50 I = 1, M*M
            WORK(I) = ZERO
   50    CONTINUE
         WORK( 1 ) = ONE
         T( 1, 1 ) = ONE
         IDUM = LWORK - M*M - 2
         IF( N2.GT.1 ) THEN
            CALL <a name="SLAGV2.487"></a><a href="slagv2.f.html#SLAGV2.1">SLAGV2</a>( A( J1, J1 ), LDA, B( J1, J1 ), LDB, AR, AI, BE,
     $                   WORK( 1 ), WORK( 2 ), T( 1, 1 ), T( 2, 1 ) )
            WORK( M+1 ) = -WORK( 2 )
            WORK( M+2 ) = WORK( 1 )
            T( N2, N2 ) = T( 1, 1 )
            T( 1, 2 ) = -T( 2, 1 )
         END IF
         WORK( M*M ) = ONE
         T( M, M ) = ONE
<span class="comment">*</span><span class="comment">
</span>         IF( N1.GT.1 ) THEN
            CALL <a name="SLAGV2.498"></a><a href="slagv2.f.html#SLAGV2.1">SLAGV2</a>( A( J1+N2, J1+N2 ), LDA, B( J1+N2, J1+N2 ), LDB,
     $                   TAUR, TAUL, WORK( M*M+1 ), WORK( N2*M+N2+1 ),
     $                   WORK( N2*M+N2+2 ), T( N2+1, N2+1 ),
     $                   T( M, M-1 ) )
            WORK( M*M ) = WORK( N2*M+N2+1 )
            WORK( M*M-1 ) = -WORK( N2*M+N2+2 )
            T( M, M ) = T( N2+1, N2+1 )
            T( M-1, M ) = -T( M, M-1 )
         END IF
         CALL SGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, N2, N1, N2, ONE, WORK, M, A( J1, J1+N2 ),
     $               LDA, ZERO, WORK( M*M+1 ), N2 )
         CALL <a name="SLACPY.509"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'Full'</span>, N2, N1, WORK( M*M+1 ), N2, A( J1, J1+N2 ),
     $                LDA )
         CALL SGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, N2, N1, N2, ONE, WORK, M, B( J1, J1+N2 ),
     $               LDB, ZERO, WORK( M*M+1 ), N2 )
         CALL <a name="SLACPY.513"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'Full'</span>, N2, N1, WORK( M*M+1 ), N2, B( J1, J1+N2 ),
     $                LDB )
         CALL SGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, M, M, M, ONE, LI, LDST, WORK, M, ZERO,
     $               WORK( M*M+1 ), M )
         CALL <a name="SLACPY.517"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'Full'</span>, M, M, WORK( M*M+1 ), M, LI, LDST )
         CALL SGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, N2, N1, N1, ONE, A( J1, J1+N2 ), LDA,
     $               T( N2+1, N2+1 ), LDST, ZERO, WORK, N2 )
         CALL <a name="SLACPY.520"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'Full'</span>, N2, N1, WORK, N2, A( J1, J1+N2 ), LDA )
         CALL SGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, N2, N1, N1, ONE, B( J1, J1+N2 ), LDB,
     $               T( N2+1, N2+1 ), LDST, ZERO, WORK, N2 )
         CALL <a name="SLACPY.523"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'Full'</span>, N2, N1, WORK, N2, B( J1, J1+N2 ), LDB )
         CALL SGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, M, M, M, ONE, IR, LDST, T, LDST, ZERO,
     $               WORK, M )
         CALL <a name="SLACPY.526"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'Full'</span>, M, M, WORK, M, IR, LDST )
<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( WANTQ ) THEN
            CALL SGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, N, M, M, ONE, Q( 1, J1 ), LDQ, LI,
     $                  LDST, ZERO, WORK, N )
            CALL <a name="SLACPY.533"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'Full'</span>, N, M, WORK, N, Q( 1, J1 ), LDQ )
<span class="comment">*</span><span class="comment">
</span>         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( WANTZ ) THEN
            CALL SGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, N, M, M, ONE, Z( 1, J1 ), LDZ, IR,
     $                  LDST, ZERO, WORK, N )
            CALL <a name="SLACPY.540"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'Full'</span>, N, M, WORK, N, Z( 1, J1 ), LDZ )
<span class="comment">*</span><span class="comment">
</span>         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Update (A(J1:J1+M-1, M+J1:N), B(J1:J1+M-1, M+J1:N)) and
</span><span class="comment">*</span><span class="comment">                (A(1:J1-1, J1:J1+M), B(1:J1-1, J1:J1+M)).
</span><span class="comment">*</span><span class="comment">
</span>         I = J1 + M
         IF( I.LE.N ) THEN
            CALL SGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, M, N-I+1, M, ONE, LI, LDST,
     $                  A( J1, I ), LDA, ZERO, WORK, M )
            CALL <a name="SLACPY.551"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'Full'</span>, M, N-I+1, WORK, M, A( J1, I ), LDA )
            CALL SGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, M, N-I+1, M, ONE, LI, LDST,
     $                  B( J1, I ), LDB, ZERO, WORK, M )
            CALL <a name="SLACPY.554"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'Full'</span>, M, N-I+1, WORK, M, B( J1, I ), LDB )
         END IF
         I = J1 - 1
         IF( I.GT.0 ) THEN
            CALL SGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, I, M, M, ONE, A( 1, J1 ), LDA, IR,
     $                  LDST, ZERO, WORK, I )
            CALL <a name="SLACPY.560"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'Full'</span>, I, M, WORK, I, A( 1, J1 ), LDA )
            CALL SGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, I, M, M, ONE, B( 1, J1 ), LDB, IR,
     $                  LDST, ZERO, WORK, I )
            CALL <a name="SLACPY.563"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'Full'</span>, I, M, WORK, I, B( 1, J1 ), LDB )
         END IF
<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>      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Exit with INFO = 1 if swap was rejected.
</span><span class="comment">*</span><span class="comment">
</span>   70 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      INFO = 1
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="STGEX2.579"></a><a href="stgex2.f.html#STGEX2.1">STGEX2</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>