1 SUBROUTINE CGET52( LEFT, N, A, LDA, B, LDB, E, LDE, ALPHA, BETA,
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 LOGICAL LEFT
10 INTEGER LDA, LDB, LDE, N
11 * ..
12 * .. Array Arguments ..
13 REAL RESULT( 2 ), RWORK( * )
14 COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ),
15 $ BETA( * ), E( LDE, * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * CGET52 does an eigenvector check for the generalized eigenvalue
22 * problem.
23 *
24 * The basic test for right eigenvectors is:
25 *
26 * | b(i) A E(i) - a(i) B E(i) |
27 * RESULT(1) = max -------------------------------
28 * i n ulp max( |b(i) A|, |a(i) B| )
29 *
30 * using the 1-norm. Here, a(i)/b(i) = w is the i-th generalized
31 * eigenvalue of A - w B, or, equivalently, b(i)/a(i) = m is the i-th
32 * generalized eigenvalue of m A - B.
33 *
34 * H H _ _
35 * For left eigenvectors, A , B , a, and b are used.
36 *
37 * CGET52 also tests the normalization of E. Each eigenvector is
38 * supposed to be normalized so that the maximum "absolute value"
39 * of its elements is 1, where in this case, "absolute value"
40 * of a complex value x is |Re(x)| + |Im(x)| ; let us call this
41 * maximum "absolute value" norm of a vector v M(v).
42 * if a(i)=b(i)=0, then the eigenvector is set to be the jth coordinate
43 * vector. The normalization test is:
44 *
45 * RESULT(2) = max | M(v(i)) - 1 | / ( n ulp )
46 * eigenvectors v(i)
47 *
48 * Arguments
49 * =========
50 *
51 * LEFT (input) LOGICAL
52 * =.TRUE.: The eigenvectors in the columns of E are assumed
53 * to be *left* eigenvectors.
54 * =.FALSE.: The eigenvectors in the columns of E are assumed
55 * to be *right* eigenvectors.
56 *
57 * N (input) INTEGER
58 * The size of the matrices. If it is zero, CGET52 does
59 * nothing. It must be at least zero.
60 *
61 * A (input) COMPLEX array, dimension (LDA, N)
62 * The matrix A.
63 *
64 * LDA (input) INTEGER
65 * The leading dimension of A. It must be at least 1
66 * and at least N.
67 *
68 * B (input) COMPLEX array, dimension (LDB, N)
69 * The matrix B.
70 *
71 * LDB (input) INTEGER
72 * The leading dimension of B. It must be at least 1
73 * and at least N.
74 *
75 * E (input) COMPLEX array, dimension (LDE, N)
76 * The matrix of eigenvectors. It must be O( 1 ).
77 *
78 * LDE (input) INTEGER
79 * The leading dimension of E. It must be at least 1 and at
80 * least N.
81 *
82 * ALPHA (input) COMPLEX array, dimension (N)
83 * The values a(i) as described above, which, along with b(i),
84 * define the generalized eigenvalues.
85 *
86 * BETA (input) COMPLEX array, dimension (N)
87 * The values b(i) as described above, which, along with a(i),
88 * define the generalized eigenvalues.
89 *
90 * WORK (workspace) COMPLEX array, dimension (N**2)
91 *
92 * RWORK (workspace) REAL array, dimension (N)
93 *
94 * RESULT (output) REAL array, dimension (2)
95 * The values computed by the test described above. If A E or
96 * B E is likely to overflow, then RESULT(1:2) is set to
97 * 10 / ulp.
98 *
99 * =====================================================================
100 *
101 * .. Parameters ..
102 REAL ZERO, ONE
103 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
104 COMPLEX CZERO, CONE
105 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ),
106 $ CONE = ( 1.0E+0, 0.0E+0 ) )
107 * ..
108 * .. Local Scalars ..
109 CHARACTER NORMAB, TRANS
110 INTEGER J, JVEC
111 REAL ABMAX, ALFMAX, ANORM, BETMAX, BNORM, ENORM,
112 $ ENRMER, ERRNRM, SAFMAX, SAFMIN, SCALE, TEMP1,
113 $ ULP
114 COMPLEX ACOEFF, ALPHAI, BCOEFF, BETAI, X
115 * ..
116 * .. External Functions ..
117 REAL CLANGE, SLAMCH
118 EXTERNAL CLANGE, SLAMCH
119 * ..
120 * .. External Subroutines ..
121 EXTERNAL CGEMV
122 * ..
123 * .. Intrinsic Functions ..
124 INTRINSIC ABS, AIMAG, CONJG, MAX, REAL
125 * ..
126 * .. Statement Functions ..
127 REAL ABS1
128 * ..
129 * .. Statement Function definitions ..
130 ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) )
131 * ..
132 * .. Executable Statements ..
133 *
134 RESULT( 1 ) = ZERO
135 RESULT( 2 ) = ZERO
136 IF( N.LE.0 )
137 $ RETURN
138 *
139 SAFMIN = SLAMCH( 'Safe minimum' )
140 SAFMAX = ONE / SAFMIN
141 ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
142 *
143 IF( LEFT ) THEN
144 TRANS = 'C'
145 NORMAB = 'I'
146 ELSE
147 TRANS = 'N'
148 NORMAB = 'O'
149 END IF
150 *
151 * Norm of A, B, and E:
152 *
153 ANORM = MAX( CLANGE( NORMAB, N, N, A, LDA, RWORK ), SAFMIN )
154 BNORM = MAX( CLANGE( NORMAB, N, N, B, LDB, RWORK ), SAFMIN )
155 ENORM = MAX( CLANGE( 'O', N, N, E, LDE, RWORK ), ULP )
156 ALFMAX = SAFMAX / MAX( ONE, BNORM )
157 BETMAX = SAFMAX / MAX( ONE, ANORM )
158 *
159 * Compute error matrix.
160 * Column i = ( b(i) A - a(i) B ) E(i) / max( |a(i) B| |b(i) A| )
161 *
162 DO 10 JVEC = 1, N
163 ALPHAI = ALPHA( JVEC )
164 BETAI = BETA( JVEC )
165 ABMAX = MAX( ABS1( ALPHAI ), ABS1( BETAI ) )
166 IF( ABS1( ALPHAI ).GT.ALFMAX .OR. ABS1( BETAI ).GT.BETMAX .OR.
167 $ ABMAX.LT.ONE ) THEN
168 SCALE = ONE / MAX( ABMAX, SAFMIN )
169 ALPHAI = SCALE*ALPHAI
170 BETAI = SCALE*BETAI
171 END IF
172 SCALE = ONE / MAX( ABS1( ALPHAI )*BNORM, ABS1( BETAI )*ANORM,
173 $ SAFMIN )
174 ACOEFF = SCALE*BETAI
175 BCOEFF = SCALE*ALPHAI
176 IF( LEFT ) THEN
177 ACOEFF = CONJG( ACOEFF )
178 BCOEFF = CONJG( BCOEFF )
179 END IF
180 CALL CGEMV( TRANS, N, N, ACOEFF, A, LDA, E( 1, JVEC ), 1,
181 $ CZERO, WORK( N*( JVEC-1 )+1 ), 1 )
182 CALL CGEMV( TRANS, N, N, -BCOEFF, B, LDA, E( 1, JVEC ), 1,
183 $ CONE, WORK( N*( JVEC-1 )+1 ), 1 )
184 10 CONTINUE
185 *
186 ERRNRM = CLANGE( 'One', N, N, WORK, N, RWORK ) / ENORM
187 *
188 * Compute RESULT(1)
189 *
190 RESULT( 1 ) = ERRNRM / ULP
191 *
192 * Normalization of E:
193 *
194 ENRMER = ZERO
195 DO 30 JVEC = 1, N
196 TEMP1 = ZERO
197 DO 20 J = 1, N
198 TEMP1 = MAX( TEMP1, ABS1( E( J, JVEC ) ) )
199 20 CONTINUE
200 ENRMER = MAX( ENRMER, TEMP1-ONE )
201 30 CONTINUE
202 *
203 * Compute RESULT(2) : the normalization error in E.
204 *
205 RESULT( 2 ) = ENRMER / ( REAL( N )*ULP )
206 *
207 RETURN
208 *
209 * End of CGET52
210 *
211 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 LOGICAL LEFT
10 INTEGER LDA, LDB, LDE, N
11 * ..
12 * .. Array Arguments ..
13 REAL RESULT( 2 ), RWORK( * )
14 COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ),
15 $ BETA( * ), E( LDE, * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * CGET52 does an eigenvector check for the generalized eigenvalue
22 * problem.
23 *
24 * The basic test for right eigenvectors is:
25 *
26 * | b(i) A E(i) - a(i) B E(i) |
27 * RESULT(1) = max -------------------------------
28 * i n ulp max( |b(i) A|, |a(i) B| )
29 *
30 * using the 1-norm. Here, a(i)/b(i) = w is the i-th generalized
31 * eigenvalue of A - w B, or, equivalently, b(i)/a(i) = m is the i-th
32 * generalized eigenvalue of m A - B.
33 *
34 * H H _ _
35 * For left eigenvectors, A , B , a, and b are used.
36 *
37 * CGET52 also tests the normalization of E. Each eigenvector is
38 * supposed to be normalized so that the maximum "absolute value"
39 * of its elements is 1, where in this case, "absolute value"
40 * of a complex value x is |Re(x)| + |Im(x)| ; let us call this
41 * maximum "absolute value" norm of a vector v M(v).
42 * if a(i)=b(i)=0, then the eigenvector is set to be the jth coordinate
43 * vector. The normalization test is:
44 *
45 * RESULT(2) = max | M(v(i)) - 1 | / ( n ulp )
46 * eigenvectors v(i)
47 *
48 * Arguments
49 * =========
50 *
51 * LEFT (input) LOGICAL
52 * =.TRUE.: The eigenvectors in the columns of E are assumed
53 * to be *left* eigenvectors.
54 * =.FALSE.: The eigenvectors in the columns of E are assumed
55 * to be *right* eigenvectors.
56 *
57 * N (input) INTEGER
58 * The size of the matrices. If it is zero, CGET52 does
59 * nothing. It must be at least zero.
60 *
61 * A (input) COMPLEX array, dimension (LDA, N)
62 * The matrix A.
63 *
64 * LDA (input) INTEGER
65 * The leading dimension of A. It must be at least 1
66 * and at least N.
67 *
68 * B (input) COMPLEX array, dimension (LDB, N)
69 * The matrix B.
70 *
71 * LDB (input) INTEGER
72 * The leading dimension of B. It must be at least 1
73 * and at least N.
74 *
75 * E (input) COMPLEX array, dimension (LDE, N)
76 * The matrix of eigenvectors. It must be O( 1 ).
77 *
78 * LDE (input) INTEGER
79 * The leading dimension of E. It must be at least 1 and at
80 * least N.
81 *
82 * ALPHA (input) COMPLEX array, dimension (N)
83 * The values a(i) as described above, which, along with b(i),
84 * define the generalized eigenvalues.
85 *
86 * BETA (input) COMPLEX array, dimension (N)
87 * The values b(i) as described above, which, along with a(i),
88 * define the generalized eigenvalues.
89 *
90 * WORK (workspace) COMPLEX array, dimension (N**2)
91 *
92 * RWORK (workspace) REAL array, dimension (N)
93 *
94 * RESULT (output) REAL array, dimension (2)
95 * The values computed by the test described above. If A E or
96 * B E is likely to overflow, then RESULT(1:2) is set to
97 * 10 / ulp.
98 *
99 * =====================================================================
100 *
101 * .. Parameters ..
102 REAL ZERO, ONE
103 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
104 COMPLEX CZERO, CONE
105 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ),
106 $ CONE = ( 1.0E+0, 0.0E+0 ) )
107 * ..
108 * .. Local Scalars ..
109 CHARACTER NORMAB, TRANS
110 INTEGER J, JVEC
111 REAL ABMAX, ALFMAX, ANORM, BETMAX, BNORM, ENORM,
112 $ ENRMER, ERRNRM, SAFMAX, SAFMIN, SCALE, TEMP1,
113 $ ULP
114 COMPLEX ACOEFF, ALPHAI, BCOEFF, BETAI, X
115 * ..
116 * .. External Functions ..
117 REAL CLANGE, SLAMCH
118 EXTERNAL CLANGE, SLAMCH
119 * ..
120 * .. External Subroutines ..
121 EXTERNAL CGEMV
122 * ..
123 * .. Intrinsic Functions ..
124 INTRINSIC ABS, AIMAG, CONJG, MAX, REAL
125 * ..
126 * .. Statement Functions ..
127 REAL ABS1
128 * ..
129 * .. Statement Function definitions ..
130 ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) )
131 * ..
132 * .. Executable Statements ..
133 *
134 RESULT( 1 ) = ZERO
135 RESULT( 2 ) = ZERO
136 IF( N.LE.0 )
137 $ RETURN
138 *
139 SAFMIN = SLAMCH( 'Safe minimum' )
140 SAFMAX = ONE / SAFMIN
141 ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
142 *
143 IF( LEFT ) THEN
144 TRANS = 'C'
145 NORMAB = 'I'
146 ELSE
147 TRANS = 'N'
148 NORMAB = 'O'
149 END IF
150 *
151 * Norm of A, B, and E:
152 *
153 ANORM = MAX( CLANGE( NORMAB, N, N, A, LDA, RWORK ), SAFMIN )
154 BNORM = MAX( CLANGE( NORMAB, N, N, B, LDB, RWORK ), SAFMIN )
155 ENORM = MAX( CLANGE( 'O', N, N, E, LDE, RWORK ), ULP )
156 ALFMAX = SAFMAX / MAX( ONE, BNORM )
157 BETMAX = SAFMAX / MAX( ONE, ANORM )
158 *
159 * Compute error matrix.
160 * Column i = ( b(i) A - a(i) B ) E(i) / max( |a(i) B| |b(i) A| )
161 *
162 DO 10 JVEC = 1, N
163 ALPHAI = ALPHA( JVEC )
164 BETAI = BETA( JVEC )
165 ABMAX = MAX( ABS1( ALPHAI ), ABS1( BETAI ) )
166 IF( ABS1( ALPHAI ).GT.ALFMAX .OR. ABS1( BETAI ).GT.BETMAX .OR.
167 $ ABMAX.LT.ONE ) THEN
168 SCALE = ONE / MAX( ABMAX, SAFMIN )
169 ALPHAI = SCALE*ALPHAI
170 BETAI = SCALE*BETAI
171 END IF
172 SCALE = ONE / MAX( ABS1( ALPHAI )*BNORM, ABS1( BETAI )*ANORM,
173 $ SAFMIN )
174 ACOEFF = SCALE*BETAI
175 BCOEFF = SCALE*ALPHAI
176 IF( LEFT ) THEN
177 ACOEFF = CONJG( ACOEFF )
178 BCOEFF = CONJG( BCOEFF )
179 END IF
180 CALL CGEMV( TRANS, N, N, ACOEFF, A, LDA, E( 1, JVEC ), 1,
181 $ CZERO, WORK( N*( JVEC-1 )+1 ), 1 )
182 CALL CGEMV( TRANS, N, N, -BCOEFF, B, LDA, E( 1, JVEC ), 1,
183 $ CONE, WORK( N*( JVEC-1 )+1 ), 1 )
184 10 CONTINUE
185 *
186 ERRNRM = CLANGE( 'One', N, N, WORK, N, RWORK ) / ENORM
187 *
188 * Compute RESULT(1)
189 *
190 RESULT( 1 ) = ERRNRM / ULP
191 *
192 * Normalization of E:
193 *
194 ENRMER = ZERO
195 DO 30 JVEC = 1, N
196 TEMP1 = ZERO
197 DO 20 J = 1, N
198 TEMP1 = MAX( TEMP1, ABS1( E( J, JVEC ) ) )
199 20 CONTINUE
200 ENRMER = MAX( ENRMER, TEMP1-ONE )
201 30 CONTINUE
202 *
203 * Compute RESULT(2) : the normalization error in E.
204 *
205 RESULT( 2 ) = ENRMER / ( REAL( N )*ULP )
206 *
207 RETURN
208 *
209 * End of CGET52
210 *
211 END