1 SUBROUTINE DGET22( TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, WR,
2 $ WI, WORK, RESULT )
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 CHARACTER TRANSA, TRANSE, TRANSW
10 INTEGER LDA, LDE, N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION A( LDA, * ), E( LDE, * ), RESULT( 2 ), WI( * ),
14 $ WORK( * ), WR( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DGET22 does an eigenvector check.
21 *
22 * The basic test is:
23 *
24 * RESULT(1) = | A E - E W | / ( |A| |E| ulp )
25 *
26 * using the 1-norm. It also tests the normalization of E:
27 *
28 * RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp )
29 * j
30 *
31 * where E(j) is the j-th eigenvector, and m-norm is the max-norm of a
32 * vector. If an eigenvector is complex, as determined from WI(j)
33 * nonzero, then the max-norm of the vector ( er + i*ei ) is the maximum
34 * of
35 * |er(1)| + |ei(1)|, ... , |er(n)| + |ei(n)|
36 *
37 * W is a block diagonal matrix, with a 1 by 1 block for each real
38 * eigenvalue and a 2 by 2 block for each complex conjugate pair.
39 * If eigenvalues j and j+1 are a complex conjugate pair, so that
40 * WR(j) = WR(j+1) = wr and WI(j) = - WI(j+1) = wi, then the 2 by 2
41 * block corresponding to the pair will be:
42 *
43 * ( wr wi )
44 * ( -wi wr )
45 *
46 * Such a block multiplying an n by 2 matrix ( ur ui ) on the right
47 * will be the same as multiplying ur + i*ui by wr + i*wi.
48 *
49 * To handle various schemes for storage of left eigenvectors, there are
50 * options to use A-transpose instead of A, E-transpose instead of E,
51 * and/or W-transpose instead of W.
52 *
53 * Arguments
54 * ==========
55 *
56 * TRANSA (input) CHARACTER*1
57 * Specifies whether or not A is transposed.
58 * = 'N': No transpose
59 * = 'T': Transpose
60 * = 'C': Conjugate transpose (= Transpose)
61 *
62 * TRANSE (input) CHARACTER*1
63 * Specifies whether or not E is transposed.
64 * = 'N': No transpose, eigenvectors are in columns of E
65 * = 'T': Transpose, eigenvectors are in rows of E
66 * = 'C': Conjugate transpose (= Transpose)
67 *
68 * TRANSW (input) CHARACTER*1
69 * Specifies whether or not W is transposed.
70 * = 'N': No transpose
71 * = 'T': Transpose, use -WI(j) instead of WI(j)
72 * = 'C': Conjugate transpose, use -WI(j) instead of WI(j)
73 *
74 * N (input) INTEGER
75 * The order of the matrix A. N >= 0.
76 *
77 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
78 * The matrix whose eigenvectors are in E.
79 *
80 * LDA (input) INTEGER
81 * The leading dimension of the array A. LDA >= max(1,N).
82 *
83 * E (input) DOUBLE PRECISION array, dimension (LDE,N)
84 * The matrix of eigenvectors. If TRANSE = 'N', the eigenvectors
85 * are stored in the columns of E, if TRANSE = 'T' or 'C', the
86 * eigenvectors are stored in the rows of E.
87 *
88 * LDE (input) INTEGER
89 * The leading dimension of the array E. LDE >= max(1,N).
90 *
91 * WR (input) DOUBLE PRECISION array, dimension (N)
92 * WI (input) DOUBLE PRECISION array, dimension (N)
93 * The real and imaginary parts of the eigenvalues of A.
94 * Purely real eigenvalues are indicated by WI(j) = 0.
95 * Complex conjugate pairs are indicated by WR(j)=WR(j+1) and
96 * WI(j) = - WI(j+1) non-zero; the real part is assumed to be
97 * stored in the j-th row/column and the imaginary part in
98 * the (j+1)-th row/column.
99 *
100 * WORK (workspace) DOUBLE PRECISION array, dimension (N*(N+1))
101 *
102 * RESULT (output) DOUBLE PRECISION array, dimension (2)
103 * RESULT(1) = | A E - E W | / ( |A| |E| ulp )
104 * RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp )
105 *
106 * =====================================================================
107 *
108 * .. Parameters ..
109 DOUBLE PRECISION ZERO, ONE
110 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
111 * ..
112 * .. Local Scalars ..
113 CHARACTER NORMA, NORME
114 INTEGER IECOL, IEROW, INCE, IPAIR, ITRNSE, J, JCOL,
115 $ JVEC
116 DOUBLE PRECISION ANORM, ENORM, ENRMAX, ENRMIN, ERRNRM, TEMP1,
117 $ ULP, UNFL
118 * ..
119 * .. Local Arrays ..
120 DOUBLE PRECISION WMAT( 2, 2 )
121 * ..
122 * .. External Functions ..
123 LOGICAL LSAME
124 DOUBLE PRECISION DLAMCH, DLANGE
125 EXTERNAL LSAME, DLAMCH, DLANGE
126 * ..
127 * .. External Subroutines ..
128 EXTERNAL DAXPY, DGEMM, DLASET
129 * ..
130 * .. Intrinsic Functions ..
131 INTRINSIC ABS, DBLE, MAX, MIN
132 * ..
133 * .. Executable Statements ..
134 *
135 * Initialize RESULT (in case N=0)
136 *
137 RESULT( 1 ) = ZERO
138 RESULT( 2 ) = ZERO
139 IF( N.LE.0 )
140 $ RETURN
141 *
142 UNFL = DLAMCH( 'Safe minimum' )
143 ULP = DLAMCH( 'Precision' )
144 *
145 ITRNSE = 0
146 INCE = 1
147 NORMA = 'O'
148 NORME = 'O'
149 *
150 IF( LSAME( TRANSA, 'T' ) .OR. LSAME( TRANSA, 'C' ) ) THEN
151 NORMA = 'I'
152 END IF
153 IF( LSAME( TRANSE, 'T' ) .OR. LSAME( TRANSE, 'C' ) ) THEN
154 NORME = 'I'
155 ITRNSE = 1
156 INCE = LDE
157 END IF
158 *
159 * Check normalization of E
160 *
161 ENRMIN = ONE / ULP
162 ENRMAX = ZERO
163 IF( ITRNSE.EQ.0 ) THEN
164 *
165 * Eigenvectors are column vectors.
166 *
167 IPAIR = 0
168 DO 30 JVEC = 1, N
169 TEMP1 = ZERO
170 IF( IPAIR.EQ.0 .AND. JVEC.LT.N .AND. WI( JVEC ).NE.ZERO )
171 $ IPAIR = 1
172 IF( IPAIR.EQ.1 ) THEN
173 *
174 * Complex eigenvector
175 *
176 DO 10 J = 1, N
177 TEMP1 = MAX( TEMP1, ABS( E( J, JVEC ) )+
178 $ ABS( E( J, JVEC+1 ) ) )
179 10 CONTINUE
180 ENRMIN = MIN( ENRMIN, TEMP1 )
181 ENRMAX = MAX( ENRMAX, TEMP1 )
182 IPAIR = 2
183 ELSE IF( IPAIR.EQ.2 ) THEN
184 IPAIR = 0
185 ELSE
186 *
187 * Real eigenvector
188 *
189 DO 20 J = 1, N
190 TEMP1 = MAX( TEMP1, ABS( E( J, JVEC ) ) )
191 20 CONTINUE
192 ENRMIN = MIN( ENRMIN, TEMP1 )
193 ENRMAX = MAX( ENRMAX, TEMP1 )
194 IPAIR = 0
195 END IF
196 30 CONTINUE
197 *
198 ELSE
199 *
200 * Eigenvectors are row vectors.
201 *
202 DO 40 JVEC = 1, N
203 WORK( JVEC ) = ZERO
204 40 CONTINUE
205 *
206 DO 60 J = 1, N
207 IPAIR = 0
208 DO 50 JVEC = 1, N
209 IF( IPAIR.EQ.0 .AND. JVEC.LT.N .AND. WI( JVEC ).NE.ZERO )
210 $ IPAIR = 1
211 IF( IPAIR.EQ.1 ) THEN
212 WORK( JVEC ) = MAX( WORK( JVEC ),
213 $ ABS( E( J, JVEC ) )+ABS( E( J,
214 $ JVEC+1 ) ) )
215 WORK( JVEC+1 ) = WORK( JVEC )
216 ELSE IF( IPAIR.EQ.2 ) THEN
217 IPAIR = 0
218 ELSE
219 WORK( JVEC ) = MAX( WORK( JVEC ),
220 $ ABS( E( J, JVEC ) ) )
221 IPAIR = 0
222 END IF
223 50 CONTINUE
224 60 CONTINUE
225 *
226 DO 70 JVEC = 1, N
227 ENRMIN = MIN( ENRMIN, WORK( JVEC ) )
228 ENRMAX = MAX( ENRMAX, WORK( JVEC ) )
229 70 CONTINUE
230 END IF
231 *
232 * Norm of A:
233 *
234 ANORM = MAX( DLANGE( NORMA, N, N, A, LDA, WORK ), UNFL )
235 *
236 * Norm of E:
237 *
238 ENORM = MAX( DLANGE( NORME, N, N, E, LDE, WORK ), ULP )
239 *
240 * Norm of error:
241 *
242 * Error = AE - EW
243 *
244 CALL DLASET( 'Full', N, N, ZERO, ZERO, WORK, N )
245 *
246 IPAIR = 0
247 IEROW = 1
248 IECOL = 1
249 *
250 DO 80 JCOL = 1, N
251 IF( ITRNSE.EQ.1 ) THEN
252 IEROW = JCOL
253 ELSE
254 IECOL = JCOL
255 END IF
256 *
257 IF( IPAIR.EQ.0 .AND. WI( JCOL ).NE.ZERO )
258 $ IPAIR = 1
259 *
260 IF( IPAIR.EQ.1 ) THEN
261 WMAT( 1, 1 ) = WR( JCOL )
262 WMAT( 2, 1 ) = -WI( JCOL )
263 WMAT( 1, 2 ) = WI( JCOL )
264 WMAT( 2, 2 ) = WR( JCOL )
265 CALL DGEMM( TRANSE, TRANSW, N, 2, 2, ONE, E( IEROW, IECOL ),
266 $ LDE, WMAT, 2, ZERO, WORK( N*( JCOL-1 )+1 ), N )
267 IPAIR = 2
268 ELSE IF( IPAIR.EQ.2 ) THEN
269 IPAIR = 0
270 *
271 ELSE
272 *
273 CALL DAXPY( N, WR( JCOL ), E( IEROW, IECOL ), INCE,
274 $ WORK( N*( JCOL-1 )+1 ), 1 )
275 IPAIR = 0
276 END IF
277 *
278 80 CONTINUE
279 *
280 CALL DGEMM( TRANSA, TRANSE, N, N, N, ONE, A, LDA, E, LDE, -ONE,
281 $ WORK, N )
282 *
283 ERRNRM = DLANGE( 'One', N, N, WORK, N, WORK( N*N+1 ) ) / ENORM
284 *
285 * Compute RESULT(1) (avoiding under/overflow)
286 *
287 IF( ANORM.GT.ERRNRM ) THEN
288 RESULT( 1 ) = ( ERRNRM / ANORM ) / ULP
289 ELSE
290 IF( ANORM.LT.ONE ) THEN
291 RESULT( 1 ) = ( MIN( ERRNRM, ANORM ) / ANORM ) / ULP
292 ELSE
293 RESULT( 1 ) = MIN( ERRNRM / ANORM, ONE ) / ULP
294 END IF
295 END IF
296 *
297 * Compute RESULT(2) : the normalization error in E.
298 *
299 RESULT( 2 ) = MAX( ABS( ENRMAX-ONE ), ABS( ENRMIN-ONE ) ) /
300 $ ( DBLE( N )*ULP )
301 *
302 RETURN
303 *
304 * End of DGET22
305 *
306 END
2 $ WI, WORK, RESULT )
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 CHARACTER TRANSA, TRANSE, TRANSW
10 INTEGER LDA, LDE, N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION A( LDA, * ), E( LDE, * ), RESULT( 2 ), WI( * ),
14 $ WORK( * ), WR( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DGET22 does an eigenvector check.
21 *
22 * The basic test is:
23 *
24 * RESULT(1) = | A E - E W | / ( |A| |E| ulp )
25 *
26 * using the 1-norm. It also tests the normalization of E:
27 *
28 * RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp )
29 * j
30 *
31 * where E(j) is the j-th eigenvector, and m-norm is the max-norm of a
32 * vector. If an eigenvector is complex, as determined from WI(j)
33 * nonzero, then the max-norm of the vector ( er + i*ei ) is the maximum
34 * of
35 * |er(1)| + |ei(1)|, ... , |er(n)| + |ei(n)|
36 *
37 * W is a block diagonal matrix, with a 1 by 1 block for each real
38 * eigenvalue and a 2 by 2 block for each complex conjugate pair.
39 * If eigenvalues j and j+1 are a complex conjugate pair, so that
40 * WR(j) = WR(j+1) = wr and WI(j) = - WI(j+1) = wi, then the 2 by 2
41 * block corresponding to the pair will be:
42 *
43 * ( wr wi )
44 * ( -wi wr )
45 *
46 * Such a block multiplying an n by 2 matrix ( ur ui ) on the right
47 * will be the same as multiplying ur + i*ui by wr + i*wi.
48 *
49 * To handle various schemes for storage of left eigenvectors, there are
50 * options to use A-transpose instead of A, E-transpose instead of E,
51 * and/or W-transpose instead of W.
52 *
53 * Arguments
54 * ==========
55 *
56 * TRANSA (input) CHARACTER*1
57 * Specifies whether or not A is transposed.
58 * = 'N': No transpose
59 * = 'T': Transpose
60 * = 'C': Conjugate transpose (= Transpose)
61 *
62 * TRANSE (input) CHARACTER*1
63 * Specifies whether or not E is transposed.
64 * = 'N': No transpose, eigenvectors are in columns of E
65 * = 'T': Transpose, eigenvectors are in rows of E
66 * = 'C': Conjugate transpose (= Transpose)
67 *
68 * TRANSW (input) CHARACTER*1
69 * Specifies whether or not W is transposed.
70 * = 'N': No transpose
71 * = 'T': Transpose, use -WI(j) instead of WI(j)
72 * = 'C': Conjugate transpose, use -WI(j) instead of WI(j)
73 *
74 * N (input) INTEGER
75 * The order of the matrix A. N >= 0.
76 *
77 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
78 * The matrix whose eigenvectors are in E.
79 *
80 * LDA (input) INTEGER
81 * The leading dimension of the array A. LDA >= max(1,N).
82 *
83 * E (input) DOUBLE PRECISION array, dimension (LDE,N)
84 * The matrix of eigenvectors. If TRANSE = 'N', the eigenvectors
85 * are stored in the columns of E, if TRANSE = 'T' or 'C', the
86 * eigenvectors are stored in the rows of E.
87 *
88 * LDE (input) INTEGER
89 * The leading dimension of the array E. LDE >= max(1,N).
90 *
91 * WR (input) DOUBLE PRECISION array, dimension (N)
92 * WI (input) DOUBLE PRECISION array, dimension (N)
93 * The real and imaginary parts of the eigenvalues of A.
94 * Purely real eigenvalues are indicated by WI(j) = 0.
95 * Complex conjugate pairs are indicated by WR(j)=WR(j+1) and
96 * WI(j) = - WI(j+1) non-zero; the real part is assumed to be
97 * stored in the j-th row/column and the imaginary part in
98 * the (j+1)-th row/column.
99 *
100 * WORK (workspace) DOUBLE PRECISION array, dimension (N*(N+1))
101 *
102 * RESULT (output) DOUBLE PRECISION array, dimension (2)
103 * RESULT(1) = | A E - E W | / ( |A| |E| ulp )
104 * RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp )
105 *
106 * =====================================================================
107 *
108 * .. Parameters ..
109 DOUBLE PRECISION ZERO, ONE
110 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
111 * ..
112 * .. Local Scalars ..
113 CHARACTER NORMA, NORME
114 INTEGER IECOL, IEROW, INCE, IPAIR, ITRNSE, J, JCOL,
115 $ JVEC
116 DOUBLE PRECISION ANORM, ENORM, ENRMAX, ENRMIN, ERRNRM, TEMP1,
117 $ ULP, UNFL
118 * ..
119 * .. Local Arrays ..
120 DOUBLE PRECISION WMAT( 2, 2 )
121 * ..
122 * .. External Functions ..
123 LOGICAL LSAME
124 DOUBLE PRECISION DLAMCH, DLANGE
125 EXTERNAL LSAME, DLAMCH, DLANGE
126 * ..
127 * .. External Subroutines ..
128 EXTERNAL DAXPY, DGEMM, DLASET
129 * ..
130 * .. Intrinsic Functions ..
131 INTRINSIC ABS, DBLE, MAX, MIN
132 * ..
133 * .. Executable Statements ..
134 *
135 * Initialize RESULT (in case N=0)
136 *
137 RESULT( 1 ) = ZERO
138 RESULT( 2 ) = ZERO
139 IF( N.LE.0 )
140 $ RETURN
141 *
142 UNFL = DLAMCH( 'Safe minimum' )
143 ULP = DLAMCH( 'Precision' )
144 *
145 ITRNSE = 0
146 INCE = 1
147 NORMA = 'O'
148 NORME = 'O'
149 *
150 IF( LSAME( TRANSA, 'T' ) .OR. LSAME( TRANSA, 'C' ) ) THEN
151 NORMA = 'I'
152 END IF
153 IF( LSAME( TRANSE, 'T' ) .OR. LSAME( TRANSE, 'C' ) ) THEN
154 NORME = 'I'
155 ITRNSE = 1
156 INCE = LDE
157 END IF
158 *
159 * Check normalization of E
160 *
161 ENRMIN = ONE / ULP
162 ENRMAX = ZERO
163 IF( ITRNSE.EQ.0 ) THEN
164 *
165 * Eigenvectors are column vectors.
166 *
167 IPAIR = 0
168 DO 30 JVEC = 1, N
169 TEMP1 = ZERO
170 IF( IPAIR.EQ.0 .AND. JVEC.LT.N .AND. WI( JVEC ).NE.ZERO )
171 $ IPAIR = 1
172 IF( IPAIR.EQ.1 ) THEN
173 *
174 * Complex eigenvector
175 *
176 DO 10 J = 1, N
177 TEMP1 = MAX( TEMP1, ABS( E( J, JVEC ) )+
178 $ ABS( E( J, JVEC+1 ) ) )
179 10 CONTINUE
180 ENRMIN = MIN( ENRMIN, TEMP1 )
181 ENRMAX = MAX( ENRMAX, TEMP1 )
182 IPAIR = 2
183 ELSE IF( IPAIR.EQ.2 ) THEN
184 IPAIR = 0
185 ELSE
186 *
187 * Real eigenvector
188 *
189 DO 20 J = 1, N
190 TEMP1 = MAX( TEMP1, ABS( E( J, JVEC ) ) )
191 20 CONTINUE
192 ENRMIN = MIN( ENRMIN, TEMP1 )
193 ENRMAX = MAX( ENRMAX, TEMP1 )
194 IPAIR = 0
195 END IF
196 30 CONTINUE
197 *
198 ELSE
199 *
200 * Eigenvectors are row vectors.
201 *
202 DO 40 JVEC = 1, N
203 WORK( JVEC ) = ZERO
204 40 CONTINUE
205 *
206 DO 60 J = 1, N
207 IPAIR = 0
208 DO 50 JVEC = 1, N
209 IF( IPAIR.EQ.0 .AND. JVEC.LT.N .AND. WI( JVEC ).NE.ZERO )
210 $ IPAIR = 1
211 IF( IPAIR.EQ.1 ) THEN
212 WORK( JVEC ) = MAX( WORK( JVEC ),
213 $ ABS( E( J, JVEC ) )+ABS( E( J,
214 $ JVEC+1 ) ) )
215 WORK( JVEC+1 ) = WORK( JVEC )
216 ELSE IF( IPAIR.EQ.2 ) THEN
217 IPAIR = 0
218 ELSE
219 WORK( JVEC ) = MAX( WORK( JVEC ),
220 $ ABS( E( J, JVEC ) ) )
221 IPAIR = 0
222 END IF
223 50 CONTINUE
224 60 CONTINUE
225 *
226 DO 70 JVEC = 1, N
227 ENRMIN = MIN( ENRMIN, WORK( JVEC ) )
228 ENRMAX = MAX( ENRMAX, WORK( JVEC ) )
229 70 CONTINUE
230 END IF
231 *
232 * Norm of A:
233 *
234 ANORM = MAX( DLANGE( NORMA, N, N, A, LDA, WORK ), UNFL )
235 *
236 * Norm of E:
237 *
238 ENORM = MAX( DLANGE( NORME, N, N, E, LDE, WORK ), ULP )
239 *
240 * Norm of error:
241 *
242 * Error = AE - EW
243 *
244 CALL DLASET( 'Full', N, N, ZERO, ZERO, WORK, N )
245 *
246 IPAIR = 0
247 IEROW = 1
248 IECOL = 1
249 *
250 DO 80 JCOL = 1, N
251 IF( ITRNSE.EQ.1 ) THEN
252 IEROW = JCOL
253 ELSE
254 IECOL = JCOL
255 END IF
256 *
257 IF( IPAIR.EQ.0 .AND. WI( JCOL ).NE.ZERO )
258 $ IPAIR = 1
259 *
260 IF( IPAIR.EQ.1 ) THEN
261 WMAT( 1, 1 ) = WR( JCOL )
262 WMAT( 2, 1 ) = -WI( JCOL )
263 WMAT( 1, 2 ) = WI( JCOL )
264 WMAT( 2, 2 ) = WR( JCOL )
265 CALL DGEMM( TRANSE, TRANSW, N, 2, 2, ONE, E( IEROW, IECOL ),
266 $ LDE, WMAT, 2, ZERO, WORK( N*( JCOL-1 )+1 ), N )
267 IPAIR = 2
268 ELSE IF( IPAIR.EQ.2 ) THEN
269 IPAIR = 0
270 *
271 ELSE
272 *
273 CALL DAXPY( N, WR( JCOL ), E( IEROW, IECOL ), INCE,
274 $ WORK( N*( JCOL-1 )+1 ), 1 )
275 IPAIR = 0
276 END IF
277 *
278 80 CONTINUE
279 *
280 CALL DGEMM( TRANSA, TRANSE, N, N, N, ONE, A, LDA, E, LDE, -ONE,
281 $ WORK, N )
282 *
283 ERRNRM = DLANGE( 'One', N, N, WORK, N, WORK( N*N+1 ) ) / ENORM
284 *
285 * Compute RESULT(1) (avoiding under/overflow)
286 *
287 IF( ANORM.GT.ERRNRM ) THEN
288 RESULT( 1 ) = ( ERRNRM / ANORM ) / ULP
289 ELSE
290 IF( ANORM.LT.ONE ) THEN
291 RESULT( 1 ) = ( MIN( ERRNRM, ANORM ) / ANORM ) / ULP
292 ELSE
293 RESULT( 1 ) = MIN( ERRNRM / ANORM, ONE ) / ULP
294 END IF
295 END IF
296 *
297 * Compute RESULT(2) : the normalization error in E.
298 *
299 RESULT( 2 ) = MAX( ABS( ENRMAX-ONE ), ABS( ENRMIN-ONE ) ) /
300 $ ( DBLE( N )*ULP )
301 *
302 RETURN
303 *
304 * End of DGET22
305 *
306 END