 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE ZCHKGE( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, 00002 $ NSVAL, THRESH, TSTERR, NMAX, A, AFAC, AINV, B, 00003 $ X, XACT, WORK, RWORK, IWORK, NOUT ) 00004 * 00005 * -- LAPACK test routine (version 3.1.1) -- 00006 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00007 * January 2007 00008 * 00009 * .. Scalar Arguments .. 00010 LOGICAL TSTERR 00011 INTEGER NM, NMAX, NN, NNB, NNS, NOUT 00012 DOUBLE PRECISION THRESH 00013 * .. 00014 * .. Array Arguments .. 00015 LOGICAL DOTYPE( * ) 00016 INTEGER IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ), 00017 $ NVAL( * ) 00018 DOUBLE PRECISION RWORK( * ) 00019 COMPLEX*16 A( * ), AFAC( * ), AINV( * ), B( * ), 00020 $ WORK( * ), X( * ), XACT( * ) 00021 * .. 00022 * 00023 * Purpose 00024 * ======= 00025 * 00026 * ZCHKGE tests ZGETRF, -TRI, -TRS, -RFS, and -CON. 00027 * 00028 * Arguments 00029 * ========= 00030 * 00031 * DOTYPE (input) LOGICAL array, dimension (NTYPES) 00032 * The matrix types to be used for testing. Matrices of type j 00033 * (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = 00034 * .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. 00035 * 00036 * NM (input) INTEGER 00037 * The number of values of M contained in the vector MVAL. 00038 * 00039 * MVAL (input) INTEGER array, dimension (NM) 00040 * The values of the matrix row dimension M. 00041 * 00042 * NN (input) INTEGER 00043 * The number of values of N contained in the vector NVAL. 00044 * 00045 * NVAL (input) INTEGER array, dimension (NN) 00046 * The values of the matrix column dimension N. 00047 * 00048 * NNB (input) INTEGER 00049 * The number of values of NB contained in the vector NBVAL. 00050 * 00051 * NBVAL (input) INTEGER array, dimension (NBVAL) 00052 * The values of the blocksize NB. 00053 * 00054 * NNS (input) INTEGER 00055 * The number of values of NRHS contained in the vector NSVAL. 00056 * 00057 * NSVAL (input) INTEGER array, dimension (NNS) 00058 * The values of the number of right hand sides NRHS. 00059 * 00060 * NRHS (input) INTEGER 00061 * The number of right hand side vectors to be generated for 00062 * each linear system. 00063 * 00064 * THRESH (input) DOUBLE PRECISION 00065 * The threshold value for the test ratios. A result is 00066 * included in the output file if RESULT >= THRESH. To have 00067 * every test ratio printed, use THRESH = 0. 00068 * 00069 * TSTERR (input) LOGICAL 00070 * Flag that indicates whether error exits are to be tested. 00071 * 00072 * NMAX (input) INTEGER 00073 * The maximum value permitted for M or N, used in dimensioning 00074 * the work arrays. 00075 * 00076 * A (workspace) COMPLEX*16 array, dimension (NMAX*NMAX) 00077 * 00078 * AFAC (workspace) COMPLEX*16 array, dimension (NMAX*NMAX) 00079 * 00080 * AINV (workspace) COMPLEX*16 array, dimension (NMAX*NMAX) 00081 * 00082 * B (workspace) COMPLEX*16 array, dimension (NMAX*NSMAX) 00083 * where NSMAX is the largest entry in NSVAL. 00084 * 00085 * X (workspace) COMPLEX*16 array, dimension (NMAX*NSMAX) 00086 * 00087 * XACT (workspace) COMPLEX*16 array, dimension (NMAX*NSMAX) 00088 * 00089 * WORK (workspace) COMPLEX*16 array, dimension 00090 * (NMAX*max(3,NSMAX)) 00091 * 00092 * RWORK (workspace) DOUBLE PRECISION array, dimension 00093 * (max(2*NMAX,2*NSMAX+NWORK)) 00094 * 00095 * IWORK (workspace) INTEGER array, dimension (NMAX) 00096 * 00097 * NOUT (input) INTEGER 00098 * The unit number for output. 00099 * 00100 * ===================================================================== 00101 * 00102 * .. Parameters .. 00103 DOUBLE PRECISION ONE, ZERO 00104 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00105 INTEGER NTYPES 00106 PARAMETER ( NTYPES = 11 ) 00107 INTEGER NTESTS 00108 PARAMETER ( NTESTS = 8 ) 00109 INTEGER NTRAN 00110 PARAMETER ( NTRAN = 3 ) 00111 * .. 00112 * .. Local Scalars .. 00113 LOGICAL TRFCON, ZEROT 00114 CHARACTER DIST, NORM, TRANS, TYPE, XTYPE 00115 CHARACTER*3 PATH 00116 INTEGER I, IM, IMAT, IN, INB, INFO, IOFF, IRHS, ITRAN, 00117 $ IZERO, K, KL, KU, LDA, LWORK, M, MODE, N, NB, 00118 $ NERRS, NFAIL, NIMAT, NRHS, NRUN, NT 00119 DOUBLE PRECISION AINVNM, ANORM, ANORMI, ANORMO, CNDNUM, DUMMY, 00120 $ RCOND, RCONDC, RCONDI, RCONDO 00121 * .. 00122 * .. Local Arrays .. 00123 CHARACTER TRANSS( NTRAN ) 00124 INTEGER ISEED( 4 ), ISEEDY( 4 ) 00125 DOUBLE PRECISION RESULT( NTESTS ) 00126 * .. 00127 * .. External Functions .. 00128 DOUBLE PRECISION DGET06, ZLANGE 00129 EXTERNAL DGET06, ZLANGE 00130 * .. 00131 * .. External Subroutines .. 00132 EXTERNAL ALAERH, ALAHD, ALASUM, XLAENV, ZERRGE, ZGECON, 00133 $ ZGERFS, ZGET01, ZGET02, ZGET03, ZGET04, ZGET07, 00134 $ ZGETRF, ZGETRI, ZGETRS, ZLACPY, ZLARHS, ZLASET, 00135 $ ZLATB4, ZLATMS 00136 * .. 00137 * .. Intrinsic Functions .. 00138 INTRINSIC DCMPLX, MAX, MIN 00139 * .. 00140 * .. Scalars in Common .. 00141 LOGICAL LERR, OK 00142 CHARACTER*32 SRNAMT 00143 INTEGER INFOT, NUNIT 00144 * .. 00145 * .. Common blocks .. 00146 COMMON / INFOC / INFOT, NUNIT, OK, LERR 00147 COMMON / SRNAMC / SRNAMT 00148 * .. 00149 * .. Data statements .. 00150 DATA ISEEDY / 1988, 1989, 1990, 1991 / , 00151 $ TRANSS / 'N', 'T', 'C' / 00152 * .. 00153 * .. Executable Statements .. 00154 * 00155 * Initialize constants and the random number seed. 00156 * 00157 PATH( 1: 1 ) = 'Zomplex precision' 00158 PATH( 2: 3 ) = 'GE' 00159 NRUN = 0 00160 NFAIL = 0 00161 NERRS = 0 00162 DO 10 I = 1, 4 00163 ISEED( I ) = ISEEDY( I ) 00164 10 CONTINUE 00165 * 00166 * Test the error exits 00167 * 00168 CALL XLAENV( 1, 1 ) 00169 IF( TSTERR ) 00170 $ CALL ZERRGE( PATH, NOUT ) 00171 INFOT = 0 00172 CALL XLAENV( 2, 2 ) 00173 * 00174 * Do for each value of M in MVAL 00175 * 00176 DO 120 IM = 1, NM 00177 M = MVAL( IM ) 00178 LDA = MAX( 1, M ) 00179 * 00180 * Do for each value of N in NVAL 00181 * 00182 DO 110 IN = 1, NN 00183 N = NVAL( IN ) 00184 XTYPE = 'N' 00185 NIMAT = NTYPES 00186 IF( M.LE.0 .OR. N.LE.0 ) 00187 $ NIMAT = 1 00188 * 00189 DO 100 IMAT = 1, NIMAT 00190 * 00191 * Do the tests only if DOTYPE( IMAT ) is true. 00192 * 00193 IF( .NOT.DOTYPE( IMAT ) ) 00194 $ GO TO 100 00195 * 00196 * Skip types 5, 6, or 7 if the matrix size is too small. 00197 * 00198 ZEROT = IMAT.GE.5 .AND. IMAT.LE.7 00199 IF( ZEROT .AND. N.LT.IMAT-4 ) 00200 $ GO TO 100 00201 * 00202 * Set up parameters with ZLATB4 and generate a test matrix 00203 * with ZLATMS. 00204 * 00205 CALL ZLATB4( PATH, IMAT, M, N, TYPE, KL, KU, ANORM, MODE, 00206 $ CNDNUM, DIST ) 00207 * 00208 SRNAMT = 'ZLATMS' 00209 CALL ZLATMS( M, N, DIST, ISEED, TYPE, RWORK, MODE, 00210 $ CNDNUM, ANORM, KL, KU, 'No packing', A, LDA, 00211 $ WORK, INFO ) 00212 * 00213 * Check error code from ZLATMS. 00214 * 00215 IF( INFO.NE.0 ) THEN 00216 CALL ALAERH( PATH, 'ZLATMS', INFO, 0, ' ', M, N, -1, 00217 $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) 00218 GO TO 100 00219 END IF 00220 * 00221 * For types 5-7, zero one or more columns of the matrix to 00222 * test that INFO is returned correctly. 00223 * 00224 IF( ZEROT ) THEN 00225 IF( IMAT.EQ.5 ) THEN 00226 IZERO = 1 00227 ELSE IF( IMAT.EQ.6 ) THEN 00228 IZERO = MIN( M, N ) 00229 ELSE 00230 IZERO = MIN( M, N ) / 2 + 1 00231 END IF 00232 IOFF = ( IZERO-1 )*LDA 00233 IF( IMAT.LT.7 ) THEN 00234 DO 20 I = 1, M 00235 A( IOFF+I ) = ZERO 00236 20 CONTINUE 00237 ELSE 00238 CALL ZLASET( 'Full', M, N-IZERO+1, DCMPLX( ZERO ), 00239 $ DCMPLX( ZERO ), A( IOFF+1 ), LDA ) 00240 END IF 00241 ELSE 00242 IZERO = 0 00243 END IF 00244 * 00245 * These lines, if used in place of the calls in the DO 60 00246 * loop, cause the code to bomb on a Sun SPARCstation. 00247 * 00248 * ANORMO = ZLANGE( 'O', M, N, A, LDA, RWORK ) 00249 * ANORMI = ZLANGE( 'I', M, N, A, LDA, RWORK ) 00250 * 00251 * Do for each blocksize in NBVAL 00252 * 00253 DO 90 INB = 1, NNB 00254 NB = NBVAL( INB ) 00255 CALL XLAENV( 1, NB ) 00256 * 00257 * Compute the LU factorization of the matrix. 00258 * 00259 CALL ZLACPY( 'Full', M, N, A, LDA, AFAC, LDA ) 00260 SRNAMT = 'ZGETRF' 00261 CALL ZGETRF( M, N, AFAC, LDA, IWORK, INFO ) 00262 * 00263 * Check error code from ZGETRF. 00264 * 00265 IF( INFO.NE.IZERO ) 00266 $ CALL ALAERH( PATH, 'ZGETRF', INFO, IZERO, ' ', M, 00267 $ N, -1, -1, NB, IMAT, NFAIL, NERRS, 00268 $ NOUT ) 00269 TRFCON = .FALSE. 00270 * 00271 *+ TEST 1 00272 * Reconstruct matrix from factors and compute residual. 00273 * 00274 CALL ZLACPY( 'Full', M, N, AFAC, LDA, AINV, LDA ) 00275 CALL ZGET01( M, N, A, LDA, AINV, LDA, IWORK, RWORK, 00276 $ RESULT( 1 ) ) 00277 NT = 1 00278 * 00279 *+ TEST 2 00280 * Form the inverse if the factorization was successful 00281 * and compute the residual. 00282 * 00283 IF( M.EQ.N .AND. INFO.EQ.0 ) THEN 00284 CALL ZLACPY( 'Full', N, N, AFAC, LDA, AINV, LDA ) 00285 SRNAMT = 'ZGETRI' 00286 NRHS = NSVAL( 1 ) 00287 LWORK = NMAX*MAX( 3, NRHS ) 00288 CALL ZGETRI( N, AINV, LDA, IWORK, WORK, LWORK, 00289 $ INFO ) 00290 * 00291 * Check error code from ZGETRI. 00292 * 00293 IF( INFO.NE.0 ) 00294 $ CALL ALAERH( PATH, 'ZGETRI', INFO, 0, ' ', N, N, 00295 $ -1, -1, NB, IMAT, NFAIL, NERRS, 00296 $ NOUT ) 00297 * 00298 * Compute the residual for the matrix times its 00299 * inverse. Also compute the 1-norm condition number 00300 * of A. 00301 * 00302 CALL ZGET03( N, A, LDA, AINV, LDA, WORK, LDA, 00303 $ RWORK, RCONDO, RESULT( 2 ) ) 00304 ANORMO = ZLANGE( 'O', M, N, A, LDA, RWORK ) 00305 * 00306 * Compute the infinity-norm condition number of A. 00307 * 00308 ANORMI = ZLANGE( 'I', M, N, A, LDA, RWORK ) 00309 AINVNM = ZLANGE( 'I', N, N, AINV, LDA, RWORK ) 00310 IF( ANORMI.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 00311 RCONDI = ONE 00312 ELSE 00313 RCONDI = ( ONE / ANORMI ) / AINVNM 00314 END IF 00315 NT = 2 00316 ELSE 00317 * 00318 * Do only the condition estimate if INFO > 0. 00319 * 00320 TRFCON = .TRUE. 00321 ANORMO = ZLANGE( 'O', M, N, A, LDA, RWORK ) 00322 ANORMI = ZLANGE( 'I', M, N, A, LDA, RWORK ) 00323 RCONDO = ZERO 00324 RCONDI = ZERO 00325 END IF 00326 * 00327 * Print information about the tests so far that did not 00328 * pass the threshold. 00329 * 00330 DO 30 K = 1, NT 00331 IF( RESULT( K ).GE.THRESH ) THEN 00332 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00333 $ CALL ALAHD( NOUT, PATH ) 00334 WRITE( NOUT, FMT = 9999 )M, N, NB, IMAT, K, 00335 $ RESULT( K ) 00336 NFAIL = NFAIL + 1 00337 END IF 00338 30 CONTINUE 00339 NRUN = NRUN + NT 00340 * 00341 * Skip the remaining tests if this is not the first 00342 * block size or if M .ne. N. Skip the solve tests if 00343 * the matrix is singular. 00344 * 00345 IF( INB.GT.1 .OR. M.NE.N ) 00346 $ GO TO 90 00347 IF( TRFCON ) 00348 $ GO TO 70 00349 * 00350 DO 60 IRHS = 1, NNS 00351 NRHS = NSVAL( IRHS ) 00352 XTYPE = 'N' 00353 * 00354 DO 50 ITRAN = 1, NTRAN 00355 TRANS = TRANSS( ITRAN ) 00356 IF( ITRAN.EQ.1 ) THEN 00357 RCONDC = RCONDO 00358 ELSE 00359 RCONDC = RCONDI 00360 END IF 00361 * 00362 *+ TEST 3 00363 * Solve and compute residual for A * X = B. 00364 * 00365 SRNAMT = 'ZLARHS' 00366 CALL ZLARHS( PATH, XTYPE, ' ', TRANS, N, N, KL, 00367 $ KU, NRHS, A, LDA, XACT, LDA, B, 00368 $ LDA, ISEED, INFO ) 00369 XTYPE = 'C' 00370 * 00371 CALL ZLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) 00372 SRNAMT = 'ZGETRS' 00373 CALL ZGETRS( TRANS, N, NRHS, AFAC, LDA, IWORK, 00374 $ X, LDA, INFO ) 00375 * 00376 * Check error code from ZGETRS. 00377 * 00378 IF( INFO.NE.0 ) 00379 $ CALL ALAERH( PATH, 'ZGETRS', INFO, 0, TRANS, 00380 $ N, N, -1, -1, NRHS, IMAT, NFAIL, 00381 $ NERRS, NOUT ) 00382 * 00383 CALL ZLACPY( 'Full', N, NRHS, B, LDA, WORK, 00384 $ LDA ) 00385 CALL ZGET02( TRANS, N, N, NRHS, A, LDA, X, LDA, 00386 $ WORK, LDA, RWORK, RESULT( 3 ) ) 00387 * 00388 *+ TEST 4 00389 * Check solution from generated exact solution. 00390 * 00391 CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00392 $ RESULT( 4 ) ) 00393 * 00394 *+ TESTS 5, 6, and 7 00395 * Use iterative refinement to improve the 00396 * solution. 00397 * 00398 SRNAMT = 'ZGERFS' 00399 CALL ZGERFS( TRANS, N, NRHS, A, LDA, AFAC, LDA, 00400 $ IWORK, B, LDA, X, LDA, RWORK, 00401 $ RWORK( NRHS+1 ), WORK, 00402 $ RWORK( 2*NRHS+1 ), INFO ) 00403 * 00404 * Check error code from ZGERFS. 00405 * 00406 IF( INFO.NE.0 ) 00407 $ CALL ALAERH( PATH, 'ZGERFS', INFO, 0, TRANS, 00408 $ N, N, -1, -1, NRHS, IMAT, NFAIL, 00409 $ NERRS, NOUT ) 00410 * 00411 CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00412 $ RESULT( 5 ) ) 00413 CALL ZGET07( TRANS, N, NRHS, A, LDA, B, LDA, X, 00414 $ LDA, XACT, LDA, RWORK, .TRUE., 00415 $ RWORK( NRHS+1 ), RESULT( 6 ) ) 00416 * 00417 * Print information about the tests that did not 00418 * pass the threshold. 00419 * 00420 DO 40 K = 3, 7 00421 IF( RESULT( K ).GE.THRESH ) THEN 00422 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00423 $ CALL ALAHD( NOUT, PATH ) 00424 WRITE( NOUT, FMT = 9998 )TRANS, N, NRHS, 00425 $ IMAT, K, RESULT( K ) 00426 NFAIL = NFAIL + 1 00427 END IF 00428 40 CONTINUE 00429 NRUN = NRUN + 5 00430 50 CONTINUE 00431 60 CONTINUE 00432 * 00433 *+ TEST 8 00434 * Get an estimate of RCOND = 1/CNDNUM. 00435 * 00436 70 CONTINUE 00437 DO 80 ITRAN = 1, 2 00438 IF( ITRAN.EQ.1 ) THEN 00439 ANORM = ANORMO 00440 RCONDC = RCONDO 00441 NORM = 'O' 00442 ELSE 00443 ANORM = ANORMI 00444 RCONDC = RCONDI 00445 NORM = 'I' 00446 END IF 00447 SRNAMT = 'ZGECON' 00448 CALL ZGECON( NORM, N, AFAC, LDA, ANORM, RCOND, 00449 $ WORK, RWORK, INFO ) 00450 * 00451 * Check error code from ZGECON. 00452 * 00453 IF( INFO.NE.0 ) 00454 $ CALL ALAERH( PATH, 'ZGECON', INFO, 0, NORM, N, 00455 $ N, -1, -1, -1, IMAT, NFAIL, NERRS, 00456 $ NOUT ) 00457 * 00458 * This line is needed on a Sun SPARCstation. 00459 * 00460 DUMMY = RCOND 00461 * 00462 RESULT( 8 ) = DGET06( RCOND, RCONDC ) 00463 * 00464 * Print information about the tests that did not pass 00465 * the threshold. 00466 * 00467 IF( RESULT( 8 ).GE.THRESH ) THEN 00468 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00469 $ CALL ALAHD( NOUT, PATH ) 00470 WRITE( NOUT, FMT = 9997 )NORM, N, IMAT, 8, 00471 $ RESULT( 8 ) 00472 NFAIL = NFAIL + 1 00473 END IF 00474 NRUN = NRUN + 1 00475 80 CONTINUE 00476 90 CONTINUE 00477 100 CONTINUE 00478 * 00479 110 CONTINUE 00480 120 CONTINUE 00481 * 00482 * Print a summary of the results. 00483 * 00484 CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS ) 00485 * 00486 9999 FORMAT( ' M = ', I5, ', N =', I5, ', NB =', I4, ', type ', I2, 00487 $ ', test(', I2, ') =', G12.5 ) 00488 9998 FORMAT( ' TRANS=''', A1, ''', N =', I5, ', NRHS=', I3, ', type ', 00489 $ I2, ', test(', I2, ') =', G12.5 ) 00490 9997 FORMAT( ' NORM =''', A1, ''', N =', I5, ',', 10X, ' type ', I2, 00491 $ ', test(', I2, ') =', G12.5 ) 00492 RETURN 00493 * 00494 * End of ZCHKGE 00495 * 00496 END
1.7.3 
