Annotation of rpl/lapack/lapack/zlaesy.f, revision 1.3
1.1 bertrand 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
CVSweb interface <joel.bertrand@systella.fr>