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*> \ingroup complex16SYauxiliary
112*
113*  =====================================================================
114      SUBROUTINE ZLAESY( 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*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
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 ZLAESY
217*
218      END
219