[BACK]Return to dtgex2.f CVS log [TXT] [DIR] Up to [local] / rpl / lapack / lapack

Diff for /rpl/lapack/lapack/dtgex2.f between versions 1.20 and 1.21

version 1.20, 2018/05/29 07:18:10 version 1.21, 2023/08/07 08:39:12
Line 181 Line 181
 *> \author Univ. of Colorado Denver  *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.  *> \author NAG Ltd.
 *  *
 *> \date December 2016  
 *  
 *> \ingroup doubleGEauxiliary  *> \ingroup doubleGEauxiliary
 *  *
 *> \par Further Details:  *> \par Further Details:
Line 221 Line 219
       SUBROUTINE DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,        SUBROUTINE DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
      $                   LDZ, J1, N1, N2, WORK, LWORK, INFO )       $                   LDZ, J1, N1, N2, WORK, LWORK, INFO )
 *  *
 *  -- LAPACK auxiliary routine (version 3.7.0) --  *  -- LAPACK auxiliary routine --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 *     December 2016  
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       LOGICAL            WANTQ, WANTZ        LOGICAL            WANTQ, WANTZ
Line 250 Line 247
       PARAMETER          ( WANDS = .TRUE. )        PARAMETER          ( WANDS = .TRUE. )
 *     ..  *     ..
 *     .. Local Scalars ..  *     .. Local Scalars ..
       LOGICAL            DTRONG, WEAK        LOGICAL            STRONG, WEAK
       INTEGER            I, IDUM, LINFO, M        INTEGER            I, IDUM, LINFO, M
       DOUBLE PRECISION   BQRA21, BRQA21, DDUM, DNORM, DSCALE, DSUM, EPS,        DOUBLE PRECISION   BQRA21, BRQA21, DDUM, DNORMA, DNORMB, DSCALE,
      $                   F, G, SA, SB, SCALE, SMLNUM, SS, THRESH, WS       $                   DSUM, EPS, F, G, SA, SB, SCALE, SMLNUM,
        $                   THRESHA, THRESHB
 *     ..  *     ..
 *     .. Local Arrays ..  *     .. Local Arrays ..
       INTEGER            IWORK( LDST )        INTEGER            IWORK( LDST )
Line 293 Line 291
       END IF        END IF
 *  *
       WEAK = .FALSE.        WEAK = .FALSE.
       DTRONG = .FALSE.        STRONG = .FALSE.
 *  *
 *     Make a local copy of selected block  *     Make a local copy of selected block
 *  *
Line 310 Line 308
       DSUM = ONE        DSUM = ONE
       CALL DLACPY( 'Full', M, M, S, LDST, WORK, M )        CALL DLACPY( 'Full', M, M, S, LDST, WORK, M )
       CALL DLASSQ( M*M, WORK, 1, DSCALE, DSUM )        CALL DLASSQ( M*M, WORK, 1, DSCALE, DSUM )
         DNORMA = DSCALE*SQRT( DSUM )
         DSCALE = ZERO
         DSUM = ONE
       CALL DLACPY( 'Full', M, M, T, LDST, WORK, M )        CALL DLACPY( 'Full', M, M, T, LDST, WORK, M )
       CALL DLASSQ( M*M, WORK, 1, DSCALE, DSUM )        CALL DLASSQ( M*M, WORK, 1, DSCALE, DSUM )
       DNORM = DSCALE*SQRT( DSUM )        DNORMB = DSCALE*SQRT( DSUM )
 *  *
 *     THRES has been changed from  *     THRES has been changed from
 *        THRESH = MAX( TEN*EPS*SA, SMLNUM )  *        THRESH = MAX( TEN*EPS*SA, SMLNUM )
Line 322 Line 323
 *     "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by  *     "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
 *     Jim Demmel and Guillaume Revy. See forum post 1783.  *     Jim Demmel and Guillaume Revy. See forum post 1783.
 *  *
       THRESH = MAX( TWENTY*EPS*DNORM, SMLNUM )        THRESHA = MAX( TWENTY*EPS*DNORMA, SMLNUM )
         THRESHB = MAX( TWENTY*EPS*DNORMB, SMLNUM )
 *  *
       IF( M.EQ.2 ) THEN        IF( M.EQ.2 ) THEN
 *  *
