Annotation of rpl/lapack/lapack/dlaev2.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DLAEV2( A, B, C, RT1, RT2, CS1, SN1 )
! 2: *
! 3: * -- LAPACK auxiliary routine (version 3.2) --
! 4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 6: * November 2006
! 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaev2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack