Return to dlags2.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

    1:       SUBROUTINE DLAGS2( UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, CSV,
    2:      $                   SNV, CSQ, SNQ )
    3: *
    4: *  -- LAPACK auxiliary routine (version 3.2) --
    5: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    6: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    7: *     November 2006
    8: *
    9: *     .. Scalar Arguments ..
   10:       LOGICAL            UPPER
   11:       DOUBLE PRECISION   A1, A2, A3, B1, B2, B3, CSQ, CSU, CSV, SNQ,
   12:      $                   SNU, SNV
   13: *     ..
   14: *
   15: *  Purpose
   16: *  =======
   17: *
   18: *  DLAGS2 computes 2-by-2 orthogonal matrices U, V and Q, such
   19: *  that if ( UPPER ) then
   20: *
   21: *            U'*A*Q = U'*( A1 A2 )*Q = ( x  0  )
   22: *                        ( 0  A3 )     ( x  x  )
   23: *  and
   24: *            V'*B*Q = V'*( B1 B2 )*Q = ( x  0  )
   25: *                        ( 0  B3 )     ( x  x  )
   26: *
   27: *  or if ( .NOT.UPPER ) then
   28: *
   29: *            U'*A*Q = U'*( A1 0  )*Q = ( x  x  )
   30: *                        ( A2 A3 )     ( 0  x  )
   31: *  and
   32: *            V'*B*Q = V'*( B1 0  )*Q = ( x  x  )
   33: *                        ( B2 B3 )     ( 0  x  )
   34: *
   35: *  The rows of the transformed A and B are parallel, where
   36: *
   37: *    U = (  CSU  SNU ), V = (  CSV SNV ), Q = (  CSQ   SNQ )
   38: *        ( -SNU  CSU )      ( -SNV CSV )      ( -SNQ   CSQ )
   39: *
   40: *  Z' denotes the transpose of Z.
   41: *
   42: *
   43: *  Arguments
   44: *  =========
   45: *
   46: *  UPPER   (input) LOGICAL
   47: *          = .TRUE.: the input matrices A and B are upper triangular.
   48: *          = .FALSE.: the input matrices A and B are lower triangular.
   49: *
   50: *  A1      (input) DOUBLE PRECISION
   51: *  A2      (input) DOUBLE PRECISION
   52: *  A3      (input) DOUBLE PRECISION
   53: *          On entry, A1, A2 and A3 are elements of the input 2-by-2
   54: *          upper (lower) triangular matrix A.
   55: *
   56: *  B1      (input) DOUBLE PRECISION
   57: *  B2      (input) DOUBLE PRECISION
   58: *  B3      (input) DOUBLE PRECISION
   59: *          On entry, B1, B2 and B3 are elements of the input 2-by-2
   60: *          upper (lower) triangular matrix B.
   61: *
   62: *  CSU     (output) DOUBLE PRECISION
   63: *  SNU     (output) DOUBLE PRECISION
   64: *          The desired orthogonal matrix U.
   65: *
   66: *  CSV     (output) DOUBLE PRECISION
   67: *  SNV     (output) DOUBLE PRECISION
   68: *          The desired orthogonal matrix V.
   69: *
   70: *  CSQ     (output) DOUBLE PRECISION
   71: *  SNQ     (output) DOUBLE PRECISION
   72: *          The desired orthogonal matrix Q.
   73: *
   74: *  =====================================================================
   75: *
   76: *     .. Parameters ..
   77:       DOUBLE PRECISION   ZERO
   78:       PARAMETER          ( ZERO = 0.0D+0 )
   79: *     ..
   80: *     .. Local Scalars ..
   81:       DOUBLE PRECISION   A, AUA11, AUA12, AUA21, AUA22, AVB11, AVB12,
   82:      $                   AVB21, AVB22, B, C, CSL, CSR, D, R, S1, S2,
   83:      $                   SNL, SNR, UA11, UA11R, UA12, UA21, UA22, UA22R,
   84:      $                   VB11, VB11R, VB12, VB21, VB22, VB22R
   85: *     ..
   86: *     .. External Subroutines ..
   87:       EXTERNAL           DLARTG, DLASV2
   88: *     ..
   89: *     .. Intrinsic Functions ..
   90:       INTRINSIC          ABS
   91: *     ..
   92: *     .. Executable Statements ..
   93: *
   94:       IF( UPPER ) THEN
   95: *
   96: *        Input matrices A and B are upper triangular matrices
   97: *
   98: *        Form matrix C = A*adj(B) = ( a b )
   99: *                                   ( 0 d )
  100: *
  101:          A = A1*B3
  102:          D = A3*B1
  103:          B = A2*B1 - A1*B2
  104: *
  105: *        The SVD of real 2-by-2 triangular C
  106: *
  107: *         ( CSL -SNL )*( A B )*(  CSR  SNR ) = ( R 0 )
  108: *         ( SNL  CSL ) ( 0 D ) ( -SNR  CSR )   ( 0 T )
  109: *
  110:          CALL DLASV2( A, B, D, S1, S2, SNR, CSR, SNL, CSL )
  111: *
  112:          IF( ABS( CSL ).GE.ABS( SNL ) .OR. ABS( CSR ).GE.ABS( SNR ) )
  113:      $        THEN
  114: *
  115: *           Compute the (1,1) and (1,2) elements of U'*A and V'*B,
  116: *           and (1,2) element of |U|'*|A| and |V|'*|B|.
  117: *
  118:             UA11R = CSL*A1
  119:             UA12 = CSL*A2 + SNL*A3
  120: *
  121:             VB11R = CSR*B1
  122:             VB12 = CSR*B2 + SNR*B3
  123: *
  124:             AUA12 = ABS( CSL )*ABS( A2 ) + ABS( SNL )*ABS( A3 )
  125:             AVB12 = ABS( CSR )*ABS( B2 ) + ABS( SNR )*ABS( B3 )
  126: *
  127: *           zero (1,2) elements of U'*A and V'*B
  128: *
  129:             IF( ( ABS( UA11R )+ABS( UA12 ) ).NE.ZERO ) THEN
  130:                IF( AUA12 / ( ABS( UA11R )+ABS( UA12 ) ).LE.AVB12 /
  131:      $             ( ABS( VB11R )+ABS( VB12 ) ) ) THEN
  132:                   CALL DLARTG( -UA11R, UA12, CSQ, SNQ, R )
  133:                ELSE
  134:                   CALL DLARTG( -VB11R, VB12, CSQ, SNQ, R )
  135:                END IF
  136:             ELSE
  137:                CALL DLARTG( -VB11R, VB12, CSQ, SNQ, R )
  138:             END IF
  139: *
  140:             CSU = CSL
  141:             SNU = -SNL
  142:             CSV = CSR
  143:             SNV = -SNR
  144: *
  145:          ELSE
  146: *
  147: *           Compute the (2,1) and (2,2) elements of U'*A and V'*B,
  148: *           and (2,2) element of |U|'*|A| and |V|'*|B|.
  149: *
  150:             UA21 = -SNL*A1
  151:             UA22 = -SNL*A2 + CSL*A3
  152: *
  153:             VB21 = -SNR*B1
  154:             VB22 = -SNR*B2 + CSR*B3
  155: *
  156:             AUA22 = ABS( SNL )*ABS( A2 ) + ABS( CSL )*ABS( A3 )
  157:             AVB22 = ABS( SNR )*ABS( B2 ) + ABS( CSR )*ABS( B3 )
  158: *
  159: *           zero (2,2) elements of U'*A and V'*B, and then swap.
  160: *
  161:             IF( ( ABS( UA21 )+ABS( UA22 ) ).NE.ZERO ) THEN
  162:                IF( AUA22 / ( ABS( UA21 )+ABS( UA22 ) ).LE.AVB22 /
  163:      $             ( ABS( VB21 )+ABS( VB22 ) ) ) THEN
  164:                   CALL DLARTG( -UA21, UA22, CSQ, SNQ, R )
  165:                ELSE
  166:                   CALL DLARTG( -VB21, VB22, CSQ, SNQ, R )
  167:                END IF
  168:             ELSE
  169:                CALL DLARTG( -VB21, VB22, CSQ, SNQ, R )
  170:             END IF
  171: *
  172:             CSU = SNL
  173:             SNU = CSL
  174:             CSV = SNR
  175:             SNV = CSR
  176: *
  177:          END IF
  178: *
  179:       ELSE
  180: *
  181: *        Input matrices A and B are lower triangular matrices
  182: *
  183: *        Form matrix C = A*adj(B) = ( a 0 )
  184: *                                   ( c d )
  185: *
  186:          A = A1*B3
  187:          D = A3*B1
  188:          C = A2*B3 - A3*B2
  189: *
  190: *        The SVD of real 2-by-2 triangular C
  191: *
  192: *         ( CSL -SNL )*( A 0 )*(  CSR  SNR ) = ( R 0 )
  193: *         ( SNL  CSL ) ( C D ) ( -SNR  CSR )   ( 0 T )
  194: *
  195:          CALL DLASV2( A, C, D, S1, S2, SNR, CSR, SNL, CSL )
  196: *
  197:          IF( ABS( CSR ).GE.ABS( SNR ) .OR. ABS( CSL ).GE.ABS( SNL ) )
  198:      $        THEN
  199: *
  200: *           Compute the (2,1) and (2,2) elements of U'*A and V'*B,
  201: *           and (2,1) element of |U|'*|A| and |V|'*|B|.
  202: *
  203:             UA21 = -SNR*A1 + CSR*A2
  204:             UA22R = CSR*A3
  205: *
  206:             VB21 = -SNL*B1 + CSL*B2
  207:             VB22R = CSL*B3
  208: *
  209:             AUA21 = ABS( SNR )*ABS( A1 ) + ABS( CSR )*ABS( A2 )
  210:             AVB21 = ABS( SNL )*ABS( B1 ) + ABS( CSL )*ABS( B2 )
  211: *
  212: *           zero (2,1) elements of U'*A and V'*B.
  213: *
  214:             IF( ( ABS( UA21 )+ABS( UA22R ) ).NE.ZERO ) THEN
  215:                IF( AUA21 / ( ABS( UA21 )+ABS( UA22R ) ).LE.AVB21 /
  216:      $             ( ABS( VB21 )+ABS( VB22R ) ) ) THEN
  217:                   CALL DLARTG( UA22R, UA21, CSQ, SNQ, R )
  218:                ELSE
  219:                   CALL DLARTG( VB22R, VB21, CSQ, SNQ, R )
  220:                END IF
  221:             ELSE
  222:                CALL DLARTG( VB22R, VB21, CSQ, SNQ, R )
  223:             END IF
  224: *
  225:             CSU = CSR
  226:             SNU = -SNR
  227:             CSV = CSL
  228:             SNV = -SNL
  229: *
  230:          ELSE
  231: *
  232: *           Compute the (1,1) and (1,2) elements of U'*A and V'*B,
  233: *           and (1,1) element of |U|'*|A| and |V|'*|B|.
  234: *
  235:             UA11 = CSR*A1 + SNR*A2
  236:             UA12 = SNR*A3
  237: *
  238:             VB11 = CSL*B1 + SNL*B2
  239:             VB12 = SNL*B3
  240: *
  241:             AUA11 = ABS( CSR )*ABS( A1 ) + ABS( SNR )*ABS( A2 )
  242:             AVB11 = ABS( CSL )*ABS( B1 ) + ABS( SNL )*ABS( B2 )
  243: *
  244: *           zero (1,1) elements of U'*A and V'*B, and then swap.
  245: *
  246:             IF( ( ABS( UA11 )+ABS( UA12 ) ).NE.ZERO ) THEN
  247:                IF( AUA11 / ( ABS( UA11 )+ABS( UA12 ) ).LE.AVB11 /
  248:      $             ( ABS( VB11 )+ABS( VB12 ) ) ) THEN
  249:                   CALL DLARTG( UA12, UA11, CSQ, SNQ, R )
  250:                ELSE
  251:                   CALL DLARTG( VB12, VB11, CSQ, SNQ, R )
  252:                END IF
  253:             ELSE
  254:                CALL DLARTG( VB12, VB11, CSQ, SNQ, R )
  255:             END IF
  256: *
  257:             CSU = SNR
  258:             SNU = CSR
  259:             CSV = SNL
  260:             SNV = CSL
  261: *
  262:          END IF
  263: *
  264:       END IF
  265: *
  266:       RETURN
  267: *
  268: *     End of DLAGS2
  269: *
  270:       END