1 SUBROUTINE ZGET52( 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 DOUBLE PRECISION RESULT( 2 ), RWORK( * )
14 COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
15 $ BETA( * ), E( LDE, * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZGET52 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 * ZGET52 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 *
49 * Arguments
50 * =========
51 *
52 * LEFT (input) LOGICAL
53 * =.TRUE.: The eigenvectors in the columns of E are assumed
54 * to be *left* eigenvectors.
55 * =.FALSE.: The eigenvectors in the columns of E are assumed
56 * to be *right* eigenvectors.
57 *
58 * N (input) INTEGER
59 * The size of the matrices. If it is zero, ZGET52 does
60 * nothing. It must be at least zero.
61 *
62 * A (input) COMPLEX*16 array, dimension (LDA, N)
63 * The matrix A.
64 *
65 * LDA (input) INTEGER
66 * The leading dimension of A. It must be at least 1
67 * and at least N.
68 *
69 * B (input) COMPLEX*16 array, dimension (LDB, N)
70 * The matrix B.
71 *
72 * LDB (input) INTEGER
73 * The leading dimension of B. It must be at least 1
74 * and at least N.
75 *
76 * E (input) COMPLEX*16 array, dimension (LDE, N)
77 * The matrix of eigenvectors. It must be O( 1 ).
78 *
79 * LDE (input) INTEGER
80 * The leading dimension of E. It must be at least 1 and at
81 * least N.
82 *
83 * ALPHA (input) COMPLEX*16 array, dimension (N)
84 * The values a(i) as described above, which, along with b(i),
85 * define the generalized eigenvalues.
86 *
87 * BETA (input) COMPLEX*16 array, dimension (N)
88 * The values b(i) as described above, which, along with a(i),
89 * define the generalized eigenvalues.
90 *
91 * WORK (workspace) COMPLEX*16 array, dimension (N**2)
92 *
93 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
94 *
95 * RESULT (output) DOUBLE PRECISION array, dimension (2)
96 * The values computed by the test described above. If A E or
97 * B E is likely to overflow, then RESULT(1:2) is set to
98 * 10 / ulp.
99 *
100 * =====================================================================
101 *
102 * .. Parameters ..
103 DOUBLE PRECISION ZERO, ONE
104 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
105 COMPLEX*16 CZERO, CONE
106 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
107 $ CONE = ( 1.0D+0, 0.0D+0 ) )
108 * ..
109 * .. Local Scalars ..
110 CHARACTER NORMAB, TRANS
111 INTEGER J, JVEC
112 DOUBLE PRECISION ABMAX, ALFMAX, ANORM, BETMAX, BNORM, ENORM,
113 $ ENRMER, ERRNRM, SAFMAX, SAFMIN, SCALE, TEMP1,
114 $ ULP
115 COMPLEX*16 ACOEFF, ALPHAI, BCOEFF, BETAI, X
116 * ..
117 * .. External Functions ..
118 DOUBLE PRECISION DLAMCH, ZLANGE
119 EXTERNAL DLAMCH, ZLANGE
120 * ..
121 * .. External Subroutines ..
122 EXTERNAL ZGEMV
123 * ..
124 * .. Intrinsic Functions ..
125 INTRINSIC ABS, DBLE, DCONJG, DIMAG, MAX
126 * ..
127 * .. Statement Functions ..
128 DOUBLE PRECISION ABS1
129 * ..
130 * .. Statement Function definitions ..
131 ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
132 * ..
133 * .. Executable Statements ..
134 *
135 RESULT( 1 ) = ZERO
136 RESULT( 2 ) = ZERO
137 IF( N.LE.0 )
138 $ RETURN
139 *
140 SAFMIN = DLAMCH( 'Safe minimum' )
141 SAFMAX = ONE / SAFMIN
142 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
143 *
144 IF( LEFT ) THEN
145 TRANS = 'C'
146 NORMAB = 'I'
147 ELSE
148 TRANS = 'N'
149 NORMAB = 'O'
150 END IF
151 *
152 * Norm of A, B, and E:
153 *
154 ANORM = MAX( ZLANGE( NORMAB, N, N, A, LDA, RWORK ), SAFMIN )
155 BNORM = MAX( ZLANGE( NORMAB, N, N, B, LDB, RWORK ), SAFMIN )
156 ENORM = MAX( ZLANGE( 'O', N, N, E, LDE, RWORK ), ULP )
157 ALFMAX = SAFMAX / MAX( ONE, BNORM )
158 BETMAX = SAFMAX / MAX( ONE, ANORM )
159 *
160 * Compute error matrix.
161 * Column i = ( b(i) A - a(i) B ) E(i) / max( |a(i) B| |b(i) A| )
162 *
163 DO 10 JVEC = 1, N
164 ALPHAI = ALPHA( JVEC )
165 BETAI = BETA( JVEC )
166 ABMAX = MAX( ABS1( ALPHAI ), ABS1( BETAI ) )
167 IF( ABS1( ALPHAI ).GT.ALFMAX .OR. ABS1( BETAI ).GT.BETMAX .OR.
168 $ ABMAX.LT.ONE ) THEN
169 SCALE = ONE / MAX( ABMAX, SAFMIN )
170 ALPHAI = SCALE*ALPHAI
171 BETAI = SCALE*BETAI
172 END IF
173 SCALE = ONE / MAX( ABS1( ALPHAI )*BNORM, ABS1( BETAI )*ANORM,
174 $ SAFMIN )
175 ACOEFF = SCALE*BETAI
176 BCOEFF = SCALE*ALPHAI
177 IF( LEFT ) THEN
178 ACOEFF = DCONJG( ACOEFF )
179 BCOEFF = DCONJG( BCOEFF )
180 END IF
181 CALL ZGEMV( TRANS, N, N, ACOEFF, A, LDA, E( 1, JVEC ), 1,
182 $ CZERO, WORK( N*( JVEC-1 )+1 ), 1 )
183 CALL ZGEMV( TRANS, N, N, -BCOEFF, B, LDA, E( 1, JVEC ), 1,
184 $ CONE, WORK( N*( JVEC-1 )+1 ), 1 )
185 10 CONTINUE
186 *
187 ERRNRM = ZLANGE( 'One', N, N, WORK, N, RWORK ) / ENORM
188 *
189 * Compute RESULT(1)
190 *
191 RESULT( 1 ) = ERRNRM / ULP
192 *
193 * Normalization of E:
194 *
195 ENRMER = ZERO
196 DO 30 JVEC = 1, N
197 TEMP1 = ZERO
198 DO 20 J = 1, N
199 TEMP1 = MAX( TEMP1, ABS1( E( J, JVEC ) ) )
200 20 CONTINUE
201 ENRMER = MAX( ENRMER, TEMP1-ONE )
202 30 CONTINUE
203 *
204 * Compute RESULT(2) : the normalization error in E.
205 *
206 RESULT( 2 ) = ENRMER / ( DBLE( N )*ULP )
207 *
208 RETURN
209 *
210 * End of ZGET52
211 *
212 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 DOUBLE PRECISION RESULT( 2 ), RWORK( * )
14 COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
15 $ BETA( * ), E( LDE, * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZGET52 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 * ZGET52 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 *
49 * Arguments
50 * =========
51 *
52 * LEFT (input) LOGICAL
53 * =.TRUE.: The eigenvectors in the columns of E are assumed
54 * to be *left* eigenvectors.
55 * =.FALSE.: The eigenvectors in the columns of E are assumed
56 * to be *right* eigenvectors.
57 *
58 * N (input) INTEGER
59 * The size of the matrices. If it is zero, ZGET52 does
60 * nothing. It must be at least zero.
61 *
62 * A (input) COMPLEX*16 array, dimension (LDA, N)
63 * The matrix A.
64 *
65 * LDA (input) INTEGER
66 * The leading dimension of A. It must be at least 1
67 * and at least N.
68 *
69 * B (input) COMPLEX*16 array, dimension (LDB, N)
70 * The matrix B.
71 *
72 * LDB (input) INTEGER
73 * The leading dimension of B. It must be at least 1
74 * and at least N.
75 *
76 * E (input) COMPLEX*16 array, dimension (LDE, N)
77 * The matrix of eigenvectors. It must be O( 1 ).
78 *
79 * LDE (input) INTEGER
80 * The leading dimension of E. It must be at least 1 and at
81 * least N.
82 *
83 * ALPHA (input) COMPLEX*16 array, dimension (N)
84 * The values a(i) as described above, which, along with b(i),
85 * define the generalized eigenvalues.
86 *
87 * BETA (input) COMPLEX*16 array, dimension (N)
88 * The values b(i) as described above, which, along with a(i),
89 * define the generalized eigenvalues.
90 *
91 * WORK (workspace) COMPLEX*16 array, dimension (N**2)
92 *
93 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
94 *
95 * RESULT (output) DOUBLE PRECISION array, dimension (2)
96 * The values computed by the test described above. If A E or
97 * B E is likely to overflow, then RESULT(1:2) is set to
98 * 10 / ulp.
99 *
100 * =====================================================================
101 *
102 * .. Parameters ..
103 DOUBLE PRECISION ZERO, ONE
104 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
105 COMPLEX*16 CZERO, CONE
106 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
107 $ CONE = ( 1.0D+0, 0.0D+0 ) )
108 * ..
109 * .. Local Scalars ..
110 CHARACTER NORMAB, TRANS
111 INTEGER J, JVEC
112 DOUBLE PRECISION ABMAX, ALFMAX, ANORM, BETMAX, BNORM, ENORM,
113 $ ENRMER, ERRNRM, SAFMAX, SAFMIN, SCALE, TEMP1,
114 $ ULP
115 COMPLEX*16 ACOEFF, ALPHAI, BCOEFF, BETAI, X
116 * ..
117 * .. External Functions ..
118 DOUBLE PRECISION DLAMCH, ZLANGE
119 EXTERNAL DLAMCH, ZLANGE
120 * ..
121 * .. External Subroutines ..
122 EXTERNAL ZGEMV
123 * ..
124 * .. Intrinsic Functions ..
125 INTRINSIC ABS, DBLE, DCONJG, DIMAG, MAX
126 * ..
127 * .. Statement Functions ..
128 DOUBLE PRECISION ABS1
129 * ..
130 * .. Statement Function definitions ..
131 ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
132 * ..
133 * .. Executable Statements ..
134 *
135 RESULT( 1 ) = ZERO
136 RESULT( 2 ) = ZERO
137 IF( N.LE.0 )
138 $ RETURN
139 *
140 SAFMIN = DLAMCH( 'Safe minimum' )
141 SAFMAX = ONE / SAFMIN
142 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
143 *
144 IF( LEFT ) THEN
145 TRANS = 'C'
146 NORMAB = 'I'
147 ELSE
148 TRANS = 'N'
149 NORMAB = 'O'
150 END IF
151 *
152 * Norm of A, B, and E:
153 *
154 ANORM = MAX( ZLANGE( NORMAB, N, N, A, LDA, RWORK ), SAFMIN )
155 BNORM = MAX( ZLANGE( NORMAB, N, N, B, LDB, RWORK ), SAFMIN )
156 ENORM = MAX( ZLANGE( 'O', N, N, E, LDE, RWORK ), ULP )
157 ALFMAX = SAFMAX / MAX( ONE, BNORM )
158 BETMAX = SAFMAX / MAX( ONE, ANORM )
159 *
160 * Compute error matrix.
161 * Column i = ( b(i) A - a(i) B ) E(i) / max( |a(i) B| |b(i) A| )
162 *
163 DO 10 JVEC = 1, N
164 ALPHAI = ALPHA( JVEC )
165 BETAI = BETA( JVEC )
166 ABMAX = MAX( ABS1( ALPHAI ), ABS1( BETAI ) )
167 IF( ABS1( ALPHAI ).GT.ALFMAX .OR. ABS1( BETAI ).GT.BETMAX .OR.
168 $ ABMAX.LT.ONE ) THEN
169 SCALE = ONE / MAX( ABMAX, SAFMIN )
170 ALPHAI = SCALE*ALPHAI
171 BETAI = SCALE*BETAI
172 END IF
173 SCALE = ONE / MAX( ABS1( ALPHAI )*BNORM, ABS1( BETAI )*ANORM,
174 $ SAFMIN )
175 ACOEFF = SCALE*BETAI
176 BCOEFF = SCALE*ALPHAI
177 IF( LEFT ) THEN
178 ACOEFF = DCONJG( ACOEFF )
179 BCOEFF = DCONJG( BCOEFF )
180 END IF
181 CALL ZGEMV( TRANS, N, N, ACOEFF, A, LDA, E( 1, JVEC ), 1,
182 $ CZERO, WORK( N*( JVEC-1 )+1 ), 1 )
183 CALL ZGEMV( TRANS, N, N, -BCOEFF, B, LDA, E( 1, JVEC ), 1,
184 $ CONE, WORK( N*( JVEC-1 )+1 ), 1 )
185 10 CONTINUE
186 *
187 ERRNRM = ZLANGE( 'One', N, N, WORK, N, RWORK ) / ENORM
188 *
189 * Compute RESULT(1)
190 *
191 RESULT( 1 ) = ERRNRM / ULP
192 *
193 * Normalization of E:
194 *
195 ENRMER = ZERO
196 DO 30 JVEC = 1, N
197 TEMP1 = ZERO
198 DO 20 J = 1, N
199 TEMP1 = MAX( TEMP1, ABS1( E( J, JVEC ) ) )
200 20 CONTINUE
201 ENRMER = MAX( ENRMER, TEMP1-ONE )
202 30 CONTINUE
203 *
204 * Compute RESULT(2) : the normalization error in E.
205 *
206 RESULT( 2 ) = ENRMER / ( DBLE( N )*ULP )
207 *
208 RETURN
209 *
210 * End of ZGET52
211 *
212 END