SRC/dlags2.f

   1 *> \brief \b DLAGS2 computes 2-by-2 orthogonal matrices U, V, and Q, and applies them to matrices A and B such that the rows of the transformed A and B are parallel.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DLAGS2 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlags2.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlags2.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlags2.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLAGS2( UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, CSV,
  22 *                          SNV, CSQ, SNQ )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       LOGICAL            UPPER
  26 *       DOUBLE PRECISION   A1, A2, A3, B1, B2, B3, CSQ, CSU, CSV, SNQ,
  27 *      $                   SNU, SNV
  28 *       ..
  29 *
  30 *
  31 *> \par Purpose:
  32 *  =============
  33 *>
  34 *> \verbatim
  35 *>
  36 *> DLAGS2 computes 2-by-2 orthogonal matrices U, V and Q, such
  37 *> that if ( UPPER ) then
  38 *>
  39 *>           U**T *A*Q = U**T *( A1 A2 )*Q = ( x  0  )
  40 *>                             ( 0  A3 )     ( x  x  )
  41 *> and
  42 *>           V**T*B*Q = V**T *( B1 B2 )*Q = ( x  0  )
  43 *>                            ( 0  B3 )     ( x  x  )
  44 *>
  45 *> or if ( .NOT.UPPER ) then
  46 *>
  47 *>           U**T *A*Q = U**T *( A1 0  )*Q = ( x  x  )
  48 *>                             ( A2 A3 )     ( 0  x  )
  49 *> and
  50 *>           V**T*B*Q = V**T*( B1 0  )*Q = ( x  x  )
  51 *>                           ( B2 B3 )     ( 0  x  )
  52 *>
  53 *> The rows of the transformed A and B are parallel, where
  54 *>
  55 *>   U = (  CSU  SNU ), V = (  CSV SNV ), Q = (  CSQ   SNQ )
  56 *>       ( -SNU  CSU )      ( -SNV CSV )      ( -SNQ   CSQ )
  57 *>
  58 *> Z**T denotes the transpose of Z.
  59 *>
  60 *> \endverbatim
  61 *
  62 *  Arguments:
  63 *  ==========
  64 *
  65 *> \param[in] UPPER
  66 *> \verbatim
  67 *>          UPPER is LOGICAL
  68 *>          = .TRUE.: the input matrices A and B are upper triangular.
  69 *>          = .FALSE.: the input matrices A and B are lower triangular.
  70 *> \endverbatim
  71 *>
  72 *> \param[in] A1
  73 *> \verbatim
  74 *>          A1 is DOUBLE PRECISION
  75 *> \endverbatim
  76 *>
  77 *> \param[in] A2
  78 *> \verbatim
  79 *>          A2 is DOUBLE PRECISION
  80 *> \endverbatim
  81 *>
  82 *> \param[in] A3
  83 *> \verbatim
  84 *>          A3 is DOUBLE PRECISION
  85 *>          On entry, A1, A2 and A3 are elements of the input 2-by-2
  86 *>          upper (lower) triangular matrix A.
  87 *> \endverbatim
  88 *>
  89 *> \param[in] B1
  90 *> \verbatim
  91 *>          B1 is DOUBLE PRECISION
  92 *> \endverbatim
  93 *>
  94 *> \param[in] B2
  95 *> \verbatim
  96 *>          B2 is DOUBLE PRECISION
  97 *> \endverbatim
  98 *>
  99 *> \param[in] B3
 100 *> \verbatim
 101 *>          B3 is DOUBLE PRECISION
 102 *>          On entry, B1, B2 and B3 are elements of the input 2-by-2
 103 *>          upper (lower) triangular matrix B.
 104 *> \endverbatim
 105 *>
 106 *> \param[out] CSU
 107 *> \verbatim
 108 *>          CSU is DOUBLE PRECISION
 109 *> \endverbatim
 110 *>
 111 *> \param[out] SNU
 112 *> \verbatim
 113 *>          SNU is DOUBLE PRECISION
 114 *>          The desired orthogonal matrix U.
 115 *> \endverbatim
 116 *>
 117 *> \param[out] CSV
 118 *> \verbatim
 119 *>          CSV is DOUBLE PRECISION
 120 *> \endverbatim
 121 *>
 122 *> \param[out] SNV
 123 *> \verbatim
 124 *>          SNV is DOUBLE PRECISION
 125 *>          The desired orthogonal matrix V.
 126 *> \endverbatim
 127 *>
 128 *> \param[out] CSQ
 129 *> \verbatim
 130 *>          CSQ is DOUBLE PRECISION
 131 *> \endverbatim
 132 *>
 133 *> \param[out] SNQ
 134 *> \verbatim
 135 *>          SNQ is DOUBLE PRECISION
 136 *>          The desired orthogonal matrix Q.
 137 *> \endverbatim
 138 *
 139 *  Authors:
 140 *  ========
 141 *
 142 *> \author Univ. of Tennessee
 143 *> \author Univ. of California Berkeley
 144 *> \author Univ. of Colorado Denver
 145 *> \author NAG Ltd.
 146 *
 147 *> \date September 2012
 148 *
 149 *> \ingroup doubleOTHERauxiliary
 150 *
 151 *  =====================================================================
 152       SUBROUTINE DLAGS2( UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, CSV,
 153      $                   SNV, CSQ, SNQ )
 154 *
 155 *  -- LAPACK auxiliary routine (version 3.4.2) --
 156 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 157 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 158 *     September 2012
 159 *
 160 *     .. Scalar Arguments ..
 161       LOGICAL            UPPER
 162       DOUBLE PRECISION   A1, A2, A3, B1, B2, B3, CSQ, CSU, CSV, SNQ,
 163      $                   SNU, SNV
 164 *     ..
 165 *
 166 *  =====================================================================
 167 *
 168 *     .. Parameters ..
 169       DOUBLE PRECISION   ZERO
 170       PARAMETER          ( ZERO = 0.0D+0 )
 171 *     ..
 172 *     .. Local Scalars ..
 173       DOUBLE PRECISION   A, AUA11, AUA12, AUA21, AUA22, AVB11, AVB12,
 174      $                   AVB21, AVB22, B, C, CSL, CSR, D, R, S1, S2,
 175      $                   SNL, SNR, UA11, UA11R, UA12, UA21, UA22, UA22R,
 176      $                   VB11, VB11R, VB12, VB21, VB22, VB22R
 177 *     ..
 178 *     .. External Subroutines ..
 179       EXTERNAL           DLARTG, DLASV2
 180 *     ..
 181 *     .. Intrinsic Functions ..
 182       INTRINSIC          ABS
 183 *     ..
 184 *     .. Executable Statements ..
 185 *
 186       IF( UPPER ) THEN
 187 *
 188 *        Input matrices A and B are upper triangular matrices
 189 *
 190 *        Form matrix C = A*adj(B) = ( a b )
 191 *                                   ( 0 d )
 192 *
 193          A = A1*B3
 194          D = A3*B1
 195          B = A2*B1 - A1*B2
 196 *
 197 *        The SVD of real 2-by-2 triangular C
 198 *
 199 *         ( CSL -SNL )*( A B )*(  CSR  SNR ) = ( R 0 )
 200 *         ( SNL  CSL ) ( 0 D ) ( -SNR  CSR )   ( 0 T )
 201 *
 202          CALL DLASV2( A, B, D, S1, S2, SNR, CSR, SNL, CSL )
 203 *
 204          IF( ABS( CSL ).GE.ABS( SNL ) .OR. ABS( CSR ).GE.ABS( SNR ) )
 205      $        THEN
 206 *
 207 *           Compute the (1,1) and (1,2) elements of U**T *A and V**T *B,
 208 *           and (1,2) element of |U|**T *|A| and |V|**T *|B|.
 209 *
 210             UA11R = CSL*A1
 211             UA12 = CSL*A2 + SNL*A3
 212 *
 213             VB11R = CSR*B1
 214             VB12 = CSR*B2 + SNR*B3
 215 *
 216             AUA12 = ABS( CSL )*ABS( A2 ) + ABS( SNL )*ABS( A3 )
 217             AVB12 = ABS( CSR )*ABS( B2 ) + ABS( SNR )*ABS( B3 )
 218 *
 219 *           zero (1,2) elements of U**T *A and V**T *B
 220 *
 221             IF( ( ABS( UA11R )+ABS( UA12 ) ).NE.ZERO ) THEN
 222                IF( AUA12 / ( ABS( UA11R )+ABS( UA12 ) ).LE.AVB12 /
 223      $             ( ABS( VB11R )+ABS( VB12 ) ) ) THEN
 224                   CALL DLARTG( -UA11R, UA12, CSQ, SNQ, R )
 225                ELSE
 226                   CALL DLARTG( -VB11R, VB12, CSQ, SNQ, R )
 227                END IF
 228             ELSE
 229                CALL DLARTG( -VB11R, VB12, CSQ, SNQ, R )
 230             END IF
 231 *
 232             CSU = CSL
 233             SNU = -SNL
 234             CSV = CSR
 235             SNV = -SNR
 236 *
 237          ELSE
 238 *
 239 *           Compute the (2,1) and (2,2) elements of U**T *A and V**T *B,
 240 *           and (2,2) element of |U|**T *|A| and |V|**T *|B|.
 241 *
 242             UA21 = -SNL*A1
 243             UA22 = -SNL*A2 + CSL*A3
 244 *
 245             VB21 = -SNR*B1
 246             VB22 = -SNR*B2 + CSR*B3
 247 *
 248             AUA22 = ABS( SNL )*ABS( A2 ) + ABS( CSL )*ABS( A3 )
 249             AVB22 = ABS( SNR )*ABS( B2 ) + ABS( CSR )*ABS( B3 )
 250 *
 251 *           zero (2,2) elements of U**T*A and V**T*B, and then swap.
 252 *
 253             IF( ( ABS( UA21 )+ABS( UA22 ) ).NE.ZERO ) THEN
 254                IF( AUA22 / ( ABS( UA21 )+ABS( UA22 ) ).LE.AVB22 /
 255      $             ( ABS( VB21 )+ABS( VB22 ) ) ) THEN
 256                   CALL DLARTG( -UA21, UA22, CSQ, SNQ, R )
 257                ELSE
 258                   CALL DLARTG( -VB21, VB22, CSQ, SNQ, R )
 259                END IF
 260             ELSE
 261                CALL DLARTG( -VB21, VB22, CSQ, SNQ, R )
 262             END IF
 263 *
 264             CSU = SNL
 265             SNU = CSL
 266             CSV = SNR
 267             SNV = CSR
 268 *
 269          END IF
 270 *
 271       ELSE
 272 *
 273 *        Input matrices A and B are lower triangular matrices
 274 *
 275 *        Form matrix C = A*adj(B) = ( a 0 )
 276 *                                   ( c d )
 277 *
 278          A = A1*B3
 279          D = A3*B1
 280          C = A2*B3 - A3*B2
 281 *
 282 *        The SVD of real 2-by-2 triangular C
 283 *
 284 *         ( CSL -SNL )*( A 0 )*(  CSR  SNR ) = ( R 0 )
 285 *         ( SNL  CSL ) ( C D ) ( -SNR  CSR )   ( 0 T )
 286 *
 287          CALL DLASV2( A, C, D, S1, S2, SNR, CSR, SNL, CSL )
 288 *
 289          IF( ABS( CSR ).GE.ABS( SNR ) .OR. ABS( CSL ).GE.ABS( SNL ) )
 290      $        THEN
 291 *
 292 *           Compute the (2,1) and (2,2) elements of U**T *A and V**T *B,
 293 *           and (2,1) element of |U|**T *|A| and |V|**T *|B|.
 294 *
 295             UA21 = -SNR*A1 + CSR*A2
 296             UA22R = CSR*A3
 297 *
 298             VB21 = -SNL*B1 + CSL*B2
 299             VB22R = CSL*B3
 300 *
 301             AUA21 = ABS( SNR )*ABS( A1 ) + ABS( CSR )*ABS( A2 )
 302             AVB21 = ABS( SNL )*ABS( B1 ) + ABS( CSL )*ABS( B2 )
 303 *
 304 *           zero (2,1) elements of U**T *A and V**T *B.
 305 *
 306             IF( ( ABS( UA21 )+ABS( UA22R ) ).NE.ZERO ) THEN
 307                IF( AUA21 / ( ABS( UA21 )+ABS( UA22R ) ).LE.AVB21 /
 308      $             ( ABS( VB21 )+ABS( VB22R ) ) ) THEN
 309                   CALL DLARTG( UA22R, UA21, CSQ, SNQ, R )
 310                ELSE
 311                   CALL DLARTG( VB22R, VB21, CSQ, SNQ, R )
 312                END IF
 313             ELSE
 314                CALL DLARTG( VB22R, VB21, CSQ, SNQ, R )
 315             END IF
 316 *
 317             CSU = CSR
 318             SNU = -SNR
 319             CSV = CSL
 320             SNV = -SNL
 321 *
 322          ELSE
 323 *
 324 *           Compute the (1,1) and (1,2) elements of U**T *A and V**T *B,
 325 *           and (1,1) element of |U|**T *|A| and |V|**T *|B|.
 326 *
 327             UA11 = CSR*A1 + SNR*A2
 328             UA12 = SNR*A3
 329 *
 330             VB11 = CSL*B1 + SNL*B2
 331             VB12 = SNL*B3
 332 *
 333             AUA11 = ABS( CSR )*ABS( A1 ) + ABS( SNR )*ABS( A2 )
 334             AVB11 = ABS( CSL )*ABS( B1 ) + ABS( SNL )*ABS( B2 )
 335 *
 336 *           zero (1,1) elements of U**T*A and V**T*B, and then swap.
 337 *
 338             IF( ( ABS( UA11 )+ABS( UA12 ) ).NE.ZERO ) THEN
 339                IF( AUA11 / ( ABS( UA11 )+ABS( UA12 ) ).LE.AVB11 /
 340      $             ( ABS( VB11 )+ABS( VB12 ) ) ) THEN
 341                   CALL DLARTG( UA12, UA11, CSQ, SNQ, R )
 342                ELSE
 343                   CALL DLARTG( VB12, VB11, CSQ, SNQ, R )
 344                END IF
 345             ELSE
 346                CALL DLARTG( VB12, VB11, CSQ, SNQ, R )
 347             END IF
 348 *
 349             CSU = SNR
 350             SNU = CSR
 351             CSV = SNL
 352             SNV = CSL
 353 *
 354          END IF
 355 *
 356       END IF
 357 *
 358       RETURN
 359 *
 360 *     End of DLAGS2
 361 *
 362       END