Annotation of rpl/lapack/lapack/zlaesy.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZLAESY( A, B, C, RT1, RT2, EVSCAL, 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: COMPLEX*16 A, B, C, CS1, EVSCAL, RT1, RT2, SN1
! 10: * ..
! 11: *
! 12: * Purpose
! 13: * =======
! 14: *
! 15: * ZLAESY computes the eigendecomposition of a 2-by-2 symmetric matrix
! 16: * ( ( A, B );( B, C ) )
! 17: * provided the norm of the matrix of eigenvectors is larger than
! 18: * some threshold value.
! 19: *
! 20: * RT1 is the eigenvalue of larger absolute value, and RT2 of
! 21: * smaller absolute value. If the eigenvectors are computed, then
! 22: * on return ( CS1, SN1 ) is the unit eigenvector for RT1, hence
! 23: *
! 24: * [ CS1 SN1 ] . [ A B ] . [ CS1 -SN1 ] = [ RT1 0 ]
! 25: * [ -SN1 CS1 ] [ B C ] [ SN1 CS1 ] [ 0 RT2 ]
! 26: *
! 27: * Arguments
! 28: * =========
! 29: *
! 30: * A (input) COMPLEX*16
! 31: * The ( 1, 1 ) element of input matrix.
! 32: *
! 33: * B (input) COMPLEX*16
! 34: * The ( 1, 2 ) element of input matrix. The ( 2, 1 ) element
! 35: * is also given by B, since the 2-by-2 matrix is symmetric.
! 36: *
! 37: * C (input) COMPLEX*16
! 38: * The ( 2, 2 ) element of input matrix.
! 39: *
! 40: * RT1 (output) COMPLEX*16
! 41: * The eigenvalue of larger modulus.
! 42: *
! 43: * RT2 (output) COMPLEX*16
! 44: * The eigenvalue of smaller modulus.
! 45: *
! 46: * EVSCAL (output) COMPLEX*16
! 47: * The complex value by which the eigenvector matrix was scaled
! 48: * to make it orthonormal. If EVSCAL is zero, the eigenvectors
! 49: * were not computed. This means one of two things: the 2-by-2
! 50: * matrix could not be diagonalized, or the norm of the matrix
! 51: * of eigenvectors before scaling was larger than the threshold
! 52: * value THRESH (set below).
! 53: *
! 54: * CS1 (output) COMPLEX*16
! 55: * SN1 (output) COMPLEX*16
! 56: * If EVSCAL .NE. 0, ( CS1, SN1 ) is the unit right eigenvector
! 57: * for RT1.
! 58: *
! 59: * =====================================================================
! 60: *
! 61: * .. Parameters ..
! 62: DOUBLE PRECISION ZERO
! 63: PARAMETER ( ZERO = 0.0D0 )
! 64: DOUBLE PRECISION ONE
! 65: PARAMETER ( ONE = 1.0D0 )
! 66: COMPLEX*16 CONE
! 67: PARAMETER ( CONE = ( 1.0D0, 0.0D0 ) )
! 68: DOUBLE PRECISION HALF
! 69: PARAMETER ( HALF = 0.5D0 )
! 70: DOUBLE PRECISION THRESH
! 71: PARAMETER ( THRESH = 0.1D0 )
! 72: * ..
! 73: * .. Local Scalars ..
! 74: DOUBLE PRECISION BABS, EVNORM, TABS, Z
! 75: COMPLEX*16 S, T, TMP
! 76: * ..
! 77: * .. Intrinsic Functions ..
! 78: INTRINSIC ABS, MAX, SQRT
! 79: * ..
! 80: * .. Executable Statements ..
! 81: *
! 82: *
! 83: * Special case: The matrix is actually diagonal.
! 84: * To avoid divide by zero later, we treat this case separately.
! 85: *
! 86: IF( ABS( B ).EQ.ZERO ) THEN
! 87: RT1 = A
! 88: RT2 = C
! 89: IF( ABS( RT1 ).LT.ABS( RT2 ) ) THEN
! 90: TMP = RT1
! 91: RT1 = RT2
! 92: RT2 = TMP
! 93: CS1 = ZERO
! 94: SN1 = ONE
! 95: ELSE
! 96: CS1 = ONE
! 97: SN1 = ZERO
! 98: END IF
! 99: ELSE
! 100: *
! 101: * Compute the eigenvalues and eigenvectors.
! 102: * The characteristic equation is
! 103: * lambda **2 - (A+C) lambda + (A*C - B*B)
! 104: * and we solve it using the quadratic formula.
! 105: *
! 106: S = ( A+C )*HALF
! 107: T = ( A-C )*HALF
! 108: *
! 109: * Take the square root carefully to avoid over/under flow.
! 110: *
! 111: BABS = ABS( B )
! 112: TABS = ABS( T )
! 113: Z = MAX( BABS, TABS )
! 114: IF( Z.GT.ZERO )
! 115: $ T = Z*SQRT( ( T / Z )**2+( B / Z )**2 )
! 116: *
! 117: * Compute the two eigenvalues. RT1 and RT2 are exchanged
! 118: * if necessary so that RT1 will have the greater magnitude.
! 119: *
! 120: RT1 = S + T
! 121: RT2 = S - T
! 122: IF( ABS( RT1 ).LT.ABS( RT2 ) ) THEN
! 123: TMP = RT1
! 124: RT1 = RT2
! 125: RT2 = TMP
! 126: END IF
! 127: *
! 128: * Choose CS1 = 1 and SN1 to satisfy the first equation, then
! 129: * scale the components of this eigenvector so that the matrix
! 130: * of eigenvectors X satisfies X * X' = I . (No scaling is
! 131: * done if the norm of the eigenvalue matrix is less than THRESH.)
! 132: *
! 133: SN1 = ( RT1-A ) / B
! 134: TABS = ABS( SN1 )
! 135: IF( TABS.GT.ONE ) THEN
! 136: T = TABS*SQRT( ( ONE / TABS )**2+( SN1 / TABS )**2 )
! 137: ELSE
! 138: T = SQRT( CONE+SN1*SN1 )
! 139: END IF
! 140: EVNORM = ABS( T )
! 141: IF( EVNORM.GE.THRESH ) THEN
! 142: EVSCAL = CONE / T
! 143: CS1 = EVSCAL
! 144: SN1 = SN1*EVSCAL
! 145: ELSE
! 146: EVSCAL = ZERO
! 147: END IF
! 148: END IF
! 149: RETURN
! 150: *
! 151: * End of ZLAESY
! 152: *
! 153: END
CVSweb interface <joel.bertrand@systella.fr>