![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SDRGEV 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 SDRGEV( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 00012 * NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, QE, LDQE, 00013 * ALPHAR, ALPHAI, BETA, ALPHR1, ALPHI1, BETA1, 00014 * WORK, LWORK, RESULT, INFO ) 00015 * 00016 * .. Scalar Arguments .. 00017 * INTEGER INFO, LDA, LDQ, LDQE, LWORK, NOUNIT, NSIZES, 00018 * $ NTYPES 00019 * REAL THRESH 00020 * .. 00021 * .. Array Arguments .. 00022 * LOGICAL DOTYPE( * ) 00023 * INTEGER ISEED( 4 ), NN( * ) 00024 * REAL A( LDA, * ), ALPHAI( * ), ALPHI1( * ), 00025 * $ ALPHAR( * ), ALPHR1( * ), B( LDA, * ), 00026 * $ BETA( * ), BETA1( * ), Q( LDQ, * ), 00027 * $ QE( LDQE, * ), RESULT( * ), S( LDA, * ), 00028 * $ T( LDA, * ), WORK( * ), Z( LDQ, * ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> SDRGEV checks the nonsymmetric generalized eigenvalue problem driver 00038 *> routine SGGEV. 00039 *> 00040 *> SGGEV 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 SDRGEV 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 SGGEV: 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 *> SDRGES 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, SDRGES 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 SDRGES to continue the same random number 00233 *> sequence. 00234 *> \endverbatim 00235 *> 00236 *> \param[in] THRESH 00237 *> \verbatim 00238 *> THRESH is REAL 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 REAL array, 00257 *> dimension(LDA, max(NN)) 00258 *> Used to hold the original A matrix. Used as input only 00259 *> if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and 00260 *> DOTYPE(MAXTYP+1)=.TRUE. 00261 *> \endverbatim 00262 *> 00263 *> \param[in] LDA 00264 *> \verbatim 00265 *> LDA is INTEGER 00266 *> The leading dimension of A, B, S, and T. 00267 *> It must be at least 1 and at least max( NN ). 00268 *> \endverbatim 00269 *> 00270 *> \param[in,out] B 00271 *> \verbatim 00272 *> B is REAL array, 00273 *> dimension(LDA, max(NN)) 00274 *> Used to hold the original B matrix. Used as input only 00275 *> if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and 00276 *> DOTYPE(MAXTYP+1)=.TRUE. 00277 *> \endverbatim 00278 *> 00279 *> \param[out] S 00280 *> \verbatim 00281 *> S is REAL array, 00282 *> dimension (LDA, max(NN)) 00283 *> The Schur form matrix computed from A by SGGES. On exit, S 00284 *> contains the Schur form matrix corresponding to the matrix 00285 *> in A. 00286 *> \endverbatim 00287 *> 00288 *> \param[out] T 00289 *> \verbatim 00290 *> T is REAL array, 00291 *> dimension (LDA, max(NN)) 00292 *> The upper triangular matrix computed from B by SGGES. 00293 *> \endverbatim 00294 *> 00295 *> \param[out] Q 00296 *> \verbatim 00297 *> Q is REAL array, 00298 *> dimension (LDQ, max(NN)) 00299 *> The (left) eigenvectors matrix computed by SGGEV. 00300 *> \endverbatim 00301 *> 00302 *> \param[in] LDQ 00303 *> \verbatim 00304 *> LDQ is INTEGER 00305 *> The leading dimension of Q and Z. It must 00306 *> be at least 1 and at least max( NN ). 00307 *> \endverbatim 00308 *> 00309 *> \param[out] Z 00310 *> \verbatim 00311 *> Z is REAL array, dimension( LDQ, max(NN) ) 00312 *> The (right) orthogonal matrix computed by SGGES. 00313 *> \endverbatim 00314 *> 00315 *> \param[out] QE 00316 *> \verbatim 00317 *> QE is REAL array, dimension( LDQ, max(NN) ) 00318 *> QE holds the computed right or left eigenvectors. 00319 *> \endverbatim 00320 *> 00321 *> \param[in] LDQE 00322 *> \verbatim 00323 *> LDQE is INTEGER 00324 *> The leading dimension of QE. LDQE >= max(1,max(NN)). 00325 *> \endverbatim 00326 *> 00327 *> \param[out] ALPHAR 00328 *> \verbatim 00329 *> ALPHAR is REAL array, dimension (max(NN)) 00330 *> \endverbatim 00331 *> 00332 *> \param[out] ALPHAI 00333 *> \verbatim 00334 *> ALPHAI is REAL array, dimension (max(NN)) 00335 *> \endverbatim 00336 *> 00337 *> \param[out] BETA 00338 *> \verbatim 00339 *> BETA is REAL array, dimension (max(NN)) 00340 *> \verbatim 00341 *> The generalized eigenvalues of (A,B) computed by SGGEV. 00342 *> ( ALPHAR(k)+ALPHAI(k)*i ) / BETA(k) is the k-th 00343 *> generalized eigenvalue of A and B. 00344 *> \endverbatim 00345 *> 00346 *> \param[out] ALPHR1 00347 *> \verbatim 00348 *> ALPHR1 is REAL array, dimension (max(NN)) 00349 *> \endverbatim 00350 *> 00351 *> \param[out] ALPHI1 00352 *> \verbatim 00353 *> ALPHI1 is REAL array, dimension (max(NN)) 00354 *> \endverbatim 00355 *> 00356 *> \param[out] BETA1 00357 *> \verbatim 00358 *> BETA1 is REAL array, dimension (max(NN)) 00359 *> 00360 *> Like ALPHAR, ALPHAI, BETA, these arrays contain the 00361 *> eigenvalues of A and B, but those computed when SGGEV only 00362 *> computes a partial eigendecomposition, i.e. not the 00363 *> eigenvalues and left and right eigenvectors. 00364 *> \endverbatim 00365 *> 00366 *> \param[out] WORK 00367 *> \verbatim 00368 *> WORK is REAL array, dimension (LWORK) 00369 *> \endverbatim 00370 *> 00371 *> \param[in] LWORK 00372 *> \verbatim 00373 *> LWORK is INTEGER 00374 *> The number of entries in WORK. LWORK >= MAX( 8*N, N*(N+1) ). 00375 *> \endverbatim 00376 *> 00377 *> \param[out] RESULT 00378 *> \verbatim 00379 *> RESULT is REAL array, dimension (2) 00380 *> The values computed by the tests described above. 00381 *> The values are currently limited to 1/ulp, to avoid overflow. 00382 *> \endverbatim 00383 *> 00384 *> \param[out] INFO 00385 *> \verbatim 00386 *> INFO is INTEGER 00387 *> = 0: successful exit 00388 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00389 *> > 0: A routine returned an error code. INFO is the 00390 *> absolute value of the INFO value returned. 00391 *> \endverbatim 00392 * 00393 * Authors: 00394 * ======== 00395 * 00396 *> \author Univ. of Tennessee 00397 *> \author Univ. of California Berkeley 00398 *> \author Univ. of Colorado Denver 00399 *> \author NAG Ltd. 00400 * 00401 *> \date November 2011 00402 * 00403 *> \ingroup single_eig 00404 * 00405 * ===================================================================== 00406 SUBROUTINE SDRGEV( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 00407 $ NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, QE, LDQE, 00408 $ ALPHAR, ALPHAI, BETA, ALPHR1, ALPHI1, BETA1, 00409 $ WORK, LWORK, RESULT, INFO ) 00410 * 00411 * -- LAPACK test routine (version 3.4.0) -- 00412 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00413 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00414 * November 2011 00415 * 00416 * .. Scalar Arguments .. 00417 INTEGER INFO, LDA, LDQ, LDQE, LWORK, NOUNIT, NSIZES, 00418 $ NTYPES 00419 REAL THRESH 00420 * .. 00421 * .. Array Arguments .. 00422 LOGICAL DOTYPE( * ) 00423 INTEGER ISEED( 4 ), NN( * ) 00424 REAL A( LDA, * ), ALPHAI( * ), ALPHI1( * ), 00425 $ ALPHAR( * ), ALPHR1( * ), B( LDA, * ), 00426 $ BETA( * ), BETA1( * ), Q( LDQ, * ), 00427 $ QE( LDQE, * ), RESULT( * ), S( LDA, * ), 00428 $ T( LDA, * ), WORK( * ), Z( LDQ, * ) 00429 * .. 00430 * 00431 * ===================================================================== 00432 * 00433 * .. Parameters .. 00434 REAL ZERO, ONE 00435 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00436 INTEGER MAXTYP 00437 PARAMETER ( MAXTYP = 26 ) 00438 * .. 00439 * .. Local Scalars .. 00440 LOGICAL BADNN 00441 INTEGER I, IADD, IERR, IN, J, JC, JR, JSIZE, JTYPE, 00442 $ MAXWRK, MINWRK, MTYPES, N, N1, NERRS, NMATS, 00443 $ NMAX, NTESTT 00444 REAL SAFMAX, SAFMIN, ULP, ULPINV 00445 * .. 00446 * .. Local Arrays .. 00447 INTEGER IASIGN( MAXTYP ), IBSIGN( MAXTYP ), 00448 $ IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ), 00449 $ KATYPE( MAXTYP ), KAZERO( MAXTYP ), 00450 $ KBMAGN( MAXTYP ), KBTYPE( MAXTYP ), 00451 $ KBZERO( MAXTYP ), KCLASS( MAXTYP ), 00452 $ KTRIAN( MAXTYP ), KZ1( 6 ), KZ2( 6 ) 00453 REAL RMAGN( 0: 3 ) 00454 * .. 00455 * .. External Functions .. 00456 INTEGER ILAENV 00457 REAL SLAMCH, SLARND 00458 EXTERNAL ILAENV, SLAMCH, SLARND 00459 * .. 00460 * .. External Subroutines .. 00461 EXTERNAL ALASVM, SGET52, SGGEV, SLABAD, SLACPY, SLARFG, 00462 $ SLASET, SLATM4, SORM2R, XERBLA 00463 * .. 00464 * .. Intrinsic Functions .. 00465 INTRINSIC ABS, MAX, MIN, REAL, SIGN 00466 * .. 00467 * .. Data statements .. 00468 DATA KCLASS / 15*1, 10*2, 1*3 / 00469 DATA KZ1 / 0, 1, 2, 1, 3, 3 / 00470 DATA KZ2 / 0, 0, 1, 2, 1, 1 / 00471 DATA KADD / 0, 0, 0, 0, 3, 2 / 00472 DATA KATYPE / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4, 00473 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 / 00474 DATA KBTYPE / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4, 00475 $ 1, 1, -4, 2, -4, 8*8, 0 / 00476 DATA KAZERO / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3, 00477 $ 4*5, 4*3, 1 / 00478 DATA KBZERO / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4, 00479 $ 4*6, 4*4, 1 / 00480 DATA KAMAGN / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3, 00481 $ 2, 1 / 00482 DATA KBMAGN / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3, 00483 $ 2, 1 / 00484 DATA KTRIAN / 16*0, 10*1 / 00485 DATA IASIGN / 6*0, 2, 0, 2*2, 2*0, 3*2, 0, 2, 3*0, 00486 $ 5*2, 0 / 00487 DATA IBSIGN / 7*0, 2, 2*0, 2*2, 2*0, 2, 0, 2, 9*0 / 00488 * .. 00489 * .. Executable Statements .. 00490 * 00491 * Check for errors 00492 * 00493 INFO = 0 00494 * 00495 BADNN = .FALSE. 00496 NMAX = 1 00497 DO 10 J = 1, NSIZES 00498 NMAX = MAX( NMAX, NN( J ) ) 00499 IF( NN( J ).LT.0 ) 00500 $ BADNN = .TRUE. 00501 10 CONTINUE 00502 * 00503 IF( NSIZES.LT.0 ) THEN 00504 INFO = -1 00505 ELSE IF( BADNN ) THEN 00506 INFO = -2 00507 ELSE IF( NTYPES.LT.0 ) THEN 00508 INFO = -3 00509 ELSE IF( THRESH.LT.ZERO ) THEN 00510 INFO = -6 00511 ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN 00512 INFO = -9 00513 ELSE IF( LDQ.LE.1 .OR. LDQ.LT.NMAX ) THEN 00514 INFO = -14 00515 ELSE IF( LDQE.LE.1 .OR. LDQE.LT.NMAX ) THEN 00516 INFO = -17 00517 END IF 00518 * 00519 * Compute workspace 00520 * (Note: Comments in the code beginning "Workspace:" describe the 00521 * minimal amount of workspace needed at that point in the code, 00522 * as well as the preferred amount for good performance. 00523 * NB refers to the optimal block size for the immediately 00524 * following subroutine, as returned by ILAENV. 00525 * 00526 MINWRK = 1 00527 IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN 00528 MINWRK = MAX( 1, 8*NMAX, NMAX*( NMAX+1 ) ) 00529 MAXWRK = 7*NMAX + NMAX*ILAENV( 1, 'SGEQRF', ' ', NMAX, 1, NMAX, 00530 $ 0 ) 00531 MAXWRK = MAX( MAXWRK, NMAX*( NMAX+1 ) ) 00532 WORK( 1 ) = MAXWRK 00533 END IF 00534 * 00535 IF( LWORK.LT.MINWRK ) 00536 $ INFO = -25 00537 * 00538 IF( INFO.NE.0 ) THEN 00539 CALL XERBLA( 'SDRGEV', -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 SAFMIN = SLAMCH( 'Safe minimum' ) 00549 ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' ) 00550 SAFMIN = SAFMIN / ULP 00551 SAFMAX = ONE / SAFMIN 00552 CALL SLABAD( 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 / REAL( 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 * KCLASS: =1 means w/o rotation, =2 means w/ rotation, 00594 * =3 means random. 00595 * KATYPE: the "type" to be passed to SLATM4 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 * IASIGN: 1 if the diagonal elements of A are to be 00604 * multiplied by a random magnitude 1 number, =2 if 00605 * randomly chosen diagonal blocks are to be rotated 00606 * to form 2x2 blocks. 00607 * KBTYPE, KBZERO, KBMAGN, IBSIGN: the same, but for B. 00608 * KTRIAN: =0: don't fill in the upper triangle, =1: do. 00609 * KZ1, KZ2, KADD: used to implement KAZERO and KBZERO. 00610 * RMAGN: used to implement KAMAGN and KBMAGN. 00611 * 00612 IF( MTYPES.GT.MAXTYP ) 00613 $ GO TO 100 00614 IERR = 0 00615 IF( KCLASS( JTYPE ).LT.3 ) THEN 00616 * 00617 * Generate A (w/o rotation) 00618 * 00619 IF( ABS( KATYPE( JTYPE ) ).EQ.3 ) THEN 00620 IN = 2*( ( N-1 ) / 2 ) + 1 00621 IF( IN.NE.N ) 00622 $ CALL SLASET( 'Full', N, N, ZERO, ZERO, A, LDA ) 00623 ELSE 00624 IN = N 00625 END IF 00626 CALL SLATM4( KATYPE( JTYPE ), IN, KZ1( KAZERO( JTYPE ) ), 00627 $ KZ2( KAZERO( JTYPE ) ), IASIGN( JTYPE ), 00628 $ RMAGN( KAMAGN( JTYPE ) ), ULP, 00629 $ RMAGN( KTRIAN( JTYPE )*KAMAGN( JTYPE ) ), 2, 00630 $ ISEED, A, LDA ) 00631 IADD = KADD( KAZERO( JTYPE ) ) 00632 IF( IADD.GT.0 .AND. IADD.LE.N ) 00633 $ A( IADD, IADD ) = ONE 00634 * 00635 * Generate B (w/o rotation) 00636 * 00637 IF( ABS( KBTYPE( JTYPE ) ).EQ.3 ) THEN 00638 IN = 2*( ( N-1 ) / 2 ) + 1 00639 IF( IN.NE.N ) 00640 $ CALL SLASET( 'Full', N, N, ZERO, ZERO, B, LDA ) 00641 ELSE 00642 IN = N 00643 END IF 00644 CALL SLATM4( KBTYPE( JTYPE ), IN, KZ1( KBZERO( JTYPE ) ), 00645 $ KZ2( KBZERO( JTYPE ) ), IBSIGN( JTYPE ), 00646 $ RMAGN( KBMAGN( JTYPE ) ), ONE, 00647 $ RMAGN( KTRIAN( JTYPE )*KBMAGN( JTYPE ) ), 2, 00648 $ ISEED, B, LDA ) 00649 IADD = KADD( KBZERO( JTYPE ) ) 00650 IF( IADD.NE.0 .AND. IADD.LE.N ) 00651 $ B( IADD, IADD ) = ONE 00652 * 00653 IF( KCLASS( JTYPE ).EQ.2 .AND. N.GT.0 ) THEN 00654 * 00655 * Include rotations 00656 * 00657 * Generate Q, Z as Householder transformations times 00658 * a diagonal matrix. 00659 * 00660 DO 40 JC = 1, N - 1 00661 DO 30 JR = JC, N 00662 Q( JR, JC ) = SLARND( 3, ISEED ) 00663 Z( JR, JC ) = SLARND( 3, ISEED ) 00664 30 CONTINUE 00665 CALL SLARFG( N+1-JC, Q( JC, JC ), Q( JC+1, JC ), 1, 00666 $ WORK( JC ) ) 00667 WORK( 2*N+JC ) = SIGN( ONE, Q( JC, JC ) ) 00668 Q( JC, JC ) = ONE 00669 CALL SLARFG( N+1-JC, Z( JC, JC ), Z( JC+1, JC ), 1, 00670 $ WORK( N+JC ) ) 00671 WORK( 3*N+JC ) = SIGN( ONE, Z( JC, JC ) ) 00672 Z( JC, JC ) = ONE 00673 40 CONTINUE 00674 Q( N, N ) = ONE 00675 WORK( N ) = ZERO 00676 WORK( 3*N ) = SIGN( ONE, SLARND( 2, ISEED ) ) 00677 Z( N, N ) = ONE 00678 WORK( 2*N ) = ZERO 00679 WORK( 4*N ) = SIGN( ONE, SLARND( 2, ISEED ) ) 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 )*WORK( 3*N+JC )* 00686 $ A( JR, JC ) 00687 B( JR, JC ) = WORK( 2*N+JR )*WORK( 3*N+JC )* 00688 $ B( JR, JC ) 00689 50 CONTINUE 00690 60 CONTINUE 00691 CALL SORM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, A, 00692 $ LDA, WORK( 2*N+1 ), IERR ) 00693 IF( IERR.NE.0 ) 00694 $ GO TO 90 00695 CALL SORM2R( 'R', 'T', N, N, N-1, Z, LDQ, WORK( N+1 ), 00696 $ A, LDA, WORK( 2*N+1 ), IERR ) 00697 IF( IERR.NE.0 ) 00698 $ GO TO 90 00699 CALL SORM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, B, 00700 $ LDA, WORK( 2*N+1 ), IERR ) 00701 IF( IERR.NE.0 ) 00702 $ GO TO 90 00703 CALL SORM2R( 'R', 'T', N, N, N-1, Z, LDQ, WORK( N+1 ), 00704 $ B, LDA, WORK( 2*N+1 ), IERR ) 00705 IF( IERR.NE.0 ) 00706 $ GO TO 90 00707 END IF 00708 ELSE 00709 * 00710 * Random matrices 00711 * 00712 DO 80 JC = 1, N 00713 DO 70 JR = 1, N 00714 A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )* 00715 $ SLARND( 2, ISEED ) 00716 B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )* 00717 $ SLARND( 2, ISEED ) 00718 70 CONTINUE 00719 80 CONTINUE 00720 END IF 00721 * 00722 90 CONTINUE 00723 * 00724 IF( IERR.NE.0 ) THEN 00725 WRITE( NOUNIT, FMT = 9999 )'Generator', IERR, N, JTYPE, 00726 $ IOLDSD 00727 INFO = ABS( IERR ) 00728 RETURN 00729 END IF 00730 * 00731 100 CONTINUE 00732 * 00733 DO 110 I = 1, 7 00734 RESULT( I ) = -ONE 00735 110 CONTINUE 00736 * 00737 * Call SGGEV to compute eigenvalues and eigenvectors. 00738 * 00739 CALL SLACPY( ' ', N, N, A, LDA, S, LDA ) 00740 CALL SLACPY( ' ', N, N, B, LDA, T, LDA ) 00741 CALL SGGEV( 'V', 'V', N, S, LDA, T, LDA, ALPHAR, ALPHAI, 00742 $ BETA, Q, LDQ, Z, LDQ, WORK, LWORK, IERR ) 00743 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN 00744 RESULT( 1 ) = ULPINV 00745 WRITE( NOUNIT, FMT = 9999 )'SGGEV1', IERR, N, JTYPE, 00746 $ IOLDSD 00747 INFO = ABS( IERR ) 00748 GO TO 190 00749 END IF 00750 * 00751 * Do the tests (1) and (2) 00752 * 00753 CALL SGET52( .TRUE., N, A, LDA, B, LDA, Q, LDQ, ALPHAR, 00754 $ ALPHAI, BETA, WORK, RESULT( 1 ) ) 00755 IF( RESULT( 2 ).GT.THRESH ) THEN 00756 WRITE( NOUNIT, FMT = 9998 )'Left', 'SGGEV1', 00757 $ RESULT( 2 ), N, JTYPE, IOLDSD 00758 END IF 00759 * 00760 * Do the tests (3) and (4) 00761 * 00762 CALL SGET52( .FALSE., N, A, LDA, B, LDA, Z, LDQ, ALPHAR, 00763 $ ALPHAI, BETA, WORK, RESULT( 3 ) ) 00764 IF( RESULT( 4 ).GT.THRESH ) THEN 00765 WRITE( NOUNIT, FMT = 9998 )'Right', 'SGGEV1', 00766 $ RESULT( 4 ), N, JTYPE, IOLDSD 00767 END IF 00768 * 00769 * Do the test (5) 00770 * 00771 CALL SLACPY( ' ', N, N, A, LDA, S, LDA ) 00772 CALL SLACPY( ' ', N, N, B, LDA, T, LDA ) 00773 CALL SGGEV( 'N', 'N', N, S, LDA, T, LDA, ALPHR1, ALPHI1, 00774 $ BETA1, Q, LDQ, Z, LDQ, WORK, LWORK, IERR ) 00775 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN 00776 RESULT( 1 ) = ULPINV 00777 WRITE( NOUNIT, FMT = 9999 )'SGGEV2', IERR, N, JTYPE, 00778 $ IOLDSD 00779 INFO = ABS( IERR ) 00780 GO TO 190 00781 END IF 00782 * 00783 DO 120 J = 1, N 00784 IF( ALPHAR( J ).NE.ALPHR1( J ) .OR. ALPHAI( J ).NE. 00785 $ ALPHI1( J ) .OR. BETA( J ).NE.BETA1( J ) ) 00786 $ RESULT( 5 ) = ULPINV 00787 120 CONTINUE 00788 * 00789 * Do the test (6): Compute eigenvalues and left eigenvectors, 00790 * and test them 00791 * 00792 CALL SLACPY( ' ', N, N, A, LDA, S, LDA ) 00793 CALL SLACPY( ' ', N, N, B, LDA, T, LDA ) 00794 CALL SGGEV( 'V', 'N', N, S, LDA, T, LDA, ALPHR1, ALPHI1, 00795 $ BETA1, QE, LDQE, Z, LDQ, WORK, LWORK, IERR ) 00796 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN 00797 RESULT( 1 ) = ULPINV 00798 WRITE( NOUNIT, FMT = 9999 )'SGGEV3', IERR, N, JTYPE, 00799 $ IOLDSD 00800 INFO = ABS( IERR ) 00801 GO TO 190 00802 END IF 00803 * 00804 DO 130 J = 1, N 00805 IF( ALPHAR( J ).NE.ALPHR1( J ) .OR. ALPHAI( J ).NE. 00806 $ ALPHI1( J ) .OR. BETA( J ).NE.BETA1( J ) ) 00807 $ 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 the test (7): Compute eigenvalues and right eigenvectors, 00818 * and test them 00819 * 00820 CALL SLACPY( ' ', N, N, A, LDA, S, LDA ) 00821 CALL SLACPY( ' ', N, N, B, LDA, T, LDA ) 00822 CALL SGGEV( 'N', 'V', N, S, LDA, T, LDA, ALPHR1, ALPHI1, 00823 $ BETA1, Q, LDQ, QE, LDQE, WORK, LWORK, IERR ) 00824 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN 00825 RESULT( 1 ) = ULPINV 00826 WRITE( NOUNIT, FMT = 9999 )'SGGEV4', 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( ALPHAR( J ).NE.ALPHR1( J ) .OR. ALPHAI( J ).NE. 00834 $ ALPHI1( J ) .OR. BETA( J ).NE.BETA1( J ) ) 00835 $ RESULT( 7 ) = ULPINV 00836 160 CONTINUE 00837 * 00838 DO 180 J = 1, N 00839 DO 170 JC = 1, N 00840 IF( Z( J, JC ).NE.QE( J, JC ) ) 00841 $ RESULT( 7 ) = ULPINV 00842 170 CONTINUE 00843 180 CONTINUE 00844 * 00845 * End of Loop -- Check for RESULT(j) > THRESH 00846 * 00847 190 CONTINUE 00848 * 00849 NTESTT = NTESTT + 7 00850 * 00851 * Print out tests which fail. 00852 * 00853 DO 200 JR = 1, 7 00854 IF( RESULT( JR ).GE.THRESH ) THEN 00855 * 00856 * If this is the first test to fail, 00857 * print a header to the data file. 00858 * 00859 IF( NERRS.EQ.0 ) THEN 00860 WRITE( NOUNIT, FMT = 9997 )'SGV' 00861 * 00862 * Matrix types 00863 * 00864 WRITE( NOUNIT, FMT = 9996 ) 00865 WRITE( NOUNIT, FMT = 9995 ) 00866 WRITE( NOUNIT, FMT = 9994 )'Orthogonal' 00867 * 00868 * Tests performed 00869 * 00870 WRITE( NOUNIT, FMT = 9993 ) 00871 * 00872 END IF 00873 NERRS = NERRS + 1 00874 IF( RESULT( JR ).LT.10000.0 ) THEN 00875 WRITE( NOUNIT, FMT = 9992 )N, JTYPE, IOLDSD, JR, 00876 $ RESULT( JR ) 00877 ELSE 00878 WRITE( NOUNIT, FMT = 9991 )N, JTYPE, IOLDSD, JR, 00879 $ RESULT( JR ) 00880 END IF 00881 END IF 00882 200 CONTINUE 00883 * 00884 210 CONTINUE 00885 220 CONTINUE 00886 * 00887 * Summary 00888 * 00889 CALL ALASVM( 'SGV', NOUNIT, NERRS, NTESTT, 0 ) 00890 * 00891 WORK( 1 ) = MAXWRK 00892 * 00893 RETURN 00894 * 00895 9999 FORMAT( ' SDRGEV: ', A, ' returned INFO=', I6, '.', / 3X, 'N=', 00896 $ I6, ', JTYPE=', I6, ', ISEED=(', 4( I4, ',' ), I5, ')' ) 00897 * 00898 9998 FORMAT( ' SDRGEV: ', A, ' Eigenvectors from ', A, ' incorrectly ', 00899 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 3X, 00900 $ 'N=', I4, ', JTYPE=', I3, ', ISEED=(', 4( I4, ',' ), I5, 00901 $ ')' ) 00902 * 00903 9997 FORMAT( / 1X, A3, ' -- Real Generalized eigenvalue problem driver' 00904 $ ) 00905 * 00906 9996 FORMAT( ' Matrix types (see SDRGEV for details): ' ) 00907 * 00908 9995 FORMAT( ' Special Matrices:', 23X, 00909 $ '(J''=transposed Jordan block)', 00910 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ', 00911 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ', 00912 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I', 00913 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', / 00914 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ', 00915 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' ) 00916 9994 FORMAT( ' Matrices Rotated by Random ', A, ' Matrices U, V:', 00917 $ / ' 16=Transposed Jordan Blocks 19=geometric ', 00918 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ', 00919 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ', 00920 $ 'alpha, beta=0,1 21=random alpha, beta=0,1', 00921 $ / ' Large & Small Matrices:', / ' 22=(large, small) ', 00922 $ '23=(small,large) 24=(small,small) 25=(large,large)', 00923 $ / ' 26=random O(1) matrices.' ) 00924 * 00925 9993 FORMAT( / ' Tests performed: ', 00926 $ / ' 1 = max | ( b A - a B )''*l | / const.,', 00927 $ / ' 2 = | |VR(i)| - 1 | / ulp,', 00928 $ / ' 3 = max | ( b A - a B )*r | / const.', 00929 $ / ' 4 = | |VL(i)| - 1 | / ulp,', 00930 $ / ' 5 = 0 if W same no matter if r or l computed,', 00931 $ / ' 6 = 0 if l same no matter if l computed,', 00932 $ / ' 7 = 0 if r same no matter if r computed,', / 1X ) 00933 9992 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=', 00934 $ 4( I4, ',' ), ' result ', I2, ' is', 0P, F8.2 ) 00935 9991 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=', 00936 $ 4( I4, ',' ), ' result ', I2, ' is', 1P, E10.3 ) 00937 * 00938 * End of SDRGEV 00939 * 00940 END