![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CDRGSX 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 CDRGSX( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, 00012 * AI, BI, Z, Q, ALPHA, BETA, C, LDC, S, WORK, 00013 * LWORK, RWORK, IWORK, LIWORK, BWORK, INFO ) 00014 * 00015 * .. Scalar Arguments .. 00016 * INTEGER INFO, LDA, LDC, LIWORK, LWORK, NCMAX, NIN, 00017 * $ NOUT, NSIZE 00018 * REAL THRESH 00019 * .. 00020 * .. Array Arguments .. 00021 * LOGICAL BWORK( * ) 00022 * INTEGER IWORK( * ) 00023 * REAL RWORK( * ), S( * ) 00024 * COMPLEX 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 *> CDRGSX checks the nonsymmetric generalized eigenvalue (Schur form) 00037 *> problem expert driver CGGESX. 00038 *> 00039 *> CGGES 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 CDRGSX is called with NSIZE > 0, five (5) types of built-in 00055 *> matrix pairs are used to test the routine CGGESX. 00056 *> 00057 *> When CDRGSX is called with NSIZE = 0, it reads in test matrix data 00058 *> to test CGGESX. 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 CGGESX, 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 CGESVD) 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 SLATM5). 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 SGGESX. 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 CLAKF2 00189 *> \endverbatim 00190 *> 00191 *> \param[in] THRESH 00192 *> \verbatim 00193 *> THRESH is REAL 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 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 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 array, dimension (LDA, NSIZE) 00241 *> Copy of A, modified by CGGESX. 00242 *> \endverbatim 00243 *> 00244 *> \param[out] BI 00245 *> \verbatim 00246 *> BI is COMPLEX array, dimension (LDA, NSIZE) 00247 *> Copy of B, modified by CGGESX. 00248 *> \endverbatim 00249 *> 00250 *> \param[out] Z 00251 *> \verbatim 00252 *> Z is COMPLEX array, dimension (LDA, NSIZE) 00253 *> Z holds the left Schur vectors computed by CGGESX. 00254 *> \endverbatim 00255 *> 00256 *> \param[out] Q 00257 *> \verbatim 00258 *> Q is COMPLEX array, dimension (LDA, NSIZE) 00259 *> Q holds the right Schur vectors computed by CGGESX. 00260 *> \endverbatim 00261 *> 00262 *> \param[out] ALPHA 00263 *> \verbatim 00264 *> ALPHA is COMPLEX array, dimension (NSIZE) 00265 *> \endverbatim 00266 *> 00267 *> \param[out] BETA 00268 *> \verbatim 00269 *> BETA is COMPLEX 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 array, dimension (LDC, LDC) 00277 *> Store the matrix generated by subroutine CLAKF2, 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 REAL array, dimension (LDC) 00291 *> Singular values of C 00292 *> \endverbatim 00293 *> 00294 *> \param[out] WORK 00295 *> \verbatim 00296 *> WORK is COMPLEX 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 REAL 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 complex_eig 00346 * 00347 * ===================================================================== 00348 SUBROUTINE CDRGSX( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, 00349 $ AI, BI, Z, Q, ALPHA, BETA, C, LDC, S, WORK, 00350 $ LWORK, 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 REAL THRESH 00361 * .. 00362 * .. Array Arguments .. 00363 LOGICAL BWORK( * ) 00364 INTEGER IWORK( * ) 00365 REAL RWORK( * ), S( * ) 00366 COMPLEX 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 REAL ZERO, ONE, TEN 00376 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TEN = 1.0E+1 ) 00377 COMPLEX CZERO 00378 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+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 REAL ABNRM, BIGNUM, DIFTRU, PLTRU, SMLNUM, TEMP1, 00387 $ TEMP2, THRSH2, ULP, ULPINV, WEIGHT 00388 COMPLEX X 00389 * .. 00390 * .. Local Arrays .. 00391 REAL DIFEST( 2 ), PL( 2 ), RESULT( 10 ) 00392 * .. 00393 * .. External Functions .. 00394 LOGICAL CLCTSX 00395 INTEGER ILAENV 00396 REAL CLANGE, SLAMCH 00397 EXTERNAL CLCTSX, ILAENV, CLANGE, SLAMCH 00398 * .. 00399 * .. External Subroutines .. 00400 EXTERNAL ALASVM, CGESVD, CGET51, CGGESX, CLACPY, CLAKF2, 00401 $ CLASET, CLATM5, SLABAD, XERBLA 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, AIMAG, MAX, REAL, SQRT 00412 * .. 00413 * .. Statement Functions .. 00414 REAL ABS1 00415 * .. 00416 * .. Statement Function definitions .. 00417 ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) ) 00418 * .. 00419 * .. Executable Statements .. 00420 * 00421 * Check for errors 00422 * 00423 IF( NSIZE.LT.0 ) THEN 00424 INFO = -1 00425 ELSE IF( THRESH.LT.ZERO ) THEN 00426 INFO = -2 00427 ELSE IF( NIN.LE.0 ) THEN 00428 INFO = -3 00429 ELSE IF( NOUT.LE.0 ) THEN 00430 INFO = -4 00431 ELSE IF( LDA.LT.1 .OR. LDA.LT.NSIZE ) THEN 00432 INFO = -6 00433 ELSE IF( LDC.LT.1 .OR. LDC.LT.NSIZE*NSIZE / 2 ) THEN 00434 INFO = -15 00435 ELSE IF( LIWORK.LT.NSIZE+2 ) THEN 00436 INFO = -21 00437 END IF 00438 * 00439 * Compute workspace 00440 * (Note: Comments in the code beginning "Workspace:" describe the 00441 * minimal amount of workspace needed at that point in the code, 00442 * as well as the preferred amount for good performance. 00443 * NB refers to the optimal block size for the immediately 00444 * following subroutine, as returned by ILAENV.) 00445 * 00446 MINWRK = 1 00447 IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN 00448 MINWRK = 3*NSIZE*NSIZE / 2 00449 * 00450 * workspace for cggesx 00451 * 00452 MAXWRK = NSIZE*( 1+ILAENV( 1, 'CGEQRF', ' ', NSIZE, 1, NSIZE, 00453 $ 0 ) ) 00454 MAXWRK = MAX( MAXWRK, NSIZE*( 1+ILAENV( 1, 'CUNGQR', ' ', 00455 $ NSIZE, 1, NSIZE, -1 ) ) ) 00456 * 00457 * workspace for cgesvd 00458 * 00459 BDSPAC = 3*NSIZE*NSIZE / 2 00460 MAXWRK = MAX( MAXWRK, NSIZE*NSIZE* 00461 $ ( 1+ILAENV( 1, 'CGEBRD', ' ', NSIZE*NSIZE / 2, 00462 $ NSIZE*NSIZE / 2, -1, -1 ) ) ) 00463 MAXWRK = MAX( MAXWRK, BDSPAC ) 00464 * 00465 MAXWRK = MAX( MAXWRK, MINWRK ) 00466 * 00467 WORK( 1 ) = MAXWRK 00468 END IF 00469 * 00470 IF( LWORK.LT.MINWRK ) 00471 $ INFO = -18 00472 * 00473 IF( INFO.NE.0 ) THEN 00474 CALL XERBLA( 'CDRGSX', -INFO ) 00475 RETURN 00476 END IF 00477 * 00478 * Important constants 00479 * 00480 ULP = SLAMCH( 'P' ) 00481 ULPINV = ONE / ULP 00482 SMLNUM = SLAMCH( 'S' ) / ULP 00483 BIGNUM = ONE / SMLNUM 00484 CALL SLABAD( SMLNUM, BIGNUM ) 00485 THRSH2 = TEN*THRESH 00486 NTESTT = 0 00487 NERRS = 0 00488 * 00489 * Go to the tests for read-in matrix pairs 00490 * 00491 IFUNC = 0 00492 IF( NSIZE.EQ.0 ) 00493 $ GO TO 70 00494 * 00495 * Test the built-in matrix pairs. 00496 * Loop over different functions (IFUNC) of CGGESX, types (PRTYPE) 00497 * of test matrices, different size (M+N) 00498 * 00499 PRTYPE = 0 00500 QBA = 3 00501 QBB = 4 00502 WEIGHT = SQRT( ULP ) 00503 * 00504 DO 60 IFUNC = 0, 3 00505 DO 50 PRTYPE = 1, 5 00506 DO 40 M = 1, NSIZE - 1 00507 DO 30 N = 1, NSIZE - M 00508 * 00509 WEIGHT = ONE / WEIGHT 00510 MPLUSN = M + N 00511 * 00512 * Generate test matrices 00513 * 00514 FS = .TRUE. 00515 K = 0 00516 * 00517 CALL CLASET( 'Full', MPLUSN, MPLUSN, CZERO, CZERO, AI, 00518 $ LDA ) 00519 CALL CLASET( 'Full', MPLUSN, MPLUSN, CZERO, CZERO, BI, 00520 $ LDA ) 00521 * 00522 CALL CLATM5( PRTYPE, M, N, AI, LDA, AI( M+1, M+1 ), 00523 $ LDA, AI( 1, M+1 ), LDA, BI, LDA, 00524 $ BI( M+1, M+1 ), LDA, BI( 1, M+1 ), LDA, 00525 $ Q, LDA, Z, LDA, WEIGHT, QBA, QBB ) 00526 * 00527 * Compute the Schur factorization and swapping the 00528 * m-by-m (1,1)-blocks with n-by-n (2,2)-blocks. 00529 * Swapping is accomplished via the function CLCTSX 00530 * which is supplied below. 00531 * 00532 IF( IFUNC.EQ.0 ) THEN 00533 SENSE = 'N' 00534 ELSE IF( IFUNC.EQ.1 ) THEN 00535 SENSE = 'E' 00536 ELSE IF( IFUNC.EQ.2 ) THEN 00537 SENSE = 'V' 00538 ELSE IF( IFUNC.EQ.3 ) THEN 00539 SENSE = 'B' 00540 END IF 00541 * 00542 CALL CLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA ) 00543 CALL CLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA ) 00544 * 00545 CALL CGGESX( 'V', 'V', 'S', CLCTSX, SENSE, MPLUSN, AI, 00546 $ LDA, BI, LDA, MM, ALPHA, BETA, Q, LDA, Z, 00547 $ LDA, PL, DIFEST, WORK, LWORK, RWORK, 00548 $ IWORK, LIWORK, BWORK, LINFO ) 00549 * 00550 IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN 00551 RESULT( 1 ) = ULPINV 00552 WRITE( NOUT, FMT = 9999 )'CGGESX', LINFO, MPLUSN, 00553 $ PRTYPE 00554 INFO = LINFO 00555 GO TO 30 00556 END IF 00557 * 00558 * Compute the norm(A, B) 00559 * 00560 CALL CLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK, 00561 $ MPLUSN ) 00562 CALL CLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, 00563 $ WORK( MPLUSN*MPLUSN+1 ), MPLUSN ) 00564 ABNRM = CLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN, 00565 $ RWORK ) 00566 * 00567 * Do tests (1) to (4) 00568 * 00569 RESULT( 2 ) = ZERO 00570 CALL CGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z, 00571 $ LDA, WORK, RWORK, RESULT( 1 ) ) 00572 CALL CGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z, 00573 $ LDA, WORK, RWORK, RESULT( 2 ) ) 00574 CALL CGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q, 00575 $ LDA, WORK, RWORK, RESULT( 3 ) ) 00576 CALL CGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z, 00577 $ LDA, WORK, RWORK, RESULT( 4 ) ) 00578 NTEST = 4 00579 * 00580 * Do tests (5) and (6): check Schur form of A and 00581 * compare eigenvalues with diagonals. 00582 * 00583 TEMP1 = ZERO 00584 RESULT( 5 ) = ZERO 00585 RESULT( 6 ) = ZERO 00586 * 00587 DO 10 J = 1, MPLUSN 00588 ILABAD = .FALSE. 00589 TEMP2 = ( ABS1( ALPHA( J )-AI( J, J ) ) / 00590 $ MAX( SMLNUM, ABS1( ALPHA( J ) ), 00591 $ ABS1( AI( J, J ) ) )+ 00592 $ ABS1( BETA( J )-BI( J, J ) ) / 00593 $ MAX( SMLNUM, ABS1( BETA( J ) ), 00594 $ ABS1( BI( J, J ) ) ) ) / ULP 00595 IF( J.LT.MPLUSN ) THEN 00596 IF( AI( J+1, J ).NE.ZERO ) THEN 00597 ILABAD = .TRUE. 00598 RESULT( 5 ) = ULPINV 00599 END IF 00600 END IF 00601 IF( J.GT.1 ) THEN 00602 IF( AI( J, J-1 ).NE.ZERO ) THEN 00603 ILABAD = .TRUE. 00604 RESULT( 5 ) = ULPINV 00605 END IF 00606 END IF 00607 TEMP1 = MAX( TEMP1, TEMP2 ) 00608 IF( ILABAD ) THEN 00609 WRITE( NOUT, FMT = 9997 )J, MPLUSN, PRTYPE 00610 END IF 00611 10 CONTINUE 00612 RESULT( 6 ) = TEMP1 00613 NTEST = NTEST + 2 00614 * 00615 * Test (7) (if sorting worked) 00616 * 00617 RESULT( 7 ) = ZERO 00618 IF( LINFO.EQ.MPLUSN+3 ) THEN 00619 RESULT( 7 ) = ULPINV 00620 ELSE IF( MM.NE.N ) THEN 00621 RESULT( 7 ) = ULPINV 00622 END IF 00623 NTEST = NTEST + 1 00624 * 00625 * Test (8): compare the estimated value DIF and its 00626 * value. first, compute the exact DIF. 00627 * 00628 RESULT( 8 ) = ZERO 00629 MN2 = MM*( MPLUSN-MM )*2 00630 IF( IFUNC.GE.2 .AND. MN2.LE.NCMAX*NCMAX ) THEN 00631 * 00632 * Note: for either following two cases, there are 00633 * almost same number of test cases fail the test. 00634 * 00635 CALL CLAKF2( MM, MPLUSN-MM, AI, LDA, 00636 $ AI( MM+1, MM+1 ), BI, 00637 $ BI( MM+1, MM+1 ), C, LDC ) 00638 * 00639 CALL CGESVD( 'N', 'N', MN2, MN2, C, LDC, S, WORK, 00640 $ 1, WORK( 2 ), 1, WORK( 3 ), LWORK-2, 00641 $ RWORK, INFO ) 00642 DIFTRU = S( MN2 ) 00643 * 00644 IF( DIFEST( 2 ).EQ.ZERO ) THEN 00645 IF( DIFTRU.GT.ABNRM*ULP ) 00646 $ RESULT( 8 ) = ULPINV 00647 ELSE IF( DIFTRU.EQ.ZERO ) THEN 00648 IF( DIFEST( 2 ).GT.ABNRM*ULP ) 00649 $ RESULT( 8 ) = ULPINV 00650 ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR. 00651 $ ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN 00652 RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ), 00653 $ DIFEST( 2 ) / DIFTRU ) 00654 END IF 00655 NTEST = NTEST + 1 00656 END IF 00657 * 00658 * Test (9) 00659 * 00660 RESULT( 9 ) = ZERO 00661 IF( LINFO.EQ.( MPLUSN+2 ) ) THEN 00662 IF( DIFTRU.GT.ABNRM*ULP ) 00663 $ RESULT( 9 ) = ULPINV 00664 IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) ) 00665 $ RESULT( 9 ) = ULPINV 00666 IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) ) 00667 $ RESULT( 9 ) = ULPINV 00668 NTEST = NTEST + 1 00669 END IF 00670 * 00671 NTESTT = NTESTT + NTEST 00672 * 00673 * Print out tests which fail. 00674 * 00675 DO 20 J = 1, 9 00676 IF( RESULT( J ).GE.THRESH ) THEN 00677 * 00678 * If this is the first test to fail, 00679 * print a header to the data file. 00680 * 00681 IF( NERRS.EQ.0 ) THEN 00682 WRITE( NOUT, FMT = 9996 )'CGX' 00683 * 00684 * Matrix types 00685 * 00686 WRITE( NOUT, FMT = 9994 ) 00687 * 00688 * Tests performed 00689 * 00690 WRITE( NOUT, FMT = 9993 )'unitary', '''', 00691 $ 'transpose', ( '''', I = 1, 4 ) 00692 * 00693 END IF 00694 NERRS = NERRS + 1 00695 IF( RESULT( J ).LT.10000.0 ) THEN 00696 WRITE( NOUT, FMT = 9992 )MPLUSN, PRTYPE, 00697 $ WEIGHT, M, J, RESULT( J ) 00698 ELSE 00699 WRITE( NOUT, FMT = 9991 )MPLUSN, PRTYPE, 00700 $ WEIGHT, M, J, RESULT( J ) 00701 END IF 00702 END IF 00703 20 CONTINUE 00704 * 00705 30 CONTINUE 00706 40 CONTINUE 00707 50 CONTINUE 00708 60 CONTINUE 00709 * 00710 GO TO 150 00711 * 00712 70 CONTINUE 00713 * 00714 * Read in data from file to check accuracy of condition estimation 00715 * Read input data until N=0 00716 * 00717 NPTKNT = 0 00718 * 00719 80 CONTINUE 00720 READ( NIN, FMT = *, END = 140 )MPLUSN 00721 IF( MPLUSN.EQ.0 ) 00722 $ GO TO 140 00723 READ( NIN, FMT = *, END = 140 )N 00724 DO 90 I = 1, MPLUSN 00725 READ( NIN, FMT = * )( AI( I, J ), J = 1, MPLUSN ) 00726 90 CONTINUE 00727 DO 100 I = 1, MPLUSN 00728 READ( NIN, FMT = * )( BI( I, J ), J = 1, MPLUSN ) 00729 100 CONTINUE 00730 READ( NIN, FMT = * )PLTRU, DIFTRU 00731 * 00732 NPTKNT = NPTKNT + 1 00733 FS = .TRUE. 00734 K = 0 00735 M = MPLUSN - N 00736 * 00737 CALL CLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA ) 00738 CALL CLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA ) 00739 * 00740 * Compute the Schur factorization while swaping the 00741 * m-by-m (1,1)-blocks with n-by-n (2,2)-blocks. 00742 * 00743 CALL CGGESX( 'V', 'V', 'S', CLCTSX, 'B', MPLUSN, AI, LDA, BI, LDA, 00744 $ MM, ALPHA, BETA, Q, LDA, Z, LDA, PL, DIFEST, WORK, 00745 $ LWORK, RWORK, IWORK, LIWORK, BWORK, LINFO ) 00746 * 00747 IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN 00748 RESULT( 1 ) = ULPINV 00749 WRITE( NOUT, FMT = 9998 )'CGGESX', LINFO, MPLUSN, NPTKNT 00750 GO TO 130 00751 END IF 00752 * 00753 * Compute the norm(A, B) 00754 * (should this be norm of (A,B) or (AI,BI)?) 00755 * 00756 CALL CLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK, MPLUSN ) 00757 CALL CLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, 00758 $ WORK( MPLUSN*MPLUSN+1 ), MPLUSN ) 00759 ABNRM = CLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN, RWORK ) 00760 * 00761 * Do tests (1) to (4) 00762 * 00763 CALL CGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z, LDA, WORK, 00764 $ RWORK, RESULT( 1 ) ) 00765 CALL CGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z, LDA, WORK, 00766 $ RWORK, RESULT( 2 ) ) 00767 CALL CGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q, LDA, WORK, 00768 $ RWORK, RESULT( 3 ) ) 00769 CALL CGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z, LDA, WORK, 00770 $ RWORK, RESULT( 4 ) ) 00771 * 00772 * Do tests (5) and (6): check Schur form of A and compare 00773 * eigenvalues with diagonals. 00774 * 00775 NTEST = 6 00776 TEMP1 = ZERO 00777 RESULT( 5 ) = ZERO 00778 RESULT( 6 ) = ZERO 00779 * 00780 DO 110 J = 1, MPLUSN 00781 ILABAD = .FALSE. 00782 TEMP2 = ( ABS1( ALPHA( J )-AI( J, J ) ) / 00783 $ MAX( SMLNUM, ABS1( ALPHA( J ) ), ABS1( AI( J, J ) ) )+ 00784 $ ABS1( BETA( J )-BI( J, J ) ) / 00785 $ MAX( SMLNUM, ABS1( BETA( J ) ), ABS1( BI( J, J ) ) ) ) 00786 $ / ULP 00787 IF( J.LT.MPLUSN ) THEN 00788 IF( AI( J+1, J ).NE.ZERO ) THEN 00789 ILABAD = .TRUE. 00790 RESULT( 5 ) = ULPINV 00791 END IF 00792 END IF 00793 IF( J.GT.1 ) THEN 00794 IF( AI( J, J-1 ).NE.ZERO ) THEN 00795 ILABAD = .TRUE. 00796 RESULT( 5 ) = ULPINV 00797 END IF 00798 END IF 00799 TEMP1 = MAX( TEMP1, TEMP2 ) 00800 IF( ILABAD ) THEN 00801 WRITE( NOUT, FMT = 9997 )J, MPLUSN, NPTKNT 00802 END IF 00803 110 CONTINUE 00804 RESULT( 6 ) = TEMP1 00805 * 00806 * Test (7) (if sorting worked) <--------- need to be checked. 00807 * 00808 NTEST = 7 00809 RESULT( 7 ) = ZERO 00810 IF( LINFO.EQ.MPLUSN+3 ) 00811 $ RESULT( 7 ) = ULPINV 00812 * 00813 * Test (8): compare the estimated value of DIF and its true value. 00814 * 00815 NTEST = 8 00816 RESULT( 8 ) = ZERO 00817 IF( DIFEST( 2 ).EQ.ZERO ) THEN 00818 IF( DIFTRU.GT.ABNRM*ULP ) 00819 $ RESULT( 8 ) = ULPINV 00820 ELSE IF( DIFTRU.EQ.ZERO ) THEN 00821 IF( DIFEST( 2 ).GT.ABNRM*ULP ) 00822 $ RESULT( 8 ) = ULPINV 00823 ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR. 00824 $ ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN 00825 RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ), DIFEST( 2 ) / DIFTRU ) 00826 END IF 00827 * 00828 * Test (9) 00829 * 00830 NTEST = 9 00831 RESULT( 9 ) = ZERO 00832 IF( LINFO.EQ.( MPLUSN+2 ) ) THEN 00833 IF( DIFTRU.GT.ABNRM*ULP ) 00834 $ RESULT( 9 ) = ULPINV 00835 IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) ) 00836 $ RESULT( 9 ) = ULPINV 00837 IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) ) 00838 $ RESULT( 9 ) = ULPINV 00839 END IF 00840 * 00841 * Test (10): compare the estimated value of PL and it true value. 00842 * 00843 NTEST = 10 00844 RESULT( 10 ) = ZERO 00845 IF( PL( 1 ).EQ.ZERO ) THEN 00846 IF( PLTRU.GT.ABNRM*ULP ) 00847 $ RESULT( 10 ) = ULPINV 00848 ELSE IF( PLTRU.EQ.ZERO ) THEN 00849 IF( PL( 1 ).GT.ABNRM*ULP ) 00850 $ RESULT( 10 ) = ULPINV 00851 ELSE IF( ( PLTRU.GT.THRESH*PL( 1 ) ) .OR. 00852 $ ( PLTRU*THRESH.LT.PL( 1 ) ) ) THEN 00853 RESULT( 10 ) = ULPINV 00854 END IF 00855 * 00856 NTESTT = NTESTT + NTEST 00857 * 00858 * Print out tests which fail. 00859 * 00860 DO 120 J = 1, NTEST 00861 IF( RESULT( J ).GE.THRESH ) THEN 00862 * 00863 * If this is the first test to fail, 00864 * print a header to the data file. 00865 * 00866 IF( NERRS.EQ.0 ) THEN 00867 WRITE( NOUT, FMT = 9996 )'CGX' 00868 * 00869 * Matrix types 00870 * 00871 WRITE( NOUT, FMT = 9995 ) 00872 * 00873 * Tests performed 00874 * 00875 WRITE( NOUT, FMT = 9993 )'unitary', '''', 'transpose', 00876 $ ( '''', I = 1, 4 ) 00877 * 00878 END IF 00879 NERRS = NERRS + 1 00880 IF( RESULT( J ).LT.10000.0 ) THEN 00881 WRITE( NOUT, FMT = 9990 )NPTKNT, MPLUSN, J, RESULT( J ) 00882 ELSE 00883 WRITE( NOUT, FMT = 9989 )NPTKNT, MPLUSN, J, RESULT( J ) 00884 END IF 00885 END IF 00886 * 00887 120 CONTINUE 00888 * 00889 130 CONTINUE 00890 GO TO 80 00891 140 CONTINUE 00892 * 00893 150 CONTINUE 00894 * 00895 * Summary 00896 * 00897 CALL ALASVM( 'CGX', NOUT, NERRS, NTESTT, 0 ) 00898 * 00899 WORK( 1 ) = MAXWRK 00900 * 00901 RETURN 00902 * 00903 9999 FORMAT( ' CDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 00904 $ I6, ', JTYPE=', I6, ')' ) 00905 * 00906 9998 FORMAT( ' CDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 00907 $ I6, ', Input Example #', I2, ')' ) 00908 * 00909 9997 FORMAT( ' CDRGSX: S not in Schur form at eigenvalue ', I6, '.', 00910 $ / 9X, 'N=', I6, ', JTYPE=', I6, ')' ) 00911 * 00912 9996 FORMAT( / 1X, A3, ' -- Complex Expert Generalized Schur form', 00913 $ ' problem driver' ) 00914 * 00915 9995 FORMAT( 'Input Example' ) 00916 * 00917 9994 FORMAT( ' Matrix types: ', / 00918 $ ' 1: A is a block diagonal matrix of Jordan blocks ', 00919 $ 'and B is the identity ', / ' matrix, ', 00920 $ / ' 2: A and B are upper triangular matrices, ', 00921 $ / ' 3: A and B are as type 2, but each second diagonal ', 00922 $ 'block in A_11 and ', / 00923 $ ' each third diaongal block in A_22 are 2x2 blocks,', 00924 $ / ' 4: A and B are block diagonal matrices, ', 00925 $ / ' 5: (A,B) has potentially close or common ', 00926 $ 'eigenvalues.', / ) 00927 * 00928 9993 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ', 00929 $ 'Q and Z are ', A, ',', / 19X, 00930 $ ' a is alpha, b is beta, and ', A, ' means ', A, '.)', 00931 $ / ' 1 = | A - Q S Z', A, 00932 $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', A, 00933 $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', A, 00934 $ ' | / ( n ulp ) 4 = | I - ZZ', A, 00935 $ ' | / ( n ulp )', / ' 5 = 1/ULP if A is not in ', 00936 $ 'Schur form S', / ' 6 = difference between (alpha,beta)', 00937 $ ' and diagonals of (S,T)', / 00938 $ ' 7 = 1/ULP if SDIM is not the correct number of ', 00939 $ 'selected eigenvalues', / 00940 $ ' 8 = 1/ULP if DIFEST/DIFTRU > 10*THRESH or ', 00941 $ 'DIFTRU/DIFEST > 10*THRESH', 00942 $ / ' 9 = 1/ULP if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ', 00943 $ 'when reordering fails', / 00944 $ ' 10 = 1/ULP if PLEST/PLTRU > THRESH or ', 00945 $ 'PLTRU/PLEST > THRESH', / 00946 $ ' ( Test 10 is only for input examples )', / ) 00947 9992 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', E10.3, 00948 $ ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, F8.2 ) 00949 9991 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', E10.3, 00950 $ ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, E10.3 ) 00951 9990 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',', 00952 $ ' result ', I2, ' is', 0P, F8.2 ) 00953 9989 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',', 00954 $ ' result ', I2, ' is', 1P, E10.3 ) 00955 * 00956 * End of CDRGSX 00957 * 00958 END