![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DLAEV2 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DLAEV2 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaev2.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaev2.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaev2.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DLAEV2( A, B, C, RT1, RT2, CS1, SN1 ) 00022 * 00023 * .. Scalar Arguments .. 00024 * DOUBLE PRECISION A, B, C, CS1, RT1, RT2, SN1 00025 * .. 00026 * 00027 * 00028 *> \par Purpose: 00029 * ============= 00030 *> 00031 *> \verbatim 00032 *> 00033 *> DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix 00034 *> [ A B ] 00035 *> [ B C ]. 00036 *> On return, RT1 is the eigenvalue of larger absolute value, RT2 is the 00037 *> eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right 00038 *> eigenvector for RT1, giving the decomposition 00039 *> 00040 *> [ CS1 SN1 ] [ A B ] [ CS1 -SN1 ] = [ RT1 0 ] 00041 *> [-SN1 CS1 ] [ B C ] [ SN1 CS1 ] [ 0 RT2 ]. 00042 *> \endverbatim 00043 * 00044 * Arguments: 00045 * ========== 00046 * 00047 *> \param[in] A 00048 *> \verbatim 00049 *> A is DOUBLE PRECISION 00050 *> The (1,1) element of the 2-by-2 matrix. 00051 *> \endverbatim 00052 *> 00053 *> \param[in] B 00054 *> \verbatim 00055 *> B is DOUBLE PRECISION 00056 *> The (1,2) element and the conjugate of the (2,1) element of 00057 *> the 2-by-2 matrix. 00058 *> \endverbatim 00059 *> 00060 *> \param[in] C 00061 *> \verbatim 00062 *> C is DOUBLE PRECISION 00063 *> The (2,2) element of the 2-by-2 matrix. 00064 *> \endverbatim 00065 *> 00066 *> \param[out] RT1 00067 *> \verbatim 00068 *> RT1 is DOUBLE PRECISION 00069 *> The eigenvalue of larger absolute value. 00070 *> \endverbatim 00071 *> 00072 *> \param[out] RT2 00073 *> \verbatim 00074 *> RT2 is DOUBLE PRECISION 00075 *> The eigenvalue of smaller absolute value. 00076 *> \endverbatim 00077 *> 00078 *> \param[out] CS1 00079 *> \verbatim 00080 *> CS1 is DOUBLE PRECISION 00081 *> \endverbatim 00082 *> 00083 *> \param[out] SN1 00084 *> \verbatim 00085 *> SN1 is DOUBLE PRECISION 00086 *> The vector (CS1, SN1) is a unit right eigenvector for RT1. 00087 *> \endverbatim 00088 * 00089 * Authors: 00090 * ======== 00091 * 00092 *> \author Univ. of Tennessee 00093 *> \author Univ. of California Berkeley 00094 *> \author Univ. of Colorado Denver 00095 *> \author NAG Ltd. 00096 * 00097 *> \date November 2011 00098 * 00099 *> \ingroup auxOTHERauxiliary 00100 * 00101 *> \par Further Details: 00102 * ===================== 00103 *> 00104 *> \verbatim 00105 *> 00106 *> RT1 is accurate to a few ulps barring over/underflow. 00107 *> 00108 *> RT2 may be inaccurate if there is massive cancellation in the 00109 *> determinant A*C-B*B; higher precision or correctly rounded or 00110 *> correctly truncated arithmetic would be needed to compute RT2 00111 *> accurately in all cases. 00112 *> 00113 *> CS1 and SN1 are accurate to a few ulps barring over/underflow. 00114 *> 00115 *> Overflow is possible only if RT1 is within a factor of 5 of overflow. 00116 *> Underflow is harmless if the input data is 0 or exceeds 00117 *> underflow_threshold / macheps. 00118 *> \endverbatim 00119 *> 00120 * ===================================================================== 00121 SUBROUTINE DLAEV2( A, B, C, RT1, RT2, CS1, SN1 ) 00122 * 00123 * -- LAPACK auxiliary routine (version 3.4.0) -- 00124 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00125 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00126 * November 2011 00127 * 00128 * .. Scalar Arguments .. 00129 DOUBLE PRECISION A, B, C, CS1, RT1, RT2, SN1 00130 * .. 00131 * 00132 * ===================================================================== 00133 * 00134 * .. Parameters .. 00135 DOUBLE PRECISION ONE 00136 PARAMETER ( ONE = 1.0D0 ) 00137 DOUBLE PRECISION TWO 00138 PARAMETER ( TWO = 2.0D0 ) 00139 DOUBLE PRECISION ZERO 00140 PARAMETER ( ZERO = 0.0D0 ) 00141 DOUBLE PRECISION HALF 00142 PARAMETER ( HALF = 0.5D0 ) 00143 * .. 00144 * .. Local Scalars .. 00145 INTEGER SGN1, SGN2 00146 DOUBLE PRECISION AB, ACMN, ACMX, ACS, ADF, CS, CT, DF, RT, SM, 00147 $ TB, TN 00148 * .. 00149 * .. Intrinsic Functions .. 00150 INTRINSIC ABS, SQRT 00151 * .. 00152 * .. Executable Statements .. 00153 * 00154 * Compute the eigenvalues 00155 * 00156 SM = A + C 00157 DF = A - C 00158 ADF = ABS( DF ) 00159 TB = B + B 00160 AB = ABS( TB ) 00161 IF( ABS( A ).GT.ABS( C ) ) THEN 00162 ACMX = A 00163 ACMN = C 00164 ELSE 00165 ACMX = C 00166 ACMN = A 00167 END IF 00168 IF( ADF.GT.AB ) THEN 00169 RT = ADF*SQRT( ONE+( AB / ADF )**2 ) 00170 ELSE IF( ADF.LT.AB ) THEN 00171 RT = AB*SQRT( ONE+( ADF / AB )**2 ) 00172 ELSE 00173 * 00174 * Includes case AB=ADF=0 00175 * 00176 RT = AB*SQRT( TWO ) 00177 END IF 00178 IF( SM.LT.ZERO ) THEN 00179 RT1 = HALF*( SM-RT ) 00180 SGN1 = -1 00181 * 00182 * Order of execution important. 00183 * To get fully accurate smaller eigenvalue, 00184 * next line needs to be executed in higher precision. 00185 * 00186 RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B 00187 ELSE IF( SM.GT.ZERO ) THEN 00188 RT1 = HALF*( SM+RT ) 00189 SGN1 = 1 00190 * 00191 * Order of execution important. 00192 * To get fully accurate smaller eigenvalue, 00193 * next line needs to be executed in higher precision. 00194 * 00195 RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B 00196 ELSE 00197 * 00198 * Includes case RT1 = RT2 = 0 00199 * 00200 RT1 = HALF*RT 00201 RT2 = -HALF*RT 00202 SGN1 = 1 00203 END IF 00204 * 00205 * Compute the eigenvector 00206 * 00207 IF( DF.GE.ZERO ) THEN 00208 CS = DF + RT 00209 SGN2 = 1 00210 ELSE 00211 CS = DF - RT 00212 SGN2 = -1 00213 END IF 00214 ACS = ABS( CS ) 00215 IF( ACS.GT.AB ) THEN 00216 CT = -TB / CS 00217 SN1 = ONE / SQRT( ONE+CT*CT ) 00218 CS1 = CT*SN1 00219 ELSE 00220 IF( AB.EQ.ZERO ) THEN 00221 CS1 = ONE 00222 SN1 = ZERO 00223 ELSE 00224 TN = -CS / TB 00225 CS1 = ONE / SQRT( ONE+TN*TN ) 00226 SN1 = TN*CS1 00227 END IF 00228 END IF 00229 IF( SGN1.EQ.SGN2 ) THEN 00230 TN = CS1 00231 CS1 = -SN1 00232 SN1 = TN 00233 END IF 00234 RETURN 00235 * 00236 * End of DLAEV2 00237 * 00238 END