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