LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
zlaesy.f
Go to the documentation of this file.
1 *> \brief \b ZLAESY 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
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">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZLAESY( A, B, C, RT1, RT2, EVSCAL, CS1, SN1 )
22 *
23 * .. Scalar Arguments ..
24 * COMPLEX*16 A, B, C, CS1, EVSCAL, RT1, RT2, SN1
25 * ..
26 *
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 *
106 *> \author Univ. of Tennessee
107 *> \author Univ. of California Berkeley
108 *> \author Univ. of Colorado Denver
109 *> \author NAG Ltd.
110 *
111 *> \date September 2012
112 *
113 *> \ingroup complex16SYauxiliary
114 *
115 * =====================================================================
116  SUBROUTINE zlaesy( A, B, C, RT1, RT2, EVSCAL, CS1, SN1 )
117 *
118 * -- LAPACK auxiliary routine (version 3.4.2) --
119 * -- LAPACK is a software package provided by Univ. of Tennessee, --
120 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
121 * September 2012
122 *
123 * .. Scalar Arguments ..
124  COMPLEX*16 a, b, c, cs1, evscal, rt1, rt2, sn1
125 * ..
126 *
127 * =====================================================================
128 *
129 * .. Parameters ..
130  DOUBLE PRECISION zero
131  parameter( zero = 0.0d0 )
132  DOUBLE PRECISION one
133  parameter( one = 1.0d0 )
134  COMPLEX*16 cone
135  parameter( cone = ( 1.0d0, 0.0d0 ) )
136  DOUBLE PRECISION half
137  parameter( half = 0.5d0 )
138  DOUBLE PRECISION thresh
139  parameter( thresh = 0.1d0 )
140 * ..
141 * .. Local Scalars ..
142  DOUBLE PRECISION babs, evnorm, tabs, z
143  COMPLEX*16 s, t, tmp
144 * ..
145 * .. Intrinsic Functions ..
146  INTRINSIC abs, max, sqrt
147 * ..
148 * .. Executable Statements ..
149 *
150 *
151 * Special case: The matrix is actually diagonal.
152 * To avoid divide by zero later, we treat this case separately.
153 *
154  IF( abs( b ).EQ.zero ) THEN
155  rt1 = a
156  rt2 = c
157  IF( abs( rt1 ).LT.abs( rt2 ) ) THEN
158  tmp = rt1
159  rt1 = rt2
160  rt2 = tmp
161  cs1 = zero
162  sn1 = one
163  ELSE
164  cs1 = one
165  sn1 = zero
166  END IF
167  ELSE
168 *
169 * Compute the eigenvalues and eigenvectors.
170 * The characteristic equation is
171 * lambda **2 - (A+C) lambda + (A*C - B*B)
172 * and we solve it using the quadratic formula.
173 *
174  s = ( a+c )*half
175  t = ( a-c )*half
176 *
177 * Take the square root carefully to avoid over/under flow.
178 *
179  babs = abs( b )
180  tabs = abs( t )
181  z = max( babs, tabs )
182  IF( z.GT.zero )
183  $ t = z*sqrt( ( t / z )**2+( b / z )**2 )
184 *
185 * Compute the two eigenvalues. RT1 and RT2 are exchanged
186 * if necessary so that RT1 will have the greater magnitude.
187 *
188  rt1 = s + t
189  rt2 = s - t
190  IF( abs( rt1 ).LT.abs( rt2 ) ) THEN
191  tmp = rt1
192  rt1 = rt2
193  rt2 = tmp
194  END IF
195 *
196 * Choose CS1 = 1 and SN1 to satisfy the first equation, then
197 * scale the components of this eigenvector so that the matrix
198 * of eigenvectors X satisfies X * X**T = I . (No scaling is
199 * done if the norm of the eigenvalue matrix is less than THRESH.)
200 *
201  sn1 = ( rt1-a ) / b
202  tabs = abs( sn1 )
203  IF( tabs.GT.one ) THEN
204  t = tabs*sqrt( ( one / tabs )**2+( sn1 / tabs )**2 )
205  ELSE
206  t = sqrt( cone+sn1*sn1 )
207  END IF
208  evnorm = abs( t )
209  IF( evnorm.GE.thresh ) THEN
210  evscal = cone / t
211  cs1 = evscal
212  sn1 = sn1*evscal
213  ELSE
214  evscal = zero
215  END IF
216  END IF
217  RETURN
218 *
219 * End of ZLAESY
220 *
221  END