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