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

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

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

Purpose:
 ZLAESY 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*16
          The ( 1, 1 ) element of input matrix.
[in]B
          B is COMPLEX*16
          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*16
          The ( 2, 2 ) element of input matrix.
[out]RT1
          RT1 is COMPLEX*16
          The eigenvalue of larger modulus.
[out]RT2
          RT2 is COMPLEX*16
          The eigenvalue of smaller modulus.
[out]EVSCAL
          EVSCAL is COMPLEX*16
          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*16
[out]SN1
          SN1 is COMPLEX*16
          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.
Date
September 2012

Definition at line 117 of file zlaesy.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*16 a, b, c, cs1, evscal, rt1, rt2, sn1
125 * ..
126 *
127 * =====================================================================
128 *
129 * .. Parameters ..
130  DOUBLE PRECISION zero
131  parameter ( zero = 0.0d0 )
132  DOUBLE PRECISION one
133  parameter ( one = 1.0d0 )
134  COMPLEX*16 cone
135  parameter ( cone = ( 1.0d0, 0.0d0 ) )
136  DOUBLE PRECISION half
137  parameter ( half = 0.5d0 )
138  DOUBLE PRECISION thresh
139  parameter ( thresh = 0.1d0 )
140 * ..
141 * .. Local Scalars ..
142  DOUBLE PRECISION babs, evnorm, tabs, z
143  COMPLEX*16 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 ZLAESY
220 *