1*> \brief \b DLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix.
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DLAEV2 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaev2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaev2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaev2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DLAEV2( A, B, C, RT1, RT2, CS1, SN1 )
22*
23*       .. Scalar Arguments ..
24*       DOUBLE PRECISION   A, B, C, CS1, RT1, RT2, SN1
25*       ..
26*
27*
28*> \par Purpose:
29*  =============
30*>
31*> \verbatim
32*>
33*> DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix
34*>    [  A   B  ]
35*>    [  B   C  ].
36*> On return, RT1 is the eigenvalue of larger absolute value, RT2 is the
37*> eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right
38*> eigenvector for RT1, giving the decomposition
39*>
40*>    [ CS1  SN1 ] [  A   B  ] [ CS1 -SN1 ]  =  [ RT1  0  ]
41*>    [-SN1  CS1 ] [  B   C  ] [ SN1  CS1 ]     [  0  RT2 ].
42*> \endverbatim
43*
44*  Arguments:
45*  ==========
46*
47*> \param[in] A
48*> \verbatim
49*>          A is DOUBLE PRECISION
50*>          The (1,1) element of the 2-by-2 matrix.
51*> \endverbatim
52*>
53*> \param[in] B
54*> \verbatim
55*>          B is DOUBLE PRECISION
56*>          The (1,2) element and the conjugate of the (2,1) element of
57*>          the 2-by-2 matrix.
58*> \endverbatim
59*>
60*> \param[in] C
61*> \verbatim
62*>          C is DOUBLE PRECISION
63*>          The (2,2) element of the 2-by-2 matrix.
64*> \endverbatim
65*>
66*> \param[out] RT1
67*> \verbatim
68*>          RT1 is DOUBLE PRECISION
69*>          The eigenvalue of larger absolute value.
70*> \endverbatim
71*>
72*> \param[out] RT2
73*> \verbatim
74*>          RT2 is DOUBLE PRECISION
75*>          The eigenvalue of smaller absolute value.
76*> \endverbatim
77*>
78*> \param[out] CS1
79*> \verbatim
80*>          CS1 is DOUBLE PRECISION
81*> \endverbatim
82*>
83*> \param[out] SN1
84*> \verbatim
85*>          SN1 is DOUBLE PRECISION
86*>          The vector (CS1, SN1) is a unit right eigenvector for RT1.
87*> \endverbatim
88*
89*  Authors:
90*  ========
91*
92*> \author Univ. of Tennessee
93*> \author Univ. of California Berkeley
94*> \author Univ. of Colorado Denver
95*> \author NAG Ltd.
96*
97*> \date September 2012
98*
99*> \ingroup auxOTHERauxiliary
100*
101*> \par Further Details:
102*  =====================
103*>
104*> \verbatim
105*>
106*>  RT1 is accurate to a few ulps barring over/underflow.
107*>
108*>  RT2 may be inaccurate if there is massive cancellation in the
109*>  determinant A*C-B*B; higher precision or correctly rounded or
110*>  correctly truncated arithmetic would be needed to compute RT2
111*>  accurately in all cases.
112*>
113*>  CS1 and SN1 are accurate to a few ulps barring over/underflow.
114*>
115*>  Overflow is possible only if RT1 is within a factor of 5 of overflow.
116*>  Underflow is harmless if the input data is 0 or exceeds
117*>     underflow_threshold / macheps.
118*> \endverbatim
119*>
120*  =====================================================================
121      SUBROUTINE GAL_DLAEV2(A,B,C,RT1,RT2,CS1,SN1)
122*
123*  -- LAPACK auxiliary routine (version 3.4.2) --
124*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
125*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
126*     September 2012
127*
128*     .. Scalar Arguments ..
129      DOUBLE PRECISION A,B,C,CS1,RT1,RT2,SN1
130*     ..
131*
132* =====================================================================
133*
134*     .. Parameters ..
135      DOUBLE PRECISION ONE
136      PARAMETER (ONE=1.0D0)
137      DOUBLE PRECISION TWO
138      PARAMETER (TWO=2.0D0)
139      DOUBLE PRECISION ZERO
140      PARAMETER (ZERO=0.0D0)
141      DOUBLE PRECISION HALF
142      PARAMETER (HALF=0.5D0)
143*     ..
144*     .. Local Scalars ..
145      INTEGER SGN1,SGN2
146      DOUBLE PRECISION AB,ACMN,ACMX,ACS,ADF,CS,CT,DF,RT,SM,
147     $TB,TN
148*     ..
149*     .. Intrinsic Functions ..
150      INTRINSIC ABS,SQRT
151*     ..
152*     .. Executable Statements ..
153*
154*     Compute the eigenvalues
155*
156      SM=A+C
157      DF=A-C
158      ADF=ABS(DF)
159      TB=B+B
160      AB=ABS(TB)
161      IF(ABS(A).GT.ABS(C))THEN
162      ACMX=A
163      ACMN=C
164      ELSE
165      ACMX=C
166      ACMN=A
167      END IF
168      IF(ADF.GT.AB)THEN
169      RT=ADF*SQRT(ONE+(AB/ADF)**2)
170      ELSE IF(ADF.LT.AB)THEN
171      RT=AB*SQRT(ONE+(ADF/AB)**2)
172      ELSE
173*
174*        Includes case AB=ADF=0
175*
176      RT=AB*SQRT(TWO)
177      END IF
178      IF(SM.LT.ZERO)THEN
179      RT1=HALF*(SM-RT)
180      SGN1=-1
181*
182*        Order of execution important.
183*        To get fully accurate smaller eigenvalue,
184*        next line needs to be executed in higher precision.
185*
186      RT2=(ACMX/RT1)*ACMN-(B/RT1)*B
187      ELSE IF(SM.GT.ZERO)THEN
188      RT1=HALF*(SM+RT)
189      SGN1=1
190*
191*        Order of execution important.
192*        To get fully accurate smaller eigenvalue,
193*        next line needs to be executed in higher precision.
194*
195      RT2=(ACMX/RT1)*ACMN-(B/RT1)*B
196      ELSE
197*
198*        Includes case RT1 = RT2 = 0
199*
200      RT1=HALF*RT
201      RT2=-HALF*RT
202      SGN1=1
203      END IF
204*
205*     Compute the eigenvector
206*
207      IF(DF.GE.ZERO)THEN
208      CS=DF+RT
209      SGN2=1
210      ELSE
211      CS=DF-RT
212      SGN2=-1
213      END IF
214      ACS=ABS(CS)
215      IF(ACS.GT.AB)THEN
216      CT=-TB/CS
217      SN1=ONE/SQRT(ONE+CT*CT)
218      CS1=CT*SN1
219      ELSE
220      IF(AB.EQ.ZERO)THEN
221      CS1=ONE
222      SN1=ZERO
223      ELSE
224      TN=-CS/TB
225      CS1=ONE/SQRT(ONE+TN*TN)
226      SN1=TN*CS1
227      END IF
228      END IF
229      IF(SGN1.EQ.SGN2)THEN
230      TN=CS1
231      CS1=-SN1
232      SN1=TN
233      END IF
234      RETURN
235*
236*     End of DLAEV2
237*
238      END
239