Line 333 Line 335
 *  *
          F = S( 2, 2 )*T( 1, 1 ) - T( 2, 2 )*S( 1, 1 )           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 )           G = S( 2, 2 )*T( 1, 2 ) - T( 2, 2 )*S( 1, 2 )
          SB = ABS( T( 2, 2 ) )           SA = ABS( S( 2, 2 ) ) * ABS( T( 1, 1 ) )
          SA = ABS( S( 2, 2 ) )           SB = ABS( S( 1, 1 ) ) * ABS( T( 2, 2 ) )
          CALL DLARTG( F, G, IR( 1, 2 ), IR( 1, 1 ), DDUM )           CALL DLARTG( F, G, IR( 1, 2 ), IR( 1, 1 ), DDUM )
          IR( 2, 1 ) = -IR( 1, 2 )           IR( 2, 1 ) = -IR( 1, 2 )
          IR( 2, 2 ) = IR( 1, 1 )           IR( 2, 2 ) = IR( 1, 1 )
Line 356 Line 358
          LI( 2, 2 ) = LI( 1, 1 )           LI( 2, 2 ) = LI( 1, 1 )
          LI( 1, 2 ) = -LI( 2, 1 )           LI( 1, 2 ) = -LI( 2, 1 )
 *  *
 *        Weak stability test:  *        Weak stability test: |S21| <= O(EPS F-norm((A)))
 *           |S21| + |T21| <= O(EPS * F-norm((S, T)))  *                           and  |T21| <= O(EPS F-norm((B)))
 *  *
          WS = ABS( S( 2, 1 ) ) + ABS( T( 2, 1 ) )           WEAK = ABS( S( 2, 1 ) ) .LE. THRESHA .AND.
          WEAK = WS.LE.THRESH       $      ABS( T( 2, 1 ) ) .LE. THRESHB
          IF( .NOT.WEAK )           IF( .NOT.WEAK )
      $      GO TO 70       $      GO TO 70
 *  *
          IF( WANDS ) THEN           IF( WANDS ) THEN
 *  *
 *           Strong stability test:  *           Strong stability test:
 *             F-norm((A-QL**T*S*QR, B-QL**T*T*QR)) <= O(EPS*F-norm((A,B)))  *               F-norm((A-QL**H*S*QR)) <= O(EPS*F-norm((A)))
   *               and
   *               F-norm((B-QL**H*T*QR)) <= O(EPS*F-norm((B)))
 *  *
             CALL DLACPY( 'Full', M, M, A( J1, J1 ), LDA, WORK( M*M+1 ),              CALL DLACPY( 'Full', M, M, A( J1, J1 ), LDA, WORK( M*M+1 ),
      $                   M )       $                   M )
Line 378 Line 382
             DSCALE = ZERO              DSCALE = ZERO
             DSUM = ONE              DSUM = ONE
             CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )              CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )
               SA = DSCALE*SQRT( DSUM )
 *  *
             CALL DLACPY( 'Full', M, M, B( J1, J1 ), LDB, WORK( M*M+1 ),              CALL DLACPY( 'Full', M, M, B( J1, J1 ), LDB, WORK( M*M+1 ),
      $                   M )       $                   M )
Line 385 Line 390
      $                  WORK, M )       $                  WORK, M )
             CALL DGEMM( 'N', 'T', M, M, M, -ONE, WORK, M, IR, LDST, ONE,              CALL DGEMM( 'N', 'T', M, M, M, -ONE, WORK, M, IR, LDST, ONE,
      $                  WORK( M*M+1 ), M )       $                  WORK( M*M+1 ), M )
               DSCALE = ZERO
               DSUM = ONE
             CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )              CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )
             SS = DSCALE*SQRT( DSUM )              SB = DSCALE*SQRT( DSUM )
             DTRONG = SS.LE.THRESH              STRONG = SA.LE.THRESHA .AND. SB.LE.THRESHB
             IF( .NOT.DTRONG )              IF( .NOT.STRONG )
      $         GO TO 70       $         GO TO 70
          END IF           END IF
 *  *
