|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE CDRVPP( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, NMAX, 00002 $ A, AFAC, ASAV, B, BSAV, X, XACT, S, WORK, 00003 $ RWORK, NOUT ) 00004 * 00005 * -- LAPACK test routine (version 3.1) -- 00006 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 LOGICAL TSTERR 00011 INTEGER NMAX, NN, NOUT, NRHS 00012 REAL THRESH 00013 * .. 00014 * .. Array Arguments .. 00015 LOGICAL DOTYPE( * ) 00016 INTEGER NVAL( * ) 00017 REAL RWORK( * ), S( * ) 00018 COMPLEX A( * ), AFAC( * ), ASAV( * ), B( * ), 00019 $ BSAV( * ), WORK( * ), X( * ), XACT( * ) 00020 * .. 00021 * 00022 * Purpose 00023 * ======= 00024 * 00025 * CDRVPP tests the driver routines CPPSV and -SVX. 00026 * 00027 * Arguments 00028 * ========= 00029 * 00030 * DOTYPE (input) LOGICAL array, dimension (NTYPES) 00031 * The matrix types to be used for testing. Matrices of type j 00032 * (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = 00033 * .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. 00034 * 00035 * NN (input) INTEGER 00036 * The number of values of N contained in the vector NVAL. 00037 * 00038 * NVAL (input) INTEGER array, dimension (NN) 00039 * The values of the matrix dimension N. 00040 * 00041 * NRHS (input) INTEGER 00042 * The number of right hand side vectors to be generated for 00043 * each linear system. 00044 * 00045 * THRESH (input) REAL 00046 * The threshold value for the test ratios. A result is 00047 * included in the output file if RESULT >= THRESH. To have 00048 * every test ratio printed, use THRESH = 0. 00049 * 00050 * TSTERR (input) LOGICAL 00051 * Flag that indicates whether error exits are to be tested. 00052 * 00053 * NMAX (input) INTEGER 00054 * The maximum value permitted for N, used in dimensioning the 00055 * work arrays. 00056 * 00057 * A (workspace) COMPLEX array, dimension (NMAX*(NMAX+1)/2) 00058 * 00059 * AFAC (workspace) COMPLEX array, dimension (NMAX*(NMAX+1)/2) 00060 * 00061 * ASAV (workspace) COMPLEX array, dimension (NMAX*(NMAX+1)/2) 00062 * 00063 * B (workspace) COMPLEX array, dimension (NMAX*NRHS) 00064 * 00065 * BSAV (workspace) COMPLEX array, dimension (NMAX*NRHS) 00066 * 00067 * X (workspace) COMPLEX array, dimension (NMAX*NRHS) 00068 * 00069 * XACT (workspace) COMPLEX array, dimension (NMAX*NRHS) 00070 * 00071 * S (workspace) REAL array, dimension (NMAX) 00072 * 00073 * WORK (workspace) COMPLEX array, dimension 00074 * (NMAX*max(3,NRHS)) 00075 * 00076 * RWORK (workspace) REAL array, dimension (NMAX+2*NRHS) 00077 * 00078 * NOUT (input) INTEGER 00079 * The unit number for output. 00080 * 00081 * ===================================================================== 00082 * 00083 * .. Parameters .. 00084 REAL ONE, ZERO 00085 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00086 INTEGER NTYPES 00087 PARAMETER ( NTYPES = 9 ) 00088 INTEGER NTESTS 00089 PARAMETER ( NTESTS = 6 ) 00090 * .. 00091 * .. Local Scalars .. 00092 LOGICAL EQUIL, NOFACT, PREFAC, ZEROT 00093 CHARACTER DIST, EQUED, FACT, PACKIT, TYPE, UPLO, XTYPE 00094 CHARACTER*3 PATH 00095 INTEGER I, IEQUED, IFACT, IMAT, IN, INFO, IOFF, IUPLO, 00096 $ IZERO, K, K1, KL, KU, LDA, MODE, N, NERRS, 00097 $ NFACT, NFAIL, NIMAT, NPP, NRUN, NT 00098 REAL AINVNM, AMAX, ANORM, CNDNUM, RCOND, RCONDC, 00099 $ ROLDC, SCOND 00100 * .. 00101 * .. Local Arrays .. 00102 CHARACTER EQUEDS( 2 ), FACTS( 3 ), PACKS( 2 ), UPLOS( 2 ) 00103 INTEGER ISEED( 4 ), ISEEDY( 4 ) 00104 REAL RESULT( NTESTS ) 00105 * .. 00106 * .. External Functions .. 00107 LOGICAL LSAME 00108 REAL CLANHP, SGET06 00109 EXTERNAL LSAME, CLANHP, SGET06 00110 * .. 00111 * .. External Subroutines .. 00112 EXTERNAL ALADHD, ALAERH, ALASVM, CCOPY, CERRVX, CGET04, 00113 $ CLACPY, CLAIPD, CLAQHP, CLARHS, CLASET, CLATB4, 00114 $ CLATMS, CPPEQU, CPPSV, CPPSVX, CPPT01, CPPT02, 00115 $ CPPT05, CPPTRF, CPPTRI 00116 * .. 00117 * .. Scalars in Common .. 00118 LOGICAL LERR, OK 00119 CHARACTER*32 SRNAMT 00120 INTEGER INFOT, NUNIT 00121 * .. 00122 * .. Common blocks .. 00123 COMMON / INFOC / INFOT, NUNIT, OK, LERR 00124 COMMON / SRNAMC / SRNAMT 00125 * .. 00126 * .. Intrinsic Functions .. 00127 INTRINSIC CMPLX, MAX 00128 * .. 00129 * .. Data statements .. 00130 DATA ISEEDY / 1988, 1989, 1990, 1991 / 00131 DATA UPLOS / 'U', 'L' / , FACTS / 'F', 'N', 'E' / , 00132 $ PACKS / 'C', 'R' / , EQUEDS / 'N', 'Y' / 00133 * .. 00134 * .. Executable Statements .. 00135 * 00136 * Initialize constants and the random number seed. 00137 * 00138 PATH( 1: 1 ) = 'Complex precision' 00139 PATH( 2: 3 ) = 'PP' 00140 NRUN = 0 00141 NFAIL = 0 00142 NERRS = 0 00143 DO 10 I = 1, 4 00144 ISEED( I ) = ISEEDY( I ) 00145 10 CONTINUE 00146 * 00147 * Test the error exits 00148 * 00149 IF( TSTERR ) 00150 $ CALL CERRVX( PATH, NOUT ) 00151 INFOT = 0 00152 * 00153 * Do for each value of N in NVAL 00154 * 00155 DO 140 IN = 1, NN 00156 N = NVAL( IN ) 00157 LDA = MAX( N, 1 ) 00158 NPP = N*( N+1 ) / 2 00159 XTYPE = 'N' 00160 NIMAT = NTYPES 00161 IF( N.LE.0 ) 00162 $ NIMAT = 1 00163 * 00164 DO 130 IMAT = 1, NIMAT 00165 * 00166 * Do the tests only if DOTYPE( IMAT ) is true. 00167 * 00168 IF( .NOT.DOTYPE( IMAT ) ) 00169 $ GO TO 130 00170 * 00171 * Skip types 3, 4, or 5 if the matrix size is too small. 00172 * 00173 ZEROT = IMAT.GE.3 .AND. IMAT.LE.5 00174 IF( ZEROT .AND. N.LT.IMAT-2 ) 00175 $ GO TO 130 00176 * 00177 * Do first for UPLO = 'U', then for UPLO = 'L' 00178 * 00179 DO 120 IUPLO = 1, 2 00180 UPLO = UPLOS( IUPLO ) 00181 PACKIT = PACKS( IUPLO ) 00182 * 00183 * Set up parameters with CLATB4 and generate a test matrix 00184 * with CLATMS. 00185 * 00186 CALL CLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE, 00187 $ CNDNUM, DIST ) 00188 RCONDC = ONE / CNDNUM 00189 * 00190 SRNAMT = 'CLATMS' 00191 CALL CLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, 00192 $ CNDNUM, ANORM, KL, KU, PACKIT, A, LDA, WORK, 00193 $ INFO ) 00194 * 00195 * Check error code from CLATMS. 00196 * 00197 IF( INFO.NE.0 ) THEN 00198 CALL ALAERH( PATH, 'CLATMS', INFO, 0, UPLO, N, N, -1, 00199 $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) 00200 GO TO 120 00201 END IF 00202 * 00203 * For types 3-5, zero one row and column of the matrix to 00204 * test that INFO is returned correctly. 00205 * 00206 IF( ZEROT ) THEN 00207 IF( IMAT.EQ.3 ) THEN 00208 IZERO = 1 00209 ELSE IF( IMAT.EQ.4 ) THEN 00210 IZERO = N 00211 ELSE 00212 IZERO = N / 2 + 1 00213 END IF 00214 * 00215 * Set row and column IZERO of A to 0. 00216 * 00217 IF( IUPLO.EQ.1 ) THEN 00218 IOFF = ( IZERO-1 )*IZERO / 2 00219 DO 20 I = 1, IZERO - 1 00220 A( IOFF+I ) = ZERO 00221 20 CONTINUE 00222 IOFF = IOFF + IZERO 00223 DO 30 I = IZERO, N 00224 A( IOFF ) = ZERO 00225 IOFF = IOFF + I 00226 30 CONTINUE 00227 ELSE 00228 IOFF = IZERO 00229 DO 40 I = 1, IZERO - 1 00230 A( IOFF ) = ZERO 00231 IOFF = IOFF + N - I 00232 40 CONTINUE 00233 IOFF = IOFF - IZERO 00234 DO 50 I = IZERO, N 00235 A( IOFF+I ) = ZERO 00236 50 CONTINUE 00237 END IF 00238 ELSE 00239 IZERO = 0 00240 END IF 00241 * 00242 * Set the imaginary part of the diagonals. 00243 * 00244 IF( IUPLO.EQ.1 ) THEN 00245 CALL CLAIPD( N, A, 2, 1 ) 00246 ELSE 00247 CALL CLAIPD( N, A, N, -1 ) 00248 END IF 00249 * 00250 * Save a copy of the matrix A in ASAV. 00251 * 00252 CALL CCOPY( NPP, A, 1, ASAV, 1 ) 00253 * 00254 DO 110 IEQUED = 1, 2 00255 EQUED = EQUEDS( IEQUED ) 00256 IF( IEQUED.EQ.1 ) THEN 00257 NFACT = 3 00258 ELSE 00259 NFACT = 1 00260 END IF 00261 * 00262 DO 100 IFACT = 1, NFACT 00263 FACT = FACTS( IFACT ) 00264 PREFAC = LSAME( FACT, 'F' ) 00265 NOFACT = LSAME( FACT, 'N' ) 00266 EQUIL = LSAME( FACT, 'E' ) 00267 * 00268 IF( ZEROT ) THEN 00269 IF( PREFAC ) 00270 $ GO TO 100 00271 RCONDC = ZERO 00272 * 00273 ELSE IF( .NOT.LSAME( FACT, 'N' ) ) THEN 00274 * 00275 * Compute the condition number for comparison with 00276 * the value returned by CPPSVX (FACT = 'N' reuses 00277 * the condition number from the previous iteration 00278 * with FACT = 'F'). 00279 * 00280 CALL CCOPY( NPP, ASAV, 1, AFAC, 1 ) 00281 IF( EQUIL .OR. IEQUED.GT.1 ) THEN 00282 * 00283 * Compute row and column scale factors to 00284 * equilibrate the matrix A. 00285 * 00286 CALL CPPEQU( UPLO, N, AFAC, S, SCOND, AMAX, 00287 $ INFO ) 00288 IF( INFO.EQ.0 .AND. N.GT.0 ) THEN 00289 IF( IEQUED.GT.1 ) 00290 $ SCOND = ZERO 00291 * 00292 * Equilibrate the matrix. 00293 * 00294 CALL CLAQHP( UPLO, N, AFAC, S, SCOND, 00295 $ AMAX, EQUED ) 00296 END IF 00297 END IF 00298 * 00299 * Save the condition number of the 00300 * non-equilibrated system for use in CGET04. 00301 * 00302 IF( EQUIL ) 00303 $ ROLDC = RCONDC 00304 * 00305 * Compute the 1-norm of A. 00306 * 00307 ANORM = CLANHP( '1', UPLO, N, AFAC, RWORK ) 00308 * 00309 * Factor the matrix A. 00310 * 00311 CALL CPPTRF( UPLO, N, AFAC, INFO ) 00312 * 00313 * Form the inverse of A. 00314 * 00315 CALL CCOPY( NPP, AFAC, 1, A, 1 ) 00316 CALL CPPTRI( UPLO, N, A, INFO ) 00317 * 00318 * Compute the 1-norm condition number of A. 00319 * 00320 AINVNM = CLANHP( '1', UPLO, N, A, RWORK ) 00321 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 00322 RCONDC = ONE 00323 ELSE 00324 RCONDC = ( ONE / ANORM ) / AINVNM 00325 END IF 00326 END IF 00327 * 00328 * Restore the matrix A. 00329 * 00330 CALL CCOPY( NPP, ASAV, 1, A, 1 ) 00331 * 00332 * Form an exact solution and set the right hand side. 00333 * 00334 SRNAMT = 'CLARHS' 00335 CALL CLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU, 00336 $ NRHS, A, LDA, XACT, LDA, B, LDA, 00337 $ ISEED, INFO ) 00338 XTYPE = 'C' 00339 CALL CLACPY( 'Full', N, NRHS, B, LDA, BSAV, LDA ) 00340 * 00341 IF( NOFACT ) THEN 00342 * 00343 * --- Test CPPSV --- 00344 * 00345 * Compute the L*L' or U'*U factorization of the 00346 * matrix and solve the system. 00347 * 00348 CALL CCOPY( NPP, A, 1, AFAC, 1 ) 00349 CALL CLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) 00350 * 00351 SRNAMT = 'CPPSV ' 00352 CALL CPPSV( UPLO, N, NRHS, AFAC, X, LDA, INFO ) 00353 * 00354 * Check error code from CPPSV . 00355 * 00356 IF( INFO.NE.IZERO ) THEN 00357 CALL ALAERH( PATH, 'CPPSV ', INFO, IZERO, 00358 $ UPLO, N, N, -1, -1, NRHS, IMAT, 00359 $ NFAIL, NERRS, NOUT ) 00360 GO TO 70 00361 ELSE IF( INFO.NE.0 ) THEN 00362 GO TO 70 00363 END IF 00364 * 00365 * Reconstruct matrix from factors and compute 00366 * residual. 00367 * 00368 CALL CPPT01( UPLO, N, A, AFAC, RWORK, 00369 $ RESULT( 1 ) ) 00370 * 00371 * Compute residual of the computed solution. 00372 * 00373 CALL CLACPY( 'Full', N, NRHS, B, LDA, WORK, 00374 $ LDA ) 00375 CALL CPPT02( UPLO, N, NRHS, A, X, LDA, WORK, 00376 $ LDA, RWORK, RESULT( 2 ) ) 00377 * 00378 * Check solution from generated exact solution. 00379 * 00380 CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00381 $ RESULT( 3 ) ) 00382 NT = 3 00383 * 00384 * Print information about the tests that did not 00385 * pass the threshold. 00386 * 00387 DO 60 K = 1, NT 00388 IF( RESULT( K ).GE.THRESH ) THEN 00389 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00390 $ CALL ALADHD( NOUT, PATH ) 00391 WRITE( NOUT, FMT = 9999 )'CPPSV ', UPLO, 00392 $ N, IMAT, K, RESULT( K ) 00393 NFAIL = NFAIL + 1 00394 END IF 00395 60 CONTINUE 00396 NRUN = NRUN + NT 00397 70 CONTINUE 00398 END IF 00399 * 00400 * --- Test CPPSVX --- 00401 * 00402 IF( .NOT.PREFAC .AND. NPP.GT.0 ) 00403 $ CALL CLASET( 'Full', NPP, 1, CMPLX( ZERO ), 00404 $ CMPLX( ZERO ), AFAC, NPP ) 00405 CALL CLASET( 'Full', N, NRHS, CMPLX( ZERO ), 00406 $ CMPLX( ZERO ), X, LDA ) 00407 IF( IEQUED.GT.1 .AND. N.GT.0 ) THEN 00408 * 00409 * Equilibrate the matrix if FACT='F' and 00410 * EQUED='Y'. 00411 * 00412 CALL CLAQHP( UPLO, N, A, S, SCOND, AMAX, EQUED ) 00413 END IF 00414 * 00415 * Solve the system and compute the condition number 00416 * and error bounds using CPPSVX. 00417 * 00418 SRNAMT = 'CPPSVX' 00419 CALL CPPSVX( FACT, UPLO, N, NRHS, A, AFAC, EQUED, 00420 $ S, B, LDA, X, LDA, RCOND, RWORK, 00421 $ RWORK( NRHS+1 ), WORK, 00422 $ RWORK( 2*NRHS+1 ), INFO ) 00423 * 00424 * Check the error code from CPPSVX. 00425 * 00426 IF( INFO.NE.IZERO ) THEN 00427 CALL ALAERH( PATH, 'CPPSVX', INFO, IZERO, 00428 $ FACT // UPLO, N, N, -1, -1, NRHS, 00429 $ IMAT, NFAIL, NERRS, NOUT ) 00430 GO TO 90 00431 END IF 00432 * 00433 IF( INFO.EQ.0 ) THEN 00434 IF( .NOT.PREFAC ) THEN 00435 * 00436 * Reconstruct matrix from factors and compute 00437 * residual. 00438 * 00439 CALL CPPT01( UPLO, N, A, AFAC, 00440 $ RWORK( 2*NRHS+1 ), RESULT( 1 ) ) 00441 K1 = 1 00442 ELSE 00443 K1 = 2 00444 END IF 00445 * 00446 * Compute residual of the computed solution. 00447 * 00448 CALL CLACPY( 'Full', N, NRHS, BSAV, LDA, WORK, 00449 $ LDA ) 00450 CALL CPPT02( UPLO, N, NRHS, ASAV, X, LDA, WORK, 00451 $ LDA, RWORK( 2*NRHS+1 ), 00452 $ RESULT( 2 ) ) 00453 * 00454 * Check solution from generated exact solution. 00455 * 00456 IF( NOFACT .OR. ( PREFAC .AND. LSAME( EQUED, 00457 $ 'N' ) ) ) THEN 00458 CALL CGET04( N, NRHS, X, LDA, XACT, LDA, 00459 $ RCONDC, RESULT( 3 ) ) 00460 ELSE 00461 CALL CGET04( N, NRHS, X, LDA, XACT, LDA, 00462 $ ROLDC, RESULT( 3 ) ) 00463 END IF 00464 * 00465 * Check the error bounds from iterative 00466 * refinement. 00467 * 00468 CALL CPPT05( UPLO, N, NRHS, ASAV, B, LDA, X, 00469 $ LDA, XACT, LDA, RWORK, 00470 $ RWORK( NRHS+1 ), RESULT( 4 ) ) 00471 ELSE 00472 K1 = 6 00473 END IF 00474 * 00475 * Compare RCOND from CPPSVX with the computed value 00476 * in RCONDC. 00477 * 00478 RESULT( 6 ) = SGET06( RCOND, RCONDC ) 00479 * 00480 * Print information about the tests that did not pass 00481 * the threshold. 00482 * 00483 DO 80 K = K1, 6 00484 IF( RESULT( K ).GE.THRESH ) THEN 00485 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00486 $ CALL ALADHD( NOUT, PATH ) 00487 IF( PREFAC ) THEN 00488 WRITE( NOUT, FMT = 9997 )'CPPSVX', FACT, 00489 $ UPLO, N, EQUED, IMAT, K, RESULT( K ) 00490 ELSE 00491 WRITE( NOUT, FMT = 9998 )'CPPSVX', FACT, 00492 $ UPLO, N, IMAT, K, RESULT( K ) 00493 END IF 00494 NFAIL = NFAIL + 1 00495 END IF 00496 80 CONTINUE 00497 NRUN = NRUN + 7 - K1 00498 90 CONTINUE 00499 100 CONTINUE 00500 110 CONTINUE 00501 120 CONTINUE 00502 130 CONTINUE 00503 140 CONTINUE 00504 * 00505 * Print a summary of the results. 00506 * 00507 CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS ) 00508 * 00509 9999 FORMAT( 1X, A, ', UPLO=''', A1, ''', N =', I5, ', type ', I1, 00510 $ ', test(', I1, ')=', G12.5 ) 00511 9998 FORMAT( 1X, A, ', FACT=''', A1, ''', UPLO=''', A1, ''', N=', I5, 00512 $ ', type ', I1, ', test(', I1, ')=', G12.5 ) 00513 9997 FORMAT( 1X, A, ', FACT=''', A1, ''', UPLO=''', A1, ''', N=', I5, 00514 $ ', EQUED=''', A1, ''', type ', I1, ', test(', I1, ')=', 00515 $ G12.5 ) 00516 RETURN 00517 * 00518 * End of CDRVPP 00519 * 00520 END
1.7.3 
