![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZDRGEV 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 ZDRGEV( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 00012 * NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, QE, LDQE, 00013 * ALPHA, BETA, ALPHA1, BETA1, WORK, LWORK, RWORK, 00014 * RESULT, INFO ) 00015 * 00016 * .. Scalar Arguments .. 00017 * INTEGER INFO, LDA, LDQ, LDQE, LWORK, NOUNIT, NSIZES, 00018 * $ NTYPES 00019 * DOUBLE PRECISION THRESH 00020 * .. 00021 * .. Array Arguments .. 00022 * LOGICAL DOTYPE( * ) 00023 * INTEGER ISEED( 4 ), NN( * ) 00024 * DOUBLE PRECISION RESULT( * ), RWORK( * ) 00025 * COMPLEX*16 A( LDA, * ), ALPHA( * ), ALPHA1( * ), 00026 * $ B( LDA, * ), BETA( * ), BETA1( * ), 00027 * $ Q( LDQ, * ), QE( LDQE, * ), S( LDA, * ), 00028 * $ T( LDA, * ), WORK( * ), Z( LDQ, * ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> ZDRGEV checks the nonsymmetric generalized eigenvalue problem driver 00038 *> routine ZGGEV. 00039 *> 00040 *> ZGGEV computes for a pair of n-by-n nonsymmetric matrices (A,B) the 00041 *> generalized eigenvalues and, optionally, the left and right 00042 *> eigenvectors. 00043 *> 00044 *> A generalized eigenvalue for a pair of matrices (A,B) is a scalar w 00045 *> or a ratio alpha/beta = w, such that A - w*B is singular. It is 00046 *> usually represented as the pair (alpha,beta), as there is reasonalbe 00047 *> interpretation for beta=0, and even for both being zero. 00048 *> 00049 *> A right generalized eigenvector corresponding to a generalized 00050 *> eigenvalue w for a pair of matrices (A,B) is a vector r such that 00051 *> (A - wB) * r = 0. A left generalized eigenvector is a vector l such 00052 *> that l**H * (A - wB) = 0, where l**H is the conjugate-transpose of l. 00053 *> 00054 *> When ZDRGEV is called, a number of matrix "sizes" ("n's") and a 00055 *> number of matrix "types" are specified. For each size ("n") 00056 *> and each type of matrix, a pair of matrices (A, B) will be generated 00057 *> and used for testing. For each matrix pair, the following tests 00058 *> will be performed and compared with the threshhold THRESH. 00059 *> 00060 *> Results from ZGGEV: 00061 *> 00062 *> (1) max over all left eigenvalue/-vector pairs (alpha/beta,l) of 00063 *> 00064 *> | VL**H * (beta A - alpha B) |/( ulp max(|beta A|, |alpha B|) ) 00065 *> 00066 *> where VL**H is the conjugate-transpose of VL. 00067 *> 00068 *> (2) | |VL(i)| - 1 | / ulp and whether largest component real 00069 *> 00070 *> VL(i) denotes the i-th column of VL. 00071 *> 00072 *> (3) max over all left eigenvalue/-vector pairs (alpha/beta,r) of 00073 *> 00074 *> | (beta A - alpha B) * VR | / ( ulp max(|beta A|, |alpha B|) ) 00075 *> 00076 *> (4) | |VR(i)| - 1 | / ulp and whether largest component real 00077 *> 00078 *> VR(i) denotes the i-th column of VR. 00079 *> 00080 *> (5) W(full) = W(partial) 00081 *> W(full) denotes the eigenvalues computed when both l and r 00082 *> are also computed, and W(partial) denotes the eigenvalues 00083 *> computed when only W, only W and r, or only W and l are 00084 *> computed. 00085 *> 00086 *> (6) VL(full) = VL(partial) 00087 *> VL(full) denotes the left eigenvectors computed when both l 00088 *> and r are computed, and VL(partial) denotes the result 00089 *> when only l is computed. 00090 *> 00091 *> (7) VR(full) = VR(partial) 00092 *> VR(full) denotes the right eigenvectors computed when both l 00093 *> and r are also computed, and VR(partial) denotes the result 00094 *> when only l is computed. 00095 *> 00096 *> 00097 *> Test Matrices 00098 *> ---- -------- 00099 *> 00100 *> The sizes of the test matrices are specified by an array 00101 *> NN(1:NSIZES); the value of each element NN(j) specifies one size. 00102 *> The "types" are specified by a logical array DOTYPE( 1:NTYPES ); if 00103 *> DOTYPE(j) is .TRUE., then matrix type "j" will be generated. 00104 *> Currently, the list of possible types is: 00105 *> 00106 *> (1) ( 0, 0 ) (a pair of zero matrices) 00107 *> 00108 *> (2) ( I, 0 ) (an identity and a zero matrix) 00109 *> 00110 *> (3) ( 0, I ) (an identity and a zero matrix) 00111 *> 00112 *> (4) ( I, I ) (a pair of identity matrices) 00113 *> 00114 *> t t 00115 *> (5) ( J , J ) (a pair of transposed Jordan blocks) 00116 *> 00117 *> t ( I 0 ) 00118 *> (6) ( X, Y ) where X = ( J 0 ) and Y = ( t ) 00119 *> ( 0 I ) ( 0 J ) 00120 *> and I is a k x k identity and J a (k+1)x(k+1) 00121 *> Jordan block; k=(N-1)/2 00122 *> 00123 *> (7) ( D, I ) where D is diag( 0, 1,..., N-1 ) (a diagonal 00124 *> matrix with those diagonal entries.) 00125 *> (8) ( I, D ) 00126 *> 00127 *> (9) ( big*D, small*I ) where "big" is near overflow and small=1/big 00128 *> 00129 *> (10) ( small*D, big*I ) 00130 *> 00131 *> (11) ( big*I, small*D ) 00132 *> 00133 *> (12) ( small*I, big*D ) 00134 *> 00135 *> (13) ( big*D, big*I ) 00136 *> 00137 *> (14) ( small*D, small*I ) 00138 *> 00139 *> (15) ( D1, D2 ) where D1 is diag( 0, 0, 1, ..., N-3, 0 ) and 00140 *> D2 is diag( 0, N-3, N-4,..., 1, 0, 0 ) 00141 *> t t 00142 *> (16) Q ( J , J ) Z where Q and Z are random orthogonal matrices. 00143 *> 00144 *> (17) Q ( T1, T2 ) Z where T1 and T2 are upper triangular matrices 00145 *> with random O(1) entries above the diagonal 00146 *> and diagonal entries diag(T1) = 00147 *> ( 0, 0, 1, ..., N-3, 0 ) and diag(T2) = 00148 *> ( 0, N-3, N-4,..., 1, 0, 0 ) 00149 *> 00150 *> (18) Q ( T1, T2 ) Z diag(T1) = ( 0, 0, 1, 1, s, ..., s, 0 ) 00151 *> diag(T2) = ( 0, 1, 0, 1,..., 1, 0 ) 00152 *> s = machine precision. 00153 *> 00154 *> (19) Q ( T1, T2 ) Z diag(T1)=( 0,0,1,1, 1-d, ..., 1-(N-5)*d=s, 0 ) 00155 *> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0 ) 00156 *> 00157 *> N-5 00158 *> (20) Q ( T1, T2 ) Z diag(T1)=( 0, 0, 1, 1, a, ..., a =s, 0 ) 00159 *> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 ) 00160 *> 00161 *> (21) Q ( T1, T2 ) Z diag(T1)=( 0, 0, 1, r1, r2, ..., r(N-4), 0 ) 00162 *> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 ) 00163 *> where r1,..., r(N-4) are random. 00164 *> 00165 *> (22) Q ( big*T1, small*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 ) 00166 *> diag(T2) = ( 0, 1, ..., 1, 0, 0 ) 00167 *> 00168 *> (23) Q ( small*T1, big*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 ) 00169 *> diag(T2) = ( 0, 1, ..., 1, 0, 0 ) 00170 *> 00171 *> (24) Q ( small*T1, small*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 ) 00172 *> diag(T2) = ( 0, 1, ..., 1, 0, 0 ) 00173 *> 00174 *> (25) Q ( big*T1, big*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 ) 00175 *> diag(T2) = ( 0, 1, ..., 1, 0, 0 ) 00176 *> 00177 *> (26) Q ( T1, T2 ) Z where T1 and T2 are random upper-triangular 00178 *> matrices. 00179 *> 00180 *> \endverbatim 00181 * 00182 * Arguments: 00183 * ========== 00184 * 00185 *> \param[in] NSIZES 00186 *> \verbatim 00187 *> NSIZES is INTEGER 00188 *> The number of sizes of matrices to use. If it is zero, 00189 *> ZDRGES does nothing. NSIZES >= 0. 00190 *> \endverbatim 00191 *> 00192 *> \param[in] NN 00193 *> \verbatim 00194 *> NN is INTEGER array, dimension (NSIZES) 00195 *> An array containing the sizes to be used for the matrices. 00196 *> Zero values will be skipped. NN >= 0. 00197 *> \endverbatim 00198 *> 00199 *> \param[in] NTYPES 00200 *> \verbatim 00201 *> NTYPES is INTEGER 00202 *> The number of elements in DOTYPE. If it is zero, ZDRGEV 00203 *> does nothing. It must be at least zero. If it is MAXTYP+1 00204 *> and NSIZES is 1, then an additional type, MAXTYP+1 is 00205 *> defined, which is to use whatever matrix is in A. This 00206 *> is only useful if DOTYPE(1:MAXTYP) is .FALSE. and 00207 *> DOTYPE(MAXTYP+1) is .TRUE. . 00208 *> \endverbatim 00209 *> 00210 *> \param[in] DOTYPE 00211 *> \verbatim 00212 *> DOTYPE is LOGICAL array, dimension (NTYPES) 00213 *> If DOTYPE(j) is .TRUE., then for each size in NN a 00214 *> matrix of that size and of type j will be generated. 00215 *> If NTYPES is smaller than the maximum number of types 00216 *> defined (PARAMETER MAXTYP), then types NTYPES+1 through 00217 *> MAXTYP will not be generated. If NTYPES is larger 00218 *> than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES) 00219 *> will be ignored. 00220 *> \endverbatim 00221 *> 00222 *> \param[in,out] ISEED 00223 *> \verbatim 00224 *> ISEED is INTEGER array, dimension (4) 00225 *> On entry ISEED specifies the seed of the random number 00226 *> generator. The array elements should be between 0 and 4095; 00227 *> if not they will be reduced mod 4096. Also, ISEED(4) must 00228 *> be odd. The random number generator uses a linear 00229 *> congruential sequence limited to small integers, and so 00230 *> should produce machine independent random numbers. The 00231 *> values of ISEED are changed on exit, and can be used in the 00232 *> next call to ZDRGES to continue the same random number 00233 *> sequence. 00234 *> \endverbatim 00235 *> 00236 *> \param[in] THRESH 00237 *> \verbatim 00238 *> THRESH is DOUBLE PRECISION 00239 *> A test will count as "failed" if the "error", computed as 00240 *> described above, exceeds THRESH. Note that the error is 00241 *> scaled to be O(1), so THRESH should be a reasonably small 00242 *> multiple of 1, e.g., 10 or 100. In particular, it should 00243 *> not depend on the precision (single vs. double) or the size 00244 *> of the matrix. It must be at least zero. 00245 *> \endverbatim 00246 *> 00247 *> \param[in] NOUNIT 00248 *> \verbatim 00249 *> NOUNIT is INTEGER 00250 *> The FORTRAN unit number for printing out error messages 00251 *> (e.g., if a routine returns IERR not equal to 0.) 00252 *> \endverbatim 00253 *> 00254 *> \param[in,out] A 00255 *> \verbatim 00256 *> A is COMPLEX*16 array, dimension(LDA, max(NN)) 00257 *> Used to hold the original A matrix. Used as input only 00258 *> if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and 00259 *> DOTYPE(MAXTYP+1)=.TRUE. 00260 *> \endverbatim 00261 *> 00262 *> \param[in] LDA 00263 *> \verbatim 00264 *> LDA is INTEGER 00265 *> The leading dimension of A, B, S, and T. 00266 *> It must be at least 1 and at least max( NN ). 00267 *> \endverbatim 00268 *> 00269 *> \param[in,out] B 00270 *> \verbatim 00271 *> B is COMPLEX*16 array, dimension(LDA, max(NN)) 00272 *> Used to hold the original B matrix. Used as input only 00273 *> if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and 00274 *> DOTYPE(MAXTYP+1)=.TRUE. 00275 *> \endverbatim 00276 *> 00277 *> \param[out] S 00278 *> \verbatim 00279 *> S is COMPLEX*16 array, dimension (LDA, max(NN)) 00280 *> The Schur form matrix computed from A by ZGGEV. On exit, S 00281 *> contains the Schur form matrix corresponding to the matrix 00282 *> in A. 00283 *> \endverbatim 00284 *> 00285 *> \param[out] T 00286 *> \verbatim 00287 *> T is COMPLEX*16 array, dimension (LDA, max(NN)) 00288 *> The upper triangular matrix computed from B by ZGGEV. 00289 *> \endverbatim 00290 *> 00291 *> \param[out] Q 00292 *> \verbatim 00293 *> Q is COMPLEX*16 array, dimension (LDQ, max(NN)) 00294 *> The (left) eigenvectors matrix computed by ZGGEV. 00295 *> \endverbatim 00296 *> 00297 *> \param[in] LDQ 00298 *> \verbatim 00299 *> LDQ is INTEGER 00300 *> The leading dimension of Q and Z. It must 00301 *> be at least 1 and at least max( NN ). 00302 *> \endverbatim 00303 *> 00304 *> \param[out] Z 00305 *> \verbatim 00306 *> Z is COMPLEX*16 array, dimension( LDQ, max(NN) ) 00307 *> The (right) orthogonal matrix computed by ZGGEV. 00308 *> \endverbatim 00309 *> 00310 *> \param[out] QE 00311 *> \verbatim 00312 *> QE is COMPLEX*16 array, dimension( LDQ, max(NN) ) 00313 *> QE holds the computed right or left eigenvectors. 00314 *> \endverbatim 00315 *> 00316 *> \param[in] LDQE 00317 *> \verbatim 00318 *> LDQE is INTEGER 00319 *> The leading dimension of QE. LDQE >= max(1,max(NN)). 00320 *> \endverbatim 00321 *> 00322 *> \param[out] ALPHA 00323 *> \verbatim 00324 *> ALPHA is COMPLEX*16 array, dimension (max(NN)) 00325 *> \endverbatim 00326 *> 00327 *> \param[out] BETA 00328 *> \verbatim 00329 *> BETA is COMPLEX*16 array, dimension (max(NN)) 00330 *> 00331 *> The generalized eigenvalues of (A,B) computed by ZGGEV. 00332 *> ( ALPHAR(k)+ALPHAI(k)*i ) / BETA(k) is the k-th 00333 *> generalized eigenvalue of A and B. 00334 *> \endverbatim 00335 *> 00336 *> \param[out] ALPHA1 00337 *> \verbatim 00338 *> ALPHA1 is COMPLEX*16 array, dimension (max(NN)) 00339 *> \endverbatim 00340 *> 00341 *> \param[out] BETA1 00342 *> \verbatim 00343 *> BETA1 is COMPLEX*16 array, dimension (max(NN)) 00344 *> 00345 *> Like ALPHAR, ALPHAI, BETA, these arrays contain the 00346 *> eigenvalues of A and B, but those computed when ZGGEV only 00347 *> computes a partial eigendecomposition, i.e. not the 00348 *> eigenvalues and left and right eigenvectors. 00349 *> \endverbatim 00350 *> 00351 *> \param[out] WORK 00352 *> \verbatim 00353 *> WORK is COMPLEX*16 array, dimension (LWORK) 00354 *> \endverbatim 00355 *> 00356 *> \param[in] LWORK 00357 *> \verbatim 00358 *> LWORK is INTEGER 00359 *> The number of entries in WORK. LWORK >= N*(N+1) 00360 *> \endverbatim 00361 *> 00362 *> \param[out] RWORK 00363 *> \verbatim 00364 *> RWORK is DOUBLE PRECISION array, dimension (8*N) 00365 *> Real workspace. 00366 *> \endverbatim 00367 *> 00368 *> \param[out] RESULT 00369 *> \verbatim 00370 *> RESULT is DOUBLE PRECISION array, dimension (2) 00371 *> The values computed by the tests described above. 00372 *> The values are currently limited to 1/ulp, to avoid overflow. 00373 *> \endverbatim 00374 *> 00375 *> \param[out] INFO 00376 *> \verbatim 00377 *> INFO is INTEGER 00378 *> = 0: successful exit 00379 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00380 *> > 0: A routine returned an error code. INFO is the 00381 *> absolute value of the INFO value returned. 00382 *> \endverbatim 00383 * 00384 * Authors: 00385 * ======== 00386 * 00387 *> \author Univ. of Tennessee 00388 *> \author Univ. of California Berkeley 00389 *> \author Univ. of Colorado Denver 00390 *> \author NAG Ltd. 00391 * 00392 *> \date November 2011 00393 * 00394 *> \ingroup complex16_eig 00395 * 00396 * ===================================================================== 00397 SUBROUTINE ZDRGEV( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 00398 $ NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, QE, LDQE, 00399 $ ALPHA, BETA, ALPHA1, BETA1, WORK, LWORK, RWORK, 00400 $ RESULT, INFO ) 00401 * 00402 * -- LAPACK test routine (version 3.4.0) -- 00403 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00404 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00405 * November 2011 00406 * 00407 * .. Scalar Arguments .. 00408 INTEGER INFO, LDA, LDQ, LDQE, LWORK, NOUNIT, NSIZES, 00409 $ NTYPES 00410 DOUBLE PRECISION THRESH 00411 * .. 00412 * .. Array Arguments .. 00413 LOGICAL DOTYPE( * ) 00414 INTEGER ISEED( 4 ), NN( * ) 00415 DOUBLE PRECISION RESULT( * ), RWORK( * ) 00416 COMPLEX*16 A( LDA, * ), ALPHA( * ), ALPHA1( * ), 00417 $ B( LDA, * ), BETA( * ), BETA1( * ), 00418 $ Q( LDQ, * ), QE( LDQE, * ), S( LDA, * ), 00419 $ T( LDA, * ), WORK( * ), Z( LDQ, * ) 00420 * .. 00421 * 00422 * ===================================================================== 00423 * 00424 * .. Parameters .. 00425 DOUBLE PRECISION ZERO, ONE 00426 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00427 COMPLEX*16 CZERO, CONE 00428 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 00429 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 00430 INTEGER MAXTYP 00431 PARAMETER ( MAXTYP = 26 ) 00432 * .. 00433 * .. Local Scalars .. 00434 LOGICAL BADNN 00435 INTEGER I, IADD, IERR, IN, J, JC, JR, JSIZE, JTYPE, 00436 $ MAXWRK, MINWRK, MTYPES, N, N1, NB, NERRS, 00437 $ NMATS, NMAX, NTESTT 00438 DOUBLE PRECISION SAFMAX, SAFMIN, ULP, ULPINV 00439 COMPLEX*16 CTEMP 00440 * .. 00441 * .. Local Arrays .. 00442 LOGICAL LASIGN( MAXTYP ), LBSIGN( MAXTYP ) 00443 INTEGER IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ), 00444 $ KATYPE( MAXTYP ), KAZERO( MAXTYP ), 00445 $ KBMAGN( MAXTYP ), KBTYPE( MAXTYP ), 00446 $ KBZERO( MAXTYP ), KCLASS( MAXTYP ), 00447 $ KTRIAN( MAXTYP ), KZ1( 6 ), KZ2( 6 ) 00448 DOUBLE PRECISION RMAGN( 0: 3 ) 00449 * .. 00450 * .. External Functions .. 00451 INTEGER ILAENV 00452 DOUBLE PRECISION DLAMCH 00453 COMPLEX*16 ZLARND 00454 EXTERNAL ILAENV, DLAMCH, ZLARND 00455 * .. 00456 * .. External Subroutines .. 00457 EXTERNAL ALASVM, DLABAD, XERBLA, ZGET52, ZGGEV, ZLACPY, 00458 $ ZLARFG, ZLASET, ZLATM4, ZUNM2R 00459 * .. 00460 * .. Intrinsic Functions .. 00461 INTRINSIC ABS, DBLE, DCONJG, MAX, MIN, SIGN 00462 * .. 00463 * .. Data statements .. 00464 DATA KCLASS / 15*1, 10*2, 1*3 / 00465 DATA KZ1 / 0, 1, 2, 1, 3, 3 / 00466 DATA KZ2 / 0, 0, 1, 2, 1, 1 / 00467 DATA KADD / 0, 0, 0, 0, 3, 2 / 00468 DATA KATYPE / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4, 00469 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 / 00470 DATA KBTYPE / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4, 00471 $ 1, 1, -4, 2, -4, 8*8, 0 / 00472 DATA KAZERO / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3, 00473 $ 4*5, 4*3, 1 / 00474 DATA KBZERO / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4, 00475 $ 4*6, 4*4, 1 / 00476 DATA KAMAGN / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3, 00477 $ 2, 1 / 00478 DATA KBMAGN / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3, 00479 $ 2, 1 / 00480 DATA KTRIAN / 16*0, 10*1 / 00481 DATA LASIGN / 6*.FALSE., .TRUE., .FALSE., 2*.TRUE., 00482 $ 2*.FALSE., 3*.TRUE., .FALSE., .TRUE., 00483 $ 3*.FALSE., 5*.TRUE., .FALSE. / 00484 DATA LBSIGN / 7*.FALSE., .TRUE., 2*.FALSE., 00485 $ 2*.TRUE., 2*.FALSE., .TRUE., .FALSE., .TRUE., 00486 $ 9*.FALSE. / 00487 * .. 00488 * .. Executable Statements .. 00489 * 00490 * Check for errors 00491 * 00492 INFO = 0 00493 * 00494 BADNN = .FALSE. 00495 NMAX = 1 00496 DO 10 J = 1, NSIZES 00497 NMAX = MAX( NMAX, NN( J ) ) 00498 IF( NN( J ).LT.0 ) 00499 $ BADNN = .TRUE. 00500 10 CONTINUE 00501 * 00502 IF( NSIZES.LT.0 ) THEN 00503 INFO = -1 00504 ELSE IF( BADNN ) THEN 00505 INFO = -2 00506 ELSE IF( NTYPES.LT.0 ) THEN 00507 INFO = -3 00508 ELSE IF( THRESH.LT.ZERO ) THEN 00509 INFO = -6 00510 ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN 00511 INFO = -9 00512 ELSE IF( LDQ.LE.1 .OR. LDQ.LT.NMAX ) THEN 00513 INFO = -14 00514 ELSE IF( LDQE.LE.1 .OR. LDQE.LT.NMAX ) THEN 00515 INFO = -17 00516 END IF 00517 * 00518 * Compute workspace 00519 * (Note: Comments in the code beginning "Workspace:" describe the 00520 * minimal amount of workspace needed at that point in the code, 00521 * as well as the preferred amount for good performance. 00522 * NB refers to the optimal block size for the immediately 00523 * following subroutine, as returned by ILAENV. 00524 * 00525 MINWRK = 1 00526 IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN 00527 MINWRK = NMAX*( NMAX+1 ) 00528 NB = MAX( 1, ILAENV( 1, 'ZGEQRF', ' ', NMAX, NMAX, -1, -1 ), 00529 $ ILAENV( 1, 'ZUNMQR', 'LC', NMAX, NMAX, NMAX, -1 ), 00530 $ ILAENV( 1, 'ZUNGQR', ' ', NMAX, NMAX, NMAX, -1 ) ) 00531 MAXWRK = MAX( 2*NMAX, NMAX*( NB+1 ), NMAX*( NMAX+1 ) ) 00532 WORK( 1 ) = MAXWRK 00533 END IF 00534 * 00535 IF( LWORK.LT.MINWRK ) 00536 $ INFO = -23 00537 * 00538 IF( INFO.NE.0 ) THEN 00539 CALL XERBLA( 'ZDRGEV', -INFO ) 00540 RETURN 00541 END IF 00542 * 00543 * Quick return if possible 00544 * 00545 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 ) 00546 $ RETURN 00547 * 00548 ULP = DLAMCH( 'Precision' ) 00549 SAFMIN = DLAMCH( 'Safe minimum' ) 00550 SAFMIN = SAFMIN / ULP 00551 SAFMAX = ONE / SAFMIN 00552 CALL DLABAD( SAFMIN, SAFMAX ) 00553 ULPINV = ONE / ULP 00554 * 00555 * The values RMAGN(2:3) depend on N, see below. 00556 * 00557 RMAGN( 0 ) = ZERO 00558 RMAGN( 1 ) = ONE 00559 * 00560 * Loop over sizes, types 00561 * 00562 NTESTT = 0 00563 NERRS = 0 00564 NMATS = 0 00565 * 00566 DO 220 JSIZE = 1, NSIZES 00567 N = NN( JSIZE ) 00568 N1 = MAX( 1, N ) 00569 RMAGN( 2 ) = SAFMAX*ULP / DBLE( N1 ) 00570 RMAGN( 3 ) = SAFMIN*ULPINV*N1 00571 * 00572 IF( NSIZES.NE.1 ) THEN 00573 MTYPES = MIN( MAXTYP, NTYPES ) 00574 ELSE 00575 MTYPES = MIN( MAXTYP+1, NTYPES ) 00576 END IF 00577 * 00578 DO 210 JTYPE = 1, MTYPES 00579 IF( .NOT.DOTYPE( JTYPE ) ) 00580 $ GO TO 210 00581 NMATS = NMATS + 1 00582 * 00583 * Save ISEED in case of an error. 00584 * 00585 DO 20 J = 1, 4 00586 IOLDSD( J ) = ISEED( J ) 00587 20 CONTINUE 00588 * 00589 * Generate test matrices A and B 00590 * 00591 * Description of control parameters: 00592 * 00593 * KZLASS: =1 means w/o rotation, =2 means w/ rotation, 00594 * =3 means random. 00595 * KATYPE: the "type" to be passed to ZLATM4 for computing A. 00596 * KAZERO: the pattern of zeros on the diagonal for A: 00597 * =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ), 00598 * =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ), 00599 * =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of 00600 * non-zero entries.) 00601 * KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1), 00602 * =2: large, =3: small. 00603 * LASIGN: .TRUE. if the diagonal elements of A are to be 00604 * multiplied by a random magnitude 1 number. 00605 * KBTYPE, KBZERO, KBMAGN, LBSIGN: the same, but for B. 00606 * KTRIAN: =0: don't fill in the upper triangle, =1: do. 00607 * KZ1, KZ2, KADD: used to implement KAZERO and KBZERO. 00608 * RMAGN: used to implement KAMAGN and KBMAGN. 00609 * 00610 IF( MTYPES.GT.MAXTYP ) 00611 $ GO TO 100 00612 IERR = 0 00613 IF( KCLASS( JTYPE ).LT.3 ) THEN 00614 * 00615 * Generate A (w/o rotation) 00616 * 00617 IF( ABS( KATYPE( JTYPE ) ).EQ.3 ) THEN 00618 IN = 2*( ( N-1 ) / 2 ) + 1 00619 IF( IN.NE.N ) 00620 $ CALL ZLASET( 'Full', N, N, CZERO, CZERO, A, LDA ) 00621 ELSE 00622 IN = N 00623 END IF 00624 CALL ZLATM4( KATYPE( JTYPE ), IN, KZ1( KAZERO( JTYPE ) ), 00625 $ KZ2( KAZERO( JTYPE ) ), LASIGN( JTYPE ), 00626 $ RMAGN( KAMAGN( JTYPE ) ), ULP, 00627 $ RMAGN( KTRIAN( JTYPE )*KAMAGN( JTYPE ) ), 2, 00628 $ ISEED, A, LDA ) 00629 IADD = KADD( KAZERO( JTYPE ) ) 00630 IF( IADD.GT.0 .AND. IADD.LE.N ) 00631 $ A( IADD, IADD ) = RMAGN( KAMAGN( JTYPE ) ) 00632 * 00633 * Generate B (w/o rotation) 00634 * 00635 IF( ABS( KBTYPE( JTYPE ) ).EQ.3 ) THEN 00636 IN = 2*( ( N-1 ) / 2 ) + 1 00637 IF( IN.NE.N ) 00638 $ CALL ZLASET( 'Full', N, N, CZERO, CZERO, B, LDA ) 00639 ELSE 00640 IN = N 00641 END IF 00642 CALL ZLATM4( KBTYPE( JTYPE ), IN, KZ1( KBZERO( JTYPE ) ), 00643 $ KZ2( KBZERO( JTYPE ) ), LBSIGN( JTYPE ), 00644 $ RMAGN( KBMAGN( JTYPE ) ), ONE, 00645 $ RMAGN( KTRIAN( JTYPE )*KBMAGN( JTYPE ) ), 2, 00646 $ ISEED, B, LDA ) 00647 IADD = KADD( KBZERO( JTYPE ) ) 00648 IF( IADD.NE.0 .AND. IADD.LE.N ) 00649 $ B( IADD, IADD ) = RMAGN( KBMAGN( JTYPE ) ) 00650 * 00651 IF( KCLASS( JTYPE ).EQ.2 .AND. N.GT.0 ) THEN 00652 * 00653 * Include rotations 00654 * 00655 * Generate Q, Z as Householder transformations times 00656 * a diagonal matrix. 00657 * 00658 DO 40 JC = 1, N - 1 00659 DO 30 JR = JC, N 00660 Q( JR, JC ) = ZLARND( 3, ISEED ) 00661 Z( JR, JC ) = ZLARND( 3, ISEED ) 00662 30 CONTINUE 00663 CALL ZLARFG( N+1-JC, Q( JC, JC ), Q( JC+1, JC ), 1, 00664 $ WORK( JC ) ) 00665 WORK( 2*N+JC ) = SIGN( ONE, DBLE( Q( JC, JC ) ) ) 00666 Q( JC, JC ) = CONE 00667 CALL ZLARFG( N+1-JC, Z( JC, JC ), Z( JC+1, JC ), 1, 00668 $ WORK( N+JC ) ) 00669 WORK( 3*N+JC ) = SIGN( ONE, DBLE( Z( JC, JC ) ) ) 00670 Z( JC, JC ) = CONE 00671 40 CONTINUE 00672 CTEMP = ZLARND( 3, ISEED ) 00673 Q( N, N ) = CONE 00674 WORK( N ) = CZERO 00675 WORK( 3*N ) = CTEMP / ABS( CTEMP ) 00676 CTEMP = ZLARND( 3, ISEED ) 00677 Z( N, N ) = CONE 00678 WORK( 2*N ) = CZERO 00679 WORK( 4*N ) = CTEMP / ABS( CTEMP ) 00680 * 00681 * Apply the diagonal matrices 00682 * 00683 DO 60 JC = 1, N 00684 DO 50 JR = 1, N 00685 A( JR, JC ) = WORK( 2*N+JR )* 00686 $ DCONJG( WORK( 3*N+JC ) )* 00687 $ A( JR, JC ) 00688 B( JR, JC ) = WORK( 2*N+JR )* 00689 $ DCONJG( WORK( 3*N+JC ) )* 00690 $ B( JR, JC ) 00691 50 CONTINUE 00692 60 CONTINUE 00693 CALL ZUNM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, A, 00694 $ LDA, WORK( 2*N+1 ), IERR ) 00695 IF( IERR.NE.0 ) 00696 $ GO TO 90 00697 CALL ZUNM2R( 'R', 'C', N, N, N-1, Z, LDQ, WORK( N+1 ), 00698 $ A, LDA, WORK( 2*N+1 ), IERR ) 00699 IF( IERR.NE.0 ) 00700 $ GO TO 90 00701 CALL ZUNM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, B, 00702 $ LDA, WORK( 2*N+1 ), IERR ) 00703 IF( IERR.NE.0 ) 00704 $ GO TO 90 00705 CALL ZUNM2R( 'R', 'C', N, N, N-1, Z, LDQ, WORK( N+1 ), 00706 $ B, LDA, WORK( 2*N+1 ), IERR ) 00707 IF( IERR.NE.0 ) 00708 $ GO TO 90 00709 END IF 00710 ELSE 00711 * 00712 * Random matrices 00713 * 00714 DO 80 JC = 1, N 00715 DO 70 JR = 1, N 00716 A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )* 00717 $ ZLARND( 4, ISEED ) 00718 B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )* 00719 $ ZLARND( 4, ISEED ) 00720 70 CONTINUE 00721 80 CONTINUE 00722 END IF 00723 * 00724 90 CONTINUE 00725 * 00726 IF( IERR.NE.0 ) THEN 00727 WRITE( NOUNIT, FMT = 9999 )'Generator', IERR, N, JTYPE, 00728 $ IOLDSD 00729 INFO = ABS( IERR ) 00730 RETURN 00731 END IF 00732 * 00733 100 CONTINUE 00734 * 00735 DO 110 I = 1, 7 00736 RESULT( I ) = -ONE 00737 110 CONTINUE 00738 * 00739 * Call ZGGEV to compute eigenvalues and eigenvectors. 00740 * 00741 CALL ZLACPY( ' ', N, N, A, LDA, S, LDA ) 00742 CALL ZLACPY( ' ', N, N, B, LDA, T, LDA ) 00743 CALL ZGGEV( 'V', 'V', N, S, LDA, T, LDA, ALPHA, BETA, Q, 00744 $ LDQ, Z, LDQ, WORK, LWORK, RWORK, IERR ) 00745 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN 00746 RESULT( 1 ) = ULPINV 00747 WRITE( NOUNIT, FMT = 9999 )'ZGGEV1', IERR, N, JTYPE, 00748 $ IOLDSD 00749 INFO = ABS( IERR ) 00750 GO TO 190 00751 END IF 00752 * 00753 * Do the tests (1) and (2) 00754 * 00755 CALL ZGET52( .TRUE., N, A, LDA, B, LDA, Q, LDQ, ALPHA, BETA, 00756 $ WORK, RWORK, RESULT( 1 ) ) 00757 IF( RESULT( 2 ).GT.THRESH ) THEN 00758 WRITE( NOUNIT, FMT = 9998 )'Left', 'ZGGEV1', 00759 $ RESULT( 2 ), N, JTYPE, IOLDSD 00760 END IF 00761 * 00762 * Do the tests (3) and (4) 00763 * 00764 CALL ZGET52( .FALSE., N, A, LDA, B, LDA, Z, LDQ, ALPHA, 00765 $ BETA, WORK, RWORK, RESULT( 3 ) ) 00766 IF( RESULT( 4 ).GT.THRESH ) THEN 00767 WRITE( NOUNIT, FMT = 9998 )'Right', 'ZGGEV1', 00768 $ RESULT( 4 ), N, JTYPE, IOLDSD 00769 END IF 00770 * 00771 * Do test (5) 00772 * 00773 CALL ZLACPY( ' ', N, N, A, LDA, S, LDA ) 00774 CALL ZLACPY( ' ', N, N, B, LDA, T, LDA ) 00775 CALL ZGGEV( 'N', 'N', N, S, LDA, T, LDA, ALPHA1, BETA1, Q, 00776 $ LDQ, Z, LDQ, WORK, LWORK, RWORK, IERR ) 00777 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN 00778 RESULT( 1 ) = ULPINV 00779 WRITE( NOUNIT, FMT = 9999 )'ZGGEV2', IERR, N, JTYPE, 00780 $ IOLDSD 00781 INFO = ABS( IERR ) 00782 GO TO 190 00783 END IF 00784 * 00785 DO 120 J = 1, N 00786 IF( ALPHA( J ).NE.ALPHA1( J ) .OR. BETA( J ).NE. 00787 $ BETA1( J ) )RESULT( 5 ) = ULPINV 00788 120 CONTINUE 00789 * 00790 * Do test (6): Compute eigenvalues and left eigenvectors, 00791 * and test them 00792 * 00793 CALL ZLACPY( ' ', N, N, A, LDA, S, LDA ) 00794 CALL ZLACPY( ' ', N, N, B, LDA, T, LDA ) 00795 CALL ZGGEV( 'V', 'N', N, S, LDA, T, LDA, ALPHA1, BETA1, QE, 00796 $ LDQE, Z, LDQ, WORK, LWORK, RWORK, IERR ) 00797 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN 00798 RESULT( 1 ) = ULPINV 00799 WRITE( NOUNIT, FMT = 9999 )'ZGGEV3', IERR, N, JTYPE, 00800 $ IOLDSD 00801 INFO = ABS( IERR ) 00802 GO TO 190 00803 END IF 00804 * 00805 DO 130 J = 1, N 00806 IF( ALPHA( J ).NE.ALPHA1( J ) .OR. BETA( J ).NE. 00807 $ BETA1( J ) )RESULT( 6 ) = ULPINV 00808 130 CONTINUE 00809 * 00810 DO 150 J = 1, N 00811 DO 140 JC = 1, N 00812 IF( Q( J, JC ).NE.QE( J, JC ) ) 00813 $ RESULT( 6 ) = ULPINV 00814 140 CONTINUE 00815 150 CONTINUE 00816 * 00817 * Do test (7): Compute eigenvalues and right eigenvectors, 00818 * and test them 00819 * 00820 CALL ZLACPY( ' ', N, N, A, LDA, S, LDA ) 00821 CALL ZLACPY( ' ', N, N, B, LDA, T, LDA ) 00822 CALL ZGGEV( 'N', 'V', N, S, LDA, T, LDA, ALPHA1, BETA1, Q, 00823 $ LDQ, QE, LDQE, WORK, LWORK, RWORK, IERR ) 00824 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN 00825 RESULT( 1 ) = ULPINV 00826 WRITE( NOUNIT, FMT = 9999 )'ZGGEV4', IERR, N, JTYPE, 00827 $ IOLDSD 00828 INFO = ABS( IERR ) 00829 GO TO 190 00830 END IF 00831 * 00832 DO 160 J = 1, N 00833 IF( ALPHA( J ).NE.ALPHA1( J ) .OR. BETA( J ).NE. 00834 $ BETA1( J ) )RESULT( 7 ) = ULPINV 00835 160 CONTINUE 00836 * 00837 DO 180 J = 1, N 00838 DO 170 JC = 1, N 00839 IF( Z( J, JC ).NE.QE( J, JC ) ) 00840 $ RESULT( 7 ) = ULPINV 00841 170 CONTINUE 00842 180 CONTINUE 00843 * 00844 * End of Loop -- Check for RESULT(j) > THRESH 00845 * 00846 190 CONTINUE 00847 * 00848 NTESTT = NTESTT + 7 00849 * 00850 * Print out tests which fail. 00851 * 00852 DO 200 JR = 1, 7 00853 IF( RESULT( JR ).GE.THRESH ) THEN 00854 * 00855 * If this is the first test to fail, 00856 * print a header to the data file. 00857 * 00858 IF( NERRS.EQ.0 ) THEN 00859 WRITE( NOUNIT, FMT = 9997 )'ZGV' 00860 * 00861 * Matrix types 00862 * 00863 WRITE( NOUNIT, FMT = 9996 ) 00864 WRITE( NOUNIT, FMT = 9995 ) 00865 WRITE( NOUNIT, FMT = 9994 )'Orthogonal' 00866 * 00867 * Tests performed 00868 * 00869 WRITE( NOUNIT, FMT = 9993 ) 00870 * 00871 END IF 00872 NERRS = NERRS + 1 00873 IF( RESULT( JR ).LT.10000.0D0 ) THEN 00874 WRITE( NOUNIT, FMT = 9992 )N, JTYPE, IOLDSD, JR, 00875 $ RESULT( JR ) 00876 ELSE 00877 WRITE( NOUNIT, FMT = 9991 )N, JTYPE, IOLDSD, JR, 00878 $ RESULT( JR ) 00879 END IF 00880 END IF 00881 200 CONTINUE 00882 * 00883 210 CONTINUE 00884 220 CONTINUE 00885 * 00886 * Summary 00887 * 00888 CALL ALASVM( 'ZGV', NOUNIT, NERRS, NTESTT, 0 ) 00889 * 00890 WORK( 1 ) = MAXWRK 00891 * 00892 RETURN 00893 * 00894 9999 FORMAT( ' ZDRGEV: ', A, ' returned INFO=', I6, '.', / 3X, 'N=', 00895 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 00896 * 00897 9998 FORMAT( ' ZDRGEV: ', A, ' Eigenvectors from ', A, ' incorrectly ', 00898 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 3X, 00899 $ 'N=', I4, ', JTYPE=', I3, ', ISEED=(', 3( I4, ',' ), I5, 00900 $ ')' ) 00901 * 00902 9997 FORMAT( / 1X, A3, ' -- Complex Generalized eigenvalue problem ', 00903 $ 'driver' ) 00904 * 00905 9996 FORMAT( ' Matrix types (see ZDRGEV for details): ' ) 00906 * 00907 9995 FORMAT( ' Special Matrices:', 23X, 00908 $ '(J''=transposed Jordan block)', 00909 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ', 00910 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ', 00911 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I', 00912 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', / 00913 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ', 00914 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' ) 00915 9994 FORMAT( ' Matrices Rotated by Random ', A, ' Matrices U, V:', 00916 $ / ' 16=Transposed Jordan Blocks 19=geometric ', 00917 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ', 00918 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ', 00919 $ 'alpha, beta=0,1 21=random alpha, beta=0,1', 00920 $ / ' Large & Small Matrices:', / ' 22=(large, small) ', 00921 $ '23=(small,large) 24=(small,small) 25=(large,large)', 00922 $ / ' 26=random O(1) matrices.' ) 00923 * 00924 9993 FORMAT( / ' Tests performed: ', 00925 $ / ' 1 = max | ( b A - a B )''*l | / const.,', 00926 $ / ' 2 = | |VR(i)| - 1 | / ulp,', 00927 $ / ' 3 = max | ( b A - a B )*r | / const.', 00928 $ / ' 4 = | |VL(i)| - 1 | / ulp,', 00929 $ / ' 5 = 0 if W same no matter if r or l computed,', 00930 $ / ' 6 = 0 if l same no matter if l computed,', 00931 $ / ' 7 = 0 if r same no matter if r computed,', / 1X ) 00932 9992 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=', 00933 $ 4( I4, ',' ), ' result ', I2, ' is', 0P, F8.2 ) 00934 9991 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=', 00935 $ 4( I4, ',' ), ' result ', I2, ' is', 1P, D10.3 ) 00936 * 00937 * End of ZDRGEV 00938 * 00939 END