 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CCHKSY( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, 00002 $ THRESH, TSTERR, NMAX, A, AFAC, AINV, B, X, 00003 $ XACT, WORK, RWORK, IWORK, NOUT ) 00004 * 00005 * -- LAPACK test routine (version 3.1) -- 00006 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00007 * June 2010 00008 * 00009 * .. Scalar Arguments .. 00010 LOGICAL TSTERR 00011 INTEGER NMAX, NN, NNB, NNS, NOUT 00012 REAL THRESH 00013 * .. 00014 * .. Array Arguments .. 00015 LOGICAL DOTYPE( * ) 00016 INTEGER IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * ) 00017 REAL RWORK( * ) 00018 COMPLEX A( * ), AFAC( * ), AINV( * ), B( * ), 00019 $ WORK( * ), X( * ), XACT( * ) 00020 * .. 00021 * 00022 * Purpose 00023 * ======= 00024 * 00025 * CCHKSY tests CSYTRF, -TRI2, -TRS, -TRS2, -RFS, and -CON. 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 * NNB (input) INTEGER 00042 * The number of values of NB contained in the vector NBVAL. 00043 * 00044 * NBVAL (input) INTEGER array, dimension (NBVAL) 00045 * The values of the blocksize NB. 00046 * 00047 * NNS (input) INTEGER 00048 * The number of values of NRHS contained in the vector NSVAL. 00049 * 00050 * NSVAL (input) INTEGER array, dimension (NNS) 00051 * The values of the number of right hand sides NRHS. 00052 * 00053 * THRESH (input) REAL 00054 * The threshold value for the test ratios. A result is 00055 * included in the output file if RESULT >= THRESH. To have 00056 * every test ratio printed, use THRESH = 0. 00057 * 00058 * TSTERR (input) LOGICAL 00059 * Flag that indicates whether error exits are to be tested. 00060 * 00061 * NMAX (input) INTEGER 00062 * The maximum value permitted for N, used in dimensioning the 00063 * work arrays. 00064 * 00065 * A (workspace) COMPLEX array, dimension (NMAX*NMAX) 00066 * 00067 * AFAC (workspace) COMPLEX array, dimension (NMAX*NMAX) 00068 * 00069 * AINV (workspace) COMPLEX array, dimension (NMAX*NMAX) 00070 * 00071 * B (workspace) COMPLEX array, dimension (NMAX*NSMAX) 00072 * where NSMAX is the largest entry in NSVAL. 00073 * 00074 * X (workspace) COMPLEX array, dimension (NMAX*NSMAX) 00075 * 00076 * XACT (workspace) COMPLEX array, dimension (NMAX*NSMAX) 00077 * 00078 * WORK (workspace) COMPLEX array, dimension 00079 * (NMAX*max(2,NSMAX)) 00080 * 00081 * RWORK (workspace) REAL array, 00082 * dimension (NMAX+2*NSMAX) 00083 * 00084 * IWORK (workspace) INTEGER array, dimension (NMAX) 00085 * 00086 * NOUT (input) INTEGER 00087 * The unit number for output. 00088 * 00089 * ===================================================================== 00090 * 00091 * .. Parameters .. 00092 REAL ZERO 00093 PARAMETER ( ZERO = 0.0E+0 ) 00094 INTEGER NTYPES 00095 PARAMETER ( NTYPES = 11 ) 00096 INTEGER NTESTS 00097 PARAMETER ( NTESTS = 9 ) 00098 * .. 00099 * .. Local Scalars .. 00100 LOGICAL TRFCON, ZEROT 00101 CHARACTER DIST, TYPE, UPLO, XTYPE 00102 CHARACTER*3 PATH 00103 INTEGER I, I1, I2, IMAT, IN, INB, INFO, IOFF, IRHS, 00104 $ IUPLO, IZERO, J, K, KL, KU, LDA, LWORK, MODE, 00105 $ N, NB, NERRS, NFAIL, NIMAT, NRHS, NRUN, NT 00106 REAL ANORM, CNDNUM, RCOND, RCONDC 00107 * .. 00108 * .. Local Arrays .. 00109 CHARACTER UPLOS( 2 ) 00110 INTEGER ISEED( 4 ), ISEEDY( 4 ) 00111 REAL RESULT( NTESTS ) 00112 * .. 00113 * .. External Functions .. 00114 REAL CLANSY, SGET06 00115 EXTERNAL CLANSY, SGET06 00116 * .. 00117 * .. External Subroutines .. 00118 EXTERNAL ALAERH, ALAHD, ALASUM, CERRSY, CGET04, CLACPY, 00119 $ CLARHS, CLATB4, CLATMS, CLATSY, CPOT05, CSYCON, 00120 $ CSYRFS, CSYT01, CSYT02, CSYT03, CSYTRF, 00121 $ CSYTRI2, CSYTRS, XLAENV 00122 * .. 00123 * .. Intrinsic Functions .. 00124 INTRINSIC MAX, MIN 00125 * .. 00126 * .. Scalars in Common .. 00127 LOGICAL LERR, OK 00128 CHARACTER*32 SRNAMT 00129 INTEGER INFOT, NUNIT 00130 * .. 00131 * .. Common blocks .. 00132 COMMON / INFOC / INFOT, NUNIT, OK, LERR 00133 COMMON / SRNAMC / SRNAMT 00134 * .. 00135 * .. Data statements .. 00136 DATA ISEEDY / 1988, 1989, 1990, 1991 / 00137 DATA UPLOS / 'U', 'L' / 00138 * .. 00139 * .. Executable Statements .. 00140 * 00141 * Initialize constants and the random number seed. 00142 * 00143 PATH( 1: 1 ) = 'Complex precision' 00144 PATH( 2: 3 ) = 'SY' 00145 NRUN = 0 00146 NFAIL = 0 00147 NERRS = 0 00148 DO 10 I = 1, 4 00149 ISEED( I ) = ISEEDY( I ) 00150 10 CONTINUE 00151 * 00152 * Test the error exits 00153 * 00154 IF( TSTERR ) 00155 $ CALL CERRSY( PATH, NOUT ) 00156 INFOT = 0 00157 * 00158 * Do for each value of N in NVAL 00159 * 00160 DO 180 IN = 1, NN 00161 N = NVAL( IN ) 00162 LDA = MAX( N, 1 ) 00163 XTYPE = 'N' 00164 NIMAT = NTYPES 00165 IF( N.LE.0 ) 00166 $ NIMAT = 1 00167 * 00168 IZERO = 0 00169 DO 170 IMAT = 1, NIMAT 00170 * 00171 * Do the tests only if DOTYPE( IMAT ) is true. 00172 * 00173 IF( .NOT.DOTYPE( IMAT ) ) 00174 $ GO TO 170 00175 * 00176 * Skip types 3, 4, 5, or 6 if the matrix size is too small. 00177 * 00178 ZEROT = IMAT.GE.3 .AND. IMAT.LE.6 00179 IF( ZEROT .AND. N.LT.IMAT-2 ) 00180 $ GO TO 170 00181 * 00182 * Do first for UPLO = 'U', then for UPLO = 'L' 00183 * 00184 DO 160 IUPLO = 1, 2 00185 UPLO = UPLOS( IUPLO ) 00186 * 00187 IF( IMAT.NE.NTYPES ) THEN 00188 * 00189 * Set up parameters with CLATB4 and generate a test 00190 * matrix with CLATMS. 00191 * 00192 CALL CLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, 00193 $ MODE, CNDNUM, DIST ) 00194 * 00195 SRNAMT = 'CLATMS' 00196 CALL CLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, 00197 $ CNDNUM, ANORM, KL, KU, 'N', 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, 00204 $ -1, -1, -1, IMAT, NFAIL, NERRS, NOUT ) 00205 GO TO 160 00206 END IF 00207 * 00208 * For types 3-6, zero one or more rows and columns of 00209 * the matrix to 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 * 00220 IF( IMAT.LT.6 ) THEN 00221 * 00222 * Set row and column IZERO to zero. 00223 * 00224 IF( IUPLO.EQ.1 ) THEN 00225 IOFF = ( IZERO-1 )*LDA 00226 DO 20 I = 1, IZERO - 1 00227 A( IOFF+I ) = ZERO 00228 20 CONTINUE 00229 IOFF = IOFF + IZERO 00230 DO 30 I = IZERO, N 00231 A( IOFF ) = ZERO 00232 IOFF = IOFF + LDA 00233 30 CONTINUE 00234 ELSE 00235 IOFF = IZERO 00236 DO 40 I = 1, IZERO - 1 00237 A( IOFF ) = ZERO 00238 IOFF = IOFF + LDA 00239 40 CONTINUE 00240 IOFF = IOFF - IZERO 00241 DO 50 I = IZERO, N 00242 A( IOFF+I ) = ZERO 00243 50 CONTINUE 00244 END IF 00245 ELSE 00246 IF( IUPLO.EQ.1 ) THEN 00247 * 00248 * Set the first IZERO rows to zero. 00249 * 00250 IOFF = 0 00251 DO 70 J = 1, N 00252 I2 = MIN( J, IZERO ) 00253 DO 60 I = 1, I2 00254 A( IOFF+I ) = ZERO 00255 60 CONTINUE 00256 IOFF = IOFF + LDA 00257 70 CONTINUE 00258 ELSE 00259 * 00260 * Set the last IZERO rows to zero. 00261 * 00262 IOFF = 0 00263 DO 90 J = 1, N 00264 I1 = MAX( J, IZERO ) 00265 DO 80 I = I1, N 00266 A( IOFF+I ) = ZERO 00267 80 CONTINUE 00268 IOFF = IOFF + LDA 00269 90 CONTINUE 00270 END IF 00271 END IF 00272 ELSE 00273 IZERO = 0 00274 END IF 00275 ELSE 00276 * 00277 * Use a special block diagonal matrix to test alternate 00278 * code for the 2 x 2 blocks. 00279 * 00280 CALL CLATSY( UPLO, N, A, LDA, ISEED ) 00281 END IF 00282 * 00283 * Do for each value of NB in NBVAL 00284 * 00285 DO 150 INB = 1, NNB 00286 NB = NBVAL( INB ) 00287 CALL XLAENV( 1, NB ) 00288 * 00289 * Compute the L*D*L' or U*D*U' factorization of the 00290 * matrix. 00291 * 00292 CALL CLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) 00293 LWORK = MAX( 2, NB )*LDA 00294 SRNAMT = 'CSYTRF' 00295 CALL CSYTRF( UPLO, N, AFAC, LDA, IWORK, AINV, LWORK, 00296 $ INFO ) 00297 * 00298 * Adjust the expected value of INFO to account for 00299 * pivoting. 00300 * 00301 K = IZERO 00302 IF( K.GT.0 ) THEN 00303 100 CONTINUE 00304 IF( IWORK( K ).LT.0 ) THEN 00305 IF( IWORK( K ).NE.-K ) THEN 00306 K = -IWORK( K ) 00307 GO TO 100 00308 END IF 00309 ELSE IF( IWORK( K ).NE.K ) THEN 00310 K = IWORK( K ) 00311 GO TO 100 00312 END IF 00313 END IF 00314 * 00315 * Check error code from CSYTRF. 00316 * 00317 IF( INFO.NE.K ) 00318 $ CALL ALAERH( PATH, 'CSYTRF', INFO, K, UPLO, N, N, 00319 $ -1, -1, NB, IMAT, NFAIL, NERRS, NOUT ) 00320 IF( INFO.NE.0 ) THEN 00321 TRFCON = .TRUE. 00322 ELSE 00323 TRFCON = .FALSE. 00324 END IF 00325 * 00326 *+ TEST 1 00327 * Reconstruct matrix from factors and compute residual. 00328 * 00329 CALL CSYT01( UPLO, N, A, LDA, AFAC, LDA, IWORK, AINV, 00330 $ LDA, RWORK, RESULT( 1 ) ) 00331 NT = 1 00332 * 00333 *+ TEST 2 00334 * Form the inverse and compute the residual. 00335 * 00336 IF( INB.EQ.1 .AND. .NOT.TRFCON ) THEN 00337 CALL CLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA ) 00338 SRNAMT = 'CSYTRI2' 00339 LWORK = (N+NB+1)*(NB+3) 00340 CALL CSYTRI2( UPLO, N, AINV, LDA, IWORK, WORK, 00341 $ LWORK, INFO ) 00342 * 00343 * Check error code from CSYTRI. 00344 * 00345 IF( INFO.NE.0 ) 00346 $ CALL ALAERH( PATH, 'CSYTRI', INFO, 0, UPLO, N, 00347 $ N, -1, -1, -1, IMAT, NFAIL, NERRS, 00348 $ NOUT ) 00349 * 00350 CALL CSYT03( UPLO, N, A, LDA, AINV, LDA, WORK, LDA, 00351 $ RWORK, RCONDC, RESULT( 2 ) ) 00352 NT = 2 00353 END IF 00354 * 00355 * Print information about the tests that did not pass 00356 * the threshold. 00357 * 00358 DO 110 K = 1, NT 00359 IF( RESULT( K ).GE.THRESH ) THEN 00360 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00361 $ CALL ALAHD( NOUT, PATH ) 00362 WRITE( NOUT, FMT = 9999 )UPLO, N, NB, IMAT, K, 00363 $ RESULT( K ) 00364 NFAIL = NFAIL + 1 00365 END IF 00366 110 CONTINUE 00367 NRUN = NRUN + NT 00368 * 00369 * Skip the other tests if this is not the first block 00370 * size. 00371 * 00372 IF( INB.GT.1 ) 00373 $ GO TO 150 00374 * 00375 * Do only the condition estimate if INFO is not 0. 00376 * 00377 IF( TRFCON ) THEN 00378 RCONDC = ZERO 00379 GO TO 140 00380 END IF 00381 * 00382 DO 130 IRHS = 1, NNS 00383 NRHS = NSVAL( IRHS ) 00384 * 00385 *+ TEST 3 00386 * Solve and compute residual for A * X = B. 00387 * 00388 SRNAMT = 'CLARHS' 00389 CALL CLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU, 00390 $ NRHS, A, LDA, XACT, LDA, B, LDA, 00391 $ ISEED, INFO ) 00392 CALL CLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) 00393 * 00394 SRNAMT = 'CSYTRS' 00395 CALL CSYTRS( UPLO, N, NRHS, AFAC, LDA, IWORK, X, 00396 $ LDA, INFO ) 00397 * 00398 * Check error code from CSYTRS. 00399 * 00400 IF( INFO.NE.0 ) 00401 $ CALL ALAERH( PATH, 'CSYTRS', INFO, 0, UPLO, N, 00402 $ N, -1, -1, NRHS, IMAT, NFAIL, 00403 $ NERRS, NOUT ) 00404 * 00405 CALL CLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) 00406 CALL CSYT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK, 00407 $ LDA, RWORK, RESULT( 3 ) ) 00408 * 00409 *+ TEST 4 00410 * Solve and compute residual for A * X = B. 00411 * 00412 SRNAMT = 'CLARHS' 00413 CALL CLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU, 00414 $ NRHS, A, LDA, XACT, LDA, B, LDA, 00415 $ ISEED, INFO ) 00416 CALL CLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) 00417 * 00418 SRNAMT = 'CSYTRS2' 00419 CALL CSYTRS2( UPLO, N, NRHS, AFAC, LDA, IWORK, X, 00420 $ LDA, WORK, INFO ) 00421 * 00422 * Check error code from CSYTRS2. 00423 * 00424 IF( INFO.NE.0 ) 00425 $ CALL ALAERH( PATH, 'CSYTRS2', INFO, 0, UPLO, N, 00426 $ N, -1, -1, NRHS, IMAT, NFAIL, 00427 $ NERRS, NOUT ) 00428 * 00429 CALL CLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) 00430 CALL CSYT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK, 00431 $ LDA, RWORK, RESULT( 4 ) ) 00432 * 00433 *+ TEST 5 00434 * Check solution from generated exact solution. 00435 * 00436 CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00437 $ RESULT( 5 ) ) 00438 * 00439 *+ TESTS 6, 7, and 8 00440 * Use iterative refinement to improve the solution. 00441 * 00442 SRNAMT = 'CSYRFS' 00443 CALL CSYRFS( UPLO, N, NRHS, A, LDA, AFAC, LDA, 00444 $ IWORK, B, LDA, X, LDA, RWORK, 00445 $ RWORK( NRHS+1 ), WORK, 00446 $ RWORK( 2*NRHS+1 ), INFO ) 00447 * 00448 * Check error code from CSYRFS. 00449 * 00450 IF( INFO.NE.0 ) 00451 $ CALL ALAERH( PATH, 'CSYRFS', INFO, 0, UPLO, N, 00452 $ N, -1, -1, NRHS, IMAT, NFAIL, 00453 $ NERRS, NOUT ) 00454 * 00455 CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00456 $ RESULT( 6 ) ) 00457 CALL CPOT05( UPLO, N, NRHS, A, LDA, B, LDA, X, LDA, 00458 $ XACT, LDA, RWORK, RWORK( NRHS+1 ), 00459 $ RESULT( 7 ) ) 00460 * 00461 * Print information about the tests that did not pass 00462 * the threshold. 00463 * 00464 DO 120 K = 3, 8 00465 IF( RESULT( K ).GE.THRESH ) THEN 00466 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00467 $ CALL ALAHD( NOUT, PATH ) 00468 WRITE( NOUT, FMT = 9998 )UPLO, N, NRHS, 00469 $ IMAT, K, RESULT( K ) 00470 NFAIL = NFAIL + 1 00471 END IF 00472 120 CONTINUE 00473 NRUN = NRUN + 5 00474 130 CONTINUE 00475 * 00476 *+ TEST 9 00477 * Get an estimate of RCOND = 1/CNDNUM. 00478 * 00479 140 CONTINUE 00480 ANORM = CLANSY( '1', UPLO, N, A, LDA, RWORK ) 00481 SRNAMT = 'CSYCON' 00482 CALL CSYCON( UPLO, N, AFAC, LDA, IWORK, ANORM, RCOND, 00483 $ WORK, INFO ) 00484 * 00485 * Check error code from CSYCON. 00486 * 00487 IF( INFO.NE.0 ) 00488 $ CALL ALAERH( PATH, 'CSYCON', INFO, 0, UPLO, N, N, 00489 $ -1, -1, -1, IMAT, NFAIL, NERRS, NOUT ) 00490 * 00491 RESULT( 9 ) = SGET06( RCOND, RCONDC ) 00492 * 00493 * Print information about the tests that did not pass 00494 * the threshold. 00495 * 00496 IF( RESULT( 9 ).GE.THRESH ) THEN 00497 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00498 $ CALL ALAHD( NOUT, PATH ) 00499 WRITE( NOUT, FMT = 9997 )UPLO, N, IMAT, 9, 00500 $ RESULT( 9 ) 00501 NFAIL = NFAIL + 1 00502 END IF 00503 NRUN = NRUN + 1 00504 150 CONTINUE 00505 160 CONTINUE 00506 170 CONTINUE 00507 180 CONTINUE 00508 * 00509 * Print a summary of the results. 00510 * 00511 CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS ) 00512 * 00513 9999 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NB =', I4, ', type ', 00514 $ I2, ', test ', I2, ', ratio =', G12.5 ) 00515 9998 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NRHS=', I3, ', type ', 00516 $ I2, ', test(', I2, ') =', G12.5 ) 00517 9997 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ',', 10X, ' type ', I2, 00518 $ ', test(', I2, ') =', G12.5 ) 00519 RETURN 00520 * 00521 * End of CCHKSY 00522 * 00523 END
1.7.3 
