![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZLAGS2( 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, A3, B1, B3, CSQ, CSU, CSV 12: COMPLEX*16 A2, B2, SNQ, SNU, SNV 13: * .. 14: * 15: * Purpose 16: * ======= 17: * 18: * ZLAGS2 computes 2-by-2 unitary 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: * where 35: * 36: * U = ( CSU SNU ), V = ( CSV SNV ), 37: * ( -CONJG(SNU) CSU ) ( -CONJG(SNV) CSV ) 38: * 39: * Q = ( CSQ SNQ ) 40: * ( -CONJG(SNQ) CSQ ) 41: * 42: * Z' denotes the conjugate transpose of Z. 43: * 44: * The rows of the transformed A and B are parallel. Moreover, if the 45: * input 2-by-2 matrix A is not zero, then the transformed (1,1) entry 46: * of A is not zero. If the input matrices A and B are both not zero, 47: * then the transformed (2,2) element of B is not zero, except when the 48: * first rows of input A and B are parallel and the second rows are 49: * zero. 50: * 51: * Arguments 52: * ========= 53: * 54: * UPPER (input) LOGICAL 55: * = .TRUE.: the input matrices A and B are upper triangular. 56: * = .FALSE.: the input matrices A and B are lower triangular. 57: * 58: * A1 (input) DOUBLE PRECISION 59: * A2 (input) COMPLEX*16 60: * A3 (input) DOUBLE PRECISION 61: * On entry, A1, A2 and A3 are elements of the input 2-by-2 62: * upper (lower) triangular matrix A. 63: * 64: * B1 (input) DOUBLE PRECISION 65: * B2 (input) COMPLEX*16 66: * B3 (input) DOUBLE PRECISION 67: * On entry, B1, B2 and B3 are elements of the input 2-by-2 68: * upper (lower) triangular matrix B. 69: * 70: * CSU (output) DOUBLE PRECISION 71: * SNU (output) COMPLEX*16 72: * The desired unitary matrix U. 73: * 74: * CSV (output) DOUBLE PRECISION 75: * SNV (output) COMPLEX*16 76: * The desired unitary matrix V. 77: * 78: * CSQ (output) DOUBLE PRECISION 79: * SNQ (output) COMPLEX*16 80: * The desired unitary matrix Q. 81: * 82: * ===================================================================== 83: * 84: * .. Parameters .. 85: DOUBLE PRECISION ZERO, ONE 86: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 87: * .. 88: * .. Local Scalars .. 89: DOUBLE PRECISION A, AUA11, AUA12, AUA21, AUA22, AVB12, AVB11, 90: $ AVB21, AVB22, CSL, CSR, D, FB, FC, S1, S2, 91: $ SNL, SNR, UA11R, UA22R, VB11R, VB22R 92: COMPLEX*16 B, C, D1, R, T, UA11, UA12, UA21, UA22, VB11, 93: $ VB12, VB21, VB22 94: * .. 95: * .. External Subroutines .. 96: EXTERNAL DLASV2, ZLARTG 97: * .. 98: * .. Intrinsic Functions .. 99: INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG 100: * .. 101: * .. Statement Functions .. 102: DOUBLE PRECISION ABS1 103: * .. 104: * .. Statement Function definitions .. 105: ABS1( T ) = ABS( DBLE( T ) ) + ABS( DIMAG( T ) ) 106: * .. 107: * .. Executable Statements .. 108: * 109: IF( UPPER ) THEN 110: * 111: * Input matrices A and B are upper triangular matrices 112: * 113: * Form matrix C = A*adj(B) = ( a b ) 114: * ( 0 d ) 115: * 116: A = A1*B3 117: D = A3*B1 118: B = A2*B1 - A1*B2 119: FB = ABS( B ) 120: * 121: * Transform complex 2-by-2 matrix C to real matrix by unitary 122: * diagonal matrix diag(1,D1). 123: * 124: D1 = ONE 125: IF( FB.NE.ZERO ) 126: $ D1 = B / FB 127: * 128: * The SVD of real 2 by 2 triangular C 129: * 130: * ( CSL -SNL )*( A B )*( CSR SNR ) = ( R 0 ) 131: * ( SNL CSL ) ( 0 D ) ( -SNR CSR ) ( 0 T ) 132: * 133: CALL DLASV2( A, FB, D, S1, S2, SNR, CSR, SNL, CSL ) 134: * 135: IF( ABS( CSL ).GE.ABS( SNL ) .OR. ABS( CSR ).GE.ABS( SNR ) ) 136: $ THEN 137: * 138: * Compute the (1,1) and (1,2) elements of U'*A and V'*B, 139: * and (1,2) element of |U|'*|A| and |V|'*|B|. 140: * 141: UA11R = CSL*A1 142: UA12 = CSL*A2 + D1*SNL*A3 143: * 144: VB11R = CSR*B1 145: VB12 = CSR*B2 + D1*SNR*B3 146: * 147: AUA12 = ABS( CSL )*ABS1( A2 ) + ABS( SNL )*ABS( A3 ) 148: AVB12 = ABS( CSR )*ABS1( B2 ) + ABS( SNR )*ABS( B3 ) 149: * 150: * zero (1,2) elements of U'*A and V'*B 151: * 152: IF( ( ABS( UA11R )+ABS1( UA12 ) ).EQ.ZERO ) THEN 153: CALL ZLARTG( -DCMPLX( VB11R ), DCONJG( VB12 ), CSQ, SNQ, 154: $ R ) 155: ELSE IF( ( ABS( VB11R )+ABS1( VB12 ) ).EQ.ZERO ) THEN 156: CALL ZLARTG( -DCMPLX( UA11R ), DCONJG( UA12 ), CSQ, SNQ, 157: $ R ) 158: ELSE IF( AUA12 / ( ABS( UA11R )+ABS1( UA12 ) ).LE.AVB12 / 159: $ ( ABS( VB11R )+ABS1( VB12 ) ) ) THEN 160: CALL ZLARTG( -DCMPLX( UA11R ), DCONJG( UA12 ), CSQ, SNQ, 161: $ R ) 162: ELSE 163: CALL ZLARTG( -DCMPLX( VB11R ), DCONJG( VB12 ), CSQ, SNQ, 164: $ R ) 165: END IF 166: * 167: CSU = CSL 168: SNU = -D1*SNL 169: CSV = CSR 170: SNV = -D1*SNR 171: * 172: ELSE 173: * 174: * Compute the (2,1) and (2,2) elements of U'*A and V'*B, 175: * and (2,2) element of |U|'*|A| and |V|'*|B|. 176: * 177: UA21 = -DCONJG( D1 )*SNL*A1 178: UA22 = -DCONJG( D1 )*SNL*A2 + CSL*A3 179: * 180: VB21 = -DCONJG( D1 )*SNR*B1 181: VB22 = -DCONJG( D1 )*SNR*B2 + CSR*B3 182: * 183: AUA22 = ABS( SNL )*ABS1( A2 ) + ABS( CSL )*ABS( A3 ) 184: AVB22 = ABS( SNR )*ABS1( B2 ) + ABS( CSR )*ABS( B3 ) 185: * 186: * zero (2,2) elements of U'*A and V'*B, and then swap. 187: * 188: IF( ( ABS1( UA21 )+ABS1( UA22 ) ).EQ.ZERO ) THEN 189: CALL ZLARTG( -DCONJG( VB21 ), DCONJG( VB22 ), CSQ, SNQ, 190: $ R ) 191: ELSE IF( ( ABS1( VB21 )+ABS( VB22 ) ).EQ.ZERO ) THEN 192: CALL ZLARTG( -DCONJG( UA21 ), DCONJG( UA22 ), CSQ, SNQ, 193: $ R ) 194: ELSE IF( AUA22 / ( ABS1( UA21 )+ABS1( UA22 ) ).LE.AVB22 / 195: $ ( ABS1( VB21 )+ABS1( VB22 ) ) ) THEN 196: CALL ZLARTG( -DCONJG( UA21 ), DCONJG( UA22 ), CSQ, SNQ, 197: $ R ) 198: ELSE 199: CALL ZLARTG( -DCONJG( VB21 ), DCONJG( VB22 ), CSQ, SNQ, 200: $ R ) 201: END IF 202: * 203: CSU = SNL 204: SNU = D1*CSL 205: CSV = SNR 206: SNV = D1*CSR 207: * 208: END IF 209: * 210: ELSE 211: * 212: * Input matrices A and B are lower triangular matrices 213: * 214: * Form matrix C = A*adj(B) = ( a 0 ) 215: * ( c d ) 216: * 217: A = A1*B3 218: D = A3*B1 219: C = A2*B3 - A3*B2 220: FC = ABS( C ) 221: * 222: * Transform complex 2-by-2 matrix C to real matrix by unitary 223: * diagonal matrix diag(d1,1). 224: * 225: D1 = ONE 226: IF( FC.NE.ZERO ) 227: $ D1 = C / FC 228: * 229: * The SVD of real 2 by 2 triangular C 230: * 231: * ( CSL -SNL )*( A 0 )*( CSR SNR ) = ( R 0 ) 232: * ( SNL CSL ) ( C D ) ( -SNR CSR ) ( 0 T ) 233: * 234: CALL DLASV2( A, FC, D, S1, S2, SNR, CSR, SNL, CSL ) 235: * 236: IF( ABS( CSR ).GE.ABS( SNR ) .OR. ABS( CSL ).GE.ABS( SNL ) ) 237: $ THEN 238: * 239: * Compute the (2,1) and (2,2) elements of U'*A and V'*B, 240: * and (2,1) element of |U|'*|A| and |V|'*|B|. 241: * 242: UA21 = -D1*SNR*A1 + CSR*A2 243: UA22R = CSR*A3 244: * 245: VB21 = -D1*SNL*B1 + CSL*B2 246: VB22R = CSL*B3 247: * 248: AUA21 = ABS( SNR )*ABS( A1 ) + ABS( CSR )*ABS1( A2 ) 249: AVB21 = ABS( SNL )*ABS( B1 ) + ABS( CSL )*ABS1( B2 ) 250: * 251: * zero (2,1) elements of U'*A and V'*B. 252: * 253: IF( ( ABS1( UA21 )+ABS( UA22R ) ).EQ.ZERO ) THEN 254: CALL ZLARTG( DCMPLX( VB22R ), VB21, CSQ, SNQ, R ) 255: ELSE IF( ( ABS1( VB21 )+ABS( VB22R ) ).EQ.ZERO ) THEN 256: CALL ZLARTG( DCMPLX( UA22R ), UA21, CSQ, SNQ, R ) 257: ELSE IF( AUA21 / ( ABS1( UA21 )+ABS( UA22R ) ).LE.AVB21 / 258: $ ( ABS1( VB21 )+ABS( VB22R ) ) ) THEN 259: CALL ZLARTG( DCMPLX( UA22R ), UA21, CSQ, SNQ, R ) 260: ELSE 261: CALL ZLARTG( DCMPLX( VB22R ), VB21, CSQ, SNQ, R ) 262: END IF 263: * 264: CSU = CSR 265: SNU = -DCONJG( D1 )*SNR 266: CSV = CSL 267: SNV = -DCONJG( D1 )*SNL 268: * 269: ELSE 270: * 271: * Compute the (1,1) and (1,2) elements of U'*A and V'*B, 272: * and (1,1) element of |U|'*|A| and |V|'*|B|. 273: * 274: UA11 = CSR*A1 + DCONJG( D1 )*SNR*A2 275: UA12 = DCONJG( D1 )*SNR*A3 276: * 277: VB11 = CSL*B1 + DCONJG( D1 )*SNL*B2 278: VB12 = DCONJG( D1 )*SNL*B3 279: * 280: AUA11 = ABS( CSR )*ABS( A1 ) + ABS( SNR )*ABS1( A2 ) 281: AVB11 = ABS( CSL )*ABS( B1 ) + ABS( SNL )*ABS1( B2 ) 282: * 283: * zero (1,1) elements of U'*A and V'*B, and then swap. 284: * 285: IF( ( ABS1( UA11 )+ABS1( UA12 ) ).EQ.ZERO ) THEN 286: CALL ZLARTG( VB12, VB11, CSQ, SNQ, R ) 287: ELSE IF( ( ABS1( VB11 )+ABS1( VB12 ) ).EQ.ZERO ) THEN 288: CALL ZLARTG( UA12, UA11, CSQ, SNQ, R ) 289: ELSE IF( AUA11 / ( ABS1( UA11 )+ABS1( UA12 ) ).LE.AVB11 / 290: $ ( ABS1( VB11 )+ABS1( VB12 ) ) ) THEN 291: CALL ZLARTG( UA12, UA11, CSQ, SNQ, R ) 292: ELSE 293: CALL ZLARTG( VB12, VB11, CSQ, SNQ, R ) 294: END IF 295: * 296: CSU = SNR 297: SNU = DCONJG( D1 )*CSR 298: CSV = SNL 299: SNV = DCONJG( D1 )*CSL 300: * 301: END IF 302: * 303: END IF 304: * 305: RETURN 306: * 307: * End of ZLAGS2 308: * 309: END