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