Line 439 Line 446
      $                IR( N2+1, N1+1 ), LDST, T, LDST, T( N1+1, N1+1 ),       $                IR( N2+1, N1+1 ), LDST, T, LDST, T( N1+1, N1+1 ),
      $                LDST, LI, LDST, SCALE, DSUM, DSCALE, IWORK, IDUM,       $                LDST, LI, LDST, SCALE, DSUM, DSCALE, IWORK, IDUM,
      $                LINFO )       $                LINFO )
            IF( LINFO.NE.0 )
        $      GO TO 70
 *  *
 *        Compute orthogonal matrix QL:  *        Compute orthogonal matrix QL:
 *  *
Line 538 Line 547
 *  *
 *        Decide which method to use.  *        Decide which method to use.
 *          Weak stability test:  *          Weak stability test:
 *             F-norm(S21) <= O(EPS * F-norm((S, T)))  *             F-norm(S21) <= O(EPS * F-norm((S)))
 *  *
          IF( BQRA21.LE.BRQA21 .AND. BQRA21.LE.THRESH ) THEN           IF( BQRA21.LE.BRQA21 .AND. BQRA21.LE.THRESHA ) THEN
             CALL DLACPY( 'F', M, M, SCPY, LDST, S, LDST )              CALL DLACPY( 'F', M, M, SCPY, LDST, S, LDST )
             CALL DLACPY( 'F', M, M, TCPY, LDST, T, LDST )              CALL DLACPY( 'F', M, M, TCPY, LDST, T, LDST )
             CALL DLACPY( 'F', M, M, IRCOP, LDST, IR, LDST )              CALL DLACPY( 'F', M, M, IRCOP, LDST, IR, LDST )
             CALL DLACPY( 'F', M, M, LICOP, LDST, LI, LDST )              CALL DLACPY( 'F', M, M, LICOP, LDST, LI, LDST )
          ELSE IF( BRQA21.GE.THRESH ) THEN           ELSE IF( BRQA21.GE.THRESHA ) THEN
             GO TO 70              GO TO 70
          END IF           END IF
 *  *
Line 556 Line 565
          IF( WANDS ) THEN           IF( WANDS ) THEN
 *  *
 *           Strong stability test:  *           Strong stability test:
 *              F-norm((A-QL*S*QR**T, B-QL*T*QR**T)) <= O(EPS*F-norm((A,B)))  *               F-norm((A-QL**H*S*QR)) <= O(EPS*F-norm((A)))
   *               and
   *               F-norm((B-QL**H*T*QR)) <= O(EPS*F-norm((B)))
 *  *
             CALL DLACPY( 'Full', M, M, A( J1, J1 ), LDA, WORK( M*M+1 ),              CALL DLACPY( 'Full', M, M, A( J1, J1 ), LDA, WORK( M*M+1 ),
      $                   M )       $                   M )
Line 567 Line 578
             DSCALE = ZERO              DSCALE = ZERO
             DSUM = ONE              DSUM = ONE
             CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )              CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )
               SA = DSCALE*SQRT( DSUM )
 *  *
             CALL DLACPY( 'Full', M, M, B( J1, J1 ), LDB, WORK( M*M+1 ),              CALL DLACPY( 'Full', M, M, B( J1, J1 ), LDB, WORK( M*M+1 ),
      $                   M )       $                   M )
Line 574 Line 586
      $                  WORK, M )       $                  WORK, M )
             CALL DGEMM( 'N', 'N', M, M, M, -ONE, WORK, M, IR, LDST, ONE,              CALL DGEMM( 'N', 'N', M, M, M, -ONE, WORK, M, IR, LDST, ONE,
      $                  WORK( M*M+1 ), M )       $                  WORK( M*M+1 ), M )
               DSCALE = ZERO
               DSUM = ONE
             CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )              CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM )
             SS = DSCALE*SQRT( DSUM )              SB = DSCALE*SQRT( DSUM )
             DTRONG = ( SS.LE.THRESH )              STRONG = SA.LE.THRESHA .AND. SB.LE.THRESHB
             IF( .NOT.DTRONG )              IF( .NOT.STRONG )
      $         GO TO 70       $         GO TO 70
 *  *
          END IF           END IF
Removed from v.1.20  
changed lines
  Added in v.1.21

CVSweb interface <joel.bertrand@systella.fr>