![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CLAESY 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CLAESY + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/claesy.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/claesy.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/claesy.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CLAESY( A, B, C, RT1, RT2, EVSCAL, CS1, SN1 ) 00022 * 00023 * .. Scalar Arguments .. 00024 * COMPLEX A, B, C, CS1, EVSCAL, RT1, RT2, SN1 00025 * .. 00026 * 00027 * 00028 *> \par Purpose: 00029 * ============= 00030 *> 00031 *> \verbatim 00032 *> 00033 *> CLAESY computes the eigendecomposition of a 2-by-2 symmetric matrix 00034 *> ( ( A, B );( B, C ) ) 00035 *> provided the norm of the matrix of eigenvectors is larger than 00036 *> some threshold value. 00037 *> 00038 *> RT1 is the eigenvalue of larger absolute value, and RT2 of 00039 *> smaller absolute value. If the eigenvectors are computed, then 00040 *> on return ( CS1, SN1 ) is the unit eigenvector for RT1, hence 00041 *> 00042 *> [ CS1 SN1 ] . [ A B ] . [ CS1 -SN1 ] = [ RT1 0 ] 00043 *> [ -SN1 CS1 ] [ B C ] [ SN1 CS1 ] [ 0 RT2 ] 00044 *> \endverbatim 00045 * 00046 * Arguments: 00047 * ========== 00048 * 00049 *> \param[in] A 00050 *> \verbatim 00051 *> A is COMPLEX 00052 *> The ( 1, 1 ) element of input matrix. 00053 *> \endverbatim 00054 *> 00055 *> \param[in] B 00056 *> \verbatim 00057 *> B is COMPLEX 00058 *> The ( 1, 2 ) element of input matrix. The ( 2, 1 ) element 00059 *> is also given by B, since the 2-by-2 matrix is symmetric. 00060 *> \endverbatim 00061 *> 00062 *> \param[in] C 00063 *> \verbatim 00064 *> C is COMPLEX 00065 *> The ( 2, 2 ) element of input matrix. 00066 *> \endverbatim 00067 *> 00068 *> \param[out] RT1 00069 *> \verbatim 00070 *> RT1 is COMPLEX 00071 *> The eigenvalue of larger modulus. 00072 *> \endverbatim 00073 *> 00074 *> \param[out] RT2 00075 *> \verbatim 00076 *> RT2 is COMPLEX 00077 *> The eigenvalue of smaller modulus. 00078 *> \endverbatim 00079 *> 00080 *> \param[out] EVSCAL 00081 *> \verbatim 00082 *> EVSCAL is COMPLEX 00083 *> The complex value by which the eigenvector matrix was scaled 00084 *> to make it orthonormal. If EVSCAL is zero, the eigenvectors 00085 *> were not computed. This means one of two things: the 2-by-2 00086 *> matrix could not be diagonalized, or the norm of the matrix 00087 *> of eigenvectors before scaling was larger than the threshold 00088 *> value THRESH (set below). 00089 *> \endverbatim 00090 *> 00091 *> \param[out] CS1 00092 *> \verbatim 00093 *> CS1 is COMPLEX 00094 *> \endverbatim 00095 *> 00096 *> \param[out] SN1 00097 *> \verbatim 00098 *> SN1 is COMPLEX 00099 *> If EVSCAL .NE. 0, ( CS1, SN1 ) is the unit right eigenvector 00100 *> for RT1. 00101 *> \endverbatim 00102 * 00103 * Authors: 00104 * ======== 00105 * 00106 *> \author Univ. of Tennessee 00107 *> \author Univ. of California Berkeley 00108 *> \author Univ. of Colorado Denver 00109 *> \author NAG Ltd. 00110 * 00111 *> \date November 2011 00112 * 00113 *> \ingroup complexSYauxiliary 00114 * 00115 * ===================================================================== 00116 SUBROUTINE CLAESY( A, B, C, RT1, RT2, EVSCAL, CS1, SN1 ) 00117 * 00118 * -- LAPACK auxiliary routine (version 3.4.0) -- 00119 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00120 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00121 * November 2011 00122 * 00123 * .. Scalar Arguments .. 00124 COMPLEX A, B, C, CS1, EVSCAL, RT1, RT2, SN1 00125 * .. 00126 * 00127 * ===================================================================== 00128 * 00129 * .. Parameters .. 00130 REAL ZERO 00131 PARAMETER ( ZERO = 0.0E0 ) 00132 REAL ONE 00133 PARAMETER ( ONE = 1.0E0 ) 00134 COMPLEX CONE 00135 PARAMETER ( CONE = ( 1.0E0, 0.0E0 ) ) 00136 REAL HALF 00137 PARAMETER ( HALF = 0.5E0 ) 00138 REAL THRESH 00139 PARAMETER ( THRESH = 0.1E0 ) 00140 * .. 00141 * .. Local Scalars .. 00142 REAL BABS, EVNORM, TABS, Z 00143 COMPLEX S, T, TMP 00144 * .. 00145 * .. Intrinsic Functions .. 00146 INTRINSIC ABS, MAX, SQRT 00147 * .. 00148 * .. Executable Statements .. 00149 * 00150 * 00151 * Special case: The matrix is actually diagonal. 00152 * To avoid divide by zero later, we treat this case separately. 00153 * 00154 IF( ABS( B ).EQ.ZERO ) THEN 00155 RT1 = A 00156 RT2 = C 00157 IF( ABS( RT1 ).LT.ABS( RT2 ) ) THEN 00158 TMP = RT1 00159 RT1 = RT2 00160 RT2 = TMP 00161 CS1 = ZERO 00162 SN1 = ONE 00163 ELSE 00164 CS1 = ONE 00165 SN1 = ZERO 00166 END IF 00167 ELSE 00168 * 00169 * Compute the eigenvalues and eigenvectors. 00170 * The characteristic equation is 00171 * lambda **2 - (A+C) lambda + (A*C - B*B) 00172 * and we solve it using the quadratic formula. 00173 * 00174 S = ( A+C )*HALF 00175 T = ( A-C )*HALF 00176 * 00177 * Take the square root carefully to avoid over/under flow. 00178 * 00179 BABS = ABS( B ) 00180 TABS = ABS( T ) 00181 Z = MAX( BABS, TABS ) 00182 IF( Z.GT.ZERO ) 00183 $ T = Z*SQRT( ( T / Z )**2+( B / Z )**2 ) 00184 * 00185 * Compute the two eigenvalues. RT1 and RT2 are exchanged 00186 * if necessary so that RT1 will have the greater magnitude. 00187 * 00188 RT1 = S + T 00189 RT2 = S - T 00190 IF( ABS( RT1 ).LT.ABS( RT2 ) ) THEN 00191 TMP = RT1 00192 RT1 = RT2 00193 RT2 = TMP 00194 END IF 00195 * 00196 * Choose CS1 = 1 and SN1 to satisfy the first equation, then 00197 * scale the components of this eigenvector so that the matrix 00198 * of eigenvectors X satisfies X * X**T = I . (No scaling is 00199 * done if the norm of the eigenvalue matrix is less than THRESH.) 00200 * 00201 SN1 = ( RT1-A ) / B 00202 TABS = ABS( SN1 ) 00203 IF( TABS.GT.ONE ) THEN 00204 T = TABS*SQRT( ( ONE / TABS )**2+( SN1 / TABS )**2 ) 00205 ELSE 00206 T = SQRT( CONE+SN1*SN1 ) 00207 END IF 00208 EVNORM = ABS( T ) 00209 IF( EVNORM.GE.THRESH ) THEN 00210 EVSCAL = CONE / T 00211 CS1 = EVSCAL 00212 SN1 = SN1*EVSCAL 00213 ELSE 00214 EVSCAL = ZERO 00215 END IF 00216 END IF 00217 RETURN 00218 * 00219 * End of CLAESY 00220 * 00221 END