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**T = 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
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**T = 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