 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DLAE2( A, B, C, RT1, RT2 ) 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 DOUBLE PRECISION A, B, C, RT1, RT2 00010 * .. 00011 * 00012 * Purpose 00013 * ======= 00014 * 00015 * DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix 00016 * [ A B ] 00017 * [ B C ]. 00018 * On return, RT1 is the eigenvalue of larger absolute value, and RT2 00019 * is the eigenvalue of smaller absolute value. 00020 * 00021 * Arguments 00022 * ========= 00023 * 00024 * A (input) DOUBLE PRECISION 00025 * The (1,1) element of the 2-by-2 matrix. 00026 * 00027 * B (input) DOUBLE PRECISION 00028 * The (1,2) and (2,1) elements of the 2-by-2 matrix. 00029 * 00030 * C (input) DOUBLE PRECISION 00031 * The (2,2) element of the 2-by-2 matrix. 00032 * 00033 * RT1 (output) DOUBLE PRECISION 00034 * The eigenvalue of larger absolute value. 00035 * 00036 * RT2 (output) DOUBLE PRECISION 00037 * The eigenvalue of smaller absolute value. 00038 * 00039 * Further Details 00040 * =============== 00041 * 00042 * RT1 is accurate to a few ulps barring over/underflow. 00043 * 00044 * RT2 may be inaccurate if there is massive cancellation in the 00045 * determinant A*C-B*B; higher precision or correctly rounded or 00046 * correctly truncated arithmetic would be needed to compute RT2 00047 * accurately in all cases. 00048 * 00049 * Overflow is possible only if RT1 is within a factor of 5 of overflow. 00050 * Underflow is harmless if the input data is 0 or exceeds 00051 * underflow_threshold / macheps. 00052 * 00053 * ===================================================================== 00054 * 00055 * .. Parameters .. 00056 DOUBLE PRECISION ONE 00057 PARAMETER ( ONE = 1.0D0 ) 00058 DOUBLE PRECISION TWO 00059 PARAMETER ( TWO = 2.0D0 ) 00060 DOUBLE PRECISION ZERO 00061 PARAMETER ( ZERO = 0.0D0 ) 00062 DOUBLE PRECISION HALF 00063 PARAMETER ( HALF = 0.5D0 ) 00064 * .. 00065 * .. Local Scalars .. 00066 DOUBLE PRECISION AB, ACMN, ACMX, ADF, DF, RT, SM, TB 00067 * .. 00068 * .. Intrinsic Functions .. 00069 INTRINSIC ABS, SQRT 00070 * .. 00071 * .. Executable Statements .. 00072 * 00073 * Compute the eigenvalues 00074 * 00075 SM = A + C 00076 DF = A - C 00077 ADF = ABS( DF ) 00078 TB = B + B 00079 AB = ABS( TB ) 00080 IF( ABS( A ).GT.ABS( C ) ) THEN 00081 ACMX = A 00082 ACMN = C 00083 ELSE 00084 ACMX = C 00085 ACMN = A 00086 END IF 00087 IF( ADF.GT.AB ) THEN 00088 RT = ADF*SQRT( ONE+( AB / ADF )**2 ) 00089 ELSE IF( ADF.LT.AB ) THEN 00090 RT = AB*SQRT( ONE+( ADF / AB )**2 ) 00091 ELSE 00092 * 00093 * Includes case AB=ADF=0 00094 * 00095 RT = AB*SQRT( TWO ) 00096 END IF 00097 IF( SM.LT.ZERO ) THEN 00098 RT1 = HALF*( SM-RT ) 00099 * 00100 * Order of execution important. 00101 * To get fully accurate smaller eigenvalue, 00102 * next line needs to be executed in higher precision. 00103 * 00104 RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B 00105 ELSE IF( SM.GT.ZERO ) THEN 00106 RT1 = HALF*( SM+RT ) 00107 * 00108 * Order of execution important. 00109 * To get fully accurate smaller eigenvalue, 00110 * next line needs to be executed in higher precision. 00111 * 00112 RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B 00113 ELSE 00114 * 00115 * Includes case RT1 = RT2 = 0 00116 * 00117 RT1 = HALF*RT 00118 RT2 = -HALF*RT 00119 END IF 00120 RETURN 00121 * 00122 * End of DLAE2 00123 * 00124 END
1.7.3 
