1 SUBROUTINE CGET22( TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, W,
2 $ WORK, RWORK, 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 REAL RESULT( 2 ), RWORK( * )
14 COMPLEX A( LDA, * ), E( LDE, * ), W( * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * CGET22 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. The max-norm of a complex n-vector x in this case is the
33 * maximum of |re(x(i)| + |im(x(i)| over i = 1, ..., n.
34 *
35 * Arguments
36 * ==========
37 *
38 * TRANSA (input) CHARACTER*1
39 * Specifies whether or not A is transposed.
40 * = 'N': No transpose
41 * = 'T': Transpose
42 * = 'C': Conjugate transpose
43 *
44 * TRANSE (input) CHARACTER*1
45 * Specifies whether or not E is transposed.
46 * = 'N': No transpose, eigenvectors are in columns of E
47 * = 'T': Transpose, eigenvectors are in rows of E
48 * = 'C': Conjugate transpose, eigenvectors are in rows of E
49 *
50 * TRANSW (input) CHARACTER*1
51 * Specifies whether or not W is transposed.
52 * = 'N': No transpose
53 * = 'T': Transpose, same as TRANSW = 'N'
54 * = 'C': Conjugate transpose, use -WI(j) instead of WI(j)
55 *
56 * N (input) INTEGER
57 * The order of the matrix A. N >= 0.
58 *
59 * A (input) COMPLEX array, dimension (LDA,N)
60 * The matrix whose eigenvectors are in E.
61 *
62 * LDA (input) INTEGER
63 * The leading dimension of the array A. LDA >= max(1,N).
64 *
65 * E (input) COMPLEX array, dimension (LDE,N)
66 * The matrix of eigenvectors. If TRANSE = 'N', the eigenvectors
67 * are stored in the columns of E, if TRANSE = 'T' or 'C', the
68 * eigenvectors are stored in the rows of E.
69 *
70 * LDE (input) INTEGER
71 * The leading dimension of the array E. LDE >= max(1,N).
72 *
73 * W (input) COMPLEX array, dimension (N)
74 * The eigenvalues of A.
75 *
76 * WORK (workspace) COMPLEX array, dimension (N*N)
77 *
78 * RWORK (workspace) REAL array, dimension (N)
79 *
80 * RESULT (output) REAL array, dimension (2)
81 * RESULT(1) = | A E - E W | / ( |A| |E| ulp )
82 * RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp )
83 * j
84 * =====================================================================
85 *
86 * .. Parameters ..
87 REAL ZERO, ONE
88 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
89 COMPLEX CZERO, CONE
90 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ),
91 $ CONE = ( 1.0E+0, 0.0E+0 ) )
92 * ..
93 * .. Local Scalars ..
94 CHARACTER NORMA, NORME
95 INTEGER ITRNSE, ITRNSW, J, JCOL, JOFF, JROW, JVEC
96 REAL ANORM, ENORM, ENRMAX, ENRMIN, ERRNRM, TEMP1,
97 $ ULP, UNFL
98 COMPLEX WTEMP
99 * ..
100 * .. External Functions ..
101 LOGICAL LSAME
102 REAL CLANGE, SLAMCH
103 EXTERNAL LSAME, CLANGE, SLAMCH
104 * ..
105 * .. External Subroutines ..
106 EXTERNAL CGEMM, CLASET
107 * ..
108 * .. Intrinsic Functions ..
109 INTRINSIC ABS, AIMAG, CONJG, MAX, MIN, REAL
110 * ..
111 * .. Executable Statements ..
112 *
113 * Initialize RESULT (in case N=0)
114 *
115 RESULT( 1 ) = ZERO
116 RESULT( 2 ) = ZERO
117 IF( N.LE.0 )
118 $ RETURN
119 *
120 UNFL = SLAMCH( 'Safe minimum' )
121 ULP = SLAMCH( 'Precision' )
122 *
123 ITRNSE = 0
124 ITRNSW = 0
125 NORMA = 'O'
126 NORME = 'O'
127 *
128 IF( LSAME( TRANSA, 'T' ) .OR. LSAME( TRANSA, 'C' ) ) THEN
129 NORMA = 'I'
130 END IF
131 *
132 IF( LSAME( TRANSE, 'T' ) ) THEN
133 ITRNSE = 1
134 NORME = 'I'
135 ELSE IF( LSAME( TRANSE, 'C' ) ) THEN
136 ITRNSE = 2
137 NORME = 'I'
138 END IF
139 *
140 IF( LSAME( TRANSW, 'C' ) ) THEN
141 ITRNSW = 1
142 END IF
143 *
144 * Normalization of E:
145 *
146 ENRMIN = ONE / ULP
147 ENRMAX = ZERO
148 IF( ITRNSE.EQ.0 ) THEN
149 DO 20 JVEC = 1, N
150 TEMP1 = ZERO
151 DO 10 J = 1, N
152 TEMP1 = MAX( TEMP1, ABS( REAL( E( J, JVEC ) ) )+
153 $ ABS( AIMAG( E( J, JVEC ) ) ) )
154 10 CONTINUE
155 ENRMIN = MIN( ENRMIN, TEMP1 )
156 ENRMAX = MAX( ENRMAX, TEMP1 )
157 20 CONTINUE
158 ELSE
159 DO 30 JVEC = 1, N
160 RWORK( JVEC ) = ZERO
161 30 CONTINUE
162 *
163 DO 50 J = 1, N
164 DO 40 JVEC = 1, N
165 RWORK( JVEC ) = MAX( RWORK( JVEC ),
166 $ ABS( REAL( E( JVEC, J ) ) )+
167 $ ABS( AIMAG( E( JVEC, J ) ) ) )
168 40 CONTINUE
169 50 CONTINUE
170 *
171 DO 60 JVEC = 1, N
172 ENRMIN = MIN( ENRMIN, RWORK( JVEC ) )
173 ENRMAX = MAX( ENRMAX, RWORK( JVEC ) )
174 60 CONTINUE
175 END IF
176 *
177 * Norm of A:
178 *
179 ANORM = MAX( CLANGE( NORMA, N, N, A, LDA, RWORK ), UNFL )
180 *
181 * Norm of E:
182 *
183 ENORM = MAX( CLANGE( NORME, N, N, E, LDE, RWORK ), ULP )
184 *
185 * Norm of error:
186 *
187 * Error = AE - EW
188 *
189 CALL CLASET( 'Full', N, N, CZERO, CZERO, WORK, N )
190 *
191 JOFF = 0
192 DO 100 JCOL = 1, N
193 IF( ITRNSW.EQ.0 ) THEN
194 WTEMP = W( JCOL )
195 ELSE
196 WTEMP = CONJG( W( JCOL ) )
197 END IF
198 *
199 IF( ITRNSE.EQ.0 ) THEN
200 DO 70 JROW = 1, N
201 WORK( JOFF+JROW ) = E( JROW, JCOL )*WTEMP
202 70 CONTINUE
203 ELSE IF( ITRNSE.EQ.1 ) THEN
204 DO 80 JROW = 1, N
205 WORK( JOFF+JROW ) = E( JCOL, JROW )*WTEMP
206 80 CONTINUE
207 ELSE
208 DO 90 JROW = 1, N
209 WORK( JOFF+JROW ) = CONJG( E( JCOL, JROW ) )*WTEMP
210 90 CONTINUE
211 END IF
212 JOFF = JOFF + N
213 100 CONTINUE
214 *
215 CALL CGEMM( TRANSA, TRANSE, N, N, N, CONE, A, LDA, E, LDE, -CONE,
216 $ WORK, N )
217 *
218 ERRNRM = CLANGE( 'One', N, N, WORK, N, RWORK ) / ENORM
219 *
220 * Compute RESULT(1) (avoiding under/overflow)
221 *
222 IF( ANORM.GT.ERRNRM ) THEN
223 RESULT( 1 ) = ( ERRNRM / ANORM ) / ULP
224 ELSE
225 IF( ANORM.LT.ONE ) THEN
226 RESULT( 1 ) = ( MIN( ERRNRM, ANORM ) / ANORM ) / ULP
227 ELSE
228 RESULT( 1 ) = MIN( ERRNRM / ANORM, ONE ) / ULP
229 END IF
230 END IF
231 *
232 * Compute RESULT(2) : the normalization error in E.
233 *
234 RESULT( 2 ) = MAX( ABS( ENRMAX-ONE ), ABS( ENRMIN-ONE ) ) /
235 $ ( REAL( N )*ULP )
236 *
237 RETURN
238 *
239 * End of CGET22
240 *
241 END
2 $ WORK, RWORK, 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 REAL RESULT( 2 ), RWORK( * )
14 COMPLEX A( LDA, * ), E( LDE, * ), W( * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * CGET22 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. The max-norm of a complex n-vector x in this case is the
33 * maximum of |re(x(i)| + |im(x(i)| over i = 1, ..., n.
34 *
35 * Arguments
36 * ==========
37 *
38 * TRANSA (input) CHARACTER*1
39 * Specifies whether or not A is transposed.
40 * = 'N': No transpose
41 * = 'T': Transpose
42 * = 'C': Conjugate transpose
43 *
44 * TRANSE (input) CHARACTER*1
45 * Specifies whether or not E is transposed.
46 * = 'N': No transpose, eigenvectors are in columns of E
47 * = 'T': Transpose, eigenvectors are in rows of E
48 * = 'C': Conjugate transpose, eigenvectors are in rows of E
49 *
50 * TRANSW (input) CHARACTER*1
51 * Specifies whether or not W is transposed.
52 * = 'N': No transpose
53 * = 'T': Transpose, same as TRANSW = 'N'
54 * = 'C': Conjugate transpose, use -WI(j) instead of WI(j)
55 *
56 * N (input) INTEGER
57 * The order of the matrix A. N >= 0.
58 *
59 * A (input) COMPLEX array, dimension (LDA,N)
60 * The matrix whose eigenvectors are in E.
61 *
62 * LDA (input) INTEGER
63 * The leading dimension of the array A. LDA >= max(1,N).
64 *
65 * E (input) COMPLEX array, dimension (LDE,N)
66 * The matrix of eigenvectors. If TRANSE = 'N', the eigenvectors
67 * are stored in the columns of E, if TRANSE = 'T' or 'C', the
68 * eigenvectors are stored in the rows of E.
69 *
70 * LDE (input) INTEGER
71 * The leading dimension of the array E. LDE >= max(1,N).
72 *
73 * W (input) COMPLEX array, dimension (N)
74 * The eigenvalues of A.
75 *
76 * WORK (workspace) COMPLEX array, dimension (N*N)
77 *
78 * RWORK (workspace) REAL array, dimension (N)
79 *
80 * RESULT (output) REAL array, dimension (2)
81 * RESULT(1) = | A E - E W | / ( |A| |E| ulp )
82 * RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp )
83 * j
84 * =====================================================================
85 *
86 * .. Parameters ..
87 REAL ZERO, ONE
88 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
89 COMPLEX CZERO, CONE
90 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ),
91 $ CONE = ( 1.0E+0, 0.0E+0 ) )
92 * ..
93 * .. Local Scalars ..
94 CHARACTER NORMA, NORME
95 INTEGER ITRNSE, ITRNSW, J, JCOL, JOFF, JROW, JVEC
96 REAL ANORM, ENORM, ENRMAX, ENRMIN, ERRNRM, TEMP1,
97 $ ULP, UNFL
98 COMPLEX WTEMP
99 * ..
100 * .. External Functions ..
101 LOGICAL LSAME
102 REAL CLANGE, SLAMCH
103 EXTERNAL LSAME, CLANGE, SLAMCH
104 * ..
105 * .. External Subroutines ..
106 EXTERNAL CGEMM, CLASET
107 * ..
108 * .. Intrinsic Functions ..
109 INTRINSIC ABS, AIMAG, CONJG, MAX, MIN, REAL
110 * ..
111 * .. Executable Statements ..
112 *
113 * Initialize RESULT (in case N=0)
114 *
115 RESULT( 1 ) = ZERO
116 RESULT( 2 ) = ZERO
117 IF( N.LE.0 )
118 $ RETURN
119 *
120 UNFL = SLAMCH( 'Safe minimum' )
121 ULP = SLAMCH( 'Precision' )
122 *
123 ITRNSE = 0
124 ITRNSW = 0
125 NORMA = 'O'
126 NORME = 'O'
127 *
128 IF( LSAME( TRANSA, 'T' ) .OR. LSAME( TRANSA, 'C' ) ) THEN
129 NORMA = 'I'
130 END IF
131 *
132 IF( LSAME( TRANSE, 'T' ) ) THEN
133 ITRNSE = 1
134 NORME = 'I'
135 ELSE IF( LSAME( TRANSE, 'C' ) ) THEN
136 ITRNSE = 2
137 NORME = 'I'
138 END IF
139 *
140 IF( LSAME( TRANSW, 'C' ) ) THEN
141 ITRNSW = 1
142 END IF
143 *
144 * Normalization of E:
145 *
146 ENRMIN = ONE / ULP
147 ENRMAX = ZERO
148 IF( ITRNSE.EQ.0 ) THEN
149 DO 20 JVEC = 1, N
150 TEMP1 = ZERO
151 DO 10 J = 1, N
152 TEMP1 = MAX( TEMP1, ABS( REAL( E( J, JVEC ) ) )+
153 $ ABS( AIMAG( E( J, JVEC ) ) ) )
154 10 CONTINUE
155 ENRMIN = MIN( ENRMIN, TEMP1 )
156 ENRMAX = MAX( ENRMAX, TEMP1 )
157 20 CONTINUE
158 ELSE
159 DO 30 JVEC = 1, N
160 RWORK( JVEC ) = ZERO
161 30 CONTINUE
162 *
163 DO 50 J = 1, N
164 DO 40 JVEC = 1, N
165 RWORK( JVEC ) = MAX( RWORK( JVEC ),
166 $ ABS( REAL( E( JVEC, J ) ) )+
167 $ ABS( AIMAG( E( JVEC, J ) ) ) )
168 40 CONTINUE
169 50 CONTINUE
170 *
171 DO 60 JVEC = 1, N
172 ENRMIN = MIN( ENRMIN, RWORK( JVEC ) )
173 ENRMAX = MAX( ENRMAX, RWORK( JVEC ) )
174 60 CONTINUE
175 END IF
176 *
177 * Norm of A:
178 *
179 ANORM = MAX( CLANGE( NORMA, N, N, A, LDA, RWORK ), UNFL )
180 *
181 * Norm of E:
182 *
183 ENORM = MAX( CLANGE( NORME, N, N, E, LDE, RWORK ), ULP )
184 *
185 * Norm of error:
186 *
187 * Error = AE - EW
188 *
189 CALL CLASET( 'Full', N, N, CZERO, CZERO, WORK, N )
190 *
191 JOFF = 0
192 DO 100 JCOL = 1, N
193 IF( ITRNSW.EQ.0 ) THEN
194 WTEMP = W( JCOL )
195 ELSE
196 WTEMP = CONJG( W( JCOL ) )
197 END IF
198 *
199 IF( ITRNSE.EQ.0 ) THEN
200 DO 70 JROW = 1, N
201 WORK( JOFF+JROW ) = E( JROW, JCOL )*WTEMP
202 70 CONTINUE
203 ELSE IF( ITRNSE.EQ.1 ) THEN
204 DO 80 JROW = 1, N
205 WORK( JOFF+JROW ) = E( JCOL, JROW )*WTEMP
206 80 CONTINUE
207 ELSE
208 DO 90 JROW = 1, N
209 WORK( JOFF+JROW ) = CONJG( E( JCOL, JROW ) )*WTEMP
210 90 CONTINUE
211 END IF
212 JOFF = JOFF + N
213 100 CONTINUE
214 *
215 CALL CGEMM( TRANSA, TRANSE, N, N, N, CONE, A, LDA, E, LDE, -CONE,
216 $ WORK, N )
217 *
218 ERRNRM = CLANGE( 'One', N, N, WORK, N, RWORK ) / ENORM
219 *
220 * Compute RESULT(1) (avoiding under/overflow)
221 *
222 IF( ANORM.GT.ERRNRM ) THEN
223 RESULT( 1 ) = ( ERRNRM / ANORM ) / ULP
224 ELSE
225 IF( ANORM.LT.ONE ) THEN
226 RESULT( 1 ) = ( MIN( ERRNRM, ANORM ) / ANORM ) / ULP
227 ELSE
228 RESULT( 1 ) = MIN( ERRNRM / ANORM, ONE ) / ULP
229 END IF
230 END IF
231 *
232 * Compute RESULT(2) : the normalization error in E.
233 *
234 RESULT( 2 ) = MAX( ABS( ENRMAX-ONE ), ABS( ENRMIN-ONE ) ) /
235 $ ( REAL( N )*ULP )
236 *
237 RETURN
238 *
239 * End of CGET22
240 *
241 END