 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DDRVRFP( NOUT, NN, NVAL, NNS, NSVAL, NNT, NTVAL, 00002 + THRESH, A, ASAV, AFAC, AINV, B, 00003 + BSAV, XACT, X, ARF, ARFINV, 00004 + D_WORK_DLATMS, D_WORK_DPOT01, D_TEMP_DPOT02, 00005 + D_TEMP_DPOT03, D_WORK_DLANSY, 00006 + D_WORK_DPOT02, D_WORK_DPOT03 ) 00007 * 00008 * -- LAPACK test routine (version 3.2.0) -- 00009 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00010 * November 2008 00011 * 00012 * .. Scalar Arguments .. 00013 INTEGER NN, NNS, NNT, NOUT 00014 DOUBLE PRECISION THRESH 00015 * .. 00016 * .. Array Arguments .. 00017 INTEGER NVAL( NN ), NSVAL( NNS ), NTVAL( NNT ) 00018 DOUBLE PRECISION A( * ) 00019 DOUBLE PRECISION AINV( * ) 00020 DOUBLE PRECISION ASAV( * ) 00021 DOUBLE PRECISION B( * ) 00022 DOUBLE PRECISION BSAV( * ) 00023 DOUBLE PRECISION AFAC( * ) 00024 DOUBLE PRECISION ARF( * ) 00025 DOUBLE PRECISION ARFINV( * ) 00026 DOUBLE PRECISION XACT( * ) 00027 DOUBLE PRECISION X( * ) 00028 DOUBLE PRECISION D_WORK_DLATMS( * ) 00029 DOUBLE PRECISION D_WORK_DPOT01( * ) 00030 DOUBLE PRECISION D_TEMP_DPOT02( * ) 00031 DOUBLE PRECISION D_TEMP_DPOT03( * ) 00032 DOUBLE PRECISION D_WORK_DLANSY( * ) 00033 DOUBLE PRECISION D_WORK_DPOT02( * ) 00034 DOUBLE PRECISION D_WORK_DPOT03( * ) 00035 * .. 00036 * 00037 * Purpose 00038 * ======= 00039 * 00040 * DDRVRFP tests the LAPACK RFP routines: 00041 * DPFTRF, DPFTRS, and DPFTRI. 00042 * 00043 * This testing routine follow the same tests as DDRVPO (test for the full 00044 * format Symmetric Positive Definite solver). 00045 * 00046 * The tests are performed in Full Format, convertion back and forth from 00047 * full format to RFP format are performed using the routines DTRTTF and 00048 * DTFTTR. 00049 * 00050 * First, a specific matrix A of size N is created. There is nine types of 00051 * different matrixes possible. 00052 * 1. Diagonal 6. Random, CNDNUM = sqrt(0.1/EPS) 00053 * 2. Random, CNDNUM = 2 7. Random, CNDNUM = 0.1/EPS 00054 * *3. First row and column zero 8. Scaled near underflow 00055 * *4. Last row and column zero 9. Scaled near overflow 00056 * *5. Middle row and column zero 00057 * (* - tests error exits from DPFTRF, no test ratios are computed) 00058 * A solution XACT of size N-by-NRHS is created and the associated right 00059 * hand side B as well. Then DPFTRF is called to compute L (or U), the 00060 * Cholesky factor of A. Then L (or U) is used to solve the linear system 00061 * of equations AX = B. This gives X. Then L (or U) is used to compute the 00062 * inverse of A, AINV. The following four tests are then performed: 00063 * (1) norm( L*L' - A ) / ( N * norm(A) * EPS ) or 00064 * norm( U'*U - A ) / ( N * norm(A) * EPS ), 00065 * (2) norm(B - A*X) / ( norm(A) * norm(X) * EPS ), 00066 * (3) norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ), 00067 * (4) ( norm(X-XACT) * RCOND ) / ( norm(XACT) * EPS ), 00068 * where EPS is the machine precision, RCOND the condition number of A, and 00069 * norm( . ) the 1-norm for (1,2,3) and the inf-norm for (4). 00070 * Errors occur when INFO parameter is not as expected. Failures occur when 00071 * a test ratios is greater than THRES. 00072 * 00073 * Arguments 00074 * ========= 00075 * 00076 * NOUT (input) INTEGER 00077 * The unit number for output. 00078 * 00079 * NN (input) INTEGER 00080 * The number of values of N contained in the vector NVAL. 00081 * 00082 * NVAL (input) INTEGER array, dimension (NN) 00083 * The values of the matrix dimension N. 00084 * 00085 * NNS (input) INTEGER 00086 * The number of values of NRHS contained in the vector NSVAL. 00087 * 00088 * NSVAL (input) INTEGER array, dimension (NNS) 00089 * The values of the number of right-hand sides NRHS. 00090 * 00091 * NNT (input) INTEGER 00092 * The number of values of MATRIX TYPE contained in the vector NTVAL. 00093 * 00094 * NTVAL (input) INTEGER array, dimension (NNT) 00095 * The values of matrix type (between 0 and 9 for PO/PP/PF matrices). 00096 * 00097 * THRESH (input) DOUBLE PRECISION 00098 * The threshold value for the test ratios. A result is 00099 * included in the output file if RESULT >= THRESH. To have 00100 * every test ratio printed, use THRESH = 0. 00101 * 00102 * A (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX) 00103 * 00104 * ASAV (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX) 00105 * 00106 * AFAC (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX) 00107 * 00108 * AINV (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX) 00109 * 00110 * B (workspace) DOUBLE PRECISION array, dimension (NMAX*MAXRHS) 00111 * 00112 * BSAV (workspace) DOUBLE PRECISION array, dimension (NMAX*MAXRHS) 00113 * 00114 * XACT (workspace) DOUBLE PRECISION array, dimension (NMAX*MAXRHS) 00115 * 00116 * X (workspace) DOUBLE PRECISION array, dimension (NMAX*MAXRHS) 00117 * 00118 * ARF (workspace) DOUBLE PRECISION array, dimension ((NMAX*(NMAX+1))/2) 00119 * 00120 * ARFINV (workspace) DOUBLE PRECISION array, dimension ((NMAX*(NMAX+1))/2) 00121 * 00122 * D_WORK_DLATMS (workspace) DOUBLE PRECISION array, dimension ( 3*NMAX ) 00123 * 00124 * D_WORK_DPOT01 (workspace) DOUBLE PRECISION array, dimension ( NMAX ) 00125 * 00126 * D_TEMP_DPOT02 (workspace) DOUBLE PRECISION array, dimension ( NMAX*MAXRHS ) 00127 * 00128 * D_TEMP_DPOT03 (workspace) DOUBLE PRECISION array, dimension ( NMAX*NMAX ) 00129 * 00130 * D_WORK_DLATMS (workspace) DOUBLE PRECISION array, dimension ( NMAX ) 00131 * 00132 * D_WORK_DLANSY (workspace) DOUBLE PRECISION array, dimension ( NMAX ) 00133 * 00134 * D_WORK_DPOT02 (workspace) DOUBLE PRECISION array, dimension ( NMAX ) 00135 * 00136 * D_WORK_DPOT03 (workspace) DOUBLE PRECISION array, dimension ( NMAX ) 00137 * 00138 * ===================================================================== 00139 * 00140 * .. Parameters .. 00141 DOUBLE PRECISION ONE, ZERO 00142 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00143 INTEGER NTESTS 00144 PARAMETER ( NTESTS = 4 ) 00145 * .. 00146 * .. Local Scalars .. 00147 LOGICAL ZEROT 00148 INTEGER I, INFO, IUPLO, LDA, LDB, IMAT, NERRS, NFAIL, 00149 + NRHS, NRUN, IZERO, IOFF, K, NT, N, IFORM, IIN, 00150 + IIT, IIS 00151 CHARACTER DIST, CTYPE, UPLO, CFORM 00152 INTEGER KL, KU, MODE 00153 DOUBLE PRECISION ANORM, AINVNM, CNDNUM, RCONDC 00154 * .. 00155 * .. Local Arrays .. 00156 CHARACTER UPLOS( 2 ), FORMS( 2 ) 00157 INTEGER ISEED( 4 ), ISEEDY( 4 ) 00158 DOUBLE PRECISION RESULT( NTESTS ) 00159 * .. 00160 * .. External Functions .. 00161 DOUBLE PRECISION DLANSY 00162 EXTERNAL DLANSY 00163 * .. 00164 * .. External Subroutines .. 00165 EXTERNAL ALADHD, ALAERH, ALASVM, DGET04, DTFTTR, DLACPY, 00166 + DLARHS, DLATB4, DLATMS, DPFTRI, DPFTRF, DPFTRS, 00167 + DPOT01, DPOT02, DPOT03, DPOTRI, DPOTRF, DTRTTF 00168 * .. 00169 * .. Scalars in Common .. 00170 CHARACTER*32 SRNAMT 00171 * .. 00172 * .. Common blocks .. 00173 COMMON / SRNAMC / SRNAMT 00174 * .. 00175 * .. Data statements .. 00176 DATA ISEEDY / 1988, 1989, 1990, 1991 / 00177 DATA UPLOS / 'U', 'L' / 00178 DATA FORMS / 'N', 'T' / 00179 * .. 00180 * .. Executable Statements .. 00181 * 00182 * Initialize constants and the random number seed. 00183 * 00184 NRUN = 0 00185 NFAIL = 0 00186 NERRS = 0 00187 DO 10 I = 1, 4 00188 ISEED( I ) = ISEEDY( I ) 00189 10 CONTINUE 00190 * 00191 DO 130 IIN = 1, NN 00192 * 00193 N = NVAL( IIN ) 00194 LDA = MAX( N, 1 ) 00195 LDB = MAX( N, 1 ) 00196 * 00197 DO 980 IIS = 1, NNS 00198 * 00199 NRHS = NSVAL( IIS ) 00200 * 00201 DO 120 IIT = 1, NNT 00202 * 00203 IMAT = NTVAL( IIT ) 00204 * 00205 * If N.EQ.0, only consider the first type 00206 * 00207 IF( N.EQ.0 .AND. IIT.GT.1 ) GO TO 120 00208 * 00209 * Skip types 3, 4, or 5 if the matrix size is too small. 00210 * 00211 IF( IMAT.EQ.4 .AND. N.LE.1 ) GO TO 120 00212 IF( IMAT.EQ.5 .AND. N.LE.2 ) GO TO 120 00213 * 00214 * Do first for UPLO = 'U', then for UPLO = 'L' 00215 * 00216 DO 110 IUPLO = 1, 2 00217 UPLO = UPLOS( IUPLO ) 00218 * 00219 * Do first for CFORM = 'N', then for CFORM = 'C' 00220 * 00221 DO 100 IFORM = 1, 2 00222 CFORM = FORMS( IFORM ) 00223 * 00224 * Set up parameters with DLATB4 and generate a test 00225 * matrix with DLATMS. 00226 * 00227 CALL DLATB4( 'DPO', IMAT, N, N, CTYPE, KL, KU, 00228 + ANORM, MODE, CNDNUM, DIST ) 00229 * 00230 SRNAMT = 'DLATMS' 00231 CALL DLATMS( N, N, DIST, ISEED, CTYPE, 00232 + D_WORK_DLATMS, 00233 + MODE, CNDNUM, ANORM, KL, KU, UPLO, A, 00234 + LDA, D_WORK_DLATMS, INFO ) 00235 * 00236 * Check error code from DLATMS. 00237 * 00238 IF( INFO.NE.0 ) THEN 00239 CALL ALAERH( 'DPF', 'DLATMS', INFO, 0, UPLO, N, 00240 + N, -1, -1, -1, IIT, NFAIL, NERRS, 00241 + NOUT ) 00242 GO TO 100 00243 END IF 00244 * 00245 * For types 3-5, zero one row and column of the matrix to 00246 * test that INFO is returned correctly. 00247 * 00248 ZEROT = IMAT.GE.3 .AND. IMAT.LE.5 00249 IF( ZEROT ) THEN 00250 IF( IIT.EQ.3 ) THEN 00251 IZERO = 1 00252 ELSE IF( IIT.EQ.4 ) THEN 00253 IZERO = N 00254 ELSE 00255 IZERO = N / 2 + 1 00256 END IF 00257 IOFF = ( IZERO-1 )*LDA 00258 * 00259 * Set row and column IZERO of A to 0. 00260 * 00261 IF( IUPLO.EQ.1 ) THEN 00262 DO 20 I = 1, IZERO - 1 00263 A( IOFF+I ) = ZERO 00264 20 CONTINUE 00265 IOFF = IOFF + IZERO 00266 DO 30 I = IZERO, N 00267 A( IOFF ) = ZERO 00268 IOFF = IOFF + LDA 00269 30 CONTINUE 00270 ELSE 00271 IOFF = IZERO 00272 DO 40 I = 1, IZERO - 1 00273 A( IOFF ) = ZERO 00274 IOFF = IOFF + LDA 00275 40 CONTINUE 00276 IOFF = IOFF - IZERO 00277 DO 50 I = IZERO, N 00278 A( IOFF+I ) = ZERO 00279 50 CONTINUE 00280 END IF 00281 ELSE 00282 IZERO = 0 00283 END IF 00284 * 00285 * Save a copy of the matrix A in ASAV. 00286 * 00287 CALL DLACPY( UPLO, N, N, A, LDA, ASAV, LDA ) 00288 * 00289 * Compute the condition number of A (RCONDC). 00290 * 00291 IF( ZEROT ) THEN 00292 RCONDC = ZERO 00293 ELSE 00294 * 00295 * Compute the 1-norm of A. 00296 * 00297 ANORM = DLANSY( '1', UPLO, N, A, LDA, 00298 + D_WORK_DLANSY ) 00299 * 00300 * Factor the matrix A. 00301 * 00302 CALL DPOTRF( UPLO, N, A, LDA, INFO ) 00303 * 00304 * Form the inverse of A. 00305 * 00306 CALL DPOTRI( UPLO, N, A, LDA, INFO ) 00307 * 00308 * Compute the 1-norm condition number of A. 00309 * 00310 AINVNM = DLANSY( '1', UPLO, N, A, LDA, 00311 + D_WORK_DLANSY ) 00312 RCONDC = ( ONE / ANORM ) / AINVNM 00313 * 00314 * Restore the matrix A. 00315 * 00316 CALL DLACPY( UPLO, N, N, ASAV, LDA, A, LDA ) 00317 * 00318 END IF 00319 * 00320 * Form an exact solution and set the right hand side. 00321 * 00322 SRNAMT = 'DLARHS' 00323 CALL DLARHS( 'DPO', 'N', UPLO, ' ', N, N, KL, KU, 00324 + NRHS, A, LDA, XACT, LDA, B, LDA, 00325 + ISEED, INFO ) 00326 CALL DLACPY( 'Full', N, NRHS, B, LDA, BSAV, LDA ) 00327 * 00328 * Compute the L*L' or U'*U factorization of the 00329 * matrix and solve the system. 00330 * 00331 CALL DLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) 00332 CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDB ) 00333 * 00334 SRNAMT = 'DTRTTF' 00335 CALL DTRTTF( CFORM, UPLO, N, AFAC, LDA, ARF, INFO ) 00336 SRNAMT = 'DPFTRF' 00337 CALL DPFTRF( CFORM, UPLO, N, ARF, INFO ) 00338 * 00339 * Check error code from DPFTRF. 00340 * 00341 IF( INFO.NE.IZERO ) THEN 00342 * 00343 * LANGOU: there is a small hick here: IZERO should 00344 * always be INFO however if INFO is ZERO, ALAERH does not 00345 * complain. 00346 * 00347 CALL ALAERH( 'DPF', 'DPFSV ', INFO, IZERO, 00348 + UPLO, N, N, -1, -1, NRHS, IIT, 00349 + NFAIL, NERRS, NOUT ) 00350 GO TO 100 00351 END IF 00352 * 00353 * Skip the tests if INFO is not 0. 00354 * 00355 IF( INFO.NE.0 ) THEN 00356 GO TO 100 00357 END IF 00358 * 00359 SRNAMT = 'DPFTRS' 00360 CALL DPFTRS( CFORM, UPLO, N, NRHS, ARF, X, LDB, 00361 + INFO ) 00362 * 00363 SRNAMT = 'DTFTTR' 00364 CALL DTFTTR( CFORM, UPLO, N, ARF, AFAC, LDA, INFO ) 00365 * 00366 * Reconstruct matrix from factors and compute 00367 * residual. 00368 * 00369 CALL DLACPY( UPLO, N, N, AFAC, LDA, ASAV, LDA ) 00370 CALL DPOT01( UPLO, N, A, LDA, AFAC, LDA, 00371 + D_WORK_DPOT01, RESULT( 1 ) ) 00372 CALL DLACPY( UPLO, N, N, ASAV, LDA, AFAC, LDA ) 00373 * 00374 * Form the inverse and compute the residual. 00375 * 00376 IF(MOD(N,2).EQ.0)THEN 00377 CALL DLACPY( 'A', N+1, N/2, ARF, N+1, ARFINV, 00378 + N+1 ) 00379 ELSE 00380 CALL DLACPY( 'A', N, (N+1)/2, ARF, N, ARFINV, 00381 + N ) 00382 END IF 00383 * 00384 SRNAMT = 'DPFTRI' 00385 CALL DPFTRI( CFORM, UPLO, N, ARFINV , INFO ) 00386 * 00387 SRNAMT = 'DTFTTR' 00388 CALL DTFTTR( CFORM, UPLO, N, ARFINV, AINV, LDA, 00389 + INFO ) 00390 * 00391 * Check error code from DPFTRI. 00392 * 00393 IF( INFO.NE.0 ) 00394 + CALL ALAERH( 'DPO', 'DPFTRI', INFO, 0, UPLO, N, 00395 + N, -1, -1, -1, IMAT, NFAIL, NERRS, 00396 + NOUT ) 00397 * 00398 CALL DPOT03( UPLO, N, A, LDA, AINV, LDA, 00399 + D_TEMP_DPOT03, LDA, D_WORK_DPOT03, 00400 + RCONDC, RESULT( 2 ) ) 00401 * 00402 * Compute residual of the computed solution. 00403 * 00404 CALL DLACPY( 'Full', N, NRHS, B, LDA, 00405 + D_TEMP_DPOT02, LDA ) 00406 CALL DPOT02( UPLO, N, NRHS, A, LDA, X, LDA, 00407 + D_TEMP_DPOT02, LDA, D_WORK_DPOT02, 00408 + RESULT( 3 ) ) 00409 * 00410 * Check solution from generated exact solution. 00411 00412 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00413 + RESULT( 4 ) ) 00414 NT = 4 00415 * 00416 * Print information about the tests that did not 00417 * pass the threshold. 00418 * 00419 DO 60 K = 1, NT 00420 IF( RESULT( K ).GE.THRESH ) THEN 00421 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00422 + CALL ALADHD( NOUT, 'DPF' ) 00423 WRITE( NOUT, FMT = 9999 )'DPFSV ', UPLO, 00424 + N, IIT, K, RESULT( K ) 00425 NFAIL = NFAIL + 1 00426 END IF 00427 60 CONTINUE 00428 NRUN = NRUN + NT 00429 100 CONTINUE 00430 110 CONTINUE 00431 120 CONTINUE 00432 980 CONTINUE 00433 130 CONTINUE 00434 * 00435 * Print a summary of the results. 00436 * 00437 CALL ALASVM( 'DPF', NOUT, NFAIL, NRUN, NERRS ) 00438 * 00439 9999 FORMAT( 1X, A6, ', UPLO=''', A1, ''', N =', I5, ', type ', I1, 00440 + ', test(', I1, ')=', G12.5 ) 00441 * 00442 RETURN 00443 * 00444 * End of DDRVRFP 00445 * 00446 END
1.7.3 
