 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE SLAEV2( A, B, C, RT1, RT2, CS1, SN1 ) 00002 * 00003 * -- LAPACK auxiliary routine (version 3.2) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 REAL A, B, C, CS1, RT1, RT2, SN1 00010 * .. 00011 * 00012 * Purpose 00013 * ======= 00014 * 00015 * SLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix 00016 * [ A B ] 00017 * [ B C ]. 00018 * On return, RT1 is the eigenvalue of larger absolute value, RT2 is the 00019 * eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right 00020 * eigenvector for RT1, giving the decomposition 00021 * 00022 * [ CS1 SN1 ] [ A B ] [ CS1 -SN1 ] = [ RT1 0 ] 00023 * [-SN1 CS1 ] [ B C ] [ SN1 CS1 ] [ 0 RT2 ]. 00024 * 00025 * Arguments 00026 * ========= 00027 * 00028 * A (input) REAL 00029 * The (1,1) element of the 2-by-2 matrix. 00030 * 00031 * B (input) REAL 00032 * The (1,2) element and the conjugate of the (2,1) element of 00033 * the 2-by-2 matrix. 00034 * 00035 * C (input) REAL 00036 * The (2,2) element of the 2-by-2 matrix. 00037 * 00038 * RT1 (output) REAL 00039 * The eigenvalue of larger absolute value. 00040 * 00041 * RT2 (output) REAL 00042 * The eigenvalue of smaller absolute value. 00043 * 00044 * CS1 (output) REAL 00045 * SN1 (output) REAL 00046 * The vector (CS1, SN1) is a unit right eigenvector for RT1. 00047 * 00048 * Further Details 00049 * =============== 00050 * 00051 * RT1 is accurate to a few ulps barring over/underflow. 00052 * 00053 * RT2 may be inaccurate if there is massive cancellation in the 00054 * determinant A*C-B*B; higher precision or correctly rounded or 00055 * correctly truncated arithmetic would be needed to compute RT2 00056 * accurately in all cases. 00057 * 00058 * CS1 and SN1 are accurate to a few ulps barring over/underflow. 00059 * 00060 * Overflow is possible only if RT1 is within a factor of 5 of overflow. 00061 * Underflow is harmless if the input data is 0 or exceeds 00062 * underflow_threshold / macheps. 00063 * 00064 * ===================================================================== 00065 * 00066 * .. Parameters .. 00067 REAL ONE 00068 PARAMETER ( ONE = 1.0E0 ) 00069 REAL TWO 00070 PARAMETER ( TWO = 2.0E0 ) 00071 REAL ZERO 00072 PARAMETER ( ZERO = 0.0E0 ) 00073 REAL HALF 00074 PARAMETER ( HALF = 0.5E0 ) 00075 * .. 00076 * .. Local Scalars .. 00077 INTEGER SGN1, SGN2 00078 REAL AB, ACMN, ACMX, ACS, ADF, CS, CT, DF, RT, SM, 00079 $ TB, TN 00080 * .. 00081 * .. Intrinsic Functions .. 00082 INTRINSIC ABS, SQRT 00083 * .. 00084 * .. Executable Statements .. 00085 * 00086 * Compute the eigenvalues 00087 * 00088 SM = A + C 00089 DF = A - C 00090 ADF = ABS( DF ) 00091 TB = B + B 00092 AB = ABS( TB ) 00093 IF( ABS( A ).GT.ABS( C ) ) THEN 00094 ACMX = A 00095 ACMN = C 00096 ELSE 00097 ACMX = C 00098 ACMN = A 00099 END IF 00100 IF( ADF.GT.AB ) THEN 00101 RT = ADF*SQRT( ONE+( AB / ADF )**2 ) 00102 ELSE IF( ADF.LT.AB ) THEN 00103 RT = AB*SQRT( ONE+( ADF / AB )**2 ) 00104 ELSE 00105 * 00106 * Includes case AB=ADF=0 00107 * 00108 RT = AB*SQRT( TWO ) 00109 END IF 00110 IF( SM.LT.ZERO ) THEN 00111 RT1 = HALF*( SM-RT ) 00112 SGN1 = -1 00113 * 00114 * Order of execution important. 00115 * To get fully accurate smaller eigenvalue, 00116 * next line needs to be executed in higher precision. 00117 * 00118 RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B 00119 ELSE IF( SM.GT.ZERO ) THEN 00120 RT1 = HALF*( SM+RT ) 00121 SGN1 = 1 00122 * 00123 * Order of execution important. 00124 * To get fully accurate smaller eigenvalue, 00125 * next line needs to be executed in higher precision. 00126 * 00127 RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B 00128 ELSE 00129 * 00130 * Includes case RT1 = RT2 = 0 00131 * 00132 RT1 = HALF*RT 00133 RT2 = -HALF*RT 00134 SGN1 = 1 00135 END IF 00136 * 00137 * Compute the eigenvector 00138 * 00139 IF( DF.GE.ZERO ) THEN 00140 CS = DF + RT 00141 SGN2 = 1 00142 ELSE 00143 CS = DF - RT 00144 SGN2 = -1 00145 END IF 00146 ACS = ABS( CS ) 00147 IF( ACS.GT.AB ) THEN 00148 CT = -TB / CS 00149 SN1 = ONE / SQRT( ONE+CT*CT ) 00150 CS1 = CT*SN1 00151 ELSE 00152 IF( AB.EQ.ZERO ) THEN 00153 CS1 = ONE 00154 SN1 = ZERO 00155 ELSE 00156 TN = -CS / TB 00157 CS1 = ONE / SQRT( ONE+TN*TN ) 00158 SN1 = TN*CS1 00159 END IF 00160 END IF 00161 IF( SGN1.EQ.SGN2 ) THEN 00162 TN = CS1 00163 CS1 = -SN1 00164 SN1 = TN 00165 END IF 00166 RETURN 00167 * 00168 * End of SLAEV2 00169 * 00170 END
1.7.3 
