![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SDRGSX 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 SDRGSX( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, 00012 * AI, 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 * REAL THRESH 00019 * .. 00020 * .. Array Arguments .. 00021 * LOGICAL BWORK( * ) 00022 * INTEGER IWORK( * ) 00023 * REAL 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 *> SDRGSX checks the nonsymmetric generalized eigenvalue (Schur form) 00036 *> problem expert driver SGGESX. 00037 *> 00038 *> SGGESX 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 SDRGSX is called with NSIZE > 0, five (5) types of built-in 00057 *> matrix pairs are used to test the routine SGGESX. 00058 *> 00059 *> When SDRGSX is called with NSIZE = 0, it reads in test matrix data 00060 *> to test SGGESX. 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 SGGESX, 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 SGESVD) 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 SLATM5). 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 SGGESX. 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 SLAKF2 00199 *> \endverbatim 00200 *> 00201 *> \param[in] THRESH 00202 *> \verbatim 00203 *> THRESH is REAL 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 REAL 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 REAL 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 REAL array, dimension (LDA, NSIZE) 00251 *> Copy of A, modified by SGGESX. 00252 *> \endverbatim 00253 *> 00254 *> \param[out] BI 00255 *> \verbatim 00256 *> BI is REAL array, dimension (LDA, NSIZE) 00257 *> Copy of B, modified by SGGESX. 00258 *> \endverbatim 00259 *> 00260 *> \param[out] Z 00261 *> \verbatim 00262 *> Z is REAL array, dimension (LDA, NSIZE) 00263 *> Z holds the left Schur vectors computed by SGGESX. 00264 *> \endverbatim 00265 *> 00266 *> \param[out] Q 00267 *> \verbatim 00268 *> Q is REAL array, dimension (LDA, NSIZE) 00269 *> Q holds the right Schur vectors computed by SGGESX. 00270 *> \endverbatim 00271 *> 00272 *> \param[out] ALPHAR 00273 *> \verbatim 00274 *> ALPHAR is REAL array, dimension (NSIZE) 00275 *> \endverbatim 00276 *> 00277 *> \param[out] ALPHAI 00278 *> \verbatim 00279 *> ALPHAI is REAL array, dimension (NSIZE) 00280 *> \endverbatim 00281 *> 00282 *> \param[out] BETA 00283 *> \verbatim 00284 *> BETA is REAL array, dimension (NSIZE) 00285 *> \verbatim 00286 *> On exit, (ALPHAR + ALPHAI*i)/BETA are the eigenvalues. 00287 *> \endverbatim 00288 *> 00289 *> \param[out] C 00290 *> \verbatim 00291 *> C is REAL array, dimension (LDC, LDC) 00292 *> Store the matrix generated by subroutine SLAKF2, 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 REAL array, dimension (LDC) 00306 *> Singular values of C 00307 *> \endverbatim 00308 *> 00309 *> \param[out] WORK 00310 *> \verbatim 00311 *> WORK is REAL 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 single_eig 00356 * 00357 * ===================================================================== 00358 SUBROUTINE SDRGSX( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, 00359 $ AI, 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 REAL THRESH 00371 * .. 00372 * .. Array Arguments .. 00373 LOGICAL BWORK( * ) 00374 INTEGER IWORK( * ) 00375 REAL 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 REAL ZERO, ONE, TEN 00385 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TEN = 1.0E+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 REAL ABNRM, BIGNUM, DIFTRU, PLTRU, SMLNUM, TEMP1, 00394 $ TEMP2, THRSH2, ULP, ULPINV, WEIGHT 00395 * .. 00396 * .. Local Arrays .. 00397 REAL DIFEST( 2 ), PL( 2 ), RESULT( 10 ) 00398 * .. 00399 * .. External Functions .. 00400 LOGICAL SLCTSX 00401 INTEGER ILAENV 00402 REAL SLAMCH, SLANGE 00403 EXTERNAL SLCTSX, ILAENV, SLAMCH, SLANGE 00404 * .. 00405 * .. External Subroutines .. 00406 EXTERNAL ALASVM, SGESVD, SGET51, SGET53, SGGESX, SLABAD, 00407 $ SLACPY, SLAKF2, SLASET, SLATM5, 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 c MINWRK = MAX( 10*( NSIZE+1 ), 5*NSIZE*NSIZE / 2-2 ) 00449 MINWRK = MAX( 10*( NSIZE+1 ), 5*NSIZE*NSIZE / 2 ) 00450 * 00451 * workspace for sggesx 00452 * 00453 MAXWRK = 9*( NSIZE+1 ) + NSIZE* 00454 $ ILAENV( 1, 'SGEQRF', ' ', NSIZE, 1, NSIZE, 0 ) 00455 MAXWRK = MAX( MAXWRK, 9*( NSIZE+1 )+NSIZE* 00456 $ ILAENV( 1, 'SORGQR', ' ', NSIZE, 1, NSIZE, -1 ) ) 00457 * 00458 * workspace for sgesvd 00459 * 00460 BDSPAC = 5*NSIZE*NSIZE / 2 00461 MAXWRK = MAX( MAXWRK, 3*NSIZE*NSIZE / 2+NSIZE*NSIZE* 00462 $ ILAENV( 1, 'SGEBRD', ' ', 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 = -19 00473 * 00474 IF( INFO.NE.0 ) THEN 00475 CALL XERBLA( 'SDRGSX', -INFO ) 00476 RETURN 00477 END IF 00478 * 00479 * Important constants 00480 * 00481 ULP = SLAMCH( 'P' ) 00482 ULPINV = ONE / ULP 00483 SMLNUM = SLAMCH( 'S' ) / ULP 00484 BIGNUM = ONE / SMLNUM 00485 CALL SLABAD( 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 SGGESX, 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 SLASET( 'Full', MPLUSN, MPLUSN, ZERO, ZERO, AI, 00519 $ LDA ) 00520 CALL SLASET( 'Full', MPLUSN, MPLUSN, ZERO, ZERO, BI, 00521 $ LDA ) 00522 * 00523 CALL SLATM5( 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 SLCTSX 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 SLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA ) 00544 CALL SLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA ) 00545 * 00546 CALL SGGESX( 'V', 'V', 'S', SLCTSX, SENSE, MPLUSN, AI, 00547 $ LDA, BI, LDA, MM, ALPHAR, ALPHAI, BETA, 00548 $ Q, LDA, Z, LDA, PL, DIFEST, WORK, LWORK, 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 )'SGGESX', LINFO, MPLUSN, 00554 $ PRTYPE 00555 INFO = LINFO 00556 GO TO 30 00557 END IF 00558 * 00559 * Compute the norm(A, B) 00560 * 00561 CALL SLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK, 00562 $ MPLUSN ) 00563 CALL SLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, 00564 $ WORK( MPLUSN*MPLUSN+1 ), MPLUSN ) 00565 ABNRM = SLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN, 00566 $ WORK ) 00567 * 00568 * Do tests (1) to (4) 00569 * 00570 CALL SGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z, 00571 $ LDA, WORK, RESULT( 1 ) ) 00572 CALL SGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z, 00573 $ LDA, WORK, RESULT( 2 ) ) 00574 CALL SGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q, 00575 $ LDA, WORK, RESULT( 3 ) ) 00576 CALL SGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z, 00577 $ LDA, WORK, 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 IF( ALPHAI( J ).EQ.ZERO ) THEN 00590 TEMP2 = ( ABS( ALPHAR( J )-AI( J, J ) ) / 00591 $ MAX( SMLNUM, ABS( ALPHAR( J ) ), 00592 $ ABS( AI( J, J ) ) )+ 00593 $ ABS( BETA( J )-BI( J, J ) ) / 00594 $ MAX( SMLNUM, ABS( BETA( J ) ), 00595 $ ABS( 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 ELSE 00609 IF( ALPHAI( J ).GT.ZERO ) THEN 00610 I1 = J 00611 ELSE 00612 I1 = J - 1 00613 END IF 00614 IF( I1.LE.0 .OR. I1.GE.MPLUSN ) THEN 00615 ILABAD = .TRUE. 00616 ELSE IF( I1.LT.MPLUSN-1 ) THEN 00617 IF( AI( I1+2, I1+1 ).NE.ZERO ) THEN 00618 ILABAD = .TRUE. 00619 RESULT( 5 ) = ULPINV 00620 END IF 00621 ELSE IF( I1.GT.1 ) THEN 00622 IF( AI( I1, I1-1 ).NE.ZERO ) THEN 00623 ILABAD = .TRUE. 00624 RESULT( 5 ) = ULPINV 00625 END IF 00626 END IF 00627 IF( .NOT.ILABAD ) THEN 00628 CALL SGET53( AI( I1, I1 ), LDA, BI( I1, I1 ), 00629 $ LDA, BETA( J ), ALPHAR( J ), 00630 $ ALPHAI( J ), TEMP2, IINFO ) 00631 IF( IINFO.GE.3 ) THEN 00632 WRITE( NOUT, FMT = 9997 )IINFO, J, 00633 $ MPLUSN, PRTYPE 00634 INFO = ABS( IINFO ) 00635 END IF 00636 ELSE 00637 TEMP2 = ULPINV 00638 END IF 00639 END IF 00640 TEMP1 = MAX( TEMP1, TEMP2 ) 00641 IF( ILABAD ) THEN 00642 WRITE( NOUT, FMT = 9996 )J, MPLUSN, PRTYPE 00643 END IF 00644 10 CONTINUE 00645 RESULT( 6 ) = TEMP1 00646 NTEST = NTEST + 2 00647 * 00648 * Test (7) (if sorting worked) 00649 * 00650 RESULT( 7 ) = ZERO 00651 IF( LINFO.EQ.MPLUSN+3 ) THEN 00652 RESULT( 7 ) = ULPINV 00653 ELSE IF( MM.NE.N ) THEN 00654 RESULT( 7 ) = ULPINV 00655 END IF 00656 NTEST = NTEST + 1 00657 * 00658 * Test (8): compare the estimated value DIF and its 00659 * value. first, compute the exact DIF. 00660 * 00661 RESULT( 8 ) = ZERO 00662 MN2 = MM*( MPLUSN-MM )*2 00663 IF( IFUNC.GE.2 .AND. MN2.LE.NCMAX*NCMAX ) THEN 00664 * 00665 * Note: for either following two causes, there are 00666 * almost same number of test cases fail the test. 00667 * 00668 CALL SLAKF2( MM, MPLUSN-MM, AI, LDA, 00669 $ AI( MM+1, MM+1 ), BI, 00670 $ BI( MM+1, MM+1 ), C, LDC ) 00671 * 00672 CALL SGESVD( 'N', 'N', MN2, MN2, C, LDC, S, WORK, 00673 $ 1, WORK( 2 ), 1, WORK( 3 ), LWORK-2, 00674 $ INFO ) 00675 DIFTRU = S( MN2 ) 00676 * 00677 IF( DIFEST( 2 ).EQ.ZERO ) THEN 00678 IF( DIFTRU.GT.ABNRM*ULP ) 00679 $ RESULT( 8 ) = ULPINV 00680 ELSE IF( DIFTRU.EQ.ZERO ) THEN 00681 IF( DIFEST( 2 ).GT.ABNRM*ULP ) 00682 $ RESULT( 8 ) = ULPINV 00683 ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR. 00684 $ ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN 00685 RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ), 00686 $ DIFEST( 2 ) / DIFTRU ) 00687 END IF 00688 NTEST = NTEST + 1 00689 END IF 00690 * 00691 * Test (9) 00692 * 00693 RESULT( 9 ) = ZERO 00694 IF( LINFO.EQ.( MPLUSN+2 ) ) THEN 00695 IF( DIFTRU.GT.ABNRM*ULP ) 00696 $ RESULT( 9 ) = ULPINV 00697 IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) ) 00698 $ RESULT( 9 ) = ULPINV 00699 IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) ) 00700 $ RESULT( 9 ) = ULPINV 00701 NTEST = NTEST + 1 00702 END IF 00703 * 00704 NTESTT = NTESTT + NTEST 00705 * 00706 * Print out tests which fail. 00707 * 00708 DO 20 J = 1, 9 00709 IF( RESULT( J ).GE.THRESH ) THEN 00710 * 00711 * If this is the first test to fail, 00712 * print a header to the data file. 00713 * 00714 IF( NERRS.EQ.0 ) THEN 00715 WRITE( NOUT, FMT = 9995 )'SGX' 00716 * 00717 * Matrix types 00718 * 00719 WRITE( NOUT, FMT = 9993 ) 00720 * 00721 * Tests performed 00722 * 00723 WRITE( NOUT, FMT = 9992 )'orthogonal', '''', 00724 $ 'transpose', ( '''', I = 1, 4 ) 00725 * 00726 END IF 00727 NERRS = NERRS + 1 00728 IF( RESULT( J ).LT.10000.0 ) THEN 00729 WRITE( NOUT, FMT = 9991 )MPLUSN, PRTYPE, 00730 $ WEIGHT, M, J, RESULT( J ) 00731 ELSE 00732 WRITE( NOUT, FMT = 9990 )MPLUSN, PRTYPE, 00733 $ WEIGHT, M, J, RESULT( J ) 00734 END IF 00735 END IF 00736 20 CONTINUE 00737 * 00738 30 CONTINUE 00739 40 CONTINUE 00740 50 CONTINUE 00741 60 CONTINUE 00742 * 00743 GO TO 150 00744 * 00745 70 CONTINUE 00746 * 00747 * Read in data from file to check accuracy of condition estimation 00748 * Read input data until N=0 00749 * 00750 NPTKNT = 0 00751 * 00752 80 CONTINUE 00753 READ( NIN, FMT = *, END = 140 )MPLUSN 00754 IF( MPLUSN.EQ.0 ) 00755 $ GO TO 140 00756 READ( NIN, FMT = *, END = 140 )N 00757 DO 90 I = 1, MPLUSN 00758 READ( NIN, FMT = * )( AI( I, J ), J = 1, MPLUSN ) 00759 90 CONTINUE 00760 DO 100 I = 1, MPLUSN 00761 READ( NIN, FMT = * )( BI( I, J ), J = 1, MPLUSN ) 00762 100 CONTINUE 00763 READ( NIN, FMT = * )PLTRU, DIFTRU 00764 * 00765 NPTKNT = NPTKNT + 1 00766 FS = .TRUE. 00767 K = 0 00768 M = MPLUSN - N 00769 * 00770 CALL SLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA ) 00771 CALL SLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA ) 00772 * 00773 * Compute the Schur factorization while swaping the 00774 * m-by-m (1,1)-blocks with n-by-n (2,2)-blocks. 00775 * 00776 CALL SGGESX( 'V', 'V', 'S', SLCTSX, 'B', MPLUSN, AI, LDA, BI, LDA, 00777 $ MM, ALPHAR, ALPHAI, BETA, Q, LDA, Z, LDA, PL, DIFEST, 00778 $ WORK, LWORK, IWORK, LIWORK, BWORK, LINFO ) 00779 * 00780 IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN 00781 RESULT( 1 ) = ULPINV 00782 WRITE( NOUT, FMT = 9998 )'SGGESX', LINFO, MPLUSN, NPTKNT 00783 GO TO 130 00784 END IF 00785 * 00786 * Compute the norm(A, B) 00787 * (should this be norm of (A,B) or (AI,BI)?) 00788 * 00789 CALL SLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK, MPLUSN ) 00790 CALL SLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, 00791 $ WORK( MPLUSN*MPLUSN+1 ), MPLUSN ) 00792 ABNRM = SLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN, WORK ) 00793 * 00794 * Do tests (1) to (4) 00795 * 00796 CALL SGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z, LDA, WORK, 00797 $ RESULT( 1 ) ) 00798 CALL SGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z, LDA, WORK, 00799 $ RESULT( 2 ) ) 00800 CALL SGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q, LDA, WORK, 00801 $ RESULT( 3 ) ) 00802 CALL SGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z, LDA, WORK, 00803 $ RESULT( 4 ) ) 00804 * 00805 * Do tests (5) and (6): check Schur form of A and compare 00806 * eigenvalues with diagonals. 00807 * 00808 NTEST = 6 00809 TEMP1 = ZERO 00810 RESULT( 5 ) = ZERO 00811 RESULT( 6 ) = ZERO 00812 * 00813 DO 110 J = 1, MPLUSN 00814 ILABAD = .FALSE. 00815 IF( ALPHAI( J ).EQ.ZERO ) THEN 00816 TEMP2 = ( ABS( ALPHAR( J )-AI( J, J ) ) / 00817 $ MAX( SMLNUM, ABS( ALPHAR( J ) ), ABS( AI( J, 00818 $ J ) ) )+ABS( BETA( J )-BI( J, J ) ) / 00819 $ MAX( SMLNUM, ABS( BETA( J ) ), ABS( BI( J, J ) ) ) ) 00820 $ / ULP 00821 IF( J.LT.MPLUSN ) THEN 00822 IF( AI( J+1, J ).NE.ZERO ) THEN 00823 ILABAD = .TRUE. 00824 RESULT( 5 ) = ULPINV 00825 END IF 00826 END IF 00827 IF( J.GT.1 ) THEN 00828 IF( AI( J, J-1 ).NE.ZERO ) THEN 00829 ILABAD = .TRUE. 00830 RESULT( 5 ) = ULPINV 00831 END IF 00832 END IF 00833 ELSE 00834 IF( ALPHAI( J ).GT.ZERO ) THEN 00835 I1 = J 00836 ELSE 00837 I1 = J - 1 00838 END IF 00839 IF( I1.LE.0 .OR. I1.GE.MPLUSN ) THEN 00840 ILABAD = .TRUE. 00841 ELSE IF( I1.LT.MPLUSN-1 ) THEN 00842 IF( AI( I1+2, I1+1 ).NE.ZERO ) THEN 00843 ILABAD = .TRUE. 00844 RESULT( 5 ) = ULPINV 00845 END IF 00846 ELSE IF( I1.GT.1 ) THEN 00847 IF( AI( I1, I1-1 ).NE.ZERO ) THEN 00848 ILABAD = .TRUE. 00849 RESULT( 5 ) = ULPINV 00850 END IF 00851 END IF 00852 IF( .NOT.ILABAD ) THEN 00853 CALL SGET53( AI( I1, I1 ), LDA, BI( I1, I1 ), LDA, 00854 $ BETA( J ), ALPHAR( J ), ALPHAI( J ), TEMP2, 00855 $ IINFO ) 00856 IF( IINFO.GE.3 ) THEN 00857 WRITE( NOUT, FMT = 9997 )IINFO, J, MPLUSN, NPTKNT 00858 INFO = ABS( IINFO ) 00859 END IF 00860 ELSE 00861 TEMP2 = ULPINV 00862 END IF 00863 END IF 00864 TEMP1 = MAX( TEMP1, TEMP2 ) 00865 IF( ILABAD ) THEN 00866 WRITE( NOUT, FMT = 9996 )J, MPLUSN, NPTKNT 00867 END IF 00868 110 CONTINUE 00869 RESULT( 6 ) = TEMP1 00870 * 00871 * Test (7) (if sorting worked) <--------- need to be checked. 00872 * 00873 NTEST = 7 00874 RESULT( 7 ) = ZERO 00875 IF( LINFO.EQ.MPLUSN+3 ) 00876 $ RESULT( 7 ) = ULPINV 00877 * 00878 * Test (8): compare the estimated value of DIF and its true value. 00879 * 00880 NTEST = 8 00881 RESULT( 8 ) = ZERO 00882 IF( DIFEST( 2 ).EQ.ZERO ) THEN 00883 IF( DIFTRU.GT.ABNRM*ULP ) 00884 $ RESULT( 8 ) = ULPINV 00885 ELSE IF( DIFTRU.EQ.ZERO ) THEN 00886 IF( DIFEST( 2 ).GT.ABNRM*ULP ) 00887 $ RESULT( 8 ) = ULPINV 00888 ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR. 00889 $ ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN 00890 RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ), DIFEST( 2 ) / DIFTRU ) 00891 END IF 00892 * 00893 * Test (9) 00894 * 00895 NTEST = 9 00896 RESULT( 9 ) = ZERO 00897 IF( LINFO.EQ.( MPLUSN+2 ) ) THEN 00898 IF( DIFTRU.GT.ABNRM*ULP ) 00899 $ RESULT( 9 ) = ULPINV 00900 IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) ) 00901 $ RESULT( 9 ) = ULPINV 00902 IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) ) 00903 $ RESULT( 9 ) = ULPINV 00904 END IF 00905 * 00906 * Test (10): compare the estimated value of PL and it true value. 00907 * 00908 NTEST = 10 00909 RESULT( 10 ) = ZERO 00910 IF( PL( 1 ).EQ.ZERO ) THEN 00911 IF( PLTRU.GT.ABNRM*ULP ) 00912 $ RESULT( 10 ) = ULPINV 00913 ELSE IF( PLTRU.EQ.ZERO ) THEN 00914 IF( PL( 1 ).GT.ABNRM*ULP ) 00915 $ RESULT( 10 ) = ULPINV 00916 ELSE IF( ( PLTRU.GT.THRESH*PL( 1 ) ) .OR. 00917 $ ( PLTRU*THRESH.LT.PL( 1 ) ) ) THEN 00918 RESULT( 10 ) = ULPINV 00919 END IF 00920 * 00921 NTESTT = NTESTT + NTEST 00922 * 00923 * Print out tests which fail. 00924 * 00925 DO 120 J = 1, NTEST 00926 IF( RESULT( J ).GE.THRESH ) THEN 00927 * 00928 * If this is the first test to fail, 00929 * print a header to the data file. 00930 * 00931 IF( NERRS.EQ.0 ) THEN 00932 WRITE( NOUT, FMT = 9995 )'SGX' 00933 * 00934 * Matrix types 00935 * 00936 WRITE( NOUT, FMT = 9994 ) 00937 * 00938 * Tests performed 00939 * 00940 WRITE( NOUT, FMT = 9992 )'orthogonal', '''', 00941 $ 'transpose', ( '''', I = 1, 4 ) 00942 * 00943 END IF 00944 NERRS = NERRS + 1 00945 IF( RESULT( J ).LT.10000.0 ) THEN 00946 WRITE( NOUT, FMT = 9989 )NPTKNT, MPLUSN, J, RESULT( J ) 00947 ELSE 00948 WRITE( NOUT, FMT = 9988 )NPTKNT, MPLUSN, J, RESULT( J ) 00949 END IF 00950 END IF 00951 * 00952 120 CONTINUE 00953 * 00954 130 CONTINUE 00955 GO TO 80 00956 140 CONTINUE 00957 * 00958 150 CONTINUE 00959 * 00960 * Summary 00961 * 00962 CALL ALASVM( 'SGX', NOUT, NERRS, NTESTT, 0 ) 00963 * 00964 WORK( 1 ) = MAXWRK 00965 * 00966 RETURN 00967 * 00968 9999 FORMAT( ' SDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 00969 $ I6, ', JTYPE=', I6, ')' ) 00970 * 00971 9998 FORMAT( ' SDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 00972 $ I6, ', Input Example #', I2, ')' ) 00973 * 00974 9997 FORMAT( ' SDRGSX: SGET53 returned INFO=', I1, ' for eigenvalue ', 00975 $ I6, '.', / 9X, 'N=', I6, ', JTYPE=', I6, ')' ) 00976 * 00977 9996 FORMAT( ' SDRGSX: S not in Schur form at eigenvalue ', I6, '.', 00978 $ / 9X, 'N=', I6, ', JTYPE=', I6, ')' ) 00979 * 00980 9995 FORMAT( / 1X, A3, ' -- Real Expert Generalized Schur form', 00981 $ ' problem driver' ) 00982 * 00983 9994 FORMAT( 'Input Example' ) 00984 * 00985 9993 FORMAT( ' Matrix types: ', / 00986 $ ' 1: A is a block diagonal matrix of Jordan blocks ', 00987 $ 'and B is the identity ', / ' matrix, ', 00988 $ / ' 2: A and B are upper triangular matrices, ', 00989 $ / ' 3: A and B are as type 2, but each second diagonal ', 00990 $ 'block in A_11 and ', / 00991 $ ' each third diaongal block in A_22 are 2x2 blocks,', 00992 $ / ' 4: A and B are block diagonal matrices, ', 00993 $ / ' 5: (A,B) has potentially close or common ', 00994 $ 'eigenvalues.', / ) 00995 * 00996 9992 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ', 00997 $ 'Q and Z are ', A, ',', / 19X, 00998 $ ' a is alpha, b is beta, and ', A, ' means ', A, '.)', 00999 $ / ' 1 = | A - Q S Z', A, 01000 $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', A, 01001 $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', A, 01002 $ ' | / ( n ulp ) 4 = | I - ZZ', A, 01003 $ ' | / ( n ulp )', / ' 5 = 1/ULP if A is not in ', 01004 $ 'Schur form S', / ' 6 = difference between (alpha,beta)', 01005 $ ' and diagonals of (S,T)', / 01006 $ ' 7 = 1/ULP if SDIM is not the correct number of ', 01007 $ 'selected eigenvalues', / 01008 $ ' 8 = 1/ULP if DIFEST/DIFTRU > 10*THRESH or ', 01009 $ 'DIFTRU/DIFEST > 10*THRESH', 01010 $ / ' 9 = 1/ULP if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ', 01011 $ 'when reordering fails', / 01012 $ ' 10 = 1/ULP if PLEST/PLTRU > THRESH or ', 01013 $ 'PLTRU/PLEST > THRESH', / 01014 $ ' ( Test 10 is only for input examples )', / ) 01015 9991 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', E10.3, 01016 $ ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, F8.2 ) 01017 9990 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', E10.3, 01018 $ ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, E10.3 ) 01019 9989 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',', 01020 $ ' result ', I2, ' is', 0P, F8.2 ) 01021 9988 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',', 01022 $ ' result ', I2, ' is', 1P, E10.3 ) 01023 * 01024 * End of SDRGSX 01025 * 01026 END