![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.2.2.
1: SUBROUTINE DLANV2( A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN ) 2: * 3: * -- LAPACK auxiliary routine (version 3.2.2) -- 4: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 6: * June 2010 7: * 8: * .. Scalar Arguments .. 9: DOUBLE PRECISION A, B, C, CS, D, RT1I, RT1R, RT2I, RT2R, SN 10: * .. 11: * 12: * Purpose 13: * ======= 14: * 15: * DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric 16: * matrix in standard form: 17: * 18: * [ A B ] = [ CS -SN ] [ AA BB ] [ CS SN ] 19: * [ C D ] [ SN CS ] [ CC DD ] [-SN CS ] 20: * 21: * where either 22: * 1) CC = 0 so that AA and DD are real eigenvalues of the matrix, or 23: * 2) AA = DD and BB*CC < 0, so that AA + or - sqrt(BB*CC) are complex 24: * conjugate eigenvalues. 25: * 26: * Arguments 27: * ========= 28: * 29: * A (input/output) DOUBLE PRECISION 30: * B (input/output) DOUBLE PRECISION 31: * C (input/output) DOUBLE PRECISION 32: * D (input/output) DOUBLE PRECISION 33: * On entry, the elements of the input matrix. 34: * On exit, they are overwritten by the elements of the 35: * standardised Schur form. 36: * 37: * RT1R (output) DOUBLE PRECISION 38: * RT1I (output) DOUBLE PRECISION 39: * RT2R (output) DOUBLE PRECISION 40: * RT2I (output) DOUBLE PRECISION 41: * The real and imaginary parts of the eigenvalues. If the 42: * eigenvalues are a complex conjugate pair, RT1I > 0. 43: * 44: * CS (output) DOUBLE PRECISION 45: * SN (output) DOUBLE PRECISION 46: * Parameters of the rotation matrix. 47: * 48: * Further Details 49: * =============== 50: * 51: * Modified by V. Sima, Research Institute for Informatics, Bucharest, 52: * Romania, to reduce the risk of cancellation errors, 53: * when computing real eigenvalues, and to ensure, if possible, that 54: * abs(RT1R) >= abs(RT2R). 55: * 56: * ===================================================================== 57: * 58: * .. Parameters .. 59: DOUBLE PRECISION ZERO, HALF, ONE 60: PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 ) 61: DOUBLE PRECISION MULTPL 62: PARAMETER ( MULTPL = 4.0D+0 ) 63: * .. 64: * .. Local Scalars .. 65: DOUBLE PRECISION AA, BB, BCMAX, BCMIS, CC, CS1, DD, EPS, P, SAB, 66: $ SAC, SCALE, SIGMA, SN1, TAU, TEMP, Z 67: * .. 68: * .. External Functions .. 69: DOUBLE PRECISION DLAMCH, DLAPY2 70: EXTERNAL DLAMCH, DLAPY2 71: * .. 72: * .. Intrinsic Functions .. 73: INTRINSIC ABS, MAX, MIN, SIGN, SQRT 74: * .. 75: * .. Executable Statements .. 76: * 77: EPS = DLAMCH( 'P' ) 78: IF( C.EQ.ZERO ) THEN 79: CS = ONE 80: SN = ZERO 81: GO TO 10 82: * 83: ELSE IF( B.EQ.ZERO ) THEN 84: * 85: * Swap rows and columns 86: * 87: CS = ZERO 88: SN = ONE 89: TEMP = D 90: D = A 91: A = TEMP 92: B = -C 93: C = ZERO 94: GO TO 10 95: ELSE IF( ( A-D ).EQ.ZERO .AND. SIGN( ONE, B ).NE.SIGN( ONE, C ) ) 96: $ THEN 97: CS = ONE 98: SN = ZERO 99: GO TO 10 100: ELSE 101: * 102: TEMP = A - D 103: P = HALF*TEMP 104: BCMAX = MAX( ABS( B ), ABS( C ) ) 105: BCMIS = MIN( ABS( B ), ABS( C ) )*SIGN( ONE, B )*SIGN( ONE, C ) 106: SCALE = MAX( ABS( P ), BCMAX ) 107: Z = ( P / SCALE )*P + ( BCMAX / SCALE )*BCMIS 108: * 109: * If Z is of the order of the machine accuracy, postpone the 110: * decision on the nature of eigenvalues 111: * 112: IF( Z.GE.MULTPL*EPS ) THEN 113: * 114: * Real eigenvalues. Compute A and D. 115: * 116: Z = P + SIGN( SQRT( SCALE )*SQRT( Z ), P ) 117: A = D + Z 118: D = D - ( BCMAX / Z )*BCMIS 119: * 120: * Compute B and the rotation matrix 121: * 122: TAU = DLAPY2( C, Z ) 123: CS = Z / TAU 124: SN = C / TAU 125: B = B - C 126: C = ZERO 127: ELSE 128: * 129: * Complex eigenvalues, or real (almost) equal eigenvalues. 130: * Make diagonal elements equal. 131: * 132: SIGMA = B + C 133: TAU = DLAPY2( SIGMA, TEMP ) 134: CS = SQRT( HALF*( ONE+ABS( SIGMA ) / TAU ) ) 135: SN = -( P / ( TAU*CS ) )*SIGN( ONE, SIGMA ) 136: * 137: * Compute [ AA BB ] = [ A B ] [ CS -SN ] 138: * [ CC DD ] [ C D ] [ SN CS ] 139: * 140: AA = A*CS + B*SN 141: BB = -A*SN + B*CS 142: CC = C*CS + D*SN 143: DD = -C*SN + D*CS 144: * 145: * Compute [ A B ] = [ CS SN ] [ AA BB ] 146: * [ C D ] [-SN CS ] [ CC DD ] 147: * 148: A = AA*CS + CC*SN 149: B = BB*CS + DD*SN 150: C = -AA*SN + CC*CS 151: D = -BB*SN + DD*CS 152: * 153: TEMP = HALF*( A+D ) 154: A = TEMP 155: D = TEMP 156: * 157: IF( C.NE.ZERO ) THEN 158: IF( B.NE.ZERO ) THEN 159: IF( SIGN( ONE, B ).EQ.SIGN( ONE, C ) ) THEN 160: * 161: * Real eigenvalues: reduce to upper triangular form 162: * 163: SAB = SQRT( ABS( B ) ) 164: SAC = SQRT( ABS( C ) ) 165: P = SIGN( SAB*SAC, C ) 166: TAU = ONE / SQRT( ABS( B+C ) ) 167: A = TEMP + P 168: D = TEMP - P 169: B = B - C 170: C = ZERO 171: CS1 = SAB*TAU 172: SN1 = SAC*TAU 173: TEMP = CS*CS1 - SN*SN1 174: SN = CS*SN1 + SN*CS1 175: CS = TEMP 176: END IF 177: ELSE 178: B = -C 179: C = ZERO 180: TEMP = CS 181: CS = -SN 182: SN = TEMP 183: END IF 184: END IF 185: END IF 186: * 187: END IF 188: * 189: 10 CONTINUE 190: * 191: * Store eigenvalues in (RT1R,RT1I) and (RT2R,RT2I). 192: * 193: RT1R = A 194: RT2R = D 195: IF( C.EQ.ZERO ) THEN 196: RT1I = ZERO 197: RT2I = ZERO 198: ELSE 199: RT1I = SQRT( ABS( B ) )*SQRT( ABS( C ) ) 200: RT2I = -RT1I 201: END IF 202: RETURN 203: * 204: * End of DLANV2 205: * 206: END