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