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
9*> Download CLAESY + dependencies
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*> \date December 2016
112*
113*> \ingroup complexSYauxiliary
114*
115*  =====================================================================
116      SUBROUTINE CLAESY( A, B, C, RT1, RT2, EVSCAL, CS1, SN1 )
117*
118*  -- LAPACK auxiliary routine (version 3.7.0) --
119*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
120*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
121*     December 2016
122*
123*     .. Scalar Arguments ..
124      COMPLEX            A, B, C, CS1, EVSCAL, RT1, RT2, SN1
125*     ..
126*
127* =====================================================================
128*
129*     .. Parameters ..
130      REAL               ZERO
131      PARAMETER          ( ZERO = 0.0E0 )
132      REAL               ONE
133      PARAMETER          ( ONE = 1.0E0 )
134      COMPLEX            CONE
135      PARAMETER          ( CONE = ( 1.0E0, 0.0E0 ) )
136      REAL               HALF
137      PARAMETER          ( HALF = 0.5E0 )
138      REAL               THRESH
139      PARAMETER          ( THRESH = 0.1E0 )
140*     ..
141*     .. Local Scalars ..
142      REAL               BABS, EVNORM, TABS, Z
143      COMPLEX            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 CLAESY
220*
221      END
222