 LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
 subroutine dlaev2 ( double precision A, double precision B, double precision C, double precision RT1, double precision RT2, double precision CS1, double precision SN1 )

DLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix.

Purpose:
DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix
[  A   B  ]
[  B   C  ].
On return, RT1 is the eigenvalue of larger absolute value, RT2 is the
eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right
eigenvector for RT1, giving the decomposition

[ CS1  SN1 ] [  A   B  ] [ CS1 -SN1 ]  =  [ RT1  0  ]
[-SN1  CS1 ] [  B   C  ] [ SN1  CS1 ]     [  0  RT2 ].
Parameters
 [in] A A is DOUBLE PRECISION The (1,1) element of the 2-by-2 matrix. [in] B B is DOUBLE PRECISION The (1,2) element and the conjugate of the (2,1) element of the 2-by-2 matrix. [in] C C is DOUBLE PRECISION The (2,2) element of the 2-by-2 matrix. [out] RT1 RT1 is DOUBLE PRECISION The eigenvalue of larger absolute value. [out] RT2 RT2 is DOUBLE PRECISION The eigenvalue of smaller absolute value. [out] CS1 CS1 is DOUBLE PRECISION [out] SN1 SN1 is DOUBLE PRECISION The vector (CS1, SN1) is a unit right eigenvector for RT1.
Date
September 2012
Further Details:
RT1 is accurate to a few ulps barring over/underflow.

RT2 may be inaccurate if there is massive cancellation in the
determinant A*C-B*B; higher precision or correctly rounded or
correctly truncated arithmetic would be needed to compute RT2
accurately in all cases.

CS1 and SN1 are accurate to a few ulps barring over/underflow.

Overflow is possible only if RT1 is within a factor of 5 of overflow.
Underflow is harmless if the input data is 0 or exceeds
underflow_threshold / macheps.

Definition at line 122 of file dlaev2.f.

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
170  ELSE IF( adf.LT.ab ) THEN
171  rt = ab*sqrt( one+( adf / ab )**2 )
172  ELSE
173 *
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 *

Here is the caller graph for this function: