Annotation of rpl/lapack/lapack/zlaesy.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b ZLAESY
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZLAESY + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaesy.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaesy.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaesy.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZLAESY( A, B, C, RT1, RT2, EVSCAL, CS1, SN1 )
! 22: *
! 23: * .. Scalar Arguments ..
! 24: * COMPLEX*16 A, B, C, CS1, EVSCAL, RT1, RT2, SN1
! 25: * ..
! 26: *
! 27: *
! 28: *> \par Purpose:
! 29: * =============
! 30: *>
! 31: *> \verbatim
! 32: *>
! 33: *> ZLAESY computes the eigendecomposition of a 2-by-2 symmetric matrix
! 34: *> ( ( A, B );( B, C ) )
! 35: *> provided the norm of the matrix of eigenvectors is larger than
! 36: *> some threshold value.
! 37: *>
! 38: *> RT1 is the eigenvalue of larger absolute value, and RT2 of
! 39: *> smaller absolute value. If the eigenvectors are computed, then
! 40: *> on return ( CS1, SN1 ) is the unit eigenvector for RT1, hence
! 41: *>
! 42: *> [ CS1 SN1 ] . [ A B ] . [ CS1 -SN1 ] = [ RT1 0 ]
! 43: *> [ -SN1 CS1 ] [ B C ] [ SN1 CS1 ] [ 0 RT2 ]
! 44: *> \endverbatim
! 45: *
! 46: * Arguments:
! 47: * ==========
! 48: *
! 49: *> \param[in] A
! 50: *> \verbatim
! 51: *> A is COMPLEX*16
! 52: *> The ( 1, 1 ) element of input matrix.
! 53: *> \endverbatim
! 54: *>
! 55: *> \param[in] B
! 56: *> \verbatim
! 57: *> B is COMPLEX*16
! 58: *> The ( 1, 2 ) element of input matrix. The ( 2, 1 ) element
! 59: *> is also given by B, since the 2-by-2 matrix is symmetric.
! 60: *> \endverbatim
! 61: *>
! 62: *> \param[in] C
! 63: *> \verbatim
! 64: *> C is COMPLEX*16
! 65: *> The ( 2, 2 ) element of input matrix.
! 66: *> \endverbatim
! 67: *>
! 68: *> \param[out] RT1
! 69: *> \verbatim
! 70: *> RT1 is COMPLEX*16
! 71: *> The eigenvalue of larger modulus.
! 72: *> \endverbatim
! 73: *>
! 74: *> \param[out] RT2
! 75: *> \verbatim
! 76: *> RT2 is COMPLEX*16
! 77: *> The eigenvalue of smaller modulus.
! 78: *> \endverbatim
! 79: *>
! 80: *> \param[out] EVSCAL
! 81: *> \verbatim
! 82: *> EVSCAL is COMPLEX*16
! 83: *> The complex value by which the eigenvector matrix was scaled
! 84: *> to make it orthonormal. If EVSCAL is zero, the eigenvectors
! 85: *> were not computed. This means one of two things: the 2-by-2
! 86: *> matrix could not be diagonalized, or the norm of the matrix
! 87: *> of eigenvectors before scaling was larger than the threshold
! 88: *> value THRESH (set below).
! 89: *> \endverbatim
! 90: *>
! 91: *> \param[out] CS1
! 92: *> \verbatim
! 93: *> CS1 is COMPLEX*16
! 94: *> \endverbatim
! 95: *>
! 96: *> \param[out] SN1
! 97: *> \verbatim
! 98: *> SN1 is COMPLEX*16
! 99: *> If EVSCAL .NE. 0, ( CS1, SN1 ) is the unit right eigenvector
! 100: *> for RT1.
! 101: *> \endverbatim
! 102: *
! 103: * Authors:
! 104: * ========
! 105: *
! 106: *> \author Univ. of Tennessee
! 107: *> \author Univ. of California Berkeley
! 108: *> \author Univ. of Colorado Denver
! 109: *> \author NAG Ltd.
! 110: *
! 111: *> \date November 2011
! 112: *
! 113: *> \ingroup complex16SYauxiliary
! 114: *
! 115: * =====================================================================
1.1 bertrand 116: SUBROUTINE ZLAESY( A, B, C, RT1, RT2, EVSCAL, CS1, SN1 )
117: *
1.9 ! bertrand 118: * -- LAPACK auxiliary routine (version 3.4.0) --
1.1 bertrand 119: * -- LAPACK is a software package provided by Univ. of Tennessee, --
120: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 121: * November 2011
1.1 bertrand 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
1.8 bertrand 198: * of eigenvectors X satisfies X * X**T = I . (No scaling is
1.1 bertrand 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: *
221: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlaesy.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack