1 SUBROUTINE DLAEV2( A, B, C, RT1, RT2, CS1, SN1 ) 2* 3* -- LAPACK auxiliary routine (version 3.0) -- 4* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 5* Courant Institute, Argonne National Lab, and Rice University 6* October 31, 1992 7* 8* .. Scalar Arguments .. 9 DOUBLE PRECISION A, B, C, CS1, RT1, RT2, SN1 10* .. 11* 12* Purpose 13* ======= 14* 15* DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix 16* [ A B ] 17* [ B C ]. 18* On return, RT1 is the eigenvalue of larger absolute value, RT2 is the 19* eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right 20* eigenvector for RT1, giving the decomposition 21* 22* [ CS1 SN1 ] [ A B ] [ CS1 -SN1 ] = [ RT1 0 ] 23* [-SN1 CS1 ] [ B C ] [ SN1 CS1 ] [ 0 RT2 ]. 24* 25* Arguments 26* ========= 27* 28* A (input) DOUBLE PRECISION 29* The (1,1) element of the 2-by-2 matrix. 30* 31* B (input) DOUBLE PRECISION 32* The (1,2) element and the conjugate of the (2,1) element of 33* the 2-by-2 matrix. 34* 35* C (input) DOUBLE PRECISION 36* The (2,2) element of the 2-by-2 matrix. 37* 38* RT1 (output) DOUBLE PRECISION 39* The eigenvalue of larger absolute value. 40* 41* RT2 (output) DOUBLE PRECISION 42* The eigenvalue of smaller absolute value. 43* 44* CS1 (output) DOUBLE PRECISION 45* SN1 (output) DOUBLE PRECISION 46* The vector (CS1, SN1) is a unit right eigenvector for RT1. 47* 48* Further Details 49* =============== 50* 51* RT1 is accurate to a few ulps barring over/underflow. 52* 53* RT2 may be inaccurate if there is massive cancellation in the 54* determinant A*C-B*B; higher precision or correctly rounded or 55* correctly truncated arithmetic would be needed to compute RT2 56* accurately in all cases. 57* 58* CS1 and SN1 are accurate to a few ulps barring over/underflow. 59* 60* Overflow is possible only if RT1 is within a factor of 5 of overflow. 61* Underflow is harmless if the input data is 0 or exceeds 62* underflow_threshold / macheps. 63* 64* ===================================================================== 65* 66* .. Parameters .. 67 DOUBLE PRECISION ONE 68 PARAMETER ( ONE = 1.0D0 ) 69 DOUBLE PRECISION TWO 70 PARAMETER ( TWO = 2.0D0 ) 71 DOUBLE PRECISION ZERO 72 PARAMETER ( ZERO = 0.0D0 ) 73 DOUBLE PRECISION HALF 74 PARAMETER ( HALF = 0.5D0 ) 75* .. 76* .. Local Scalars .. 77 INTEGER SGN1, SGN2 78 DOUBLE PRECISION AB, ACMN, ACMX, ACS, ADF, CS, CT, DF, RT, SM, 79 $ TB, TN 80* .. 81* .. Intrinsic Functions .. 82 INTRINSIC ABS, SQRT 83* .. 84* .. Executable Statements .. 85* 86* Compute the eigenvalues 87* 88 SM = A + C 89 DF = A - C 90 ADF = ABS( DF ) 91 TB = B + B 92 AB = ABS( TB ) 93 IF( ABS( A ).GT.ABS( C ) ) THEN 94 ACMX = A 95 ACMN = C 96 ELSE 97 ACMX = C 98 ACMN = A 99 END IF 100 IF( ADF.GT.AB ) THEN 101 RT = ADF*SQRT( ONE+( AB / ADF )**2 ) 102 ELSE IF( ADF.LT.AB ) THEN 103 RT = AB*SQRT( ONE+( ADF / AB )**2 ) 104 ELSE 105* 106* Includes case AB=ADF=0 107* 108 RT = AB*SQRT( TWO ) 109 END IF 110 IF( SM.LT.ZERO ) THEN 111 RT1 = HALF*( SM-RT ) 112 SGN1 = -1 113* 114* Order of execution important. 115* To get fully accurate smaller eigenvalue, 116* next line needs to be executed in higher precision. 117* 118 RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B 119 ELSE IF( SM.GT.ZERO ) THEN 120 RT1 = HALF*( SM+RT ) 121 SGN1 = 1 122* 123* Order of execution important. 124* To get fully accurate smaller eigenvalue, 125* next line needs to be executed in higher precision. 126* 127 RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B 128 ELSE 129* 130* Includes case RT1 = RT2 = 0 131* 132 RT1 = HALF*RT 133 RT2 = -HALF*RT 134 SGN1 = 1 135 END IF 136* 137* Compute the eigenvector 138* 139 IF( DF.GE.ZERO ) THEN 140 CS = DF + RT 141 SGN2 = 1 142 ELSE 143 CS = DF - RT 144 SGN2 = -1 145 END IF 146 ACS = ABS( CS ) 147 IF( ACS.GT.AB ) THEN 148 CT = -TB / CS 149 SN1 = ONE / SQRT( ONE+CT*CT ) 150 CS1 = CT*SN1 151 ELSE 152 IF( AB.EQ.ZERO ) THEN 153 CS1 = ONE 154 SN1 = ZERO 155 ELSE 156 TN = -CS / TB 157 CS1 = ONE / SQRT( ONE+TN*TN ) 158 SN1 = TN*CS1 159 END IF 160 END IF 161 IF( SGN1.EQ.SGN2 ) THEN 162 TN = CS1 163 CS1 = -SN1 164 SN1 = TN 165 END IF 166 RETURN 167* 168* End of DLAEV2 169* 170 END 171