dtgex2.f.html

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

</span>      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
      DOUBLE PRECISION   TEN
      PARAMETER          ( TEN = 1.0D+01 )
      INTEGER            LDST
      PARAMETER          ( LDST = 4 )
      LOGICAL            WANDS
      PARAMETER          ( WANDS = .TRUE. )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            DTRONG, WEAK
      INTEGER            I, IDUM, LINFO, M
      DOUBLE PRECISION   BQRA21, BRQA21, DDUM, DNORM, DSCALE, DSUM, EPS,
     $                   F, G, SA, SB, SCALE, SMLNUM, SS, THRESH, WS
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      INTEGER            IWORK( LDST )
      DOUBLE PRECISION   AI( 2 ), AR( 2 ), BE( 2 ), IR( LDST, LDST ),
     $                   IRCOP( LDST, LDST ), LI( LDST, LDST ),
     $                   LICOP( LDST, LDST ), S( LDST, LDST ),
     $                   SCPY( LDST, LDST ), T( LDST, LDST ),
     $                   TAUL( LDST ), TAUR( LDST ), TCPY( LDST, LDST )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      DOUBLE PRECISION   <a name="DLAMCH.158"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>
      EXTERNAL           <a name="DLAMCH.159"></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           DGEMM, <a name="DGEQR2.162"></a><a href="dgeqr2.f.html#DGEQR2.1">DGEQR2</a>, <a name="DGERQ2.162"></a><a href="dgerq2.f.html#DGERQ2.1">DGERQ2</a>, <a name="DLACPY.162"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>, <a name="DLAGV2.162"></a><a href="dlagv2.f.html#DLAGV2.1">DLAGV2</a>, <a name="DLARTG.162"></a><a href="dlartg.f.html#DLARTG.1">DLARTG</a>,
     $                   <a name="DLASET.163"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>, <a name="DLASSQ.163"></a><a href="dlassq.f.html#DLASSQ.1">DLASSQ</a>, <a name="DORG2R.163"></a><a href="dorg2r.f.html#DORG2R.1">DORG2R</a>, <a name="DORGR2.163"></a><a href="dorgr2.f.html#DORGR2.1">DORGR2</a>, <a name="DORM2R.163"></a><a href="dorm2r.f.html#DORM2R.1">DORM2R</a>, <a name="DORMR2.163"></a><a href="dormr2.f.html#DORMR2.1">DORMR2</a>,
     $                   DROT, DSCAL, <a name="DTGSY2.164"></a><a href="dtgsy2.f.html#DTGSY2.1">DTGSY2</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, 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 .OR. N1.LE.0 .OR. N2.LE.0 )
     $   RETURN
      IF( N1.GT.N .OR. ( J1+N1 ).GT.N )
     $   RETURN
      M = N1 + N2
      IF( LWORK.LT.MAX( 1, N*M, M*M*2 ) ) THEN
         INFO = -16
         WORK( 1 ) = MAX( 1, N*M, M*M*2 )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      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
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="DLASET.191"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'Full'</span>, LDST, LDST, ZERO, ZERO, LI, LDST )
      CALL <a name="DLASET.192"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'Full'</span>, LDST, LDST, ZERO, ZERO, IR, LDST )
      CALL <a name="DLACPY.193"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'Full'</span>, M, M, A( J1, J1 ), LDA, S, LDST )
      CALL <a name="DLACPY.194"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</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 threshold for testing acceptance of swapping.
</span><span class="comment">*</span><span class="comment">
</span>      EPS = <a name="DLAMCH.198"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'P'</span> )
      SMLNUM = <a name="DLAMCH.199"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'S'</span> ) / EPS
      DSCALE = ZERO
      DSUM = ONE
      CALL <a name="DLACPY.202"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'Full'</span>, M, M, S, LDST, WORK, M )
      CALL <a name="DLASSQ.203"></a><a href="dlassq.f.html#DLASSQ.1">DLASSQ</a>( M*M, WORK, 1, DSCALE, DSUM )
      CALL <a name="DLACPY.204"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'Full'</span>, M, M, T, LDST, WORK, M )
      CALL <a name="DLASSQ.205"></a><a href="dlassq.f.html#DLASSQ.1">DLASSQ</a>( M*M, WORK, 1, DSCALE, DSUM )
      DNORM = DSCALE*SQRT( DSUM )
      THRESH = MAX( TEN*EPS*DNORM, SMLNUM )
