LAPACK: SRC/dlags2.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DLAGS2( UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, CSV,
00002      $                   SNV, CSQ, SNQ )
00003 *
00004 *  -- LAPACK auxiliary routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       LOGICAL            UPPER
00011       DOUBLE PRECISION   A1, A2, A3, B1, B2, B3, CSQ, CSU, CSV, SNQ,
00012      $                   SNU, SNV
00013 *     ..
00014 *
00015 *  Purpose
00016 *  =======
00017 *
00018 *  DLAGS2 computes 2-by-2 orthogonal matrices U, V and Q, such
00019 *  that if ( UPPER ) then
00020 *
00021 *            U'*A*Q = U'*( A1 A2 )*Q = ( x  0  )
00022 *                        ( 0  A3 )     ( x  x  )
00023 *  and
00024 *            V'*B*Q = V'*( B1 B2 )*Q = ( x  0  )
00025 *                        ( 0  B3 )     ( x  x  )
00026 *
00027 *  or if ( .NOT.UPPER ) then
00028 *
00029 *            U'*A*Q = U'*( A1 0  )*Q = ( x  x  )
00030 *                        ( A2 A3 )     ( 0  x  )
00031 *  and
00032 *            V'*B*Q = V'*( B1 0  )*Q = ( x  x  )
00033 *                        ( B2 B3 )     ( 0  x  )
00034 *
00035 *  The rows of the transformed A and B are parallel, where
00036 *
00037 *    U = (  CSU  SNU ), V = (  CSV SNV ), Q = (  CSQ   SNQ )
00038 *        ( -SNU  CSU )      ( -SNV CSV )      ( -SNQ   CSQ )
00039 *
00040 *  Z' denotes the transpose of Z.
00041 *
00042 *
00043 *  Arguments
00044 *  =========
00045 *
00046 *  UPPER   (input) LOGICAL
00047 *          = .TRUE.: the input matrices A and B are upper triangular.
00048 *          = .FALSE.: the input matrices A and B are lower triangular.
00049 *
00050 *  A1      (input) DOUBLE PRECISION
00051 *  A2      (input) DOUBLE PRECISION
00052 *  A3      (input) DOUBLE PRECISION
00053 *          On entry, A1, A2 and A3 are elements of the input 2-by-2
00054 *          upper (lower) triangular matrix A.
00055 *
00056 *  B1      (input) DOUBLE PRECISION
00057 *  B2      (input) DOUBLE PRECISION
00058 *  B3      (input) DOUBLE PRECISION
00059 *          On entry, B1, B2 and B3 are elements of the input 2-by-2
00060 *          upper (lower) triangular matrix B.
00061 *
00062 *  CSU     (output) DOUBLE PRECISION
00063 *  SNU     (output) DOUBLE PRECISION
00064 *          The desired orthogonal matrix U.
00065 *
00066 *  CSV     (output) DOUBLE PRECISION
00067 *  SNV     (output) DOUBLE PRECISION
00068 *          The desired orthogonal matrix V.
00069 *
00070 *  CSQ     (output) DOUBLE PRECISION
00071 *  SNQ     (output) DOUBLE PRECISION
00072 *          The desired orthogonal matrix Q.
00073 *
00074 *  =====================================================================
00075 *
00076 *     .. Parameters ..
00077       DOUBLE PRECISION   ZERO
00078       PARAMETER          ( ZERO = 0.0D+0 )
00079 *     ..
00080 *     .. Local Scalars ..
00081       DOUBLE PRECISION   A, AUA11, AUA12, AUA21, AUA22, AVB11, AVB12,
00082      $                   AVB21, AVB22, B, C, CSL, CSR, D, R, S1, S2,
00083      $                   SNL, SNR, UA11, UA11R, UA12, UA21, UA22, UA22R,
00084      $                   VB11, VB11R, VB12, VB21, VB22, VB22R
00085 *     ..
00086 *     .. External Subroutines ..
00087       EXTERNAL           DLARTG, DLASV2
00088 *     ..
00089 *     .. Intrinsic Functions ..
00090       INTRINSIC          ABS
00091 *     ..
00092 *     .. Executable Statements ..
00093 *
00094       IF( UPPER ) THEN
00095 *
00096 *        Input matrices A and B are upper triangular matrices
00097 *
00098 *        Form matrix C = A*adj(B) = ( a b )
00099 *                                   ( 0 d )
00100 *
00101          A = A1*B3
00102          D = A3*B1
00103          B = A2*B1 - A1*B2
00104 *
00105 *        The SVD of real 2-by-2 triangular C
00106 *
00107 *         ( CSL -SNL )*( A B )*(  CSR  SNR ) = ( R 0 )
00108 *         ( SNL  CSL ) ( 0 D ) ( -SNR  CSR )   ( 0 T )
00109 *
00110          CALL DLASV2( A, B, D, S1, S2, SNR, CSR, SNL, CSL )
00111 *
00112          IF( ABS( CSL ).GE.ABS( SNL ) .OR. ABS( CSR ).GE.ABS( SNR ) )
00113      $        THEN
00114 *
00115 *           Compute the (1,1) and (1,2) elements of U'*A and V'*B,
00116 *           and (1,2) element of |U|'*|A| and |V|'*|B|.
00117 *
00118             UA11R = CSL*A1
00119             UA12 = CSL*A2 + SNL*A3
00120 *
00121             VB11R = CSR*B1
00122             VB12 = CSR*B2 + SNR*B3
00123 *
00124             AUA12 = ABS( CSL )*ABS( A2 ) + ABS( SNL )*ABS( A3 )
00125             AVB12 = ABS( CSR )*ABS( B2 ) + ABS( SNR )*ABS( B3 )
00126 *
00127 *           zero (1,2) elements of U'*A and V'*B
00128 *
00129             IF( ( ABS( UA11R )+ABS( UA12 ) ).NE.ZERO ) THEN
00130                IF( AUA12 / ( ABS( UA11R )+ABS( UA12 ) ).LE.AVB12 /
00131      $             ( ABS( VB11R )+ABS( VB12 ) ) ) THEN
00132                   CALL DLARTG( -UA11R, UA12, CSQ, SNQ, R )
00133                ELSE
00134                   CALL DLARTG( -VB11R, VB12, CSQ, SNQ, R )
00135                END IF
00136             ELSE
00137                CALL DLARTG( -VB11R, VB12, CSQ, SNQ, R )
00138             END IF
00139 *
00140             CSU = CSL
00141             SNU = -SNL
00142             CSV = CSR
00143             SNV = -SNR
00144 *
00145          ELSE
00146 *
00147 *           Compute the (2,1) and (2,2) elements of U'*A and V'*B,
00148 *           and (2,2) element of |U|'*|A| and |V|'*|B|.
00149 *
00150             UA21 = -SNL*A1
00151             UA22 = -SNL*A2 + CSL*A3
00152 *
00153             VB21 = -SNR*B1
00154             VB22 = -SNR*B2 + CSR*B3
00155 *
00156             AUA22 = ABS( SNL )*ABS( A2 ) + ABS( CSL )*ABS( A3 )
00157             AVB22 = ABS( SNR )*ABS( B2 ) + ABS( CSR )*ABS( B3 )
00158 *
00159 *           zero (2,2) elements of U'*A and V'*B, and then swap.
00160 *
00161             IF( ( ABS( UA21 )+ABS( UA22 ) ).NE.ZERO ) THEN
00162                IF( AUA22 / ( ABS( UA21 )+ABS( UA22 ) ).LE.AVB22 /
00163      $             ( ABS( VB21 )+ABS( VB22 ) ) ) THEN
00164                   CALL DLARTG( -UA21, UA22, CSQ, SNQ, R )
00165                ELSE
00166                   CALL DLARTG( -VB21, VB22, CSQ, SNQ, R )
00167                END IF
00168             ELSE
00169                CALL DLARTG( -VB21, VB22, CSQ, SNQ, R )
00170             END IF
00171 *
00172             CSU = SNL
00173             SNU = CSL
00174             CSV = SNR
00175             SNV = CSR
00176 *
00177          END IF
00178 *
00179       ELSE
00180 *
00181 *        Input matrices A and B are lower triangular matrices
00182 *
00183 *        Form matrix C = A*adj(B) = ( a 0 )
00184 *                                   ( c d )
00185 *
00186          A = A1*B3
00187          D = A3*B1
00188          C = A2*B3 - A3*B2
00189 *
00190 *        The SVD of real 2-by-2 triangular C
00191 *
00192 *         ( CSL -SNL )*( A 0 )*(  CSR  SNR ) = ( R 0 )
00193 *         ( SNL  CSL ) ( C D ) ( -SNR  CSR )   ( 0 T )
00194 *
00195          CALL DLASV2( A, C, D, S1, S2, SNR, CSR, SNL, CSL )
00196 *
00197          IF( ABS( CSR ).GE.ABS( SNR ) .OR. ABS( CSL ).GE.ABS( SNL ) )
00198      $        THEN
00199 *
00200 *           Compute the (2,1) and (2,2) elements of U'*A and V'*B,
00201 *           and (2,1) element of |U|'*|A| and |V|'*|B|.
00202 *
00203             UA21 = -SNR*A1 + CSR*A2
00204             UA22R = CSR*A3
00205 *
00206             VB21 = -SNL*B1 + CSL*B2
00207             VB22R = CSL*B3
00208 *
00209             AUA21 = ABS( SNR )*ABS( A1 ) + ABS( CSR )*ABS( A2 )
00210             AVB21 = ABS( SNL )*ABS( B1 ) + ABS( CSL )*ABS( B2 )
00211 *
00212 *           zero (2,1) elements of U'*A and V'*B.
00213 *
00214             IF( ( ABS( UA21 )+ABS( UA22R ) ).NE.ZERO ) THEN
00215                IF( AUA21 / ( ABS( UA21 )+ABS( UA22R ) ).LE.AVB21 /
00216      $             ( ABS( VB21 )+ABS( VB22R ) ) ) THEN
00217                   CALL DLARTG( UA22R, UA21, CSQ, SNQ, R )
00218                ELSE
00219                   CALL DLARTG( VB22R, VB21, CSQ, SNQ, R )
00220                END IF
00221             ELSE
00222                CALL DLARTG( VB22R, VB21, CSQ, SNQ, R )
00223             END IF
00224 *
00225             CSU = CSR
00226             SNU = -SNR
00227             CSV = CSL
00228             SNV = -SNL
00229 *
00230          ELSE
00231 *
00232 *           Compute the (1,1) and (1,2) elements of U'*A and V'*B,
00233 *           and (1,1) element of |U|'*|A| and |V|'*|B|.
00234 *
00235             UA11 = CSR*A1 + SNR*A2
00236             UA12 = SNR*A3
00237 *
00238             VB11 = CSL*B1 + SNL*B2
00239             VB12 = SNL*B3
00240 *
00241             AUA11 = ABS( CSR )*ABS( A1 ) + ABS( SNR )*ABS( A2 )
00242             AVB11 = ABS( CSL )*ABS( B1 ) + ABS( SNL )*ABS( B2 )
00243 *
00244 *           zero (1,1) elements of U'*A and V'*B, and then swap.
00245 *
00246             IF( ( ABS( UA11 )+ABS( UA12 ) ).NE.ZERO ) THEN
00247                IF( AUA11 / ( ABS( UA11 )+ABS( UA12 ) ).LE.AVB11 /
00248      $             ( ABS( VB11 )+ABS( VB12 ) ) ) THEN
00249                   CALL DLARTG( UA12, UA11, CSQ, SNQ, R )
00250                ELSE
00251                   CALL DLARTG( VB12, VB11, CSQ, SNQ, R )
00252                END IF
00253             ELSE
00254                CALL DLARTG( VB12, VB11, CSQ, SNQ, R )
00255             END IF
00256 *
00257             CSU = SNR
00258             SNU = CSR
00259             CSV = SNL
00260             SNV = CSL
00261 *
00262          END IF
00263 *
00264       END IF
00265 *
00266       RETURN
00267 *
00268 *     End of DLAGS2
00269 *
00270       END