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