![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZLAESY( A, B, C, RT1, RT2, EVSCAL, CS1, SN1 ) 2: * 3: * -- LAPACK auxiliary routine (version 3.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: * November 2006 7: * 8: * .. Scalar Arguments .. 9: COMPLEX*16 A, B, C, CS1, EVSCAL, RT1, RT2, SN1 10: * .. 11: * 12: * Purpose 13: * ======= 14: * 15: * ZLAESY computes the eigendecomposition of a 2-by-2 symmetric matrix 16: * ( ( A, B );( B, C ) ) 17: * provided the norm of the matrix of eigenvectors is larger than 18: * some threshold value. 19: * 20: * RT1 is the eigenvalue of larger absolute value, and RT2 of 21: * smaller absolute value. If the eigenvectors are computed, then 22: * on return ( CS1, SN1 ) is the unit eigenvector for RT1, hence 23: * 24: * [ CS1 SN1 ] . [ A B ] . [ CS1 -SN1 ] = [ RT1 0 ] 25: * [ -SN1 CS1 ] [ B C ] [ SN1 CS1 ] [ 0 RT2 ] 26: * 27: * Arguments 28: * ========= 29: * 30: * A (input) COMPLEX*16 31: * The ( 1, 1 ) element of input matrix. 32: * 33: * B (input) COMPLEX*16 34: * The ( 1, 2 ) element of input matrix. The ( 2, 1 ) element 35: * is also given by B, since the 2-by-2 matrix is symmetric. 36: * 37: * C (input) COMPLEX*16 38: * The ( 2, 2 ) element of input matrix. 39: * 40: * RT1 (output) COMPLEX*16 41: * The eigenvalue of larger modulus. 42: * 43: * RT2 (output) COMPLEX*16 44: * The eigenvalue of smaller modulus. 45: * 46: * EVSCAL (output) COMPLEX*16 47: * The complex value by which the eigenvector matrix was scaled 48: * to make it orthonormal. If EVSCAL is zero, the eigenvectors 49: * were not computed. This means one of two things: the 2-by-2 50: * matrix could not be diagonalized, or the norm of the matrix 51: * of eigenvectors before scaling was larger than the threshold 52: * value THRESH (set below). 53: * 54: * CS1 (output) COMPLEX*16 55: * SN1 (output) COMPLEX*16 56: * If EVSCAL .NE. 0, ( CS1, SN1 ) is the unit right eigenvector 57: * for RT1. 58: * 59: * ===================================================================== 60: * 61: * .. Parameters .. 62: DOUBLE PRECISION ZERO 63: PARAMETER ( ZERO = 0.0D0 ) 64: DOUBLE PRECISION ONE 65: PARAMETER ( ONE = 1.0D0 ) 66: COMPLEX*16 CONE 67: PARAMETER ( CONE = ( 1.0D0, 0.0D0 ) ) 68: DOUBLE PRECISION HALF 69: PARAMETER ( HALF = 0.5D0 ) 70: DOUBLE PRECISION THRESH 71: PARAMETER ( THRESH = 0.1D0 ) 72: * .. 73: * .. Local Scalars .. 74: DOUBLE PRECISION BABS, EVNORM, TABS, Z 75: COMPLEX*16 S, T, TMP 76: * .. 77: * .. Intrinsic Functions .. 78: INTRINSIC ABS, MAX, SQRT 79: * .. 80: * .. Executable Statements .. 81: * 82: * 83: * Special case: The matrix is actually diagonal. 84: * To avoid divide by zero later, we treat this case separately. 85: * 86: IF( ABS( B ).EQ.ZERO ) THEN 87: RT1 = A 88: RT2 = C 89: IF( ABS( RT1 ).LT.ABS( RT2 ) ) THEN 90: TMP = RT1 91: RT1 = RT2 92: RT2 = TMP 93: CS1 = ZERO 94: SN1 = ONE 95: ELSE 96: CS1 = ONE 97: SN1 = ZERO 98: END IF 99: ELSE 100: * 101: * Compute the eigenvalues and eigenvectors. 102: * The characteristic equation is 103: * lambda **2 - (A+C) lambda + (A*C - B*B) 104: * and we solve it using the quadratic formula. 105: * 106: S = ( A+C )*HALF 107: T = ( A-C )*HALF 108: * 109: * Take the square root carefully to avoid over/under flow. 110: * 111: BABS = ABS( B ) 112: TABS = ABS( T ) 113: Z = MAX( BABS, TABS ) 114: IF( Z.GT.ZERO ) 115: $ T = Z*SQRT( ( T / Z )**2+( B / Z )**2 ) 116: * 117: * Compute the two eigenvalues. RT1 and RT2 are exchanged 118: * if necessary so that RT1 will have the greater magnitude. 119: * 120: RT1 = S + T 121: RT2 = S - T 122: IF( ABS( RT1 ).LT.ABS( RT2 ) ) THEN 123: TMP = RT1 124: RT1 = RT2 125: RT2 = TMP 126: END IF 127: * 128: * Choose CS1 = 1 and SN1 to satisfy the first equation, then 129: * scale the components of this eigenvector so that the matrix 130: * of eigenvectors X satisfies X * X' = I . (No scaling is 131: * done if the norm of the eigenvalue matrix is less than THRESH.) 132: * 133: SN1 = ( RT1-A ) / B 134: TABS = ABS( SN1 ) 135: IF( TABS.GT.ONE ) THEN 136: T = TABS*SQRT( ( ONE / TABS )**2+( SN1 / TABS )**2 ) 137: ELSE 138: T = SQRT( CONE+SN1*SN1 ) 139: END IF 140: EVNORM = ABS( T ) 141: IF( EVNORM.GE.THRESH ) THEN 142: EVSCAL = CONE / T 143: CS1 = EVSCAL 144: SN1 = SN1*EVSCAL 145: ELSE 146: EVSCAL = ZERO 147: END IF 148: END IF 149: RETURN 150: * 151: * End of ZLAESY 152: * 153: END