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