LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ claesy()

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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 114 of file claesy.f.

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 A, B, C, CS1, EVSCAL, RT1, RT2, SN1
122 * ..
123 *
124 * =====================================================================
125 *
126 * .. Parameters ..
127  REAL ZERO
128  parameter( zero = 0.0e0 )
129  REAL ONE
130  parameter( one = 1.0e0 )
131  COMPLEX CONE
132  parameter( cone = ( 1.0e0, 0.0e0 ) )
133  REAL HALF
134  parameter( half = 0.5e0 )
135  REAL THRESH
136  parameter( thresh = 0.1e0 )
137 * ..
138 * .. Local Scalars ..
139  REAL BABS, EVNORM, TABS, Z
140  COMPLEX 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 CLAESY
217 *