1 SUBROUTINE SLATM6( TYPE, N, A, LDA, B, X, LDX, Y, LDY, ALPHA,
2 $ BETA, WX, WY, S, DIF )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 INTEGER LDA, LDX, LDY, N, TYPE
10 REAL ALPHA, BETA, WX, WY
11 * ..
12 * .. Array Arguments ..
13 REAL A( LDA, * ), B( LDA, * ), DIF( * ), S( * ),
14 $ X( LDX, * ), Y( LDY, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * SLATM6 generates test matrices for the generalized eigenvalue
21 * problem, their corresponding right and left eigenvector matrices,
22 * and also reciprocal condition numbers for all eigenvalues and
23 * the reciprocal condition numbers of eigenvectors corresponding to
24 * the 1th and 5th eigenvalues.
25 *
26 * Test Matrices
27 * =============
28 *
29 * Two kinds of test matrix pairs
30 *
31 * (A, B) = inverse(YH) * (Da, Db) * inverse(X)
32 *
33 * are used in the tests:
34 *
35 * Type 1:
36 * Da = 1+a 0 0 0 0 Db = 1 0 0 0 0
37 * 0 2+a 0 0 0 0 1 0 0 0
38 * 0 0 3+a 0 0 0 0 1 0 0
39 * 0 0 0 4+a 0 0 0 0 1 0
40 * 0 0 0 0 5+a , 0 0 0 0 1 , and
41 *
42 * Type 2:
43 * Da = 1 -1 0 0 0 Db = 1 0 0 0 0
44 * 1 1 0 0 0 0 1 0 0 0
45 * 0 0 1 0 0 0 0 1 0 0
46 * 0 0 0 1+a 1+b 0 0 0 1 0
47 * 0 0 0 -1-b 1+a , 0 0 0 0 1 .
48 *
49 * In both cases the same inverse(YH) and inverse(X) are used to compute
50 * (A, B), giving the exact eigenvectors to (A,B) as (YH, X):
51 *
52 * YH: = 1 0 -y y -y X = 1 0 -x -x x
53 * 0 1 -y y -y 0 1 x -x -x
54 * 0 0 1 0 0 0 0 1 0 0
55 * 0 0 0 1 0 0 0 0 1 0
56 * 0 0 0 0 1, 0 0 0 0 1 ,
57 *
58 * where a, b, x and y will have all values independently of each other.
59 *
60 * Arguments
61 * =========
62 *
63 * TYPE (input) INTEGER
64 * Specifies the problem type (see futher details).
65 *
66 * N (input) INTEGER
67 * Size of the matrices A and B.
68 *
69 * A (output) REAL array, dimension (LDA, N).
70 * On exit A N-by-N is initialized according to TYPE.
71 *
72 * LDA (input) INTEGER
73 * The leading dimension of A and of B.
74 *
75 * B (output) REAL array, dimension (LDA, N).
76 * On exit B N-by-N is initialized according to TYPE.
77 *
78 * X (output) REAL array, dimension (LDX, N).
79 * On exit X is the N-by-N matrix of right eigenvectors.
80 *
81 * LDX (input) INTEGER
82 * The leading dimension of X.
83 *
84 * Y (output) REAL array, dimension (LDY, N).
85 * On exit Y is the N-by-N matrix of left eigenvectors.
86 *
87 * LDY (input) INTEGER
88 * The leading dimension of Y.
89 *
90 * ALPHA (input) REAL
91 * BETA (input) REAL
92 * Weighting constants for matrix A.
93 *
94 * WX (input) REAL
95 * Constant for right eigenvector matrix.
96 *
97 * WY (input) REAL
98 * Constant for left eigenvector matrix.
99 *
100 * S (output) REAL array, dimension (N)
101 * S(i) is the reciprocal condition number for eigenvalue i.
102 *
103 * DIF (output) REAL array, dimension (N)
104 * DIF(i) is the reciprocal condition number for eigenvector i.
105 *
106 * =====================================================================
107 *
108 * .. Parameters ..
109 REAL ZERO, ONE, TWO, THREE
110 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0,
111 $ THREE = 3.0E+0 )
112 * ..
113 * .. Local Scalars ..
114 INTEGER I, INFO, J
115 * ..
116 * .. Local Arrays ..
117 REAL WORK( 100 ), Z( 12, 12 )
118 * ..
119 * .. Intrinsic Functions ..
120 INTRINSIC REAL, SQRT
121 * ..
122 * .. External Subroutines ..
123 EXTERNAL SGESVD, SLACPY, SLAKF2
124 * ..
125 * .. Executable Statements ..
126 *
127 * Generate test problem ...
128 * (Da, Db) ...
129 *
130 DO 20 I = 1, N
131 DO 10 J = 1, N
132 *
133 IF( I.EQ.J ) THEN
134 A( I, I ) = REAL( I ) + ALPHA
135 B( I, I ) = ONE
136 ELSE
137 A( I, J ) = ZERO
138 B( I, J ) = ZERO
139 END IF
140 *
141 10 CONTINUE
142 20 CONTINUE
143 *
144 * Form X and Y
145 *
146 CALL SLACPY( 'F', N, N, B, LDA, Y, LDY )
147 Y( 3, 1 ) = -WY
148 Y( 4, 1 ) = WY
149 Y( 5, 1 ) = -WY
150 Y( 3, 2 ) = -WY
151 Y( 4, 2 ) = WY
152 Y( 5, 2 ) = -WY
153 *
154 CALL SLACPY( 'F', N, N, B, LDA, X, LDX )
155 X( 1, 3 ) = -WX
156 X( 1, 4 ) = -WX
157 X( 1, 5 ) = WX
158 X( 2, 3 ) = WX
159 X( 2, 4 ) = -WX
160 X( 2, 5 ) = -WX
161 *
162 * Form (A, B)
163 *
164 B( 1, 3 ) = WX + WY
165 B( 2, 3 ) = -WX + WY
166 B( 1, 4 ) = WX - WY
167 B( 2, 4 ) = WX - WY
168 B( 1, 5 ) = -WX + WY
169 B( 2, 5 ) = WX + WY
170 IF( TYPE.EQ.1 ) THEN
171 A( 1, 3 ) = WX*A( 1, 1 ) + WY*A( 3, 3 )
172 A( 2, 3 ) = -WX*A( 2, 2 ) + WY*A( 3, 3 )
173 A( 1, 4 ) = WX*A( 1, 1 ) - WY*A( 4, 4 )
174 A( 2, 4 ) = WX*A( 2, 2 ) - WY*A( 4, 4 )
175 A( 1, 5 ) = -WX*A( 1, 1 ) + WY*A( 5, 5 )
176 A( 2, 5 ) = WX*A( 2, 2 ) + WY*A( 5, 5 )
177 ELSE IF( TYPE.EQ.2 ) THEN
178 A( 1, 3 ) = TWO*WX + WY
179 A( 2, 3 ) = WY
180 A( 1, 4 ) = -WY*( TWO+ALPHA+BETA )
181 A( 2, 4 ) = TWO*WX - WY*( TWO+ALPHA+BETA )
182 A( 1, 5 ) = -TWO*WX + WY*( ALPHA-BETA )
183 A( 2, 5 ) = WY*( ALPHA-BETA )
184 A( 1, 1 ) = ONE
185 A( 1, 2 ) = -ONE
186 A( 2, 1 ) = ONE
187 A( 2, 2 ) = A( 1, 1 )
188 A( 3, 3 ) = ONE
189 A( 4, 4 ) = ONE + ALPHA
190 A( 4, 5 ) = ONE + BETA
191 A( 5, 4 ) = -A( 4, 5 )
192 A( 5, 5 ) = A( 4, 4 )
193 END IF
194 *
195 * Compute condition numbers
196 *
197 IF( TYPE.EQ.1 ) THEN
198 *
199 S( 1 ) = ONE / SQRT( ( ONE+THREE*WY*WY ) /
200 $ ( ONE+A( 1, 1 )*A( 1, 1 ) ) )
201 S( 2 ) = ONE / SQRT( ( ONE+THREE*WY*WY ) /
202 $ ( ONE+A( 2, 2 )*A( 2, 2 ) ) )
203 S( 3 ) = ONE / SQRT( ( ONE+TWO*WX*WX ) /
204 $ ( ONE+A( 3, 3 )*A( 3, 3 ) ) )
205 S( 4 ) = ONE / SQRT( ( ONE+TWO*WX*WX ) /
206 $ ( ONE+A( 4, 4 )*A( 4, 4 ) ) )
207 S( 5 ) = ONE / SQRT( ( ONE+TWO*WX*WX ) /
208 $ ( ONE+A( 5, 5 )*A( 5, 5 ) ) )
209 *
210 CALL SLAKF2( 1, 4, A, LDA, A( 2, 2 ), B, B( 2, 2 ), Z, 12 )
211 CALL SGESVD( 'N', 'N', 8, 8, Z, 12, WORK, WORK( 9 ), 1,
212 $ WORK( 10 ), 1, WORK( 11 ), 40, INFO )
213 DIF( 1 ) = WORK( 8 )
214 *
215 CALL SLAKF2( 4, 1, A, LDA, A( 5, 5 ), B, B( 5, 5 ), Z, 12 )
216 CALL SGESVD( 'N', 'N', 8, 8, Z, 12, WORK, WORK( 9 ), 1,
217 $ WORK( 10 ), 1, WORK( 11 ), 40, INFO )
218 DIF( 5 ) = WORK( 8 )
219 *
220 ELSE IF( TYPE.EQ.2 ) THEN
221 *
222 S( 1 ) = ONE / SQRT( ONE / THREE+WY*WY )
223 S( 2 ) = S( 1 )
224 S( 3 ) = ONE / SQRT( ONE / TWO+WX*WX )
225 S( 4 ) = ONE / SQRT( ( ONE+TWO*WX*WX ) /
226 $ ( ONE+( ONE+ALPHA )*( ONE+ALPHA )+( ONE+BETA )*( ONE+
227 $ BETA ) ) )
228 S( 5 ) = S( 4 )
229 *
230 CALL SLAKF2( 2, 3, A, LDA, A( 3, 3 ), B, B( 3, 3 ), Z, 12 )
231 CALL SGESVD( 'N', 'N', 12, 12, Z, 12, WORK, WORK( 13 ), 1,
232 $ WORK( 14 ), 1, WORK( 15 ), 60, INFO )
233 DIF( 1 ) = WORK( 12 )
234 *
235 CALL SLAKF2( 3, 2, A, LDA, A( 4, 4 ), B, B( 4, 4 ), Z, 12 )
236 CALL SGESVD( 'N', 'N', 12, 12, Z, 12, WORK, WORK( 13 ), 1,
237 $ WORK( 14 ), 1, WORK( 15 ), 60, INFO )
238 DIF( 5 ) = WORK( 12 )
239 *
240 END IF
241 *
242 RETURN
243 *
244 * End of SLATM6
245 *
246 END
2 $ BETA, WX, WY, S, DIF )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 INTEGER LDA, LDX, LDY, N, TYPE
10 REAL ALPHA, BETA, WX, WY
11 * ..
12 * .. Array Arguments ..
13 REAL A( LDA, * ), B( LDA, * ), DIF( * ), S( * ),
14 $ X( LDX, * ), Y( LDY, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * SLATM6 generates test matrices for the generalized eigenvalue
21 * problem, their corresponding right and left eigenvector matrices,
22 * and also reciprocal condition numbers for all eigenvalues and
23 * the reciprocal condition numbers of eigenvectors corresponding to
24 * the 1th and 5th eigenvalues.
25 *
26 * Test Matrices
27 * =============
28 *
29 * Two kinds of test matrix pairs
30 *
31 * (A, B) = inverse(YH) * (Da, Db) * inverse(X)
32 *
33 * are used in the tests:
34 *
35 * Type 1:
36 * Da = 1+a 0 0 0 0 Db = 1 0 0 0 0
37 * 0 2+a 0 0 0 0 1 0 0 0
38 * 0 0 3+a 0 0 0 0 1 0 0
39 * 0 0 0 4+a 0 0 0 0 1 0
40 * 0 0 0 0 5+a , 0 0 0 0 1 , and
41 *
42 * Type 2:
43 * Da = 1 -1 0 0 0 Db = 1 0 0 0 0
44 * 1 1 0 0 0 0 1 0 0 0
45 * 0 0 1 0 0 0 0 1 0 0
46 * 0 0 0 1+a 1+b 0 0 0 1 0
47 * 0 0 0 -1-b 1+a , 0 0 0 0 1 .
48 *
49 * In both cases the same inverse(YH) and inverse(X) are used to compute
50 * (A, B), giving the exact eigenvectors to (A,B) as (YH, X):
51 *
52 * YH: = 1 0 -y y -y X = 1 0 -x -x x
53 * 0 1 -y y -y 0 1 x -x -x
54 * 0 0 1 0 0 0 0 1 0 0
55 * 0 0 0 1 0 0 0 0 1 0
56 * 0 0 0 0 1, 0 0 0 0 1 ,
57 *
58 * where a, b, x and y will have all values independently of each other.
59 *
60 * Arguments
61 * =========
62 *
63 * TYPE (input) INTEGER
64 * Specifies the problem type (see futher details).
65 *
66 * N (input) INTEGER
67 * Size of the matrices A and B.
68 *
69 * A (output) REAL array, dimension (LDA, N).
70 * On exit A N-by-N is initialized according to TYPE.
71 *
72 * LDA (input) INTEGER
73 * The leading dimension of A and of B.
74 *
75 * B (output) REAL array, dimension (LDA, N).
76 * On exit B N-by-N is initialized according to TYPE.
77 *
78 * X (output) REAL array, dimension (LDX, N).
79 * On exit X is the N-by-N matrix of right eigenvectors.
80 *
81 * LDX (input) INTEGER
82 * The leading dimension of X.
83 *
84 * Y (output) REAL array, dimension (LDY, N).
85 * On exit Y is the N-by-N matrix of left eigenvectors.
86 *
87 * LDY (input) INTEGER
88 * The leading dimension of Y.
89 *
90 * ALPHA (input) REAL
91 * BETA (input) REAL
92 * Weighting constants for matrix A.
93 *
94 * WX (input) REAL
95 * Constant for right eigenvector matrix.
96 *
97 * WY (input) REAL
98 * Constant for left eigenvector matrix.
99 *
100 * S (output) REAL array, dimension (N)
101 * S(i) is the reciprocal condition number for eigenvalue i.
102 *
103 * DIF (output) REAL array, dimension (N)
104 * DIF(i) is the reciprocal condition number for eigenvector i.
105 *
106 * =====================================================================
107 *
108 * .. Parameters ..
109 REAL ZERO, ONE, TWO, THREE
110 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0,
111 $ THREE = 3.0E+0 )
112 * ..
113 * .. Local Scalars ..
114 INTEGER I, INFO, J
115 * ..
116 * .. Local Arrays ..
117 REAL WORK( 100 ), Z( 12, 12 )
118 * ..
119 * .. Intrinsic Functions ..
120 INTRINSIC REAL, SQRT
121 * ..
122 * .. External Subroutines ..
123 EXTERNAL SGESVD, SLACPY, SLAKF2
124 * ..
125 * .. Executable Statements ..
126 *
127 * Generate test problem ...
128 * (Da, Db) ...
129 *
130 DO 20 I = 1, N
131 DO 10 J = 1, N
132 *
133 IF( I.EQ.J ) THEN
134 A( I, I ) = REAL( I ) + ALPHA
135 B( I, I ) = ONE
136 ELSE
137 A( I, J ) = ZERO
138 B( I, J ) = ZERO
139 END IF
140 *
141 10 CONTINUE
142 20 CONTINUE
143 *
144 * Form X and Y
145 *
146 CALL SLACPY( 'F', N, N, B, LDA, Y, LDY )
147 Y( 3, 1 ) = -WY
148 Y( 4, 1 ) = WY
149 Y( 5, 1 ) = -WY
150 Y( 3, 2 ) = -WY
151 Y( 4, 2 ) = WY
152 Y( 5, 2 ) = -WY
153 *
154 CALL SLACPY( 'F', N, N, B, LDA, X, LDX )
155 X( 1, 3 ) = -WX
156 X( 1, 4 ) = -WX
157 X( 1, 5 ) = WX
158 X( 2, 3 ) = WX
159 X( 2, 4 ) = -WX
160 X( 2, 5 ) = -WX
161 *
162 * Form (A, B)
163 *
164 B( 1, 3 ) = WX + WY
165 B( 2, 3 ) = -WX + WY
166 B( 1, 4 ) = WX - WY
167 B( 2, 4 ) = WX - WY
168 B( 1, 5 ) = -WX + WY
169 B( 2, 5 ) = WX + WY
170 IF( TYPE.EQ.1 ) THEN
171 A( 1, 3 ) = WX*A( 1, 1 ) + WY*A( 3, 3 )
172 A( 2, 3 ) = -WX*A( 2, 2 ) + WY*A( 3, 3 )
173 A( 1, 4 ) = WX*A( 1, 1 ) - WY*A( 4, 4 )
174 A( 2, 4 ) = WX*A( 2, 2 ) - WY*A( 4, 4 )
175 A( 1, 5 ) = -WX*A( 1, 1 ) + WY*A( 5, 5 )
176 A( 2, 5 ) = WX*A( 2, 2 ) + WY*A( 5, 5 )
177 ELSE IF( TYPE.EQ.2 ) THEN
178 A( 1, 3 ) = TWO*WX + WY
179 A( 2, 3 ) = WY
180 A( 1, 4 ) = -WY*( TWO+ALPHA+BETA )
181 A( 2, 4 ) = TWO*WX - WY*( TWO+ALPHA+BETA )
182 A( 1, 5 ) = -TWO*WX + WY*( ALPHA-BETA )
183 A( 2, 5 ) = WY*( ALPHA-BETA )
184 A( 1, 1 ) = ONE
185 A( 1, 2 ) = -ONE
186 A( 2, 1 ) = ONE
187 A( 2, 2 ) = A( 1, 1 )
188 A( 3, 3 ) = ONE
189 A( 4, 4 ) = ONE + ALPHA
190 A( 4, 5 ) = ONE + BETA
191 A( 5, 4 ) = -A( 4, 5 )
192 A( 5, 5 ) = A( 4, 4 )
193 END IF
194 *
195 * Compute condition numbers
196 *
197 IF( TYPE.EQ.1 ) THEN
198 *
199 S( 1 ) = ONE / SQRT( ( ONE+THREE*WY*WY ) /
200 $ ( ONE+A( 1, 1 )*A( 1, 1 ) ) )
201 S( 2 ) = ONE / SQRT( ( ONE+THREE*WY*WY ) /
202 $ ( ONE+A( 2, 2 )*A( 2, 2 ) ) )
203 S( 3 ) = ONE / SQRT( ( ONE+TWO*WX*WX ) /
204 $ ( ONE+A( 3, 3 )*A( 3, 3 ) ) )
205 S( 4 ) = ONE / SQRT( ( ONE+TWO*WX*WX ) /
206 $ ( ONE+A( 4, 4 )*A( 4, 4 ) ) )
207 S( 5 ) = ONE / SQRT( ( ONE+TWO*WX*WX ) /
208 $ ( ONE+A( 5, 5 )*A( 5, 5 ) ) )
209 *
210 CALL SLAKF2( 1, 4, A, LDA, A( 2, 2 ), B, B( 2, 2 ), Z, 12 )
211 CALL SGESVD( 'N', 'N', 8, 8, Z, 12, WORK, WORK( 9 ), 1,
212 $ WORK( 10 ), 1, WORK( 11 ), 40, INFO )
213 DIF( 1 ) = WORK( 8 )
214 *
215 CALL SLAKF2( 4, 1, A, LDA, A( 5, 5 ), B, B( 5, 5 ), Z, 12 )
216 CALL SGESVD( 'N', 'N', 8, 8, Z, 12, WORK, WORK( 9 ), 1,
217 $ WORK( 10 ), 1, WORK( 11 ), 40, INFO )
218 DIF( 5 ) = WORK( 8 )
219 *
220 ELSE IF( TYPE.EQ.2 ) THEN
221 *
222 S( 1 ) = ONE / SQRT( ONE / THREE+WY*WY )
223 S( 2 ) = S( 1 )
224 S( 3 ) = ONE / SQRT( ONE / TWO+WX*WX )
225 S( 4 ) = ONE / SQRT( ( ONE+TWO*WX*WX ) /
226 $ ( ONE+( ONE+ALPHA )*( ONE+ALPHA )+( ONE+BETA )*( ONE+
227 $ BETA ) ) )
228 S( 5 ) = S( 4 )
229 *
230 CALL SLAKF2( 2, 3, A, LDA, A( 3, 3 ), B, B( 3, 3 ), Z, 12 )
231 CALL SGESVD( 'N', 'N', 12, 12, Z, 12, WORK, WORK( 13 ), 1,
232 $ WORK( 14 ), 1, WORK( 15 ), 60, INFO )
233 DIF( 1 ) = WORK( 12 )
234 *
235 CALL SLAKF2( 3, 2, A, LDA, A( 4, 4 ), B, B( 4, 4 ), Z, 12 )
236 CALL SGESVD( 'N', 'N', 12, 12, Z, 12, WORK, WORK( 13 ), 1,
237 $ WORK( 14 ), 1, WORK( 15 ), 60, INFO )
238 DIF( 5 ) = WORK( 12 )
239 *
240 END IF
241 *
242 RETURN
243 *
244 * End of SLATM6
245 *
246 END