LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
claesy.f
Go to the documentation of this file.
1*> \brief \b CLAESY computes the eigenvalues and eigenvectors of a 2-by-2 complex symmetric matrix.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/claesy.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/claesy.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/claesy.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CLAESY( A, B, C, RT1, RT2, EVSCAL, CS1, SN1 )
22*
23* .. Scalar Arguments ..
24* COMPLEX A, B, C, CS1, EVSCAL, RT1, RT2, SN1
25* ..
26*
27*
28*> \par Purpose:
29* =============
30*>
31*> \verbatim
32*>
33*> CLAESY 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
52*> The ( 1, 1 ) element of input matrix.
53*> \endverbatim
54*>
55*> \param[in] B
56*> \verbatim
57*> B is COMPLEX
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
65*> The ( 2, 2 ) element of input matrix.
66*> \endverbatim
67*>
68*> \param[out] RT1
69*> \verbatim
70*> RT1 is COMPLEX
71*> The eigenvalue of larger modulus.
72*> \endverbatim
73*>
74*> \param[out] RT2
75*> \verbatim
76*> RT2 is COMPLEX
77*> The eigenvalue of smaller modulus.
78*> \endverbatim
79*>
80*> \param[out] EVSCAL
81*> \verbatim
82*> EVSCAL is COMPLEX
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
94*> \endverbatim
95*>
96*> \param[out] SN1
97*> \verbatim
98*> SN1 is COMPLEX
99*> If EVSCAL .NE. 0, ( CS1, SN1 ) is the unit right eigenvector
100*> for RT1.
101*> \endverbatim
102*
103* Authors:
104* ========
105*
106*> \author Univ. of Tennessee
107*> \author Univ. of California Berkeley
108*> \author Univ. of Colorado Denver
109*> \author NAG Ltd.
110*
111*> \ingroup laesy
112*
113* =====================================================================
114 SUBROUTINE claesy( A, B, C, RT1, RT2, EVSCAL, CS1, SN1 )
115*
116* -- LAPACK auxiliary routine --
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 A, B, C, CS1, EVSCAL, RT1, RT2, SN1
122* ..
123*
124* =====================================================================
125*
126* .. Parameters ..
127 REAL ZERO
128 parameter( zero = 0.0e0 )
129 REAL ONE
130 parameter( one = 1.0e0 )
131 COMPLEX CONE
132 parameter( cone = ( 1.0e0, 0.0e0 ) )
133 REAL HALF
134 parameter( half = 0.5e0 )
135 REAL THRESH
136 parameter( thresh = 0.1e0 )
137* ..
138* .. Local Scalars ..
139 REAL BABS, EVNORM, TABS, Z
140 COMPLEX 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
195* of eigenvectors X satisfies X * X**T = I . (No scaling is
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 CLAESY
217*
218 END
subroutine claesy(a, b, c, rt1, rt2, evscal, cs1, sn1)
CLAESY computes the eigenvalues and eigenvectors of a 2-by-2 complex symmetric matrix.
Definition claesy.f:115