ztgex2.f.html

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

</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            DTRONG, WEAK
      INTEGER            I, M
      DOUBLE PRECISION   CQ, CZ, EPS, SA, SB, SCALE, SMLNUM, SS, SUM,
     $                   THRESH, WS
      COMPLEX*16         CDUM, F, G, SQ, SZ
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      COMPLEX*16         S( LDST, LDST ), T( LDST, LDST ), WORK( 8 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      DOUBLE PRECISION   <a name="DLAMCH.138"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>
      EXTERNAL           <a name="DLAMCH.139"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="ZLACPY.142"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>, <a name="ZLARTG.142"></a><a href="zlartg.f.html#ZLARTG.1">ZLARTG</a>, <a name="ZLASSQ.142"></a><a href="zlassq.f.html#ZLASSQ.1">ZLASSQ</a>, <a name="ZROT.142"></a><a href="zrot.f.html#ZROT.1">ZROT</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, DBLE, DCONJG, MAX, SQRT
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.LE.1 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span>      M = LDST
      WEAK = .FALSE.
      DTRONG = .FALSE.
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Make a local copy of selected block in (A, B)
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="ZLACPY.162"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>( <span class="string">'Full'</span>, M, M, A( J1, J1 ), LDA, S, LDST )
      CALL <a name="ZLACPY.163"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>( <span class="string">'Full'</span>, M, M, B( J1, J1 ), LDB, T, LDST )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute the threshold for testing the acceptance of swapping.
</span><span class="comment">*</span><span class="comment">
</span>      EPS = <a name="DLAMCH.167"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'P'</span> )
      SMLNUM = <a name="DLAMCH.168"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'S'</span> ) / EPS
      SCALE = DBLE( CZERO )
      SUM = DBLE( CONE )
      CALL <a name="ZLACPY.171"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>( <span class="string">'Full'</span>, M, M, S, LDST, WORK, M )
      CALL <a name="ZLACPY.172"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>( <span class="string">'Full'</span>, M, M, T, LDST, WORK( M*M+1 ), M )
      CALL <a name="ZLASSQ.173"></a><a href="zlassq.f.html#ZLASSQ.1">ZLASSQ</a>( 2*M*M, WORK, 1, SCALE, SUM )
      SA = SCALE*SQRT( SUM )
      THRESH = MAX( TEN*EPS*SA, SMLNUM )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute unitary QL and RQ that swap 1-by-1 and 1-by-1 blocks
</span><span class="comment">*</span><span class="comment">     using Givens rotations and perform the swap tentatively.
</span><span class="comment">*</span><span class="comment">
</span>      F = S( 2, 2 )*T( 1, 1 ) - T( 2, 2 )*S( 1, 1 )
      G = S( 2, 2 )*T( 1, 2 ) - T( 2, 2 )*S( 1, 2 )
      SA = ABS( S( 2, 2 ) )
      SB = ABS( T( 2, 2 ) )
      CALL <a name="ZLARTG.184"></a><a href="zlartg.f.html#ZLARTG.1">ZLARTG</a>( G, F, CZ, SZ, CDUM )
      SZ = -SZ
      CALL <a name="ZROT.186"></a><a href="zrot.f.html#ZROT.1">ZROT</a>( 2, S( 1, 1 ), 1, S( 1, 2 ), 1, CZ, DCONJG( SZ ) )
      CALL <a name="ZROT.187"></a><a href="zrot.f.html#ZROT.1">ZROT</a>( 2, T( 1, 1 ), 1, T( 1, 2 ), 1, CZ, DCONJG( SZ ) )
      IF( SA.GE.SB ) THEN
         CALL <a name="ZLARTG.189"></a><a href="zlartg.f.html#ZLARTG.1">ZLARTG</a>( S( 1, 1 ), S( 2, 1 ), CQ, SQ, CDUM )
      ELSE
         CALL <a name="ZLARTG.191"></a><a href="zlartg.f.html#ZLARTG.1">ZLARTG</a>( T( 1, 1 ), T( 2, 1 ), CQ, SQ, CDUM )
      END IF
      CALL <a name="ZROT.193"></a><a href="zrot.f.html#ZROT.1">ZROT</a>( 2, S( 1, 1 ), LDST, S( 2, 1 ), LDST, CQ, SQ )
      CALL <a name="ZROT.194"></a><a href="zrot.f.html#ZROT.1">ZROT</a>( 2, T( 1, 1 ), LDST, T( 2, 1 ), LDST, CQ, SQ )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Weak stability test: |S21| + |T21| &lt;= O(EPS F-norm((S, T)))
