LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
 subroutine claesy ( complex A, complex B, complex C, complex RT1, complex RT2, complex EVSCAL, complex CS1, complex SN1 )

CLAESY computes the eigenvalues and eigenvectors of a 2-by-2 complex symmetric matrix.

Download CLAESY + dependencies [TGZ] [ZIP] [TXT]

Purpose:
``` CLAESY computes the eigendecomposition of a 2-by-2 symmetric matrix
( ( A, B );( B, C ) )
provided the norm of the matrix of eigenvectors is larger than
some threshold value.

RT1 is the eigenvalue of larger absolute value, and RT2 of
smaller absolute value.  If the eigenvectors are computed, then
on return ( CS1, SN1 ) is the unit eigenvector for RT1, hence

[  CS1     SN1   ] . [ A  B ] . [ CS1    -SN1   ] = [ RT1  0  ]
[ -SN1     CS1   ]   [ B  C ]   [ SN1     CS1   ]   [  0  RT2 ]```
Parameters
 [in] A ``` A is COMPLEX The ( 1, 1 ) element of input matrix.``` [in] B ``` B is COMPLEX The ( 1, 2 ) element of input matrix. The ( 2, 1 ) element is also given by B, since the 2-by-2 matrix is symmetric.``` [in] C ``` C is COMPLEX The ( 2, 2 ) element of input matrix.``` [out] RT1 ``` RT1 is COMPLEX The eigenvalue of larger modulus.``` [out] RT2 ``` RT2 is COMPLEX The eigenvalue of smaller modulus.``` [out] EVSCAL ``` EVSCAL is COMPLEX The complex value by which the eigenvector matrix was scaled to make it orthonormal. If EVSCAL is zero, the eigenvectors were not computed. This means one of two things: the 2-by-2 matrix could not be diagonalized, or the norm of the matrix of eigenvectors before scaling was larger than the threshold value THRESH (set below).``` [out] CS1 ` CS1 is COMPLEX` [out] SN1 ``` SN1 is COMPLEX If EVSCAL .NE. 0, ( CS1, SN1 ) is the unit right eigenvector for RT1.```
Date
September 2012

Definition at line 117 of file claesy.f.

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 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 *