 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE ZLAGS2( 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 DOUBLE PRECISION A1, A3, B1, B3, CSQ, CSU, CSV 00012 COMPLEX*16 A2, B2, SNQ, SNU, SNV 00013 * .. 00014 * 00015 * Purpose 00016 * ======= 00017 * 00018 * ZLAGS2 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) DOUBLE PRECISION 00057 * A2 (input) COMPLEX*16 00058 * A3 (input) DOUBLE PRECISION 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) DOUBLE PRECISION 00063 * B2 (input) COMPLEX*16 00064 * B3 (input) DOUBLE PRECISION 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) DOUBLE PRECISION 00069 * SNU (output) COMPLEX*16 00070 * The desired unitary matrix U. 00071 * 00072 * CSV (output) DOUBLE PRECISION 00073 * SNV (output) COMPLEX*16 00074 * The desired unitary matrix V. 00075 * 00076 * CSQ (output) DOUBLE PRECISION 00077 * SNQ (output) COMPLEX*16 00078 * The desired unitary matrix Q. 00079 * 00080 * ===================================================================== 00081 * 00082 * .. Parameters .. 00083 DOUBLE PRECISION ZERO, ONE 00084 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00085 * .. 00086 * .. Local Scalars .. 00087 DOUBLE PRECISION A, AUA11, AUA12, AUA21, AUA22, AVB12, AVB11, 00088 $ AVB21, AVB22, CSL, CSR, D, FB, FC, S1, S2, 00089 $ SNL, SNR, UA11R, UA22R, VB11R, VB22R 00090 COMPLEX*16 B, C, D1, R, T, UA11, UA12, UA21, UA22, VB11, 00091 $ VB12, VB21, VB22 00092 * .. 00093 * .. External Subroutines .. 00094 EXTERNAL DLASV2, ZLARTG 00095 * .. 00096 * .. Intrinsic Functions .. 00097 INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG 00098 * .. 00099 * .. Statement Functions .. 00100 DOUBLE PRECISION ABS1 00101 * .. 00102 * .. Statement Function definitions .. 00103 ABS1( T ) = ABS( DBLE( T ) ) + ABS( DIMAG( 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 DLASV2( 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 ZLARTG( -DCMPLX( VB11R ), DCONJG( VB12 ), CSQ, SNQ, 00152 $ R ) 00153 ELSE IF( ( ABS( VB11R )+ABS1( VB12 ) ).EQ.ZERO ) THEN 00154 CALL ZLARTG( -DCMPLX( UA11R ), DCONJG( UA12 ), CSQ, SNQ, 00155 $ R ) 00156 ELSE IF( AUA12 / ( ABS( UA11R )+ABS1( UA12 ) ).LE.AVB12 / 00157 $ ( ABS( VB11R )+ABS1( VB12 ) ) ) THEN 00158 CALL ZLARTG( -DCMPLX( UA11R ), DCONJG( UA12 ), CSQ, SNQ, 00159 $ R ) 00160 ELSE 00161 CALL ZLARTG( -DCMPLX( VB11R ), DCONJG( 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 = -DCONJG( D1 )*SNL*A1 00176 UA22 = -DCONJG( D1 )*SNL*A2 + CSL*A3 00177 * 00178 VB21 = -DCONJG( D1 )*SNR*B1 00179 VB22 = -DCONJG( 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 ZLARTG( -DCONJG( VB21 ), DCONJG( VB22 ), CSQ, SNQ, 00188 $ R ) 00189 ELSE IF( ( ABS1( VB21 )+ABS( VB22 ) ).EQ.ZERO ) THEN 00190 CALL ZLARTG( -DCONJG( UA21 ), DCONJG( UA22 ), CSQ, SNQ, 00191 $ R ) 00192 ELSE IF( AUA22 / ( ABS1( UA21 )+ABS1( UA22 ) ).LE.AVB22 / 00193 $ ( ABS1( VB21 )+ABS1( VB22 ) ) ) THEN 00194 CALL ZLARTG( -DCONJG( UA21 ), DCONJG( UA22 ), CSQ, SNQ, 00195 $ R ) 00196 ELSE 00197 CALL ZLARTG( -DCONJG( VB21 ), DCONJG( VB22 ), CSQ, SNQ, 00198 $ R ) 00199 END IF 00200 * 00201 CSU = SNL 00202 SNU = D1*CSL 00203 CSV = SNR 00204 SNV = D1*CSR 00205 * 00206 END IF 00207 * 00208 ELSE 00209 * 00210 * Input matrices A and B are lower triangular matrices 00211 * 00212 * Form matrix C = A*adj(B) = ( a 0 ) 00213 * ( c d ) 00214 * 00215 A = A1*B3 00216 D = A3*B1 00217 C = A2*B3 - A3*B2 00218 FC = ABS( C ) 00219 * 00220 * Transform complex 2-by-2 matrix C to real matrix by unitary 00221 * diagonal matrix diag(d1,1). 00222 * 00223 D1 = ONE 00224 IF( FC.NE.ZERO ) 00225 $ D1 = C / FC 00226 * 00227 * The SVD of real 2 by 2 triangular C 00228 * 00229 * ( CSL -SNL )*( A 0 )*( CSR SNR ) = ( R 0 ) 00230 * ( SNL CSL ) ( C D ) ( -SNR CSR ) ( 0 T ) 00231 * 00232 CALL DLASV2( A, FC, D, S1, S2, SNR, CSR, SNL, CSL ) 00233 * 00234 IF( ABS( CSR ).GE.ABS( SNR ) .OR. ABS( CSL ).GE.ABS( SNL ) ) 00235 $ THEN 00236 * 00237 * Compute the (2,1) and (2,2) elements of U**H *A and V**H *B, 00238 * and (2,1) element of |U|**H *|A| and |V|**H *|B|. 00239 * 00240 UA21 = -D1*SNR*A1 + CSR*A2 00241 UA22R = CSR*A3 00242 * 00243 VB21 = -D1*SNL*B1 + CSL*B2 00244 VB22R = CSL*B3 00245 * 00246 AUA21 = ABS( SNR )*ABS( A1 ) + ABS( CSR )*ABS1( A2 ) 00247 AVB21 = ABS( SNL )*ABS( B1 ) + ABS( CSL )*ABS1( B2 ) 00248 * 00249 * zero (2,1) elements of U**H *A and V**H *B. 00250 * 00251 IF( ( ABS1( UA21 )+ABS( UA22R ) ).EQ.ZERO ) THEN 00252 CALL ZLARTG( DCMPLX( VB22R ), VB21, CSQ, SNQ, R ) 00253 ELSE IF( ( ABS1( VB21 )+ABS( VB22R ) ).EQ.ZERO ) THEN 00254 CALL ZLARTG( DCMPLX( UA22R ), UA21, CSQ, SNQ, R ) 00255 ELSE IF( AUA21 / ( ABS1( UA21 )+ABS( UA22R ) ).LE.AVB21 / 00256 $ ( ABS1( VB21 )+ABS( VB22R ) ) ) THEN 00257 CALL ZLARTG( DCMPLX( UA22R ), UA21, CSQ, SNQ, R ) 00258 ELSE 00259 CALL ZLARTG( DCMPLX( VB22R ), VB21, CSQ, SNQ, R ) 00260 END IF 00261 * 00262 CSU = CSR 00263 SNU = -DCONJG( D1 )*SNR 00264 CSV = CSL 00265 SNV = -DCONJG( D1 )*SNL 00266 * 00267 ELSE 00268 * 00269 * Compute the (1,1) and (1,2) elements of U**H *A and V**H *B, 00270 * and (1,1) element of |U|**H *|A| and |V|**H *|B|. 00271 * 00272 UA11 = CSR*A1 + DCONJG( D1 )*SNR*A2 00273 UA12 = DCONJG( D1 )*SNR*A3 00274 * 00275 VB11 = CSL*B1 + DCONJG( D1 )*SNL*B2 00276 VB12 = DCONJG( D1 )*SNL*B3 00277 * 00278 AUA11 = ABS( CSR )*ABS( A1 ) + ABS( SNR )*ABS1( A2 ) 00279 AVB11 = ABS( CSL )*ABS( B1 ) + ABS( SNL )*ABS1( B2 ) 00280 * 00281 * zero (1,1) elements of U**H *A and V**H *B, and then swap. 00282 * 00283 IF( ( ABS1( UA11 )+ABS1( UA12 ) ).EQ.ZERO ) THEN 00284 CALL ZLARTG( VB12, VB11, CSQ, SNQ, R ) 00285 ELSE IF( ( ABS1( VB11 )+ABS1( VB12 ) ).EQ.ZERO ) THEN 00286 CALL ZLARTG( UA12, UA11, CSQ, SNQ, R ) 00287 ELSE IF( AUA11 / ( ABS1( UA11 )+ABS1( UA12 ) ).LE.AVB11 / 00288 $ ( ABS1( VB11 )+ABS1( VB12 ) ) ) THEN 00289 CALL ZLARTG( UA12, UA11, CSQ, SNQ, R ) 00290 ELSE 00291 CALL ZLARTG( VB12, VB11, CSQ, SNQ, R ) 00292 END IF 00293 * 00294 CSU = SNR 00295 SNU = DCONJG( D1 )*CSR 00296 CSV = SNL 00297 SNV = DCONJG( D1 )*CSL 00298 * 00299 END IF 00300 * 00301 END IF 00302 * 00303 RETURN 00304 * 00305 * End of ZLAGS2 00306 * 00307 END
1.7.3 