</span><span class="comment">*</span><span class="comment">
</span>      WS = ABS( S( 2, 1 ) ) + ABS( T( 2, 1 ) )
      WEAK = WS.LE.THRESH
      IF( .NOT.WEAK )
     $   GO TO 20
<span class="comment">*</span><span class="comment">
</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="ZLACPY.208"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>( <span class="string">'Full'</span>, M, M, S, LDST, WORK, M )
         CALL <a name="ZLACPY.209"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>( <span class="string">'Full'</span>, M, M, T, LDST, WORK( M*M+1 ), M )
         CALL <a name="ZROT.210"></a><a href="zrot.f.html#ZROT.1">ZROT</a>( 2, WORK, 1, WORK( 3 ), 1, CZ, -DCONJG( SZ ) )
         CALL <a name="ZROT.211"></a><a href="zrot.f.html#ZROT.1">ZROT</a>( 2, WORK( 5 ), 1, WORK( 7 ), 1, CZ, -DCONJG( SZ ) )
         CALL <a name="ZROT.212"></a><a href="zrot.f.html#ZROT.1">ZROT</a>( 2, WORK, 2, WORK( 2 ), 2, CQ, -SQ )
         CALL <a name="ZROT.213"></a><a href="zrot.f.html#ZROT.1">ZROT</a>( 2, WORK( 5 ), 2, WORK( 6 ), 2, CQ, -SQ )
         DO 10 I = 1, 2
            WORK( I ) = WORK( I ) - A( J1+I-1, J1 )
            WORK( I+2 ) = WORK( I+2 ) - A( J1+I-1, J1+1 )
            WORK( I+4 ) = WORK( I+4 ) - B( J1+I-1, J1 )
            WORK( I+6 ) = WORK( I+6 ) - B( J1+I-1, J1+1 )
   10    CONTINUE
         SCALE = DBLE( CZERO )
         SUM = DBLE( CONE )
         CALL <a name="ZLASSQ.222"></a><a href="zlassq.f.html#ZLASSQ.1">ZLASSQ</a>( 2*M*M, WORK, 1, SCALE, SUM )
         SS = SCALE*SQRT( SUM )
         DTRONG = SS.LE.THRESH
         IF( .NOT.DTRONG )
     $      GO TO 20
      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">     equivalence transformations to the original matrix pair (A,B)
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="ZROT.232"></a><a href="zrot.f.html#ZROT.1">ZROT</a>( J1+1, A( 1, J1 ), 1, A( 1, J1+1 ), 1, CZ,
     $           DCONJG( SZ ) )
      CALL <a name="ZROT.234"></a><a href="zrot.f.html#ZROT.1">ZROT</a>( J1+1, B( 1, J1 ), 1, B( 1, J1+1 ), 1, CZ,
     $           DCONJG( SZ ) )
      CALL <a name="ZROT.236"></a><a href="zrot.f.html#ZROT.1">ZROT</a>( N-J1+1, A( J1, J1 ), LDA, A( J1+1, J1 ), LDA, CQ, SQ )
      CALL <a name="ZROT.237"></a><a href="zrot.f.html#ZROT.1">ZROT</a>( N-J1+1, B( J1, J1 ), LDB, B( J1+1, J1 ), LDB, CQ, SQ )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set  N1 by N2 (2,1) blocks to 0
</span><span class="comment">*</span><span class="comment">
</span>      A( J1+1, J1 ) = CZERO
      B( J1+1, J1 ) = CZERO
<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 <a name="ZROT.247"></a><a href="zrot.f.html#ZROT.1">ZROT</a>( N, Z( 1, J1 ), 1, Z( 1, J1+1 ), 1, CZ,
     $              DCONJG( SZ ) )
      IF( WANTQ )
     $   CALL <a name="ZROT.250"></a><a href="zrot.f.html#ZROT.1">ZROT</a>( N, Q( 1, J1 ), 1, Q( 1, J1+1 ), 1, CQ,
     $              DCONJG( SQ ) )
<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><span class="comment">*</span><span class="comment">     Exit with INFO = 1 if swap was rejected.
</span><span class="comment">*</span><span class="comment">
</span>   20 CONTINUE
      INFO = 1
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="ZTGEX2.263"></a><a href="ztgex2.f.html#ZTGEX2.1">ZTGEX2</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>