![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CCHKGG 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 * Definition: 00009 * =========== 00010 * 00011 * SUBROUTINE CCHKGG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 00012 * TSTDIF, THRSHN, NOUNIT, A, LDA, B, H, T, S1, 00013 * S2, P1, P2, U, LDU, V, Q, Z, ALPHA1, BETA1, 00014 * ALPHA3, BETA3, EVECTL, EVECTR, WORK, LWORK, 00015 * RWORK, LLWORK, RESULT, INFO ) 00016 * 00017 * .. Scalar Arguments .. 00018 * LOGICAL TSTDIF 00019 * INTEGER INFO, LDA, LDU, LWORK, NOUNIT, NSIZES, NTYPES 00020 * REAL THRESH, THRSHN 00021 * .. 00022 * .. Array Arguments .. 00023 * LOGICAL DOTYPE( * ), LLWORK( * ) 00024 * INTEGER ISEED( 4 ), NN( * ) 00025 * REAL RESULT( 15 ), RWORK( * ) 00026 * COMPLEX A( LDA, * ), ALPHA1( * ), ALPHA3( * ), 00027 * $ B( LDA, * ), BETA1( * ), BETA3( * ), 00028 * $ EVECTL( LDU, * ), EVECTR( LDU, * ), 00029 * $ H( LDA, * ), P1( LDA, * ), P2( LDA, * ), 00030 * $ Q( LDU, * ), S1( LDA, * ), S2( LDA, * ), 00031 * $ T( LDA, * ), U( LDU, * ), V( LDU, * ), 00032 * $ WORK( * ), Z( LDU, * ) 00033 * .. 00034 * 00035 * 00036 *> \par Purpose: 00037 * ============= 00038 *> 00039 *> \verbatim 00040 *> 00041 *> CCHKGG checks the nonsymmetric generalized eigenvalue problem 00042 *> routines. 00043 *> H H H 00044 *> CGGHRD factors A and B as U H V and U T V , where means conjugate 00045 *> transpose, H is hessenberg, T is triangular and U and V are unitary. 00046 *> 00047 *> H H 00048 *> CHGEQZ factors H and T as Q S Z and Q P Z , where P and S are upper 00049 *> triangular and Q and Z are unitary. It also computes the generalized 00050 *> eigenvalues (alpha(1),beta(1)),...,(alpha(n),beta(n)), where 00051 *> alpha(j)=S(j,j) and beta(j)=P(j,j) -- thus, w(j) = alpha(j)/beta(j) 00052 *> is a root of the generalized eigenvalue problem 00053 *> 00054 *> det( A - w(j) B ) = 0 00055 *> 00056 *> and m(j) = beta(j)/alpha(j) is a root of the essentially equivalent 00057 *> problem 00058 *> 00059 *> det( m(j) A - B ) = 0 00060 *> 00061 *> CTGEVC computes the matrix L of left eigenvectors and the matrix R 00062 *> of right eigenvectors for the matrix pair ( S, P ). In the 00063 *> description below, l and r are left and right eigenvectors 00064 *> corresponding to the generalized eigenvalues (alpha,beta). 00065 *> 00066 *> When CCHKGG is called, a number of matrix "sizes" ("n's") and a 00067 *> number of matrix "types" are specified. For each size ("n") 00068 *> and each type of matrix, one matrix will be generated and used 00069 *> to test the nonsymmetric eigenroutines. For each matrix, 13 00070 *> tests will be performed. The first twelve "test ratios" should be 00071 *> small -- O(1). They will be compared with the threshhold THRESH: 00072 *> 00073 *> H 00074 *> (1) | A - U H V | / ( |A| n ulp ) 00075 *> 00076 *> H 00077 *> (2) | B - U T V | / ( |B| n ulp ) 00078 *> 00079 *> H 00080 *> (3) | I - UU | / ( n ulp ) 00081 *> 00082 *> H 00083 *> (4) | I - VV | / ( n ulp ) 00084 *> 00085 *> H 00086 *> (5) | H - Q S Z | / ( |H| n ulp ) 00087 *> 00088 *> H 00089 *> (6) | T - Q P Z | / ( |T| n ulp ) 00090 *> 00091 *> H 00092 *> (7) | I - QQ | / ( n ulp ) 00093 *> 00094 *> H 00095 *> (8) | I - ZZ | / ( n ulp ) 00096 *> 00097 *> (9) max over all left eigenvalue/-vector pairs (beta/alpha,l) of 00098 *> H 00099 *> | (beta A - alpha B) l | / ( ulp max( |beta A|, |alpha B| ) ) 00100 *> 00101 *> (10) max over all left eigenvalue/-vector pairs (beta/alpha,l') of 00102 *> H 00103 *> | (beta H - alpha T) l' | / ( ulp max( |beta H|, |alpha T| ) ) 00104 *> 00105 *> where the eigenvectors l' are the result of passing Q to 00106 *> STGEVC and back transforming (JOB='B'). 00107 *> 00108 *> (11) max over all right eigenvalue/-vector pairs (beta/alpha,r) of 00109 *> 00110 *> | (beta A - alpha B) r | / ( ulp max( |beta A|, |alpha B| ) ) 00111 *> 00112 *> (12) max over all right eigenvalue/-vector pairs (beta/alpha,r') of 00113 *> 00114 *> | (beta H - alpha T) r' | / ( ulp max( |beta H|, |alpha T| ) ) 00115 *> 00116 *> where the eigenvectors r' are the result of passing Z to 00117 *> STGEVC and back transforming (JOB='B'). 00118 *> 00119 *> The last three test ratios will usually be small, but there is no 00120 *> mathematical requirement that they be so. They are therefore 00121 *> compared with THRESH only if TSTDIF is .TRUE. 00122 *> 00123 *> (13) | S(Q,Z computed) - S(Q,Z not computed) | / ( |S| ulp ) 00124 *> 00125 *> (14) | P(Q,Z computed) - P(Q,Z not computed) | / ( |P| ulp ) 00126 *> 00127 *> (15) max( |alpha(Q,Z computed) - alpha(Q,Z not computed)|/|S| , 00128 *> |beta(Q,Z computed) - beta(Q,Z not computed)|/|P| ) / ulp 00129 *> 00130 *> In addition, the normalization of L and R are checked, and compared 00131 *> with the threshhold THRSHN. 00132 *> 00133 *> Test Matrices 00134 *> ---- -------- 00135 *> 00136 *> The sizes of the test matrices are specified by an array 00137 *> NN(1:NSIZES); the value of each element NN(j) specifies one size. 00138 *> The "types" are specified by a logical array DOTYPE( 1:NTYPES ); if 00139 *> DOTYPE(j) is .TRUE., then matrix type "j" will be generated. 00140 *> Currently, the list of possible types is: 00141 *> 00142 *> (1) ( 0, 0 ) (a pair of zero matrices) 00143 *> 00144 *> (2) ( I, 0 ) (an identity and a zero matrix) 00145 *> 00146 *> (3) ( 0, I ) (an identity and a zero matrix) 00147 *> 00148 *> (4) ( I, I ) (a pair of identity matrices) 00149 *> 00150 *> t t 00151 *> (5) ( J , J ) (a pair of transposed Jordan blocks) 00152 *> 00153 *> t ( I 0 ) 00154 *> (6) ( X, Y ) where X = ( J 0 ) and Y = ( t ) 00155 *> ( 0 I ) ( 0 J ) 00156 *> and I is a k x k identity and J a (k+1)x(k+1) 00157 *> Jordan block; k=(N-1)/2 00158 *> 00159 *> (7) ( D, I ) where D is P*D1, P is a random unitary diagonal 00160 *> matrix (i.e., with random magnitude 1 entries 00161 *> on the diagonal), and D1=diag( 0, 1,..., N-1 ) 00162 *> (i.e., a diagonal matrix with D1(1,1)=0, 00163 *> D1(2,2)=1, ..., D1(N,N)=N-1.) 00164 *> (8) ( I, D ) 00165 *> 00166 *> (9) ( big*D, small*I ) where "big" is near overflow and small=1/big 00167 *> 00168 *> (10) ( small*D, big*I ) 00169 *> 00170 *> (11) ( big*I, small*D ) 00171 *> 00172 *> (12) ( small*I, big*D ) 00173 *> 00174 *> (13) ( big*D, big*I ) 00175 *> 00176 *> (14) ( small*D, small*I ) 00177 *> 00178 *> (15) ( D1, D2 ) where D1=P*diag( 0, 0, 1, ..., N-3, 0 ) and 00179 *> D2=Q*diag( 0, N-3, N-4,..., 1, 0, 0 ), and 00180 *> P and Q are random unitary diagonal matrices. 00181 *> t t 00182 *> (16) U ( J , J ) V where U and V are random unitary matrices. 00183 *> 00184 *> (17) U ( T1, T2 ) V where T1 and T2 are upper triangular matrices 00185 *> with random O(1) entries above the diagonal 00186 *> and diagonal entries diag(T1) = 00187 *> P*( 0, 0, 1, ..., N-3, 0 ) and diag(T2) = 00188 *> Q*( 0, N-3, N-4,..., 1, 0, 0 ) 00189 *> 00190 *> (18) U ( T1, T2 ) V diag(T1) = ( 0, 0, 1, 1, s, ..., s, 0 ) 00191 *> diag(T2) = ( 0, 1, 0, 1,..., 1, 0 ) 00192 *> s = machine precision. 00193 *> 00194 *> (19) U ( T1, T2 ) V diag(T1)=( 0,0,1,1, 1-d, ..., 1-(N-5)*d=s, 0 ) 00195 *> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0 ) 00196 *> 00197 *> N-5 00198 *> (20) U ( T1, T2 ) V diag(T1)=( 0, 0, 1, 1, a, ..., a =s, 0 ) 00199 *> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 ) 00200 *> 00201 *> (21) U ( T1, T2 ) V diag(T1)=( 0, 0, 1, r1, r2, ..., r(N-4), 0 ) 00202 *> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 ) 00203 *> where r1,..., r(N-4) are random. 00204 *> 00205 *> (22) U ( big*T1, small*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 ) 00206 *> diag(T2) = ( 0, 1, ..., 1, 0, 0 ) 00207 *> 00208 *> (23) U ( small*T1, big*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 ) 00209 *> diag(T2) = ( 0, 1, ..., 1, 0, 0 ) 00210 *> 00211 *> (24) U ( small*T1, small*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 ) 00212 *> diag(T2) = ( 0, 1, ..., 1, 0, 0 ) 00213 *> 00214 *> (25) U ( big*T1, big*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 ) 00215 *> diag(T2) = ( 0, 1, ..., 1, 0, 0 ) 00216 *> 00217 *> (26) U ( T1, T2 ) V where T1 and T2 are random upper-triangular 00218 *> matrices. 00219 *> \endverbatim 00220 * 00221 * Arguments: 00222 * ========== 00223 * 00224 *> \param[in] NSIZES 00225 *> \verbatim 00226 *> NSIZES is INTEGER 00227 *> The number of sizes of matrices to use. If it is zero, 00228 *> CCHKGG does nothing. It must be at least zero. 00229 *> \endverbatim 00230 *> 00231 *> \param[in] NN 00232 *> \verbatim 00233 *> NN is INTEGER array, dimension (NSIZES) 00234 *> An array containing the sizes to be used for the matrices. 00235 *> Zero values will be skipped. The values must be at least 00236 *> zero. 00237 *> \endverbatim 00238 *> 00239 *> \param[in] NTYPES 00240 *> \verbatim 00241 *> NTYPES is INTEGER 00242 *> The number of elements in DOTYPE. If it is zero, CCHKGG 00243 *> does nothing. It must be at least zero. If it is MAXTYP+1 00244 *> and NSIZES is 1, then an additional type, MAXTYP+1 is 00245 *> defined, which is to use whatever matrix is in A. This 00246 *> is only useful if DOTYPE(1:MAXTYP) is .FALSE. and 00247 *> DOTYPE(MAXTYP+1) is .TRUE. . 00248 *> \endverbatim 00249 *> 00250 *> \param[in] DOTYPE 00251 *> \verbatim 00252 *> DOTYPE is LOGICAL array, dimension (NTYPES) 00253 *> If DOTYPE(j) is .TRUE., then for each size in NN a 00254 *> matrix of that size and of type j will be generated. 00255 *> If NTYPES is smaller than the maximum number of types 00256 *> defined (PARAMETER MAXTYP), then types NTYPES+1 through 00257 *> MAXTYP will not be generated. If NTYPES is larger 00258 *> than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES) 00259 *> will be ignored. 00260 *> \endverbatim 00261 *> 00262 *> \param[in,out] ISEED 00263 *> \verbatim 00264 *> ISEED is INTEGER array, dimension (4) 00265 *> On entry ISEED specifies the seed of the random number 00266 *> generator. The array elements should be between 0 and 4095; 00267 *> if not they will be reduced mod 4096. Also, ISEED(4) must 00268 *> be odd. The random number generator uses a linear 00269 *> congruential sequence limited to small integers, and so 00270 *> should produce machine independent random numbers. The 00271 *> values of ISEED are changed on exit, and can be used in the 00272 *> next call to CCHKGG to continue the same random number 00273 *> sequence. 00274 *> \endverbatim 00275 *> 00276 *> \param[in] THRESH 00277 *> \verbatim 00278 *> THRESH is REAL 00279 *> A test will count as "failed" if the "error", computed as 00280 *> described above, exceeds THRESH. Note that the error 00281 *> is scaled to be O(1), so THRESH should be a reasonably 00282 *> small multiple of 1, e.g., 10 or 100. In particular, 00283 *> it should not depend on the precision (single vs. double) 00284 *> or the size of the matrix. It must be at least zero. 00285 *> \endverbatim 00286 *> 00287 *> \param[in] TSTDIF 00288 *> \verbatim 00289 *> TSTDIF is LOGICAL 00290 *> Specifies whether test ratios 13-15 will be computed and 00291 *> compared with THRESH. 00292 *> = .FALSE.: Only test ratios 1-12 will be computed and tested. 00293 *> Ratios 13-15 will be set to zero. 00294 *> = .TRUE.: All the test ratios 1-15 will be computed and 00295 *> tested. 00296 *> \endverbatim 00297 *> 00298 *> \param[in] THRSHN 00299 *> \verbatim 00300 *> THRSHN is REAL 00301 *> Threshhold for reporting eigenvector normalization error. 00302 *> If the normalization of any eigenvector differs from 1 by 00303 *> more than THRSHN*ulp, then a special error message will be 00304 *> printed. (This is handled separately from the other tests, 00305 *> since only a compiler or programming error should cause an 00306 *> error message, at least if THRSHN is at least 5--10.) 00307 *> \endverbatim 00308 *> 00309 *> \param[in] NOUNIT 00310 *> \verbatim 00311 *> NOUNIT is INTEGER 00312 *> The FORTRAN unit number for printing out error messages 00313 *> (e.g., if a routine returns IINFO not equal to 0.) 00314 *> \endverbatim 00315 *> 00316 *> \param[in,out] A 00317 *> \verbatim 00318 *> A is COMPLEX array, dimension (LDA, max(NN)) 00319 *> Used to hold the original A matrix. Used as input only 00320 *> if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and 00321 *> DOTYPE(MAXTYP+1)=.TRUE. 00322 *> \endverbatim 00323 *> 00324 *> \param[in] LDA 00325 *> \verbatim 00326 *> LDA is INTEGER 00327 *> The leading dimension of A, B, H, T, S1, P1, S2, and P2. 00328 *> It must be at least 1 and at least max( NN ). 00329 *> \endverbatim 00330 *> 00331 *> \param[in,out] B 00332 *> \verbatim 00333 *> B is COMPLEX array, dimension (LDA, max(NN)) 00334 *> Used to hold the original B matrix. Used as input only 00335 *> if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and 00336 *> DOTYPE(MAXTYP+1)=.TRUE. 00337 *> \endverbatim 00338 *> 00339 *> \param[out] H 00340 *> \verbatim 00341 *> H is COMPLEX array, dimension (LDA, max(NN)) 00342 *> The upper Hessenberg matrix computed from A by CGGHRD. 00343 *> \endverbatim 00344 *> 00345 *> \param[out] T 00346 *> \verbatim 00347 *> T is COMPLEX array, dimension (LDA, max(NN)) 00348 *> The upper triangular matrix computed from B by CGGHRD. 00349 *> \endverbatim 00350 *> 00351 *> \param[out] S1 00352 *> \verbatim 00353 *> S1 is COMPLEX array, dimension (LDA, max(NN)) 00354 *> The Schur (upper triangular) matrix computed from H by CHGEQZ 00355 *> when Q and Z are also computed. 00356 *> \endverbatim 00357 *> 00358 *> \param[out] S2 00359 *> \verbatim 00360 *> S2 is COMPLEX array, dimension (LDA, max(NN)) 00361 *> The Schur (upper triangular) matrix computed from H by CHGEQZ 00362 *> when Q and Z are not computed. 00363 *> \endverbatim 00364 *> 00365 *> \param[out] P1 00366 *> \verbatim 00367 *> P1 is COMPLEX array, dimension (LDA, max(NN)) 00368 *> The upper triangular matrix computed from T by CHGEQZ 00369 *> when Q and Z are also computed. 00370 *> \endverbatim 00371 *> 00372 *> \param[out] P2 00373 *> \verbatim 00374 *> P2 is COMPLEX array, dimension (LDA, max(NN)) 00375 *> The upper triangular matrix computed from T by CHGEQZ 00376 *> when Q and Z are not computed. 00377 *> \endverbatim 00378 *> 00379 *> \param[out] U 00380 *> \verbatim 00381 *> U is COMPLEX array, dimension (LDU, max(NN)) 00382 *> The (left) unitary matrix computed by CGGHRD. 00383 *> \endverbatim 00384 *> 00385 *> \param[in] LDU 00386 *> \verbatim 00387 *> LDU is INTEGER 00388 *> The leading dimension of U, V, Q, Z, EVECTL, and EVECTR. It 00389 *> must be at least 1 and at least max( NN ). 00390 *> \endverbatim 00391 *> 00392 *> \param[out] V 00393 *> \verbatim 00394 *> V is COMPLEX array, dimension (LDU, max(NN)) 00395 *> The (right) unitary matrix computed by CGGHRD. 00396 *> \endverbatim 00397 *> 00398 *> \param[out] Q 00399 *> \verbatim 00400 *> Q is COMPLEX array, dimension (LDU, max(NN)) 00401 *> The (left) unitary matrix computed by CHGEQZ. 00402 *> \endverbatim 00403 *> 00404 *> \param[out] Z 00405 *> \verbatim 00406 *> Z is COMPLEX array, dimension (LDU, max(NN)) 00407 *> The (left) unitary matrix computed by CHGEQZ. 00408 *> \endverbatim 00409 *> 00410 *> \param[out] ALPHA1 00411 *> \verbatim 00412 *> ALPHA1 is COMPLEX array, dimension (max(NN)) 00413 *> \endverbatim 00414 *> 00415 *> \param[out] BETA1 00416 *> \verbatim 00417 *> BETA1 is COMPLEX array, dimension (max(NN)) 00418 *> The generalized eigenvalues of (A,B) computed by CHGEQZ 00419 *> when Q, Z, and the full Schur matrices are computed. 00420 *> \endverbatim 00421 *> 00422 *> \param[out] ALPHA3 00423 *> \verbatim 00424 *> ALPHA3 is COMPLEX array, dimension (max(NN)) 00425 *> \endverbatim 00426 *> 00427 *> \param[out] BETA3 00428 *> \verbatim 00429 *> BETA3 is COMPLEX array, dimension (max(NN)) 00430 *> The generalized eigenvalues of (A,B) computed by CHGEQZ 00431 *> when neither Q, Z, nor the Schur matrices are computed. 00432 *> \endverbatim 00433 *> 00434 *> \param[out] EVECTL 00435 *> \verbatim 00436 *> EVECTL is COMPLEX array, dimension (LDU, max(NN)) 00437 *> The (lower triangular) left eigenvector matrix for the 00438 *> matrices in S1 and P1. 00439 *> \endverbatim 00440 *> 00441 *> \param[out] EVECTR 00442 *> \verbatim 00443 *> EVECTR is COMPLEX array, dimension (LDU, max(NN)) 00444 *> The (upper triangular) right eigenvector matrix for the 00445 *> matrices in S1 and P1. 00446 *> \endverbatim 00447 *> 00448 *> \param[out] WORK 00449 *> \verbatim 00450 *> WORK is COMPLEX array, dimension (LWORK) 00451 *> \endverbatim 00452 *> 00453 *> \param[in] LWORK 00454 *> \verbatim 00455 *> LWORK is INTEGER 00456 *> The number of entries in WORK. This must be at least 00457 *> max( 4*N, 2 * N**2, 1 ), for all N=NN(j). 00458 *> \endverbatim 00459 *> 00460 *> \param[out] RWORK 00461 *> \verbatim 00462 *> RWORK is REAL array, dimension (2*max(NN)) 00463 *> \endverbatim 00464 *> 00465 *> \param[out] LLWORK 00466 *> \verbatim 00467 *> LLWORK is LOGICAL array, dimension (max(NN)) 00468 *> \endverbatim 00469 *> 00470 *> \param[out] RESULT 00471 *> \verbatim 00472 *> RESULT is REAL array, dimension (15) 00473 *> The values computed by the tests described above. 00474 *> The values are currently limited to 1/ulp, to avoid 00475 *> overflow. 00476 *> \endverbatim 00477 *> 00478 *> \param[out] INFO 00479 *> \verbatim 00480 *> INFO is INTEGER 00481 *> = 0: successful exit. 00482 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00483 *> > 0: A routine returned an error code. INFO is the 00484 *> absolute value of the INFO value returned. 00485 *> \endverbatim 00486 * 00487 * Authors: 00488 * ======== 00489 * 00490 *> \author Univ. of Tennessee 00491 *> \author Univ. of California Berkeley 00492 *> \author Univ. of Colorado Denver 00493 *> \author NAG Ltd. 00494 * 00495 *> \date November 2011 00496 * 00497 *> \ingroup complex_eig 00498 * 00499 * ===================================================================== 00500 SUBROUTINE CCHKGG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 00501 $ TSTDIF, THRSHN, NOUNIT, A, LDA, B, H, T, S1, 00502 $ S2, P1, P2, U, LDU, V, Q, Z, ALPHA1, BETA1, 00503 $ ALPHA3, BETA3, EVECTL, EVECTR, WORK, LWORK, 00504 $ RWORK, LLWORK, RESULT, INFO ) 00505 * 00506 * -- LAPACK test routine (version 3.4.0) -- 00507 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00508 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00509 * November 2011 00510 * 00511 * .. Scalar Arguments .. 00512 LOGICAL TSTDIF 00513 INTEGER INFO, LDA, LDU, LWORK, NOUNIT, NSIZES, NTYPES 00514 REAL THRESH, THRSHN 00515 * .. 00516 * .. Array Arguments .. 00517 LOGICAL DOTYPE( * ), LLWORK( * ) 00518 INTEGER ISEED( 4 ), NN( * ) 00519 REAL RESULT( 15 ), RWORK( * ) 00520 COMPLEX A( LDA, * ), ALPHA1( * ), ALPHA3( * ), 00521 $ B( LDA, * ), BETA1( * ), BETA3( * ), 00522 $ EVECTL( LDU, * ), EVECTR( LDU, * ), 00523 $ H( LDA, * ), P1( LDA, * ), P2( LDA, * ), 00524 $ Q( LDU, * ), S1( LDA, * ), S2( LDA, * ), 00525 $ T( LDA, * ), U( LDU, * ), V( LDU, * ), 00526 $ WORK( * ), Z( LDU, * ) 00527 * .. 00528 * 00529 * ===================================================================== 00530 * 00531 * .. Parameters .. 00532 REAL ZERO, ONE 00533 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00534 COMPLEX CZERO, CONE 00535 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00536 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 00537 INTEGER MAXTYP 00538 PARAMETER ( MAXTYP = 26 ) 00539 * .. 00540 * .. Local Scalars .. 00541 LOGICAL BADNN 00542 INTEGER I1, IADD, IINFO, IN, J, JC, JR, JSIZE, JTYPE, 00543 $ LWKOPT, MTYPES, N, N1, NERRS, NMATS, NMAX, 00544 $ NTEST, NTESTT 00545 REAL ANORM, BNORM, SAFMAX, SAFMIN, TEMP1, TEMP2, 00546 $ ULP, ULPINV 00547 COMPLEX CTEMP 00548 * .. 00549 * .. Local Arrays .. 00550 LOGICAL LASIGN( MAXTYP ), LBSIGN( MAXTYP ) 00551 INTEGER IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ), 00552 $ KATYPE( MAXTYP ), KAZERO( MAXTYP ), 00553 $ KBMAGN( MAXTYP ), KBTYPE( MAXTYP ), 00554 $ KBZERO( MAXTYP ), KCLASS( MAXTYP ), 00555 $ KTRIAN( MAXTYP ), KZ1( 6 ), KZ2( 6 ) 00556 REAL DUMMA( 4 ), RMAGN( 0: 3 ) 00557 COMPLEX CDUMMA( 4 ) 00558 * .. 00559 * .. External Functions .. 00560 REAL CLANGE, SLAMCH 00561 COMPLEX CLARND 00562 EXTERNAL CLANGE, SLAMCH, CLARND 00563 * .. 00564 * .. External Subroutines .. 00565 EXTERNAL CGEQR2, CGET51, CGET52, CGGHRD, CHGEQZ, CLACPY, 00566 $ CLARFG, CLASET, CLATM4, CTGEVC, CUNM2R, SLABAD, 00567 $ SLASUM, XERBLA 00568 * .. 00569 * .. Intrinsic Functions .. 00570 INTRINSIC ABS, CONJG, MAX, MIN, REAL, SIGN 00571 * .. 00572 * .. Data statements .. 00573 DATA KCLASS / 15*1, 10*2, 1*3 / 00574 DATA KZ1 / 0, 1, 2, 1, 3, 3 / 00575 DATA KZ2 / 0, 0, 1, 2, 1, 1 / 00576 DATA KADD / 0, 0, 0, 0, 3, 2 / 00577 DATA KATYPE / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4, 00578 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 / 00579 DATA KBTYPE / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4, 00580 $ 1, 1, -4, 2, -4, 8*8, 0 / 00581 DATA KAZERO / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3, 00582 $ 4*5, 4*3, 1 / 00583 DATA KBZERO / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4, 00584 $ 4*6, 4*4, 1 / 00585 DATA KAMAGN / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3, 00586 $ 2, 1 / 00587 DATA KBMAGN / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3, 00588 $ 2, 1 / 00589 DATA KTRIAN / 16*0, 10*1 / 00590 DATA LASIGN / 6*.FALSE., .TRUE., .FALSE., 2*.TRUE., 00591 $ 2*.FALSE., 3*.TRUE., .FALSE., .TRUE., 00592 $ 3*.FALSE., 5*.TRUE., .FALSE. / 00593 DATA LBSIGN / 7*.FALSE., .TRUE., 2*.FALSE., 00594 $ 2*.TRUE., 2*.FALSE., .TRUE., .FALSE., .TRUE., 00595 $ 9*.FALSE. / 00596 * .. 00597 * .. Executable Statements .. 00598 * 00599 * Check for errors 00600 * 00601 INFO = 0 00602 * 00603 BADNN = .FALSE. 00604 NMAX = 1 00605 DO 10 J = 1, NSIZES 00606 NMAX = MAX( NMAX, NN( J ) ) 00607 IF( NN( J ).LT.0 ) 00608 $ BADNN = .TRUE. 00609 10 CONTINUE 00610 * 00611 LWKOPT = MAX( 2*NMAX*NMAX, 4*NMAX, 1 ) 00612 * 00613 * Check for errors 00614 * 00615 IF( NSIZES.LT.0 ) THEN 00616 INFO = -1 00617 ELSE IF( BADNN ) THEN 00618 INFO = -2 00619 ELSE IF( NTYPES.LT.0 ) THEN 00620 INFO = -3 00621 ELSE IF( THRESH.LT.ZERO ) THEN 00622 INFO = -6 00623 ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN 00624 INFO = -10 00625 ELSE IF( LDU.LE.1 .OR. LDU.LT.NMAX ) THEN 00626 INFO = -19 00627 ELSE IF( LWKOPT.GT.LWORK ) THEN 00628 INFO = -30 00629 END IF 00630 * 00631 IF( INFO.NE.0 ) THEN 00632 CALL XERBLA( 'CCHKGG', -INFO ) 00633 RETURN 00634 END IF 00635 * 00636 * Quick return if possible 00637 * 00638 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 ) 00639 $ RETURN 00640 * 00641 SAFMIN = SLAMCH( 'Safe minimum' ) 00642 ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' ) 00643 SAFMIN = SAFMIN / ULP 00644 SAFMAX = ONE / SAFMIN 00645 CALL SLABAD( SAFMIN, SAFMAX ) 00646 ULPINV = ONE / ULP 00647 * 00648 * The values RMAGN(2:3) depend on N, see below. 00649 * 00650 RMAGN( 0 ) = ZERO 00651 RMAGN( 1 ) = ONE 00652 * 00653 * Loop over sizes, types 00654 * 00655 NTESTT = 0 00656 NERRS = 0 00657 NMATS = 0 00658 * 00659 DO 240 JSIZE = 1, NSIZES 00660 N = NN( JSIZE ) 00661 N1 = MAX( 1, N ) 00662 RMAGN( 2 ) = SAFMAX*ULP / REAL( N1 ) 00663 RMAGN( 3 ) = SAFMIN*ULPINV*N1 00664 * 00665 IF( NSIZES.NE.1 ) THEN 00666 MTYPES = MIN( MAXTYP, NTYPES ) 00667 ELSE 00668 MTYPES = MIN( MAXTYP+1, NTYPES ) 00669 END IF 00670 * 00671 DO 230 JTYPE = 1, MTYPES 00672 IF( .NOT.DOTYPE( JTYPE ) ) 00673 $ GO TO 230 00674 NMATS = NMATS + 1 00675 NTEST = 0 00676 * 00677 * Save ISEED in case of an error. 00678 * 00679 DO 20 J = 1, 4 00680 IOLDSD( J ) = ISEED( J ) 00681 20 CONTINUE 00682 * 00683 * Initialize RESULT 00684 * 00685 DO 30 J = 1, 15 00686 RESULT( J ) = ZERO 00687 30 CONTINUE 00688 * 00689 * Compute A and B 00690 * 00691 * Description of control parameters: 00692 * 00693 * KCLASS: =1 means w/o rotation, =2 means w/ rotation, 00694 * =3 means random. 00695 * KATYPE: the "type" to be passed to CLATM4 for computing A. 00696 * KAZERO: the pattern of zeros on the diagonal for A: 00697 * =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ), 00698 * =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ), 00699 * =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of 00700 * non-zero entries.) 00701 * KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1), 00702 * =2: large, =3: small. 00703 * LASIGN: .TRUE. if the diagonal elements of A are to be 00704 * multiplied by a random magnitude 1 number. 00705 * KBTYPE, KBZERO, KBMAGN, LBSIGN: the same, but for B. 00706 * KTRIAN: =0: don't fill in the upper triangle, =1: do. 00707 * KZ1, KZ2, KADD: used to implement KAZERO and KBZERO. 00708 * RMAGN: used to implement KAMAGN and KBMAGN. 00709 * 00710 IF( MTYPES.GT.MAXTYP ) 00711 $ GO TO 110 00712 IINFO = 0 00713 IF( KCLASS( JTYPE ).LT.3 ) THEN 00714 * 00715 * Generate A (w/o rotation) 00716 * 00717 IF( ABS( KATYPE( JTYPE ) ).EQ.3 ) THEN 00718 IN = 2*( ( N-1 ) / 2 ) + 1 00719 IF( IN.NE.N ) 00720 $ CALL CLASET( 'Full', N, N, CZERO, CZERO, A, LDA ) 00721 ELSE 00722 IN = N 00723 END IF 00724 CALL CLATM4( KATYPE( JTYPE ), IN, KZ1( KAZERO( JTYPE ) ), 00725 $ KZ2( KAZERO( JTYPE ) ), LASIGN( JTYPE ), 00726 $ RMAGN( KAMAGN( JTYPE ) ), ULP, 00727 $ RMAGN( KTRIAN( JTYPE )*KAMAGN( JTYPE ) ), 4, 00728 $ ISEED, A, LDA ) 00729 IADD = KADD( KAZERO( JTYPE ) ) 00730 IF( IADD.GT.0 .AND. IADD.LE.N ) 00731 $ A( IADD, IADD ) = RMAGN( KAMAGN( JTYPE ) ) 00732 * 00733 * Generate B (w/o rotation) 00734 * 00735 IF( ABS( KBTYPE( JTYPE ) ).EQ.3 ) THEN 00736 IN = 2*( ( N-1 ) / 2 ) + 1 00737 IF( IN.NE.N ) 00738 $ CALL CLASET( 'Full', N, N, CZERO, CZERO, B, LDA ) 00739 ELSE 00740 IN = N 00741 END IF 00742 CALL CLATM4( KBTYPE( JTYPE ), IN, KZ1( KBZERO( JTYPE ) ), 00743 $ KZ2( KBZERO( JTYPE ) ), LBSIGN( JTYPE ), 00744 $ RMAGN( KBMAGN( JTYPE ) ), ONE, 00745 $ RMAGN( KTRIAN( JTYPE )*KBMAGN( JTYPE ) ), 4, 00746 $ ISEED, B, LDA ) 00747 IADD = KADD( KBZERO( JTYPE ) ) 00748 IF( IADD.NE.0 ) 00749 $ B( IADD, IADD ) = RMAGN( KBMAGN( JTYPE ) ) 00750 * 00751 IF( KCLASS( JTYPE ).EQ.2 .AND. N.GT.0 ) THEN 00752 * 00753 * Include rotations 00754 * 00755 * Generate U, V as Householder transformations times a 00756 * diagonal matrix. (Note that CLARFG makes U(j,j) and 00757 * V(j,j) real.) 00758 * 00759 DO 50 JC = 1, N - 1 00760 DO 40 JR = JC, N 00761 U( JR, JC ) = CLARND( 3, ISEED ) 00762 V( JR, JC ) = CLARND( 3, ISEED ) 00763 40 CONTINUE 00764 CALL CLARFG( N+1-JC, U( JC, JC ), U( JC+1, JC ), 1, 00765 $ WORK( JC ) ) 00766 WORK( 2*N+JC ) = SIGN( ONE, REAL( U( JC, JC ) ) ) 00767 U( JC, JC ) = CONE 00768 CALL CLARFG( N+1-JC, V( JC, JC ), V( JC+1, JC ), 1, 00769 $ WORK( N+JC ) ) 00770 WORK( 3*N+JC ) = SIGN( ONE, REAL( V( JC, JC ) ) ) 00771 V( JC, JC ) = CONE 00772 50 CONTINUE 00773 CTEMP = CLARND( 3, ISEED ) 00774 U( N, N ) = CONE 00775 WORK( N ) = CZERO 00776 WORK( 3*N ) = CTEMP / ABS( CTEMP ) 00777 CTEMP = CLARND( 3, ISEED ) 00778 V( N, N ) = CONE 00779 WORK( 2*N ) = CZERO 00780 WORK( 4*N ) = CTEMP / ABS( CTEMP ) 00781 * 00782 * Apply the diagonal matrices 00783 * 00784 DO 70 JC = 1, N 00785 DO 60 JR = 1, N 00786 A( JR, JC ) = WORK( 2*N+JR )* 00787 $ CONJG( WORK( 3*N+JC ) )* 00788 $ A( JR, JC ) 00789 B( JR, JC ) = WORK( 2*N+JR )* 00790 $ CONJG( WORK( 3*N+JC ) )* 00791 $ B( JR, JC ) 00792 60 CONTINUE 00793 70 CONTINUE 00794 CALL CUNM2R( 'L', 'N', N, N, N-1, U, LDU, WORK, A, 00795 $ LDA, WORK( 2*N+1 ), IINFO ) 00796 IF( IINFO.NE.0 ) 00797 $ GO TO 100 00798 CALL CUNM2R( 'R', 'C', N, N, N-1, V, LDU, WORK( N+1 ), 00799 $ A, LDA, WORK( 2*N+1 ), IINFO ) 00800 IF( IINFO.NE.0 ) 00801 $ GO TO 100 00802 CALL CUNM2R( 'L', 'N', N, N, N-1, U, LDU, WORK, B, 00803 $ LDA, WORK( 2*N+1 ), IINFO ) 00804 IF( IINFO.NE.0 ) 00805 $ GO TO 100 00806 CALL CUNM2R( 'R', 'C', N, N, N-1, V, LDU, WORK( N+1 ), 00807 $ B, LDA, WORK( 2*N+1 ), IINFO ) 00808 IF( IINFO.NE.0 ) 00809 $ GO TO 100 00810 END IF 00811 ELSE 00812 * 00813 * Random matrices 00814 * 00815 DO 90 JC = 1, N 00816 DO 80 JR = 1, N 00817 A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )* 00818 $ CLARND( 4, ISEED ) 00819 B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )* 00820 $ CLARND( 4, ISEED ) 00821 80 CONTINUE 00822 90 CONTINUE 00823 END IF 00824 * 00825 ANORM = CLANGE( '1', N, N, A, LDA, RWORK ) 00826 BNORM = CLANGE( '1', N, N, B, LDA, RWORK ) 00827 * 00828 100 CONTINUE 00829 * 00830 IF( IINFO.NE.0 ) THEN 00831 WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE, 00832 $ IOLDSD 00833 INFO = ABS( IINFO ) 00834 RETURN 00835 END IF 00836 * 00837 110 CONTINUE 00838 * 00839 * Call CGEQR2, CUNM2R, and CGGHRD to compute H, T, U, and V 00840 * 00841 CALL CLACPY( ' ', N, N, A, LDA, H, LDA ) 00842 CALL CLACPY( ' ', N, N, B, LDA, T, LDA ) 00843 NTEST = 1 00844 RESULT( 1 ) = ULPINV 00845 * 00846 CALL CGEQR2( N, N, T, LDA, WORK, WORK( N+1 ), IINFO ) 00847 IF( IINFO.NE.0 ) THEN 00848 WRITE( NOUNIT, FMT = 9999 )'CGEQR2', IINFO, N, JTYPE, 00849 $ IOLDSD 00850 INFO = ABS( IINFO ) 00851 GO TO 210 00852 END IF 00853 * 00854 CALL CUNM2R( 'L', 'C', N, N, N, T, LDA, WORK, H, LDA, 00855 $ WORK( N+1 ), IINFO ) 00856 IF( IINFO.NE.0 ) THEN 00857 WRITE( NOUNIT, FMT = 9999 )'CUNM2R', IINFO, N, JTYPE, 00858 $ IOLDSD 00859 INFO = ABS( IINFO ) 00860 GO TO 210 00861 END IF 00862 * 00863 CALL CLASET( 'Full', N, N, CZERO, CONE, U, LDU ) 00864 CALL CUNM2R( 'R', 'N', N, N, N, T, LDA, WORK, U, LDU, 00865 $ WORK( N+1 ), IINFO ) 00866 IF( IINFO.NE.0 ) THEN 00867 WRITE( NOUNIT, FMT = 9999 )'CUNM2R', IINFO, N, JTYPE, 00868 $ IOLDSD 00869 INFO = ABS( IINFO ) 00870 GO TO 210 00871 END IF 00872 * 00873 CALL CGGHRD( 'V', 'I', N, 1, N, H, LDA, T, LDA, U, LDU, V, 00874 $ LDU, IINFO ) 00875 IF( IINFO.NE.0 ) THEN 00876 WRITE( NOUNIT, FMT = 9999 )'CGGHRD', IINFO, N, JTYPE, 00877 $ IOLDSD 00878 INFO = ABS( IINFO ) 00879 GO TO 210 00880 END IF 00881 NTEST = 4 00882 * 00883 * Do tests 1--4 00884 * 00885 CALL CGET51( 1, N, A, LDA, H, LDA, U, LDU, V, LDU, WORK, 00886 $ RWORK, RESULT( 1 ) ) 00887 CALL CGET51( 1, N, B, LDA, T, LDA, U, LDU, V, LDU, WORK, 00888 $ RWORK, RESULT( 2 ) ) 00889 CALL CGET51( 3, N, B, LDA, T, LDA, U, LDU, U, LDU, WORK, 00890 $ RWORK, RESULT( 3 ) ) 00891 CALL CGET51( 3, N, B, LDA, T, LDA, V, LDU, V, LDU, WORK, 00892 $ RWORK, RESULT( 4 ) ) 00893 * 00894 * Call CHGEQZ to compute S1, P1, S2, P2, Q, and Z, do tests. 00895 * 00896 * Compute T1 and UZ 00897 * 00898 * Eigenvalues only 00899 * 00900 CALL CLACPY( ' ', N, N, H, LDA, S2, LDA ) 00901 CALL CLACPY( ' ', N, N, T, LDA, P2, LDA ) 00902 NTEST = 5 00903 RESULT( 5 ) = ULPINV 00904 * 00905 CALL CHGEQZ( 'E', 'N', 'N', N, 1, N, S2, LDA, P2, LDA, 00906 $ ALPHA3, BETA3, Q, LDU, Z, LDU, WORK, LWORK, 00907 $ RWORK, IINFO ) 00908 IF( IINFO.NE.0 ) THEN 00909 WRITE( NOUNIT, FMT = 9999 )'CHGEQZ(E)', IINFO, N, JTYPE, 00910 $ IOLDSD 00911 INFO = ABS( IINFO ) 00912 GO TO 210 00913 END IF 00914 * 00915 * Eigenvalues and Full Schur Form 00916 * 00917 CALL CLACPY( ' ', N, N, H, LDA, S2, LDA ) 00918 CALL CLACPY( ' ', N, N, T, LDA, P2, LDA ) 00919 * 00920 CALL CHGEQZ( 'S', 'N', 'N', N, 1, N, S2, LDA, P2, LDA, 00921 $ ALPHA1, BETA1, Q, LDU, Z, LDU, WORK, LWORK, 00922 $ RWORK, IINFO ) 00923 IF( IINFO.NE.0 ) THEN 00924 WRITE( NOUNIT, FMT = 9999 )'CHGEQZ(S)', IINFO, N, JTYPE, 00925 $ IOLDSD 00926 INFO = ABS( IINFO ) 00927 GO TO 210 00928 END IF 00929 * 00930 * Eigenvalues, Schur Form, and Schur Vectors 00931 * 00932 CALL CLACPY( ' ', N, N, H, LDA, S1, LDA ) 00933 CALL CLACPY( ' ', N, N, T, LDA, P1, LDA ) 00934 * 00935 CALL CHGEQZ( 'S', 'I', 'I', N, 1, N, S1, LDA, P1, LDA, 00936 $ ALPHA1, BETA1, Q, LDU, Z, LDU, WORK, LWORK, 00937 $ RWORK, IINFO ) 00938 IF( IINFO.NE.0 ) THEN 00939 WRITE( NOUNIT, FMT = 9999 )'CHGEQZ(V)', IINFO, N, JTYPE, 00940 $ IOLDSD 00941 INFO = ABS( IINFO ) 00942 GO TO 210 00943 END IF 00944 * 00945 NTEST = 8 00946 * 00947 * Do Tests 5--8 00948 * 00949 CALL CGET51( 1, N, H, LDA, S1, LDA, Q, LDU, Z, LDU, WORK, 00950 $ RWORK, RESULT( 5 ) ) 00951 CALL CGET51( 1, N, T, LDA, P1, LDA, Q, LDU, Z, LDU, WORK, 00952 $ RWORK, RESULT( 6 ) ) 00953 CALL CGET51( 3, N, T, LDA, P1, LDA, Q, LDU, Q, LDU, WORK, 00954 $ RWORK, RESULT( 7 ) ) 00955 CALL CGET51( 3, N, T, LDA, P1, LDA, Z, LDU, Z, LDU, WORK, 00956 $ RWORK, RESULT( 8 ) ) 00957 * 00958 * Compute the Left and Right Eigenvectors of (S1,P1) 00959 * 00960 * 9: Compute the left eigenvector Matrix without 00961 * back transforming: 00962 * 00963 NTEST = 9 00964 RESULT( 9 ) = ULPINV 00965 * 00966 * To test "SELECT" option, compute half of the eigenvectors 00967 * in one call, and half in another 00968 * 00969 I1 = N / 2 00970 DO 120 J = 1, I1 00971 LLWORK( J ) = .TRUE. 00972 120 CONTINUE 00973 DO 130 J = I1 + 1, N 00974 LLWORK( J ) = .FALSE. 00975 130 CONTINUE 00976 * 00977 CALL CTGEVC( 'L', 'S', LLWORK, N, S1, LDA, P1, LDA, EVECTL, 00978 $ LDU, CDUMMA, LDU, N, IN, WORK, RWORK, IINFO ) 00979 IF( IINFO.NE.0 ) THEN 00980 WRITE( NOUNIT, FMT = 9999 )'CTGEVC(L,S1)', IINFO, N, 00981 $ JTYPE, IOLDSD 00982 INFO = ABS( IINFO ) 00983 GO TO 210 00984 END IF 00985 * 00986 I1 = IN 00987 DO 140 J = 1, I1 00988 LLWORK( J ) = .FALSE. 00989 140 CONTINUE 00990 DO 150 J = I1 + 1, N 00991 LLWORK( J ) = .TRUE. 00992 150 CONTINUE 00993 * 00994 CALL CTGEVC( 'L', 'S', LLWORK, N, S1, LDA, P1, LDA, 00995 $ EVECTL( 1, I1+1 ), LDU, CDUMMA, LDU, N, IN, 00996 $ WORK, RWORK, IINFO ) 00997 IF( IINFO.NE.0 ) THEN 00998 WRITE( NOUNIT, FMT = 9999 )'CTGEVC(L,S2)', IINFO, N, 00999 $ JTYPE, IOLDSD 01000 INFO = ABS( IINFO ) 01001 GO TO 210 01002 END IF 01003 * 01004 CALL CGET52( .TRUE., N, S1, LDA, P1, LDA, EVECTL, LDU, 01005 $ ALPHA1, BETA1, WORK, RWORK, DUMMA( 1 ) ) 01006 RESULT( 9 ) = DUMMA( 1 ) 01007 IF( DUMMA( 2 ).GT.THRSHN ) THEN 01008 WRITE( NOUNIT, FMT = 9998 )'Left', 'CTGEVC(HOWMNY=S)', 01009 $ DUMMA( 2 ), N, JTYPE, IOLDSD 01010 END IF 01011 * 01012 * 10: Compute the left eigenvector Matrix with 01013 * back transforming: 01014 * 01015 NTEST = 10 01016 RESULT( 10 ) = ULPINV 01017 CALL CLACPY( 'F', N, N, Q, LDU, EVECTL, LDU ) 01018 CALL CTGEVC( 'L', 'B', LLWORK, N, S1, LDA, P1, LDA, EVECTL, 01019 $ LDU, CDUMMA, LDU, N, IN, WORK, RWORK, IINFO ) 01020 IF( IINFO.NE.0 ) THEN 01021 WRITE( NOUNIT, FMT = 9999 )'CTGEVC(L,B)', IINFO, N, 01022 $ JTYPE, IOLDSD 01023 INFO = ABS( IINFO ) 01024 GO TO 210 01025 END IF 01026 * 01027 CALL CGET52( .TRUE., N, H, LDA, T, LDA, EVECTL, LDU, ALPHA1, 01028 $ BETA1, WORK, RWORK, DUMMA( 1 ) ) 01029 RESULT( 10 ) = DUMMA( 1 ) 01030 IF( DUMMA( 2 ).GT.THRSHN ) THEN 01031 WRITE( NOUNIT, FMT = 9998 )'Left', 'CTGEVC(HOWMNY=B)', 01032 $ DUMMA( 2 ), N, JTYPE, IOLDSD 01033 END IF 01034 * 01035 * 11: Compute the right eigenvector Matrix without 01036 * back transforming: 01037 * 01038 NTEST = 11 01039 RESULT( 11 ) = ULPINV 01040 * 01041 * To test "SELECT" option, compute half of the eigenvectors 01042 * in one call, and half in another 01043 * 01044 I1 = N / 2 01045 DO 160 J = 1, I1 01046 LLWORK( J ) = .TRUE. 01047 160 CONTINUE 01048 DO 170 J = I1 + 1, N 01049 LLWORK( J ) = .FALSE. 01050 170 CONTINUE 01051 * 01052 CALL CTGEVC( 'R', 'S', LLWORK, N, S1, LDA, P1, LDA, CDUMMA, 01053 $ LDU, EVECTR, LDU, N, IN, WORK, RWORK, IINFO ) 01054 IF( IINFO.NE.0 ) THEN 01055 WRITE( NOUNIT, FMT = 9999 )'CTGEVC(R,S1)', IINFO, N, 01056 $ JTYPE, IOLDSD 01057 INFO = ABS( IINFO ) 01058 GO TO 210 01059 END IF 01060 * 01061 I1 = IN 01062 DO 180 J = 1, I1 01063 LLWORK( J ) = .FALSE. 01064 180 CONTINUE 01065 DO 190 J = I1 + 1, N 01066 LLWORK( J ) = .TRUE. 01067 190 CONTINUE 01068 * 01069 CALL CTGEVC( 'R', 'S', LLWORK, N, S1, LDA, P1, LDA, CDUMMA, 01070 $ LDU, EVECTR( 1, I1+1 ), LDU, N, IN, WORK, 01071 $ RWORK, IINFO ) 01072 IF( IINFO.NE.0 ) THEN 01073 WRITE( NOUNIT, FMT = 9999 )'CTGEVC(R,S2)', IINFO, N, 01074 $ JTYPE, IOLDSD 01075 INFO = ABS( IINFO ) 01076 GO TO 210 01077 END IF 01078 * 01079 CALL CGET52( .FALSE., N, S1, LDA, P1, LDA, EVECTR, LDU, 01080 $ ALPHA1, BETA1, WORK, RWORK, DUMMA( 1 ) ) 01081 RESULT( 11 ) = DUMMA( 1 ) 01082 IF( DUMMA( 2 ).GT.THRESH ) THEN 01083 WRITE( NOUNIT, FMT = 9998 )'Right', 'CTGEVC(HOWMNY=S)', 01084 $ DUMMA( 2 ), N, JTYPE, IOLDSD 01085 END IF 01086 * 01087 * 12: Compute the right eigenvector Matrix with 01088 * back transforming: 01089 * 01090 NTEST = 12 01091 RESULT( 12 ) = ULPINV 01092 CALL CLACPY( 'F', N, N, Z, LDU, EVECTR, LDU ) 01093 CALL CTGEVC( 'R', 'B', LLWORK, N, S1, LDA, P1, LDA, CDUMMA, 01094 $ LDU, EVECTR, LDU, N, IN, WORK, RWORK, IINFO ) 01095 IF( IINFO.NE.0 ) THEN 01096 WRITE( NOUNIT, FMT = 9999 )'CTGEVC(R,B)', IINFO, N, 01097 $ JTYPE, IOLDSD 01098 INFO = ABS( IINFO ) 01099 GO TO 210 01100 END IF 01101 * 01102 CALL CGET52( .FALSE., N, H, LDA, T, LDA, EVECTR, LDU, 01103 $ ALPHA1, BETA1, WORK, RWORK, DUMMA( 1 ) ) 01104 RESULT( 12 ) = DUMMA( 1 ) 01105 IF( DUMMA( 2 ).GT.THRESH ) THEN 01106 WRITE( NOUNIT, FMT = 9998 )'Right', 'CTGEVC(HOWMNY=B)', 01107 $ DUMMA( 2 ), N, JTYPE, IOLDSD 01108 END IF 01109 * 01110 * Tests 13--15 are done only on request 01111 * 01112 IF( TSTDIF ) THEN 01113 * 01114 * Do Tests 13--14 01115 * 01116 CALL CGET51( 2, N, S1, LDA, S2, LDA, Q, LDU, Z, LDU, 01117 $ WORK, RWORK, RESULT( 13 ) ) 01118 CALL CGET51( 2, N, P1, LDA, P2, LDA, Q, LDU, Z, LDU, 01119 $ WORK, RWORK, RESULT( 14 ) ) 01120 * 01121 * Do Test 15 01122 * 01123 TEMP1 = ZERO 01124 TEMP2 = ZERO 01125 DO 200 J = 1, N 01126 TEMP1 = MAX( TEMP1, ABS( ALPHA1( J )-ALPHA3( J ) ) ) 01127 TEMP2 = MAX( TEMP2, ABS( BETA1( J )-BETA3( J ) ) ) 01128 200 CONTINUE 01129 * 01130 TEMP1 = TEMP1 / MAX( SAFMIN, ULP*MAX( TEMP1, ANORM ) ) 01131 TEMP2 = TEMP2 / MAX( SAFMIN, ULP*MAX( TEMP2, BNORM ) ) 01132 RESULT( 15 ) = MAX( TEMP1, TEMP2 ) 01133 NTEST = 15 01134 ELSE 01135 RESULT( 13 ) = ZERO 01136 RESULT( 14 ) = ZERO 01137 RESULT( 15 ) = ZERO 01138 NTEST = 12 01139 END IF 01140 * 01141 * End of Loop -- Check for RESULT(j) > THRESH 01142 * 01143 210 CONTINUE 01144 * 01145 NTESTT = NTESTT + NTEST 01146 * 01147 * Print out tests which fail. 01148 * 01149 DO 220 JR = 1, NTEST 01150 IF( RESULT( JR ).GE.THRESH ) THEN 01151 * 01152 * If this is the first test to fail, 01153 * print a header to the data file. 01154 * 01155 IF( NERRS.EQ.0 ) THEN 01156 WRITE( NOUNIT, FMT = 9997 )'CGG' 01157 * 01158 * Matrix types 01159 * 01160 WRITE( NOUNIT, FMT = 9996 ) 01161 WRITE( NOUNIT, FMT = 9995 ) 01162 WRITE( NOUNIT, FMT = 9994 )'Unitary' 01163 * 01164 * Tests performed 01165 * 01166 WRITE( NOUNIT, FMT = 9993 )'unitary', '*', 01167 $ 'conjugate transpose', ( '*', J = 1, 10 ) 01168 * 01169 END IF 01170 NERRS = NERRS + 1 01171 IF( RESULT( JR ).LT.10000.0 ) THEN 01172 WRITE( NOUNIT, FMT = 9992 )N, JTYPE, IOLDSD, JR, 01173 $ RESULT( JR ) 01174 ELSE 01175 WRITE( NOUNIT, FMT = 9991 )N, JTYPE, IOLDSD, JR, 01176 $ RESULT( JR ) 01177 END IF 01178 END IF 01179 220 CONTINUE 01180 * 01181 230 CONTINUE 01182 240 CONTINUE 01183 * 01184 * Summary 01185 * 01186 CALL SLASUM( 'CGG', NOUNIT, NERRS, NTESTT ) 01187 RETURN 01188 * 01189 9999 FORMAT( ' CCHKGG: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 01190 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 01191 * 01192 9998 FORMAT( ' CCHKGG: ', A, ' Eigenvectors from ', A, ' incorrectly ', 01193 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X, 01194 $ 'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, 01195 $ ')' ) 01196 * 01197 9997 FORMAT( 1X, A3, ' -- Complex Generalized eigenvalue problem' ) 01198 * 01199 9996 FORMAT( ' Matrix types (see CCHKGG for details): ' ) 01200 * 01201 9995 FORMAT( ' Special Matrices:', 23X, 01202 $ '(J''=transposed Jordan block)', 01203 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ', 01204 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ', 01205 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I', 01206 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', / 01207 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ', 01208 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' ) 01209 9994 FORMAT( ' Matrices Rotated by Random ', A, ' Matrices U, V:', 01210 $ / ' 16=Transposed Jordan Blocks 19=geometric ', 01211 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ', 01212 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ', 01213 $ 'alpha, beta=0,1 21=random alpha, beta=0,1', 01214 $ / ' Large & Small Matrices:', / ' 22=(large, small) ', 01215 $ '23=(small,large) 24=(small,small) 25=(large,large)', 01216 $ / ' 26=random O(1) matrices.' ) 01217 * 01218 9993 FORMAT( / ' Tests performed: (H is Hessenberg, S is Schur, B, ', 01219 $ 'T, P are triangular,', / 20X, 'U, V, Q, and Z are ', A, 01220 $ ', l and r are the', / 20X, 01221 $ 'appropriate left and right eigenvectors, resp., a is', 01222 $ / 20X, 'alpha, b is beta, and ', A, ' means ', A, '.)', 01223 $ / ' 1 = | A - U H V', A, 01224 $ ' | / ( |A| n ulp ) 2 = | B - U T V', A, 01225 $ ' | / ( |B| n ulp )', / ' 3 = | I - UU', A, 01226 $ ' | / ( n ulp ) 4 = | I - VV', A, 01227 $ ' | / ( n ulp )', / ' 5 = | H - Q S Z', A, 01228 $ ' | / ( |H| n ulp )', 6X, '6 = | T - Q P Z', A, 01229 $ ' | / ( |T| n ulp )', / ' 7 = | I - QQ', A, 01230 $ ' | / ( n ulp ) 8 = | I - ZZ', A, 01231 $ ' | / ( n ulp )', / ' 9 = max | ( b S - a P )', A, 01232 $ ' l | / const. 10 = max | ( b H - a T )', A, 01233 $ ' l | / const.', / 01234 $ ' 11= max | ( b S - a P ) r | / const. 12 = max | ( b H', 01235 $ ' - a T ) r | / const.', / 1X ) 01236 * 01237 9992 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=', 01238 $ 4( I4, ',' ), ' result ', I2, ' is', 0P, F8.2 ) 01239 9991 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=', 01240 $ 4( I4, ',' ), ' result ', I2, ' is', 1P, E10.3 ) 01241 * 01242 * End of CCHKGG 01243 * 01244 END