<span class="comment">*</span><span class="comment">
</span>      IF( M.EQ.2 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        CASE 1: Swap 1-by-1 and 1-by-1 blocks.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute orthogonal 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 )
         SB = ABS( T( 2, 2 ) )
         SA = ABS( S( 2, 2 ) )
         CALL <a name="DLARTG.220"></a><a href="dlartg.f.html#DLARTG.1">DLARTG</a>( F, G, IR( 1, 2 ), IR( 1, 1 ), DDUM )
         IR( 2, 1 ) = -IR( 1, 2 )
         IR( 2, 2 ) = IR( 1, 1 )
         CALL DROT( 2, S( 1, 1 ), 1, S( 1, 2 ), 1, IR( 1, 1 ),
     $              IR( 2, 1 ) )
         CALL DROT( 2, T( 1, 1 ), 1, T( 1, 2 ), 1, IR( 1, 1 ),
     $              IR( 2, 1 ) )
         IF( SA.GE.SB ) THEN
            CALL <a name="DLARTG.228"></a><a href="dlartg.f.html#DLARTG.1">DLARTG</a>( S( 1, 1 ), S( 2, 1 ), LI( 1, 1 ), LI( 2, 1 ),
     $                   DDUM )
         ELSE
            CALL <a name="DLARTG.231"></a><a href="dlartg.f.html#DLARTG.1">DLARTG</a>( T( 1, 1 ), T( 2, 1 ), LI( 1, 1 ), LI( 2, 1 ),
     $                   DDUM )
         END IF
         CALL DROT( 2, S( 1, 1 ), LDST, S( 2, 1 ), LDST, LI( 1, 1 ),
     $              LI( 2, 1 ) )
         CALL DROT( 2, T( 1, 1 ), LDST, T( 2, 1 ), LDST, LI( 1, 1 ),
     $              LI( 2, 1 ) )
         LI( 2, 2 ) = LI( 1, 1 )
         LI( 1, 2 ) = -LI( 2, 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Weak stability test:
</span><span class="comment">*</span><span class="comment">           |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 70
<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="DLACPY.254"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'Full'</span>, M, M, A( J1, J1 ), LDA, WORK( M*M+1 ),
     $                   M )
            CALL DGEMM( <span class="string">'N'</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, ONE,
     $                  WORK( M*M+1 ), M )
            DSCALE = ZERO
            DSUM = ONE
            CALL <a name="DLASSQ.262"></a><a href="dlassq.f.html#DLASSQ.1">DLASSQ</a>( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )
<span class="comment">*</span><span class="comment">
</span>            CALL <a name="DLACPY.264"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'Full'</span>, M, M, B( J1, J1 ), LDB, WORK( M*M+1 ),
     $                   M )
            CALL DGEMM( <span class="string">'N'</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, ONE,
     $                  WORK( M*M+1 ), M )
            CALL <a name="DLASSQ.270"></a><a href="dlassq.f.html#DLASSQ.1">DLASSQ</a>( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )
            SS = DSCALE*SQRT( DSUM )
            DTRONG = SS.LE.THRESH
            IF( .NOT.DTRONG )
     $         GO TO 70
         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>         CALL DROT( J1+1, A( 1, J1 ), 1, A( 1, J1+1 ), 1, IR( 1, 1 ),
     $              IR( 2, 1 ) )
         CALL DROT( J1+1, B( 1, J1 ), 1, B( 1, J1+1 ), 1, IR( 1, 1 ),
     $              IR( 2, 1 ) )
         CALL DROT( N-J1+1, A( J1, J1 ), LDA, A( J1+1, J1 ), LDA,
     $              LI( 1, 1 ), LI( 2, 1 ) )