Annotation of rpl/lapack/lapack/zlaesy.f, revision 1.19
1.12 bertrand 1: *> \brief \b ZLAESY computes the eigenvalues and eigenvectors of a 2-by-2 complex symmetric matrix.
1.9 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
1.16 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.9 bertrand 7: *
8: *> \htmlonly
1.16 bertrand 9: *> Download ZLAESY + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaesy.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaesy.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaesy.f">
1.9 bertrand 15: *> [TXT]</a>
1.16 bertrand 16: *> \endhtmlonly
1.9 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZLAESY( A, B, C, RT1, RT2, EVSCAL, CS1, SN1 )
1.16 bertrand 22: *
1.9 bertrand 23: * .. Scalar Arguments ..
24: * COMPLEX*16 A, B, C, CS1, EVSCAL, RT1, RT2, SN1
25: * ..
1.16 bertrand 26: *
1.9 bertrand 27: *
28: *> \par Purpose:
29: * =============
30: *>
31: *> \verbatim
32: *>
33: *> ZLAESY computes the eigendecomposition of a 2-by-2 symmetric matrix
34: *> ( ( A, B );( B, C ) )
35: *> provided the norm of the matrix of eigenvectors is larger than
36: *> some threshold value.
37: *>
38: *> RT1 is the eigenvalue of larger absolute value, and RT2 of
39: *> smaller absolute value. If the eigenvectors are computed, then
40: *> on return ( CS1, SN1 ) is the unit eigenvector for RT1, hence
41: *>
42: *> [ CS1 SN1 ] . [ A B ] . [ CS1 -SN1 ] = [ RT1 0 ]
43: *> [ -SN1 CS1 ] [ B C ] [ SN1 CS1 ] [ 0 RT2 ]
44: *> \endverbatim
45: *
46: * Arguments:
47: * ==========
48: *
49: *> \param[in] A
50: *> \verbatim
51: *> A is COMPLEX*16
52: *> The ( 1, 1 ) element of input matrix.
53: *> \endverbatim
54: *>
55: *> \param[in] B
56: *> \verbatim
57: *> B is COMPLEX*16
58: *> The ( 1, 2 ) element of input matrix. The ( 2, 1 ) element
59: *> is also given by B, since the 2-by-2 matrix is symmetric.
60: *> \endverbatim
61: *>
62: *> \param[in] C
63: *> \verbatim
64: *> C is COMPLEX*16
65: *> The ( 2, 2 ) element of input matrix.
66: *> \endverbatim
67: *>
68: *> \param[out] RT1
69: *> \verbatim
70: *> RT1 is COMPLEX*16
71: *> The eigenvalue of larger modulus.
72: *> \endverbatim
73: *>
74: *> \param[out] RT2
75: *> \verbatim
76: *> RT2 is COMPLEX*16
77: *> The eigenvalue of smaller modulus.
78: *> \endverbatim
79: *>
80: *> \param[out] EVSCAL
81: *> \verbatim
82: *> EVSCAL is COMPLEX*16
83: *> The complex value by which the eigenvector matrix was scaled
84: *> to make it orthonormal. If EVSCAL is zero, the eigenvectors
85: *> were not computed. This means one of two things: the 2-by-2
86: *> matrix could not be diagonalized, or the norm of the matrix
87: *> of eigenvectors before scaling was larger than the threshold
88: *> value THRESH (set below).
89: *> \endverbatim
90: *>
91: *> \param[out] CS1
92: *> \verbatim
93: *> CS1 is COMPLEX*16
94: *> \endverbatim
95: *>
96: *> \param[out] SN1
97: *> \verbatim
98: *> SN1 is COMPLEX*16
99: *> If EVSCAL .NE. 0, ( CS1, SN1 ) is the unit right eigenvector
100: *> for RT1.
101: *> \endverbatim
102: *
103: * Authors:
104: * ========
105: *
1.16 bertrand 106: *> \author Univ. of Tennessee
107: *> \author Univ. of California Berkeley
108: *> \author Univ. of Colorado Denver
109: *> \author NAG Ltd.
1.9 bertrand 110: *
111: *> \ingroup complex16SYauxiliary
112: *
113: * =====================================================================
1.1 bertrand 114: SUBROUTINE ZLAESY( A, B, C, RT1, RT2, EVSCAL, CS1, SN1 )
115: *
1.19 ! bertrand 116: * -- LAPACK auxiliary routine --
1.1 bertrand 117: * -- LAPACK is a software package provided by Univ. of Tennessee, --
118: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
119: *
120: * .. Scalar Arguments ..
121: COMPLEX*16 A, B, C, CS1, EVSCAL, RT1, RT2, SN1
122: * ..
123: *
124: * =====================================================================
125: *
126: * .. Parameters ..
127: DOUBLE PRECISION ZERO
128: PARAMETER ( ZERO = 0.0D0 )
129: DOUBLE PRECISION ONE
130: PARAMETER ( ONE = 1.0D0 )
131: COMPLEX*16 CONE
132: PARAMETER ( CONE = ( 1.0D0, 0.0D0 ) )
133: DOUBLE PRECISION HALF
134: PARAMETER ( HALF = 0.5D0 )
135: DOUBLE PRECISION THRESH
136: PARAMETER ( THRESH = 0.1D0 )
137: * ..
138: * .. Local Scalars ..
139: DOUBLE PRECISION BABS, EVNORM, TABS, Z
140: COMPLEX*16 S, T, TMP
141: * ..
142: * .. Intrinsic Functions ..
143: INTRINSIC ABS, MAX, SQRT
144: * ..
145: * .. Executable Statements ..
146: *
147: *
148: * Special case: The matrix is actually diagonal.
149: * To avoid divide by zero later, we treat this case separately.
150: *
151: IF( ABS( B ).EQ.ZERO ) THEN
152: RT1 = A
153: RT2 = C
154: IF( ABS( RT1 ).LT.ABS( RT2 ) ) THEN
155: TMP = RT1
156: RT1 = RT2
157: RT2 = TMP
158: CS1 = ZERO
159: SN1 = ONE
160: ELSE
161: CS1 = ONE
162: SN1 = ZERO
163: END IF
164: ELSE
165: *
166: * Compute the eigenvalues and eigenvectors.
167: * The characteristic equation is
168: * lambda **2 - (A+C) lambda + (A*C - B*B)
169: * and we solve it using the quadratic formula.
170: *
171: S = ( A+C )*HALF
172: T = ( A-C )*HALF
173: *
174: * Take the square root carefully to avoid over/under flow.
175: *
176: BABS = ABS( B )
177: TABS = ABS( T )
178: Z = MAX( BABS, TABS )
179: IF( Z.GT.ZERO )
180: $ T = Z*SQRT( ( T / Z )**2+( B / Z )**2 )
181: *
182: * Compute the two eigenvalues. RT1 and RT2 are exchanged
183: * if necessary so that RT1 will have the greater magnitude.
184: *
185: RT1 = S + T
186: RT2 = S - T
187: IF( ABS( RT1 ).LT.ABS( RT2 ) ) THEN
188: TMP = RT1
189: RT1 = RT2
190: RT2 = TMP
191: END IF
192: *
193: * Choose CS1 = 1 and SN1 to satisfy the first equation, then
194: * scale the components of this eigenvector so that the matrix
1.8 bertrand 195: * of eigenvectors X satisfies X * X**T = I . (No scaling is
1.1 bertrand 196: * done if the norm of the eigenvalue matrix is less than THRESH.)
197: *
198: SN1 = ( RT1-A ) / B
199: TABS = ABS( SN1 )
200: IF( TABS.GT.ONE ) THEN
201: T = TABS*SQRT( ( ONE / TABS )**2+( SN1 / TABS )**2 )
202: ELSE
203: T = SQRT( CONE+SN1*SN1 )
204: END IF
205: EVNORM = ABS( T )
206: IF( EVNORM.GE.THRESH ) THEN
207: EVSCAL = CONE / T
208: CS1 = EVSCAL
209: SN1 = SN1*EVSCAL
210: ELSE
211: EVSCAL = ZERO
212: END IF
213: END IF
214: RETURN
215: *
216: * End of ZLAESY
217: *
218: END
CVSweb interface <joel.bertrand@systella.fr>