|
SUBROUTINE SLAEV2( A, B, C, RT1, RT2, CS1, SN1 )
*
* -- LAPACK auxiliary routine (version 3.2) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2006
*
* .. Scalar Arguments ..
REAL A, B, C, CS1, RT1, RT2, SN1
* ..
*
* Purpose
* =======
*
* SLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix
* [ A B ]
* [ B C ].
* On return, RT1 is the eigenvalue of larger absolute value, RT2 is the
* eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right
* eigenvector for RT1, giving the decomposition
*
* [ CS1 SN1 ] [ A B ] [ CS1 -SN1 ] = [ RT1 0 ]
* [-SN1 CS1 ] [ B C ] [ SN1 CS1 ] [ 0 RT2 ].
*
* Arguments
* =========
*
* A (input) REAL
* The (1,1) element of the 2-by-2 matrix.
*
* B (input) REAL
* The (1,2) element and the conjugate of the (2,1) element of
* the 2-by-2 matrix.
*
* C (input) REAL
* The (2,2) element of the 2-by-2 matrix.
*
* RT1 (output) REAL
* The eigenvalue of larger absolute value.
*
* RT2 (output) REAL
* The eigenvalue of smaller absolute value.
*
* CS1 (output) REAL
* SN1 (output) REAL
* The vector (CS1, SN1) is a unit right eigenvector for RT1.
*
* Further Details
* ===============
*
* RT1 is accurate to a few ulps barring over/underflow.
*
* RT2 may be inaccurate if there is massive cancellation in the
* determinant A*C-B*B; higher precision or correctly rounded or
* correctly truncated arithmetic would be needed to compute RT2
* accurately in all cases.
*
* CS1 and SN1 are accurate to a few ulps barring over/underflow.
*
* Overflow is possible only if RT1 is within a factor of 5 of overflow.
* Underflow is harmless if the input data is 0 or exceeds
* underflow_threshold / macheps.
*
* =====================================================================
*
* .. Parameters ..
REAL ONE
PARAMETER ( ONE = 1.0E0 )
REAL TWO
PARAMETER ( TWO = 2.0E0 )
REAL ZERO
PARAMETER ( ZERO = 0.0E0 )
REAL HALF
PARAMETER ( HALF = 0.5E0 )
* ..
* .. Local Scalars ..
INTEGER SGN1, SGN2
REAL AB, ACMN, ACMX, ACS, ADF, CS, CT, DF, RT, SM,
$ TB, TN
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, SQRT
* ..
* .. Executable Statements ..
*
* Compute the eigenvalues
*
SM = A + C
DF = A - C
ADF = ABS( DF )
TB = B + B
AB = ABS( TB )
IF( ABS( A ).GT.ABS( C ) ) THEN
ACMX = A
ACMN = C
ELSE
ACMX = C
ACMN = A
END IF
IF( ADF.GT.AB ) THEN
RT = ADF*SQRT( ONE+( AB / ADF )**2 )
ELSE IF( ADF.LT.AB ) THEN
RT = AB*SQRT( ONE+( ADF / AB )**2 )
ELSE
*
* Includes case AB=ADF=0
*
RT = AB*SQRT( TWO )
END IF
IF( SM.LT.ZERO ) THEN
RT1 = HALF*( SM-RT )
SGN1 = -1
*
* Order of execution important.
* To get fully accurate smaller eigenvalue,
* next line needs to be executed in higher precision.
*
RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
ELSE IF( SM.GT.ZERO ) THEN
RT1 = HALF*( SM+RT )
SGN1 = 1
*
* Order of execution important.
* To get fully accurate smaller eigenvalue,
* next line needs to be executed in higher precision.
*
RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
ELSE
*
* Includes case RT1 = RT2 = 0
*
RT1 = HALF*RT
RT2 = -HALF*RT
SGN1 = 1
END IF
*
* Compute the eigenvector
*
IF( DF.GE.ZERO ) THEN
CS = DF + RT
SGN2 = 1
ELSE
CS = DF - RT
SGN2 = -1
END IF
ACS = ABS( CS )
IF( ACS.GT.AB ) THEN
CT = -TB / CS
SN1 = ONE / SQRT( ONE+CT*CT )
CS1 = CT*SN1
ELSE
IF( AB.EQ.ZERO ) THEN
CS1 = ONE
SN1 = ZERO
ELSE
TN = -CS / TB
CS1 = ONE / SQRT( ONE+TN*TN )
SN1 = TN*CS1
END IF
END IF
IF( SGN1.EQ.SGN2 ) THEN
TN = CS1
CS1 = -SN1
SN1 = TN
END IF
RETURN
*
* End of SLAEV2
*
END
|