![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZDRGSX 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 ZDRGSX( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, AI, 00012 * BI, Z, Q, ALPHA, BETA, C, LDC, S, WORK, LWORK, 00013 * RWORK, IWORK, LIWORK, BWORK, INFO ) 00014 * 00015 * .. Scalar Arguments .. 00016 * INTEGER INFO, LDA, LDC, LIWORK, LWORK, NCMAX, NIN, 00017 * $ NOUT, NSIZE 00018 * DOUBLE PRECISION THRESH 00019 * .. 00020 * .. Array Arguments .. 00021 * LOGICAL BWORK( * ) 00022 * INTEGER IWORK( * ) 00023 * DOUBLE PRECISION RWORK( * ), S( * ) 00024 * COMPLEX*16 A( LDA, * ), AI( LDA, * ), ALPHA( * ), 00025 * $ B( LDA, * ), BETA( * ), BI( LDA, * ), 00026 * $ C( LDC, * ), Q( LDA, * ), WORK( * ), 00027 * $ Z( LDA, * ) 00028 * .. 00029 * 00030 * 00031 *> \par Purpose: 00032 * ============= 00033 *> 00034 *> \verbatim 00035 *> 00036 *> ZDRGSX checks the nonsymmetric generalized eigenvalue (Schur form) 00037 *> problem expert driver ZGGESX. 00038 *> 00039 *> ZGGES factors A and B as Q*S*Z' and Q*T*Z' , where ' means conjugate 00040 *> transpose, S and T are upper triangular (i.e., in generalized Schur 00041 *> form), and Q and Z are unitary. It also computes the generalized 00042 *> eigenvalues (alpha(j),beta(j)), j=1,...,n. Thus, 00043 *> w(j) = alpha(j)/beta(j) is a root of the characteristic equation 00044 *> 00045 *> det( A - w(j) B ) = 0 00046 *> 00047 *> Optionally it also reorders the eigenvalues so that a selected 00048 *> cluster of eigenvalues appears in the leading diagonal block of the 00049 *> Schur forms; computes a reciprocal condition number for the average 00050 *> of the selected eigenvalues; and computes a reciprocal condition 00051 *> number for the right and left deflating subspaces corresponding to 00052 *> the selected eigenvalues. 00053 *> 00054 *> When ZDRGSX is called with NSIZE > 0, five (5) types of built-in 00055 *> matrix pairs are used to test the routine ZGGESX. 00056 *> 00057 *> When ZDRGSX is called with NSIZE = 0, it reads in test matrix data 00058 *> to test ZGGESX. 00059 *> (need more details on what kind of read-in data are needed). 00060 *> 00061 *> For each matrix pair, the following tests will be performed and 00062 *> compared with the threshhold THRESH except for the tests (7) and (9): 00063 *> 00064 *> (1) | A - Q S Z' | / ( |A| n ulp ) 00065 *> 00066 *> (2) | B - Q T Z' | / ( |B| n ulp ) 00067 *> 00068 *> (3) | I - QQ' | / ( n ulp ) 00069 *> 00070 *> (4) | I - ZZ' | / ( n ulp ) 00071 *> 00072 *> (5) if A is in Schur form (i.e. triangular form) 00073 *> 00074 *> (6) maximum over j of D(j) where: 00075 *> 00076 *> |alpha(j) - S(j,j)| |beta(j) - T(j,j)| 00077 *> D(j) = ------------------------ + ----------------------- 00078 *> max(|alpha(j)|,|S(j,j)|) max(|beta(j)|,|T(j,j)|) 00079 *> 00080 *> (7) if sorting worked and SDIM is the number of eigenvalues 00081 *> which were selected. 00082 *> 00083 *> (8) the estimated value DIF does not differ from the true values of 00084 *> Difu and Difl more than a factor 10*THRESH. If the estimate DIF 00085 *> equals zero the corresponding true values of Difu and Difl 00086 *> should be less than EPS*norm(A, B). If the true value of Difu 00087 *> and Difl equal zero, the estimate DIF should be less than 00088 *> EPS*norm(A, B). 00089 *> 00090 *> (9) If INFO = N+3 is returned by ZGGESX, the reordering "failed" 00091 *> and we check that DIF = PL = PR = 0 and that the true value of 00092 *> Difu and Difl is < EPS*norm(A, B). We count the events when 00093 *> INFO=N+3. 00094 *> 00095 *> For read-in test matrices, the same tests are run except that the 00096 *> exact value for DIF (and PL) is input data. Additionally, there is 00097 *> one more test run for read-in test matrices: 00098 *> 00099 *> (10) the estimated value PL does not differ from the true value of 00100 *> PLTRU more than a factor THRESH. If the estimate PL equals 00101 *> zero the corresponding true value of PLTRU should be less than 00102 *> EPS*norm(A, B). If the true value of PLTRU equal zero, the 00103 *> estimate PL should be less than EPS*norm(A, B). 00104 *> 00105 *> Note that for the built-in tests, a total of 10*NSIZE*(NSIZE-1) 00106 *> matrix pairs are generated and tested. NSIZE should be kept small. 00107 *> 00108 *> SVD (routine ZGESVD) is used for computing the true value of DIF_u 00109 *> and DIF_l when testing the built-in test problems. 00110 *> 00111 *> Built-in Test Matrices 00112 *> ====================== 00113 *> 00114 *> All built-in test matrices are the 2 by 2 block of triangular 00115 *> matrices 00116 *> 00117 *> A = [ A11 A12 ] and B = [ B11 B12 ] 00118 *> [ A22 ] [ B22 ] 00119 *> 00120 *> where for different type of A11 and A22 are given as the following. 00121 *> A12 and B12 are chosen so that the generalized Sylvester equation 00122 *> 00123 *> A11*R - L*A22 = -A12 00124 *> B11*R - L*B22 = -B12 00125 *> 00126 *> have prescribed solution R and L. 00127 *> 00128 *> Type 1: A11 = J_m(1,-1) and A_22 = J_k(1-a,1). 00129 *> B11 = I_m, B22 = I_k 00130 *> where J_k(a,b) is the k-by-k Jordan block with ``a'' on 00131 *> diagonal and ``b'' on superdiagonal. 00132 *> 00133 *> Type 2: A11 = (a_ij) = ( 2(.5-sin(i)) ) and 00134 *> B11 = (b_ij) = ( 2(.5-sin(ij)) ) for i=1,...,m, j=i,...,m 00135 *> A22 = (a_ij) = ( 2(.5-sin(i+j)) ) and 00136 *> B22 = (b_ij) = ( 2(.5-sin(ij)) ) for i=m+1,...,k, j=i,...,k 00137 *> 00138 *> Type 3: A11, A22 and B11, B22 are chosen as for Type 2, but each 00139 *> second diagonal block in A_11 and each third diagonal block 00140 *> in A_22 are made as 2 by 2 blocks. 00141 *> 00142 *> Type 4: A11 = ( 20(.5 - sin(ij)) ) and B22 = ( 2(.5 - sin(i+j)) ) 00143 *> for i=1,...,m, j=1,...,m and 00144 *> A22 = ( 20(.5 - sin(i+j)) ) and B22 = ( 2(.5 - sin(ij)) ) 00145 *> for i=m+1,...,k, j=m+1,...,k 00146 *> 00147 *> Type 5: (A,B) and have potentially close or common eigenvalues and 00148 *> very large departure from block diagonality A_11 is chosen 00149 *> as the m x m leading submatrix of A_1: 00150 *> | 1 b | 00151 *> | -b 1 | 00152 *> | 1+d b | 00153 *> | -b 1+d | 00154 *> A_1 = | d 1 | 00155 *> | -1 d | 00156 *> | -d 1 | 00157 *> | -1 -d | 00158 *> | 1 | 00159 *> and A_22 is chosen as the k x k leading submatrix of A_2: 00160 *> | -1 b | 00161 *> | -b -1 | 00162 *> | 1-d b | 00163 *> | -b 1-d | 00164 *> A_2 = | d 1+b | 00165 *> | -1-b d | 00166 *> | -d 1+b | 00167 *> | -1+b -d | 00168 *> | 1-d | 00169 *> and matrix B are chosen as identity matrices (see DLATM5). 00170 *> 00171 *> \endverbatim 00172 * 00173 * Arguments: 00174 * ========== 00175 * 00176 *> \param[in] NSIZE 00177 *> \verbatim 00178 *> NSIZE is INTEGER 00179 *> The maximum size of the matrices to use. NSIZE >= 0. 00180 *> If NSIZE = 0, no built-in tests matrices are used, but 00181 *> read-in test matrices are used to test DGGESX. 00182 *> \endverbatim 00183 *> 00184 *> \param[in] NCMAX 00185 *> \verbatim 00186 *> NCMAX is INTEGER 00187 *> Maximum allowable NMAX for generating Kroneker matrix 00188 *> in call to ZLAKF2 00189 *> \endverbatim 00190 *> 00191 *> \param[in] THRESH 00192 *> \verbatim 00193 *> THRESH is DOUBLE PRECISION 00194 *> A test will count as "failed" if the "error", computed as 00195 *> described above, exceeds THRESH. Note that the error 00196 *> is scaled to be O(1), so THRESH should be a reasonably 00197 *> small multiple of 1, e.g., 10 or 100. In particular, 00198 *> it should not depend on the precision (single vs. double) 00199 *> or the size of the matrix. THRESH >= 0. 00200 *> \endverbatim 00201 *> 00202 *> \param[in] NIN 00203 *> \verbatim 00204 *> NIN is INTEGER 00205 *> The FORTRAN unit number for reading in the data file of 00206 *> problems to solve. 00207 *> \endverbatim 00208 *> 00209 *> \param[in] NOUT 00210 *> \verbatim 00211 *> NOUT is INTEGER 00212 *> The FORTRAN unit number for printing out error messages 00213 *> (e.g., if a routine returns INFO not equal to 0.) 00214 *> \endverbatim 00215 *> 00216 *> \param[out] A 00217 *> \verbatim 00218 *> A is COMPLEX*16 array, dimension (LDA, NSIZE) 00219 *> Used to store the matrix whose eigenvalues are to be 00220 *> computed. On exit, A contains the last matrix actually used. 00221 *> \endverbatim 00222 *> 00223 *> \param[in] LDA 00224 *> \verbatim 00225 *> LDA is INTEGER 00226 *> The leading dimension of A, B, AI, BI, Z and Q, 00227 *> LDA >= max( 1, NSIZE ). For the read-in test, 00228 *> LDA >= max( 1, N ), N is the size of the test matrices. 00229 *> \endverbatim 00230 *> 00231 *> \param[out] B 00232 *> \verbatim 00233 *> B is COMPLEX*16 array, dimension (LDA, NSIZE) 00234 *> Used to store the matrix whose eigenvalues are to be 00235 *> computed. On exit, B contains the last matrix actually used. 00236 *> \endverbatim 00237 *> 00238 *> \param[out] AI 00239 *> \verbatim 00240 *> AI is COMPLEX*16 array, dimension (LDA, NSIZE) 00241 *> Copy of A, modified by ZGGESX. 00242 *> \endverbatim 00243 *> 00244 *> \param[out] BI 00245 *> \verbatim 00246 *> BI is COMPLEX*16 array, dimension (LDA, NSIZE) 00247 *> Copy of B, modified by ZGGESX. 00248 *> \endverbatim 00249 *> 00250 *> \param[out] Z 00251 *> \verbatim 00252 *> Z is COMPLEX*16 array, dimension (LDA, NSIZE) 00253 *> Z holds the left Schur vectors computed by ZGGESX. 00254 *> \endverbatim 00255 *> 00256 *> \param[out] Q 00257 *> \verbatim 00258 *> Q is COMPLEX*16 array, dimension (LDA, NSIZE) 00259 *> Q holds the right Schur vectors computed by ZGGESX. 00260 *> \endverbatim 00261 *> 00262 *> \param[out] ALPHA 00263 *> \verbatim 00264 *> ALPHA is COMPLEX*16 array, dimension (NSIZE) 00265 *> \endverbatim 00266 *> 00267 *> \param[out] BETA 00268 *> \verbatim 00269 *> BETA is COMPLEX*16 array, dimension (NSIZE) 00270 *> 00271 *> On exit, ALPHA/BETA are the eigenvalues. 00272 *> \endverbatim 00273 *> 00274 *> \param[out] C 00275 *> \verbatim 00276 *> C is COMPLEX*16 array, dimension (LDC, LDC) 00277 *> Store the matrix generated by subroutine ZLAKF2, this is the 00278 *> matrix formed by Kronecker products used for estimating 00279 *> DIF. 00280 *> \endverbatim 00281 *> 00282 *> \param[in] LDC 00283 *> \verbatim 00284 *> LDC is INTEGER 00285 *> The leading dimension of C. LDC >= max(1, LDA*LDA/2 ). 00286 *> \endverbatim 00287 *> 00288 *> \param[out] S 00289 *> \verbatim 00290 *> S is DOUBLE PRECISION array, dimension (LDC) 00291 *> Singular values of C 00292 *> \endverbatim 00293 *> 00294 *> \param[out] WORK 00295 *> \verbatim 00296 *> WORK is COMPLEX*16 array, dimension (LWORK) 00297 *> \endverbatim 00298 *> 00299 *> \param[in] LWORK 00300 *> \verbatim 00301 *> LWORK is INTEGER 00302 *> The dimension of the array WORK. LWORK >= 3*NSIZE*NSIZE/2 00303 *> \endverbatim 00304 *> 00305 *> \param[out] RWORK 00306 *> \verbatim 00307 *> RWORK is DOUBLE PRECISION array, 00308 *> dimension (5*NSIZE*NSIZE/2 - 4) 00309 *> \endverbatim 00310 *> 00311 *> \param[out] IWORK 00312 *> \verbatim 00313 *> IWORK is INTEGER array, dimension (LIWORK) 00314 *> \endverbatim 00315 *> 00316 *> \param[in] LIWORK 00317 *> \verbatim 00318 *> LIWORK is INTEGER 00319 *> The dimension of the array IWORK. LIWORK >= NSIZE + 2. 00320 *> \endverbatim 00321 *> 00322 *> \param[out] BWORK 00323 *> \verbatim 00324 *> BWORK is LOGICAL array, dimension (NSIZE) 00325 *> \endverbatim 00326 *> 00327 *> \param[out] INFO 00328 *> \verbatim 00329 *> INFO is INTEGER 00330 *> = 0: successful exit 00331 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00332 *> > 0: A routine returned an error code. 00333 *> \endverbatim 00334 * 00335 * Authors: 00336 * ======== 00337 * 00338 *> \author Univ. of Tennessee 00339 *> \author Univ. of California Berkeley 00340 *> \author Univ. of Colorado Denver 00341 *> \author NAG Ltd. 00342 * 00343 *> \date November 2011 00344 * 00345 *> \ingroup complex16_eig 00346 * 00347 * ===================================================================== 00348 SUBROUTINE ZDRGSX( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, AI, 00349 $ BI, Z, Q, ALPHA, BETA, C, LDC, S, WORK, LWORK, 00350 $ RWORK, IWORK, LIWORK, BWORK, INFO ) 00351 * 00352 * -- LAPACK test routine (version 3.4.0) -- 00353 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00354 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00355 * November 2011 00356 * 00357 * .. Scalar Arguments .. 00358 INTEGER INFO, LDA, LDC, LIWORK, LWORK, NCMAX, NIN, 00359 $ NOUT, NSIZE 00360 DOUBLE PRECISION THRESH 00361 * .. 00362 * .. Array Arguments .. 00363 LOGICAL BWORK( * ) 00364 INTEGER IWORK( * ) 00365 DOUBLE PRECISION RWORK( * ), S( * ) 00366 COMPLEX*16 A( LDA, * ), AI( LDA, * ), ALPHA( * ), 00367 $ B( LDA, * ), BETA( * ), BI( LDA, * ), 00368 $ C( LDC, * ), Q( LDA, * ), WORK( * ), 00369 $ Z( LDA, * ) 00370 * .. 00371 * 00372 * ===================================================================== 00373 * 00374 * .. Parameters .. 00375 DOUBLE PRECISION ZERO, ONE, TEN 00376 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 1.0D+1 ) 00377 COMPLEX*16 CZERO 00378 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) ) 00379 * .. 00380 * .. Local Scalars .. 00381 LOGICAL ILABAD 00382 CHARACTER SENSE 00383 INTEGER BDSPAC, I, IFUNC, J, LINFO, MAXWRK, MINWRK, MM, 00384 $ MN2, NERRS, NPTKNT, NTEST, NTESTT, PRTYPE, QBA, 00385 $ QBB 00386 DOUBLE PRECISION ABNRM, BIGNUM, DIFTRU, PLTRU, SMLNUM, TEMP1, 00387 $ TEMP2, THRSH2, ULP, ULPINV, WEIGHT 00388 COMPLEX*16 X 00389 * .. 00390 * .. Local Arrays .. 00391 DOUBLE PRECISION DIFEST( 2 ), PL( 2 ), RESULT( 10 ) 00392 * .. 00393 * .. External Functions .. 00394 LOGICAL ZLCTSX 00395 INTEGER ILAENV 00396 DOUBLE PRECISION DLAMCH, ZLANGE 00397 EXTERNAL ZLCTSX, ILAENV, DLAMCH, ZLANGE 00398 * .. 00399 * .. External Subroutines .. 00400 EXTERNAL ALASVM, DLABAD, XERBLA, ZGESVD, ZGET51, ZGGESX, 00401 $ ZLACPY, ZLAKF2, ZLASET, ZLATM5 00402 * .. 00403 * .. Scalars in Common .. 00404 LOGICAL FS 00405 INTEGER K, M, MPLUSN, N 00406 * .. 00407 * .. Common blocks .. 00408 COMMON / MN / M, N, MPLUSN, K, FS 00409 * .. 00410 * .. Intrinsic Functions .. 00411 INTRINSIC ABS, DBLE, DIMAG, MAX, SQRT 00412 * .. 00413 * .. Statement Functions .. 00414 DOUBLE PRECISION ABS1 00415 * .. 00416 * .. Statement Function definitions .. 00417 ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) ) 00418 * .. 00419 * .. Executable Statements .. 00420 * 00421 * Check for errors 00422 * 00423 INFO = 0 00424 IF( NSIZE.LT.0 ) THEN 00425 INFO = -1 00426 ELSE IF( THRESH.LT.ZERO ) THEN 00427 INFO = -2 00428 ELSE IF( NIN.LE.0 ) THEN 00429 INFO = -3 00430 ELSE IF( NOUT.LE.0 ) THEN 00431 INFO = -4 00432 ELSE IF( LDA.LT.1 .OR. LDA.LT.NSIZE ) THEN 00433 INFO = -6 00434 ELSE IF( LDC.LT.1 .OR. LDC.LT.NSIZE*NSIZE / 2 ) THEN 00435 INFO = -15 00436 ELSE IF( LIWORK.LT.NSIZE+2 ) THEN 00437 INFO = -21 00438 END IF 00439 * 00440 * Compute workspace 00441 * (Note: Comments in the code beginning "Workspace:" describe the 00442 * minimal amount of workspace needed at that point in the code, 00443 * as well as the preferred amount for good performance. 00444 * NB refers to the optimal block size for the immediately 00445 * following subroutine, as returned by ILAENV.) 00446 * 00447 MINWRK = 1 00448 IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN 00449 MINWRK = 3*NSIZE*NSIZE / 2 00450 * 00451 * workspace for cggesx 00452 * 00453 MAXWRK = NSIZE*( 1+ILAENV( 1, 'ZGEQRF', ' ', NSIZE, 1, NSIZE, 00454 $ 0 ) ) 00455 MAXWRK = MAX( MAXWRK, NSIZE*( 1+ILAENV( 1, 'ZUNGQR', ' ', 00456 $ NSIZE, 1, NSIZE, -1 ) ) ) 00457 * 00458 * workspace for zgesvd 00459 * 00460 BDSPAC = 3*NSIZE*NSIZE / 2 00461 MAXWRK = MAX( MAXWRK, NSIZE*NSIZE* 00462 $ ( 1+ILAENV( 1, 'ZGEBRD', ' ', NSIZE*NSIZE / 2, 00463 $ NSIZE*NSIZE / 2, -1, -1 ) ) ) 00464 MAXWRK = MAX( MAXWRK, BDSPAC ) 00465 * 00466 MAXWRK = MAX( MAXWRK, MINWRK ) 00467 * 00468 WORK( 1 ) = MAXWRK 00469 END IF 00470 * 00471 IF( LWORK.LT.MINWRK ) 00472 $ INFO = -18 00473 * 00474 IF( INFO.NE.0 ) THEN 00475 CALL XERBLA( 'ZDRGSX', -INFO ) 00476 RETURN 00477 END IF 00478 * 00479 * Important constants 00480 * 00481 ULP = DLAMCH( 'P' ) 00482 ULPINV = ONE / ULP 00483 SMLNUM = DLAMCH( 'S' ) / ULP 00484 BIGNUM = ONE / SMLNUM 00485 CALL DLABAD( SMLNUM, BIGNUM ) 00486 THRSH2 = TEN*THRESH 00487 NTESTT = 0 00488 NERRS = 0 00489 * 00490 * Go to the tests for read-in matrix pairs 00491 * 00492 IFUNC = 0 00493 IF( NSIZE.EQ.0 ) 00494 $ GO TO 70 00495 * 00496 * Test the built-in matrix pairs. 00497 * Loop over different functions (IFUNC) of ZGGESX, types (PRTYPE) 00498 * of test matrices, different size (M+N) 00499 * 00500 PRTYPE = 0 00501 QBA = 3 00502 QBB = 4 00503 WEIGHT = SQRT( ULP ) 00504 * 00505 DO 60 IFUNC = 0, 3 00506 DO 50 PRTYPE = 1, 5 00507 DO 40 M = 1, NSIZE - 1 00508 DO 30 N = 1, NSIZE - M 00509 * 00510 WEIGHT = ONE / WEIGHT 00511 MPLUSN = M + N 00512 * 00513 * Generate test matrices 00514 * 00515 FS = .TRUE. 00516 K = 0 00517 * 00518 CALL ZLASET( 'Full', MPLUSN, MPLUSN, CZERO, CZERO, AI, 00519 $ LDA ) 00520 CALL ZLASET( 'Full', MPLUSN, MPLUSN, CZERO, CZERO, BI, 00521 $ LDA ) 00522 * 00523 CALL ZLATM5( PRTYPE, M, N, AI, LDA, AI( M+1, M+1 ), 00524 $ LDA, AI( 1, M+1 ), LDA, BI, LDA, 00525 $ BI( M+1, M+1 ), LDA, BI( 1, M+1 ), LDA, 00526 $ Q, LDA, Z, LDA, WEIGHT, QBA, QBB ) 00527 * 00528 * Compute the Schur factorization and swapping the 00529 * m-by-m (1,1)-blocks with n-by-n (2,2)-blocks. 00530 * Swapping is accomplished via the function ZLCTSX 00531 * which is supplied below. 00532 * 00533 IF( IFUNC.EQ.0 ) THEN 00534 SENSE = 'N' 00535 ELSE IF( IFUNC.EQ.1 ) THEN 00536 SENSE = 'E' 00537 ELSE IF( IFUNC.EQ.2 ) THEN 00538 SENSE = 'V' 00539 ELSE IF( IFUNC.EQ.3 ) THEN 00540 SENSE = 'B' 00541 END IF 00542 * 00543 CALL ZLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA ) 00544 CALL ZLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA ) 00545 * 00546 CALL ZGGESX( 'V', 'V', 'S', ZLCTSX, SENSE, MPLUSN, AI, 00547 $ LDA, BI, LDA, MM, ALPHA, BETA, Q, LDA, Z, 00548 $ LDA, PL, DIFEST, WORK, LWORK, RWORK, 00549 $ IWORK, LIWORK, BWORK, LINFO ) 00550 * 00551 IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN 00552 RESULT( 1 ) = ULPINV 00553 WRITE( NOUT, FMT = 9999 )'ZGGESX', LINFO, MPLUSN, 00554 $ PRTYPE 00555 INFO = LINFO 00556 GO TO 30 00557 END IF 00558 * 00559 * Compute the norm(A, B) 00560 * 00561 CALL ZLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK, 00562 $ MPLUSN ) 00563 CALL ZLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, 00564 $ WORK( MPLUSN*MPLUSN+1 ), MPLUSN ) 00565 ABNRM = ZLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN, 00566 $ RWORK ) 00567 * 00568 * Do tests (1) to (4) 00569 * 00570 RESULT( 2 ) = ZERO 00571 CALL ZGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z, 00572 $ LDA, WORK, RWORK, RESULT( 1 ) ) 00573 CALL ZGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z, 00574 $ LDA, WORK, RWORK, RESULT( 2 ) ) 00575 CALL ZGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q, 00576 $ LDA, WORK, RWORK, RESULT( 3 ) ) 00577 CALL ZGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z, 00578 $ LDA, WORK, RWORK, RESULT( 4 ) ) 00579 NTEST = 4 00580 * 00581 * Do tests (5) and (6): check Schur form of A and 00582 * compare eigenvalues with diagonals. 00583 * 00584 TEMP1 = ZERO 00585 RESULT( 5 ) = ZERO 00586 RESULT( 6 ) = ZERO 00587 * 00588 DO 10 J = 1, MPLUSN 00589 ILABAD = .FALSE. 00590 TEMP2 = ( ABS1( ALPHA( J )-AI( J, J ) ) / 00591 $ MAX( SMLNUM, ABS1( ALPHA( J ) ), 00592 $ ABS1( AI( J, J ) ) )+ 00593 $ ABS1( BETA( J )-BI( J, J ) ) / 00594 $ MAX( SMLNUM, ABS1( BETA( J ) ), 00595 $ ABS1( BI( J, J ) ) ) ) / ULP 00596 IF( J.LT.MPLUSN ) THEN 00597 IF( AI( J+1, J ).NE.ZERO ) THEN 00598 ILABAD = .TRUE. 00599 RESULT( 5 ) = ULPINV 00600 END IF 00601 END IF 00602 IF( J.GT.1 ) THEN 00603 IF( AI( J, J-1 ).NE.ZERO ) THEN 00604 ILABAD = .TRUE. 00605 RESULT( 5 ) = ULPINV 00606 END IF 00607 END IF 00608 TEMP1 = MAX( TEMP1, TEMP2 ) 00609 IF( ILABAD ) THEN 00610 WRITE( NOUT, FMT = 9997 )J, MPLUSN, PRTYPE 00611 END IF 00612 10 CONTINUE 00613 RESULT( 6 ) = TEMP1 00614 NTEST = NTEST + 2 00615 * 00616 * Test (7) (if sorting worked) 00617 * 00618 RESULT( 7 ) = ZERO 00619 IF( LINFO.EQ.MPLUSN+3 ) THEN 00620 RESULT( 7 ) = ULPINV 00621 ELSE IF( MM.NE.N ) THEN 00622 RESULT( 7 ) = ULPINV 00623 END IF 00624 NTEST = NTEST + 1 00625 * 00626 * Test (8): compare the estimated value DIF and its 00627 * value. first, compute the exact DIF. 00628 * 00629 RESULT( 8 ) = ZERO 00630 MN2 = MM*( MPLUSN-MM )*2 00631 IF( IFUNC.GE.2 .AND. MN2.LE.NCMAX*NCMAX ) THEN 00632 * 00633 * Note: for either following two cases, there are 00634 * almost same number of test cases fail the test. 00635 * 00636 CALL ZLAKF2( MM, MPLUSN-MM, AI, LDA, 00637 $ AI( MM+1, MM+1 ), BI, 00638 $ BI( MM+1, MM+1 ), C, LDC ) 00639 * 00640 CALL ZGESVD( 'N', 'N', MN2, MN2, C, LDC, S, WORK, 00641 $ 1, WORK( 2 ), 1, WORK( 3 ), LWORK-2, 00642 $ RWORK, INFO ) 00643 DIFTRU = S( MN2 ) 00644 * 00645 IF( DIFEST( 2 ).EQ.ZERO ) THEN 00646 IF( DIFTRU.GT.ABNRM*ULP ) 00647 $ RESULT( 8 ) = ULPINV 00648 ELSE IF( DIFTRU.EQ.ZERO ) THEN 00649 IF( DIFEST( 2 ).GT.ABNRM*ULP ) 00650 $ RESULT( 8 ) = ULPINV 00651 ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR. 00652 $ ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN 00653 RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ), 00654 $ DIFEST( 2 ) / DIFTRU ) 00655 END IF 00656 NTEST = NTEST + 1 00657 END IF 00658 * 00659 * Test (9) 00660 * 00661 RESULT( 9 ) = ZERO 00662 IF( LINFO.EQ.( MPLUSN+2 ) ) THEN 00663 IF( DIFTRU.GT.ABNRM*ULP ) 00664 $ RESULT( 9 ) = ULPINV 00665 IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) ) 00666 $ RESULT( 9 ) = ULPINV 00667 IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) ) 00668 $ RESULT( 9 ) = ULPINV 00669 NTEST = NTEST + 1 00670 END IF 00671 * 00672 NTESTT = NTESTT + NTEST 00673 * 00674 * Print out tests which fail. 00675 * 00676 DO 20 J = 1, 9 00677 IF( RESULT( J ).GE.THRESH ) THEN 00678 * 00679 * If this is the first test to fail, 00680 * print a header to the data file. 00681 * 00682 IF( NERRS.EQ.0 ) THEN 00683 WRITE( NOUT, FMT = 9996 )'ZGX' 00684 * 00685 * Matrix types 00686 * 00687 WRITE( NOUT, FMT = 9994 ) 00688 * 00689 * Tests performed 00690 * 00691 WRITE( NOUT, FMT = 9993 )'unitary', '''', 00692 $ 'transpose', ( '''', I = 1, 4 ) 00693 * 00694 END IF 00695 NERRS = NERRS + 1 00696 IF( RESULT( J ).LT.10000.0D0 ) THEN 00697 WRITE( NOUT, FMT = 9992 )MPLUSN, PRTYPE, 00698 $ WEIGHT, M, J, RESULT( J ) 00699 ELSE 00700 WRITE( NOUT, FMT = 9991 )MPLUSN, PRTYPE, 00701 $ WEIGHT, M, J, RESULT( J ) 00702 END IF 00703 END IF 00704 20 CONTINUE 00705 * 00706 30 CONTINUE 00707 40 CONTINUE 00708 50 CONTINUE 00709 60 CONTINUE 00710 * 00711 GO TO 150 00712 * 00713 70 CONTINUE 00714 * 00715 * Read in data from file to check accuracy of condition estimation 00716 * Read input data until N=0 00717 * 00718 NPTKNT = 0 00719 * 00720 80 CONTINUE 00721 READ( NIN, FMT = *, END = 140 )MPLUSN 00722 IF( MPLUSN.EQ.0 ) 00723 $ GO TO 140 00724 READ( NIN, FMT = *, END = 140 )N 00725 DO 90 I = 1, MPLUSN 00726 READ( NIN, FMT = * )( AI( I, J ), J = 1, MPLUSN ) 00727 90 CONTINUE 00728 DO 100 I = 1, MPLUSN 00729 READ( NIN, FMT = * )( BI( I, J ), J = 1, MPLUSN ) 00730 100 CONTINUE 00731 READ( NIN, FMT = * )PLTRU, DIFTRU 00732 * 00733 NPTKNT = NPTKNT + 1 00734 FS = .TRUE. 00735 K = 0 00736 M = MPLUSN - N 00737 * 00738 CALL ZLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA ) 00739 CALL ZLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA ) 00740 * 00741 * Compute the Schur factorization while swaping the 00742 * m-by-m (1,1)-blocks with n-by-n (2,2)-blocks. 00743 * 00744 CALL ZGGESX( 'V', 'V', 'S', ZLCTSX, 'B', MPLUSN, AI, LDA, BI, LDA, 00745 $ MM, ALPHA, BETA, Q, LDA, Z, LDA, PL, DIFEST, WORK, 00746 $ LWORK, RWORK, IWORK, LIWORK, BWORK, LINFO ) 00747 * 00748 IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN 00749 RESULT( 1 ) = ULPINV 00750 WRITE( NOUT, FMT = 9998 )'ZGGESX', LINFO, MPLUSN, NPTKNT 00751 GO TO 130 00752 END IF 00753 * 00754 * Compute the norm(A, B) 00755 * (should this be norm of (A,B) or (AI,BI)?) 00756 * 00757 CALL ZLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK, MPLUSN ) 00758 CALL ZLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, 00759 $ WORK( MPLUSN*MPLUSN+1 ), MPLUSN ) 00760 ABNRM = ZLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN, RWORK ) 00761 * 00762 * Do tests (1) to (4) 00763 * 00764 CALL ZGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z, LDA, WORK, 00765 $ RWORK, RESULT( 1 ) ) 00766 CALL ZGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z, LDA, WORK, 00767 $ RWORK, RESULT( 2 ) ) 00768 CALL ZGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q, LDA, WORK, 00769 $ RWORK, RESULT( 3 ) ) 00770 CALL ZGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z, LDA, WORK, 00771 $ RWORK, RESULT( 4 ) ) 00772 * 00773 * Do tests (5) and (6): check Schur form of A and compare 00774 * eigenvalues with diagonals. 00775 * 00776 NTEST = 6 00777 TEMP1 = ZERO 00778 RESULT( 5 ) = ZERO 00779 RESULT( 6 ) = ZERO 00780 * 00781 DO 110 J = 1, MPLUSN 00782 ILABAD = .FALSE. 00783 TEMP2 = ( ABS1( ALPHA( J )-AI( J, J ) ) / 00784 $ MAX( SMLNUM, ABS1( ALPHA( J ) ), ABS1( AI( J, J ) ) )+ 00785 $ ABS1( BETA( J )-BI( J, J ) ) / 00786 $ MAX( SMLNUM, ABS1( BETA( J ) ), ABS1( BI( J, J ) ) ) ) 00787 $ / ULP 00788 IF( J.LT.MPLUSN ) THEN 00789 IF( AI( J+1, J ).NE.ZERO ) THEN 00790 ILABAD = .TRUE. 00791 RESULT( 5 ) = ULPINV 00792 END IF 00793 END IF 00794 IF( J.GT.1 ) THEN 00795 IF( AI( J, J-1 ).NE.ZERO ) THEN 00796 ILABAD = .TRUE. 00797 RESULT( 5 ) = ULPINV 00798 END IF 00799 END IF 00800 TEMP1 = MAX( TEMP1, TEMP2 ) 00801 IF( ILABAD ) THEN 00802 WRITE( NOUT, FMT = 9997 )J, MPLUSN, NPTKNT 00803 END IF 00804 110 CONTINUE 00805 RESULT( 6 ) = TEMP1 00806 * 00807 * Test (7) (if sorting worked) <--------- need to be checked. 00808 * 00809 NTEST = 7 00810 RESULT( 7 ) = ZERO 00811 IF( LINFO.EQ.MPLUSN+3 ) 00812 $ RESULT( 7 ) = ULPINV 00813 * 00814 * Test (8): compare the estimated value of DIF and its true value. 00815 * 00816 NTEST = 8 00817 RESULT( 8 ) = ZERO 00818 IF( DIFEST( 2 ).EQ.ZERO ) THEN 00819 IF( DIFTRU.GT.ABNRM*ULP ) 00820 $ RESULT( 8 ) = ULPINV 00821 ELSE IF( DIFTRU.EQ.ZERO ) THEN 00822 IF( DIFEST( 2 ).GT.ABNRM*ULP ) 00823 $ RESULT( 8 ) = ULPINV 00824 ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR. 00825 $ ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN 00826 RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ), DIFEST( 2 ) / DIFTRU ) 00827 END IF 00828 * 00829 * Test (9) 00830 * 00831 NTEST = 9 00832 RESULT( 9 ) = ZERO 00833 IF( LINFO.EQ.( MPLUSN+2 ) ) THEN 00834 IF( DIFTRU.GT.ABNRM*ULP ) 00835 $ RESULT( 9 ) = ULPINV 00836 IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) ) 00837 $ RESULT( 9 ) = ULPINV 00838 IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) ) 00839 $ RESULT( 9 ) = ULPINV 00840 END IF 00841 * 00842 * Test (10): compare the estimated value of PL and it true value. 00843 * 00844 NTEST = 10 00845 RESULT( 10 ) = ZERO 00846 IF( PL( 1 ).EQ.ZERO ) THEN 00847 IF( PLTRU.GT.ABNRM*ULP ) 00848 $ RESULT( 10 ) = ULPINV 00849 ELSE IF( PLTRU.EQ.ZERO ) THEN 00850 IF( PL( 1 ).GT.ABNRM*ULP ) 00851 $ RESULT( 10 ) = ULPINV 00852 ELSE IF( ( PLTRU.GT.THRESH*PL( 1 ) ) .OR. 00853 $ ( PLTRU*THRESH.LT.PL( 1 ) ) ) THEN 00854 RESULT( 10 ) = ULPINV 00855 END IF 00856 * 00857 NTESTT = NTESTT + NTEST 00858 * 00859 * Print out tests which fail. 00860 * 00861 DO 120 J = 1, NTEST 00862 IF( RESULT( J ).GE.THRESH ) THEN 00863 * 00864 * If this is the first test to fail, 00865 * print a header to the data file. 00866 * 00867 IF( NERRS.EQ.0 ) THEN 00868 WRITE( NOUT, FMT = 9996 )'ZGX' 00869 * 00870 * Matrix types 00871 * 00872 WRITE( NOUT, FMT = 9995 ) 00873 * 00874 * Tests performed 00875 * 00876 WRITE( NOUT, FMT = 9993 )'unitary', '''', 'transpose', 00877 $ ( '''', I = 1, 4 ) 00878 * 00879 END IF 00880 NERRS = NERRS + 1 00881 IF( RESULT( J ).LT.10000.0D0 ) THEN 00882 WRITE( NOUT, FMT = 9990 )NPTKNT, MPLUSN, J, RESULT( J ) 00883 ELSE 00884 WRITE( NOUT, FMT = 9989 )NPTKNT, MPLUSN, J, RESULT( J ) 00885 END IF 00886 END IF 00887 * 00888 120 CONTINUE 00889 * 00890 130 CONTINUE 00891 GO TO 80 00892 140 CONTINUE 00893 * 00894 150 CONTINUE 00895 * 00896 * Summary 00897 * 00898 CALL ALASVM( 'ZGX', NOUT, NERRS, NTESTT, 0 ) 00899 * 00900 WORK( 1 ) = MAXWRK 00901 * 00902 RETURN 00903 * 00904 9999 FORMAT( ' ZDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 00905 $ I6, ', JTYPE=', I6, ')' ) 00906 * 00907 9998 FORMAT( ' ZDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 00908 $ I6, ', Input Example #', I2, ')' ) 00909 * 00910 9997 FORMAT( ' ZDRGSX: S not in Schur form at eigenvalue ', I6, '.', 00911 $ / 9X, 'N=', I6, ', JTYPE=', I6, ')' ) 00912 * 00913 9996 FORMAT( / 1X, A3, ' -- Complex Expert Generalized Schur form', 00914 $ ' problem driver' ) 00915 * 00916 9995 FORMAT( 'Input Example' ) 00917 * 00918 9994 FORMAT( ' Matrix types: ', / 00919 $ ' 1: A is a block diagonal matrix of Jordan blocks ', 00920 $ 'and B is the identity ', / ' matrix, ', 00921 $ / ' 2: A and B are upper triangular matrices, ', 00922 $ / ' 3: A and B are as type 2, but each second diagonal ', 00923 $ 'block in A_11 and ', / 00924 $ ' each third diaongal block in A_22 are 2x2 blocks,', 00925 $ / ' 4: A and B are block diagonal matrices, ', 00926 $ / ' 5: (A,B) has potentially close or common ', 00927 $ 'eigenvalues.', / ) 00928 * 00929 9993 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ', 00930 $ 'Q and Z are ', A, ',', / 19X, 00931 $ ' a is alpha, b is beta, and ', A, ' means ', A, '.)', 00932 $ / ' 1 = | A - Q S Z', A, 00933 $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', A, 00934 $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', A, 00935 $ ' | / ( n ulp ) 4 = | I - ZZ', A, 00936 $ ' | / ( n ulp )', / ' 5 = 1/ULP if A is not in ', 00937 $ 'Schur form S', / ' 6 = difference between (alpha,beta)', 00938 $ ' and diagonals of (S,T)', / 00939 $ ' 7 = 1/ULP if SDIM is not the correct number of ', 00940 $ 'selected eigenvalues', / 00941 $ ' 8 = 1/ULP if DIFEST/DIFTRU > 10*THRESH or ', 00942 $ 'DIFTRU/DIFEST > 10*THRESH', 00943 $ / ' 9 = 1/ULP if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ', 00944 $ 'when reordering fails', / 00945 $ ' 10 = 1/ULP if PLEST/PLTRU > THRESH or ', 00946 $ 'PLTRU/PLEST > THRESH', / 00947 $ ' ( Test 10 is only for input examples )', / ) 00948 9992 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', D10.3, 00949 $ ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, F8.2 ) 00950 9991 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', D10.3, 00951 $ ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, D10.3 ) 00952 9990 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',', 00953 $ ' result ', I2, ' is', 0P, F8.2 ) 00954 9989 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',', 00955 $ ' result ', I2, ' is', 1P, D10.3 ) 00956 * 00957 * End of ZDRGSX 00958 * 00959 END