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