 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CDRVPO( 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 * CDRVPO tests the driver routines CPOSV 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) 00058 * 00059 * AFAC (workspace) COMPLEX array, dimension (NMAX*NMAX) 00060 * 00061 * ASAV (workspace) COMPLEX array, dimension (NMAX*NMAX) 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, 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, NB, NBMIN, 00097 $ NERRS, NFACT, NFAIL, NIMAT, NRUN, NT 00098 REAL AINVNM, AMAX, ANORM, CNDNUM, RCOND, RCONDC, 00099 $ ROLDC, SCOND 00100 * .. 00101 * .. Local Arrays .. 00102 CHARACTER EQUEDS( 2 ), FACTS( 3 ), UPLOS( 2 ) 00103 INTEGER ISEED( 4 ), ISEEDY( 4 ) 00104 REAL RESULT( NTESTS ) 00105 * .. 00106 * .. External Functions .. 00107 LOGICAL LSAME 00108 REAL CLANHE, SGET06 00109 EXTERNAL LSAME, CLANHE, SGET06 00110 * .. 00111 * .. External Subroutines .. 00112 EXTERNAL ALADHD, ALAERH, ALASVM, CERRVX, CGET04, CLACPY, 00113 $ CLAIPD, CLAQHE, CLARHS, CLASET, CLATB4, CLATMS, 00114 $ CPOEQU, CPOSV, CPOSVX, CPOT01, CPOT02, CPOT05, 00115 $ CPOTRF, CPOTRI, XLAENV 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' / 00132 DATA FACTS / 'F', 'N', 'E' / 00133 DATA EQUEDS / 'N', 'Y' / 00134 * .. 00135 * .. Executable Statements .. 00136 * 00137 * Initialize constants and the random number seed. 00138 * 00139 PATH( 1: 1 ) = 'Complex precision' 00140 PATH( 2: 3 ) = 'PO' 00141 NRUN = 0 00142 NFAIL = 0 00143 NERRS = 0 00144 DO 10 I = 1, 4 00145 ISEED( I ) = ISEEDY( I ) 00146 10 CONTINUE 00147 * 00148 * Test the error exits 00149 * 00150 IF( TSTERR ) 00151 $ CALL CERRVX( PATH, NOUT ) 00152 INFOT = 0 00153 * 00154 * Set the block size and minimum block size for testing. 00155 * 00156 NB = 1 00157 NBMIN = 2 00158 CALL XLAENV( 1, NB ) 00159 CALL XLAENV( 2, NBMIN ) 00160 * 00161 * Do for each value of N in NVAL 00162 * 00163 DO 130 IN = 1, NN 00164 N = NVAL( IN ) 00165 LDA = MAX( N, 1 ) 00166 XTYPE = 'N' 00167 NIMAT = NTYPES 00168 IF( N.LE.0 ) 00169 $ NIMAT = 1 00170 * 00171 DO 120 IMAT = 1, NIMAT 00172 * 00173 * Do the tests only if DOTYPE( IMAT ) is true. 00174 * 00175 IF( .NOT.DOTYPE( IMAT ) ) 00176 $ GO TO 120 00177 * 00178 * Skip types 3, 4, or 5 if the matrix size is too small. 00179 * 00180 ZEROT = IMAT.GE.3 .AND. IMAT.LE.5 00181 IF( ZEROT .AND. N.LT.IMAT-2 ) 00182 $ GO TO 120 00183 * 00184 * Do first for UPLO = 'U', then for UPLO = 'L' 00185 * 00186 DO 110 IUPLO = 1, 2 00187 UPLO = UPLOS( IUPLO ) 00188 * 00189 * Set up parameters with CLATB4 and generate a test matrix 00190 * with CLATMS. 00191 * 00192 CALL CLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE, 00193 $ CNDNUM, DIST ) 00194 * 00195 SRNAMT = 'CLATMS' 00196 CALL CLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, 00197 $ CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK, 00198 $ INFO ) 00199 * 00200 * Check error code from CLATMS. 00201 * 00202 IF( INFO.NE.0 ) THEN 00203 CALL ALAERH( PATH, 'CLATMS', INFO, 0, UPLO, N, N, -1, 00204 $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) 00205 GO TO 110 00206 END IF 00207 * 00208 * For types 3-5, zero one row and column of the matrix to 00209 * test that INFO is returned correctly. 00210 * 00211 IF( ZEROT ) THEN 00212 IF( IMAT.EQ.3 ) THEN 00213 IZERO = 1 00214 ELSE IF( IMAT.EQ.4 ) THEN 00215 IZERO = N 00216 ELSE 00217 IZERO = N / 2 + 1 00218 END IF 00219 IOFF = ( IZERO-1 )*LDA 00220 * 00221 * Set row and column IZERO of A to 0. 00222 * 00223 IF( IUPLO.EQ.1 ) THEN 00224 DO 20 I = 1, IZERO - 1 00225 A( IOFF+I ) = ZERO 00226 20 CONTINUE 00227 IOFF = IOFF + IZERO 00228 DO 30 I = IZERO, N 00229 A( IOFF ) = ZERO 00230 IOFF = IOFF + LDA 00231 30 CONTINUE 00232 ELSE 00233 IOFF = IZERO 00234 DO 40 I = 1, IZERO - 1 00235 A( IOFF ) = ZERO 00236 IOFF = IOFF + LDA 00237 40 CONTINUE 00238 IOFF = IOFF - IZERO 00239 DO 50 I = IZERO, N 00240 A( IOFF+I ) = ZERO 00241 50 CONTINUE 00242 END IF 00243 ELSE 00244 IZERO = 0 00245 END IF 00246 * 00247 * Set the imaginary part of the diagonals. 00248 * 00249 CALL CLAIPD( N, A, LDA+1, 0 ) 00250 * 00251 * Save a copy of the matrix A in ASAV. 00252 * 00253 CALL CLACPY( UPLO, N, N, A, LDA, ASAV, LDA ) 00254 * 00255 DO 100 IEQUED = 1, 2 00256 EQUED = EQUEDS( IEQUED ) 00257 IF( IEQUED.EQ.1 ) THEN 00258 NFACT = 3 00259 ELSE 00260 NFACT = 1 00261 END IF 00262 * 00263 DO 90 IFACT = 1, NFACT 00264 FACT = FACTS( IFACT ) 00265 PREFAC = LSAME( FACT, 'F' ) 00266 NOFACT = LSAME( FACT, 'N' ) 00267 EQUIL = LSAME( FACT, 'E' ) 00268 * 00269 IF( ZEROT ) THEN 00270 IF( PREFAC ) 00271 $ GO TO 90 00272 RCONDC = ZERO 00273 * 00274 ELSE IF( .NOT.LSAME( FACT, 'N' ) ) THEN 00275 * 00276 * Compute the condition number for comparison with 00277 * the value returned by CPOSVX (FACT = 'N' reuses 00278 * the condition number from the previous iteration 00279 * with FACT = 'F'). 00280 * 00281 CALL CLACPY( UPLO, N, N, ASAV, LDA, AFAC, LDA ) 00282 IF( EQUIL .OR. IEQUED.GT.1 ) THEN 00283 * 00284 * Compute row and column scale factors to 00285 * equilibrate the matrix A. 00286 * 00287 CALL CPOEQU( N, AFAC, LDA, S, SCOND, AMAX, 00288 $ INFO ) 00289 IF( INFO.EQ.0 .AND. N.GT.0 ) THEN 00290 IF( IEQUED.GT.1 ) 00291 $ SCOND = ZERO 00292 * 00293 * Equilibrate the matrix. 00294 * 00295 CALL CLAQHE( UPLO, N, AFAC, LDA, S, SCOND, 00296 $ AMAX, EQUED ) 00297 END IF 00298 END IF 00299 * 00300 * Save the condition number of the 00301 * non-equilibrated system for use in CGET04. 00302 * 00303 IF( EQUIL ) 00304 $ ROLDC = RCONDC 00305 * 00306 * Compute the 1-norm of A. 00307 * 00308 ANORM = CLANHE( '1', UPLO, N, AFAC, LDA, RWORK ) 00309 * 00310 * Factor the matrix A. 00311 * 00312 CALL CPOTRF( UPLO, N, AFAC, LDA, INFO ) 00313 * 00314 * Form the inverse of A. 00315 * 00316 CALL CLACPY( UPLO, N, N, AFAC, LDA, A, LDA ) 00317 CALL CPOTRI( UPLO, N, A, LDA, INFO ) 00318 * 00319 * Compute the 1-norm condition number of A. 00320 * 00321 AINVNM = CLANHE( '1', UPLO, N, A, LDA, RWORK ) 00322 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 00323 RCONDC = ONE 00324 ELSE 00325 RCONDC = ( ONE / ANORM ) / AINVNM 00326 END IF 00327 END IF 00328 * 00329 * Restore the matrix A. 00330 * 00331 CALL CLACPY( UPLO, N, N, ASAV, LDA, A, LDA ) 00332 * 00333 * Form an exact solution and set the right hand side. 00334 * 00335 SRNAMT = 'CLARHS' 00336 CALL CLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU, 00337 $ NRHS, A, LDA, XACT, LDA, B, LDA, 00338 $ ISEED, INFO ) 00339 XTYPE = 'C' 00340 CALL CLACPY( 'Full', N, NRHS, B, LDA, BSAV, LDA ) 00341 * 00342 IF( NOFACT ) THEN 00343 * 00344 * --- Test CPOSV --- 00345 * 00346 * Compute the L*L' or U'*U factorization of the 00347 * matrix and solve the system. 00348 * 00349 CALL CLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) 00350 CALL CLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) 00351 * 00352 SRNAMT = 'CPOSV ' 00353 CALL CPOSV( UPLO, N, NRHS, AFAC, LDA, X, LDA, 00354 $ INFO ) 00355 * 00356 * Check error code from CPOSV . 00357 * 00358 IF( INFO.NE.IZERO ) THEN 00359 CALL ALAERH( PATH, 'CPOSV ', INFO, IZERO, 00360 $ UPLO, N, N, -1, -1, NRHS, IMAT, 00361 $ NFAIL, NERRS, NOUT ) 00362 GO TO 70 00363 ELSE IF( INFO.NE.0 ) THEN 00364 GO TO 70 00365 END IF 00366 * 00367 * Reconstruct matrix from factors and compute 00368 * residual. 00369 * 00370 CALL CPOT01( UPLO, N, A, LDA, AFAC, LDA, RWORK, 00371 $ RESULT( 1 ) ) 00372 * 00373 * Compute residual of the computed solution. 00374 * 00375 CALL CLACPY( 'Full', N, NRHS, B, LDA, WORK, 00376 $ LDA ) 00377 CALL CPOT02( UPLO, N, NRHS, A, LDA, X, LDA, 00378 $ WORK, LDA, RWORK, RESULT( 2 ) ) 00379 * 00380 * Check solution from generated exact solution. 00381 * 00382 CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00383 $ RESULT( 3 ) ) 00384 NT = 3 00385 * 00386 * Print information about the tests that did not 00387 * pass the threshold. 00388 * 00389 DO 60 K = 1, NT 00390 IF( RESULT( K ).GE.THRESH ) THEN 00391 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00392 $ CALL ALADHD( NOUT, PATH ) 00393 WRITE( NOUT, FMT = 9999 )'CPOSV ', UPLO, 00394 $ N, IMAT, K, RESULT( K ) 00395 NFAIL = NFAIL + 1 00396 END IF 00397 60 CONTINUE 00398 NRUN = NRUN + NT 00399 70 CONTINUE 00400 END IF 00401 * 00402 * --- Test CPOSVX --- 00403 * 00404 IF( .NOT.PREFAC ) 00405 $ CALL CLASET( UPLO, N, N, CMPLX( ZERO ), 00406 $ CMPLX( ZERO ), AFAC, LDA ) 00407 CALL CLASET( 'Full', N, NRHS, CMPLX( ZERO ), 00408 $ CMPLX( ZERO ), X, LDA ) 00409 IF( IEQUED.GT.1 .AND. N.GT.0 ) THEN 00410 * 00411 * Equilibrate the matrix if FACT='F' and 00412 * EQUED='Y'. 00413 * 00414 CALL CLAQHE( UPLO, N, A, LDA, S, SCOND, AMAX, 00415 $ EQUED ) 00416 END IF 00417 * 00418 * Solve the system and compute the condition number 00419 * and error bounds using CPOSVX. 00420 * 00421 SRNAMT = 'CPOSVX' 00422 CALL CPOSVX( FACT, UPLO, N, NRHS, A, LDA, AFAC, 00423 $ LDA, EQUED, S, B, LDA, X, LDA, RCOND, 00424 $ RWORK, RWORK( NRHS+1 ), WORK, 00425 $ RWORK( 2*NRHS+1 ), INFO ) 00426 * 00427 * Check the error code from CPOSVX. 00428 * 00429 IF( INFO.NE.IZERO ) THEN 00430 CALL ALAERH( PATH, 'CPOSVX', INFO, IZERO, 00431 $ FACT // UPLO, N, N, -1, -1, NRHS, 00432 $ IMAT, NFAIL, NERRS, NOUT ) 00433 GO TO 90 00434 END IF 00435 * 00436 IF( INFO.EQ.0 ) THEN 00437 IF( .NOT.PREFAC ) THEN 00438 * 00439 * Reconstruct matrix from factors and compute 00440 * residual. 00441 * 00442 CALL CPOT01( UPLO, N, A, LDA, AFAC, LDA, 00443 $ RWORK( 2*NRHS+1 ), RESULT( 1 ) ) 00444 K1 = 1 00445 ELSE 00446 K1 = 2 00447 END IF 00448 * 00449 * Compute residual of the computed solution. 00450 * 00451 CALL CLACPY( 'Full', N, NRHS, BSAV, LDA, WORK, 00452 $ LDA ) 00453 CALL CPOT02( UPLO, N, NRHS, ASAV, LDA, X, LDA, 00454 $ WORK, LDA, RWORK( 2*NRHS+1 ), 00455 $ RESULT( 2 ) ) 00456 * 00457 * Check solution from generated exact solution. 00458 * 00459 IF( NOFACT .OR. ( PREFAC .AND. LSAME( EQUED, 00460 $ 'N' ) ) ) THEN 00461 CALL CGET04( N, NRHS, X, LDA, XACT, LDA, 00462 $ RCONDC, RESULT( 3 ) ) 00463 ELSE 00464 CALL CGET04( N, NRHS, X, LDA, XACT, LDA, 00465 $ ROLDC, RESULT( 3 ) ) 00466 END IF 00467 * 00468 * Check the error bounds from iterative 00469 * refinement. 00470 * 00471 CALL CPOT05( UPLO, N, NRHS, ASAV, LDA, B, LDA, 00472 $ X, LDA, XACT, LDA, RWORK, 00473 $ RWORK( NRHS+1 ), RESULT( 4 ) ) 00474 ELSE 00475 K1 = 6 00476 END IF 00477 * 00478 * Compare RCOND from CPOSVX with the computed value 00479 * in RCONDC. 00480 * 00481 RESULT( 6 ) = SGET06( RCOND, RCONDC ) 00482 * 00483 * Print information about the tests that did not pass 00484 * the threshold. 00485 * 00486 DO 80 K = K1, 6 00487 IF( RESULT( K ).GE.THRESH ) THEN 00488 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00489 $ CALL ALADHD( NOUT, PATH ) 00490 IF( PREFAC ) THEN 00491 WRITE( NOUT, FMT = 9997 )'CPOSVX', FACT, 00492 $ UPLO, N, EQUED, IMAT, K, RESULT( K ) 00493 ELSE 00494 WRITE( NOUT, FMT = 9998 )'CPOSVX', FACT, 00495 $ UPLO, N, IMAT, K, RESULT( K ) 00496 END IF 00497 NFAIL = NFAIL + 1 00498 END IF 00499 80 CONTINUE 00500 NRUN = NRUN + 7 - K1 00501 90 CONTINUE 00502 100 CONTINUE 00503 110 CONTINUE 00504 120 CONTINUE 00505 130 CONTINUE 00506 * 00507 * Print a summary of the results. 00508 * 00509 CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS ) 00510 * 00511 9999 FORMAT( 1X, A, ', UPLO=''', A1, ''', N =', I5, ', type ', I1, 00512 $ ', test(', I1, ')=', G12.5 ) 00513 9998 FORMAT( 1X, A, ', FACT=''', A1, ''', UPLO=''', A1, ''', N=', I5, 00514 $ ', type ', I1, ', test(', I1, ')=', G12.5 ) 00515 9997 FORMAT( 1X, A, ', FACT=''', A1, ''', UPLO=''', A1, ''', N=', I5, 00516 $ ', EQUED=''', A1, ''', type ', I1, ', test(', I1, ') =', 00517 $ G12.5 ) 00518 RETURN 00519 * 00520 * End of CDRVPO 00521 * 00522 END
1.7.3 
