1*> \brief \b SCHKPB 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8* Definition: 9* =========== 10* 11* SUBROUTINE SCHKPB( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, 12* THRESH, TSTERR, NMAX, A, AFAC, AINV, B, X, 13* XACT, WORK, RWORK, IWORK, NOUT ) 14* 15* .. Scalar Arguments .. 16* LOGICAL TSTERR 17* INTEGER NMAX, NN, NNB, NNS, NOUT 18* REAL THRESH 19* .. 20* .. Array Arguments .. 21* LOGICAL DOTYPE( * ) 22* INTEGER IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * ) 23* REAL A( * ), AFAC( * ), AINV( * ), B( * ), 24* $ RWORK( * ), WORK( * ), X( * ), XACT( * ) 25* .. 26* 27* 28*> \par Purpose: 29* ============= 30*> 31*> \verbatim 32*> 33*> SCHKPB tests SPBTRF, -TRS, -RFS, and -CON. 34*> \endverbatim 35* 36* Arguments: 37* ========== 38* 39*> \param[in] DOTYPE 40*> \verbatim 41*> DOTYPE is LOGICAL array, dimension (NTYPES) 42*> The matrix types to be used for testing. Matrices of type j 43*> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = 44*> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. 45*> \endverbatim 46*> 47*> \param[in] NN 48*> \verbatim 49*> NN is INTEGER 50*> The number of values of N contained in the vector NVAL. 51*> \endverbatim 52*> 53*> \param[in] NVAL 54*> \verbatim 55*> NVAL is INTEGER array, dimension (NN) 56*> The values of the matrix dimension N. 57*> \endverbatim 58*> 59*> \param[in] NNB 60*> \verbatim 61*> NNB is INTEGER 62*> The number of values of NB contained in the vector NBVAL. 63*> \endverbatim 64*> 65*> \param[in] NBVAL 66*> \verbatim 67*> NBVAL is INTEGER array, dimension (NNB) 68*> The values of the blocksize NB. 69*> \endverbatim 70*> 71*> \param[in] NNS 72*> \verbatim 73*> NNS is INTEGER 74*> The number of values of NRHS contained in the vector NSVAL. 75*> \endverbatim 76*> 77*> \param[in] NSVAL 78*> \verbatim 79*> NSVAL is INTEGER array, dimension (NNS) 80*> The values of the number of right hand sides NRHS. 81*> \endverbatim 82*> 83*> \param[in] THRESH 84*> \verbatim 85*> THRESH is REAL 86*> The threshold value for the test ratios. A result is 87*> included in the output file if RESULT >= THRESH. To have 88*> every test ratio printed, use THRESH = 0. 89*> \endverbatim 90*> 91*> \param[in] TSTERR 92*> \verbatim 93*> TSTERR is LOGICAL 94*> Flag that indicates whether error exits are to be tested. 95*> \endverbatim 96*> 97*> \param[in] NMAX 98*> \verbatim 99*> NMAX is INTEGER 100*> The maximum value permitted for N, used in dimensioning the 101*> work arrays. 102*> \endverbatim 103*> 104*> \param[out] A 105*> \verbatim 106*> A is REAL array, dimension (NMAX*NMAX) 107*> \endverbatim 108*> 109*> \param[out] AFAC 110*> \verbatim 111*> AFAC is REAL array, dimension (NMAX*NMAX) 112*> \endverbatim 113*> 114*> \param[out] AINV 115*> \verbatim 116*> AINV is REAL array, dimension (NMAX*NMAX) 117*> \endverbatim 118*> 119*> \param[out] B 120*> \verbatim 121*> B is REAL array, dimension (NMAX*NSMAX) 122*> where NSMAX is the largest entry in NSVAL. 123*> \endverbatim 124*> 125*> \param[out] X 126*> \verbatim 127*> X is REAL array, dimension (NMAX*NSMAX) 128*> \endverbatim 129*> 130*> \param[out] XACT 131*> \verbatim 132*> XACT is REAL array, dimension (NMAX*NSMAX) 133*> \endverbatim 134*> 135*> \param[out] WORK 136*> \verbatim 137*> WORK is REAL array, dimension 138*> (NMAX*max(3,NSMAX)) 139*> \endverbatim 140*> 141*> \param[out] RWORK 142*> \verbatim 143*> RWORK is REAL array, dimension 144*> (max(NMAX,2*NSMAX)) 145*> \endverbatim 146*> 147*> \param[out] IWORK 148*> \verbatim 149*> IWORK is INTEGER array, dimension (NMAX) 150*> \endverbatim 151*> 152*> \param[in] NOUT 153*> \verbatim 154*> NOUT is INTEGER 155*> The unit number for output. 156*> \endverbatim 157* 158* Authors: 159* ======== 160* 161*> \author Univ. of Tennessee 162*> \author Univ. of California Berkeley 163*> \author Univ. of Colorado Denver 164*> \author NAG Ltd. 165* 166*> \ingroup single_lin 167* 168* ===================================================================== 169 SUBROUTINE SCHKPB( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, 170 $ THRESH, TSTERR, NMAX, A, AFAC, AINV, B, X, 171 $ XACT, WORK, RWORK, IWORK, NOUT ) 172* 173* -- LAPACK test routine -- 174* -- LAPACK is a software package provided by Univ. of Tennessee, -- 175* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 176* 177* .. Scalar Arguments .. 178 LOGICAL TSTERR 179 INTEGER NMAX, NN, NNB, NNS, NOUT 180 REAL THRESH 181* .. 182* .. Array Arguments .. 183 LOGICAL DOTYPE( * ) 184 INTEGER IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * ) 185 REAL A( * ), AFAC( * ), AINV( * ), B( * ), 186 $ RWORK( * ), WORK( * ), X( * ), XACT( * ) 187* .. 188* 189* ===================================================================== 190* 191* .. Parameters .. 192 REAL ONE, ZERO 193 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 194 INTEGER NTYPES, NTESTS 195 PARAMETER ( NTYPES = 8, NTESTS = 7 ) 196 INTEGER NBW 197 PARAMETER ( NBW = 4 ) 198* .. 199* .. Local Scalars .. 200 LOGICAL ZEROT 201 CHARACTER DIST, PACKIT, TYPE, UPLO, XTYPE 202 CHARACTER*3 PATH 203 INTEGER I, I1, I2, IKD, IMAT, IN, INB, INFO, IOFF, 204 $ IRHS, IUPLO, IW, IZERO, K, KD, KL, KOFF, KU, 205 $ LDA, LDAB, MODE, N, NB, NERRS, NFAIL, NIMAT, 206 $ NKD, NRHS, NRUN 207 REAL AINVNM, ANORM, CNDNUM, RCOND, RCONDC 208* .. 209* .. Local Arrays .. 210 INTEGER ISEED( 4 ), ISEEDY( 4 ), KDVAL( NBW ) 211 REAL RESULT( NTESTS ) 212* .. 213* .. External Functions .. 214 REAL SGET06, SLANGE, SLANSB 215 EXTERNAL SGET06, SLANGE, SLANSB 216* .. 217* .. External Subroutines .. 218 EXTERNAL ALAERH, ALAHD, ALASUM, SCOPY, SERRPO, SGET04, 219 $ SLACPY, SLARHS, SLASET, SLATB4, SLATMS, SPBCON, 220 $ SPBRFS, SPBT01, SPBT02, SPBT05, SPBTRF, SPBTRS, 221 $ SSWAP, XLAENV 222* .. 223* .. Intrinsic Functions .. 224 INTRINSIC MAX, MIN 225* .. 226* .. Scalars in Common .. 227 LOGICAL LERR, OK 228 CHARACTER*32 SRNAMT 229 INTEGER INFOT, NUNIT 230* .. 231* .. Common blocks .. 232 COMMON / INFOC / INFOT, NUNIT, OK, LERR 233 COMMON / SRNAMC / SRNAMT 234* .. 235* .. Data statements .. 236 DATA ISEEDY / 1988, 1989, 1990, 1991 / 237* .. 238* .. Executable Statements .. 239* 240* Initialize constants and the random number seed. 241* 242 PATH( 1: 1 ) = 'Single precision' 243 PATH( 2: 3 ) = 'PB' 244 NRUN = 0 245 NFAIL = 0 246 NERRS = 0 247 DO 10 I = 1, 4 248 ISEED( I ) = ISEEDY( I ) 249 10 CONTINUE 250* 251* Test the error exits 252* 253 IF( TSTERR ) 254 $ CALL SERRPO( PATH, NOUT ) 255 INFOT = 0 256 CALL XLAENV( 2, 2 ) 257 KDVAL( 1 ) = 0 258* 259* Do for each value of N in NVAL 260* 261 DO 90 IN = 1, NN 262 N = NVAL( IN ) 263 LDA = MAX( N, 1 ) 264 XTYPE = 'N' 265* 266* Set limits on the number of loop iterations. 267* 268 NKD = MAX( 1, MIN( N, 4 ) ) 269 NIMAT = NTYPES 270 IF( N.EQ.0 ) 271 $ NIMAT = 1 272* 273 KDVAL( 2 ) = N + ( N+1 ) / 4 274 KDVAL( 3 ) = ( 3*N-1 ) / 4 275 KDVAL( 4 ) = ( N+1 ) / 4 276* 277 DO 80 IKD = 1, NKD 278* 279* Do for KD = 0, (5*N+1)/4, (3N-1)/4, and (N+1)/4. This order 280* makes it easier to skip redundant values for small values 281* of N. 282* 283 KD = KDVAL( IKD ) 284 LDAB = KD + 1 285* 286* Do first for UPLO = 'U', then for UPLO = 'L' 287* 288 DO 70 IUPLO = 1, 2 289 KOFF = 1 290 IF( IUPLO.EQ.1 ) THEN 291 UPLO = 'U' 292 KOFF = MAX( 1, KD+2-N ) 293 PACKIT = 'Q' 294 ELSE 295 UPLO = 'L' 296 PACKIT = 'B' 297 END IF 298* 299 DO 60 IMAT = 1, NIMAT 300* 301* Do the tests only if DOTYPE( IMAT ) is true. 302* 303 IF( .NOT.DOTYPE( IMAT ) ) 304 $ GO TO 60 305* 306* Skip types 2, 3, or 4 if the matrix size is too small. 307* 308 ZEROT = IMAT.GE.2 .AND. IMAT.LE.4 309 IF( ZEROT .AND. N.LT.IMAT-1 ) 310 $ GO TO 60 311* 312 IF( .NOT.ZEROT .OR. .NOT.DOTYPE( 1 ) ) THEN 313* 314* Set up parameters with SLATB4 and generate a test 315* matrix with SLATMS. 316* 317 CALL SLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, 318 $ MODE, CNDNUM, DIST ) 319* 320 SRNAMT = 'SLATMS' 321 CALL SLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, 322 $ CNDNUM, ANORM, KD, KD, PACKIT, 323 $ A( KOFF ), LDAB, WORK, INFO ) 324* 325* Check error code from SLATMS. 326* 327 IF( INFO.NE.0 ) THEN 328 CALL ALAERH( PATH, 'SLATMS', INFO, 0, UPLO, N, 329 $ N, KD, KD, -1, IMAT, NFAIL, NERRS, 330 $ NOUT ) 331 GO TO 60 332 END IF 333 ELSE IF( IZERO.GT.0 ) THEN 334* 335* Use the same matrix for types 3 and 4 as for type 336* 2 by copying back the zeroed out column, 337* 338 IW = 2*LDA + 1 339 IF( IUPLO.EQ.1 ) THEN 340 IOFF = ( IZERO-1 )*LDAB + KD + 1 341 CALL SCOPY( IZERO-I1, WORK( IW ), 1, 342 $ A( IOFF-IZERO+I1 ), 1 ) 343 IW = IW + IZERO - I1 344 CALL SCOPY( I2-IZERO+1, WORK( IW ), 1, 345 $ A( IOFF ), MAX( LDAB-1, 1 ) ) 346 ELSE 347 IOFF = ( I1-1 )*LDAB + 1 348 CALL SCOPY( IZERO-I1, WORK( IW ), 1, 349 $ A( IOFF+IZERO-I1 ), 350 $ MAX( LDAB-1, 1 ) ) 351 IOFF = ( IZERO-1 )*LDAB + 1 352 IW = IW + IZERO - I1 353 CALL SCOPY( I2-IZERO+1, WORK( IW ), 1, 354 $ A( IOFF ), 1 ) 355 END IF 356 END IF 357* 358* For types 2-4, zero one row and column of the matrix 359* to test that INFO is returned correctly. 360* 361 IZERO = 0 362 IF( ZEROT ) THEN 363 IF( IMAT.EQ.2 ) THEN 364 IZERO = 1 365 ELSE IF( IMAT.EQ.3 ) THEN 366 IZERO = N 367 ELSE 368 IZERO = N / 2 + 1 369 END IF 370* 371* Save the zeroed out row and column in WORK(*,3) 372* 373 IW = 2*LDA 374 DO 20 I = 1, MIN( 2*KD+1, N ) 375 WORK( IW+I ) = ZERO 376 20 CONTINUE 377 IW = IW + 1 378 I1 = MAX( IZERO-KD, 1 ) 379 I2 = MIN( IZERO+KD, N ) 380* 381 IF( IUPLO.EQ.1 ) THEN 382 IOFF = ( IZERO-1 )*LDAB + KD + 1 383 CALL SSWAP( IZERO-I1, A( IOFF-IZERO+I1 ), 1, 384 $ WORK( IW ), 1 ) 385 IW = IW + IZERO - I1 386 CALL SSWAP( I2-IZERO+1, A( IOFF ), 387 $ MAX( LDAB-1, 1 ), WORK( IW ), 1 ) 388 ELSE 389 IOFF = ( I1-1 )*LDAB + 1 390 CALL SSWAP( IZERO-I1, A( IOFF+IZERO-I1 ), 391 $ MAX( LDAB-1, 1 ), WORK( IW ), 1 ) 392 IOFF = ( IZERO-1 )*LDAB + 1 393 IW = IW + IZERO - I1 394 CALL SSWAP( I2-IZERO+1, A( IOFF ), 1, 395 $ WORK( IW ), 1 ) 396 END IF 397 END IF 398* 399* Do for each value of NB in NBVAL 400* 401 DO 50 INB = 1, NNB 402 NB = NBVAL( INB ) 403 CALL XLAENV( 1, NB ) 404* 405* Compute the L*L' or U'*U factorization of the band 406* matrix. 407* 408 CALL SLACPY( 'Full', KD+1, N, A, LDAB, AFAC, LDAB ) 409 SRNAMT = 'SPBTRF' 410 CALL SPBTRF( UPLO, N, KD, AFAC, LDAB, INFO ) 411* 412* Check error code from SPBTRF. 413* 414 IF( INFO.NE.IZERO ) THEN 415 CALL ALAERH( PATH, 'SPBTRF', INFO, IZERO, UPLO, 416 $ N, N, KD, KD, NB, IMAT, NFAIL, 417 $ NERRS, NOUT ) 418 GO TO 50 419 END IF 420* 421* Skip the tests if INFO is not 0. 422* 423 IF( INFO.NE.0 ) 424 $ GO TO 50 425* 426*+ TEST 1 427* Reconstruct matrix from factors and compute 428* residual. 429* 430 CALL SLACPY( 'Full', KD+1, N, AFAC, LDAB, AINV, 431 $ LDAB ) 432 CALL SPBT01( UPLO, N, KD, A, LDAB, AINV, LDAB, 433 $ RWORK, RESULT( 1 ) ) 434* 435* Print the test ratio if it is .GE. THRESH. 436* 437 IF( RESULT( 1 ).GE.THRESH ) THEN 438 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 439 $ CALL ALAHD( NOUT, PATH ) 440 WRITE( NOUT, FMT = 9999 )UPLO, N, KD, NB, IMAT, 441 $ 1, RESULT( 1 ) 442 NFAIL = NFAIL + 1 443 END IF 444 NRUN = NRUN + 1 445* 446* Only do other tests if this is the first blocksize. 447* 448 IF( INB.GT.1 ) 449 $ GO TO 50 450* 451* Form the inverse of A so we can get a good estimate 452* of RCONDC = 1/(norm(A) * norm(inv(A))). 453* 454 CALL SLASET( 'Full', N, N, ZERO, ONE, AINV, LDA ) 455 SRNAMT = 'SPBTRS' 456 CALL SPBTRS( UPLO, N, KD, N, AFAC, LDAB, AINV, LDA, 457 $ INFO ) 458* 459* Compute RCONDC = 1/(norm(A) * norm(inv(A))). 460* 461 ANORM = SLANSB( '1', UPLO, N, KD, A, LDAB, RWORK ) 462 AINVNM = SLANGE( '1', N, N, AINV, LDA, RWORK ) 463 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 464 RCONDC = ONE 465 ELSE 466 RCONDC = ( ONE / ANORM ) / AINVNM 467 END IF 468* 469 DO 40 IRHS = 1, NNS 470 NRHS = NSVAL( IRHS ) 471* 472*+ TEST 2 473* Solve and compute residual for A * X = B. 474* 475 SRNAMT = 'SLARHS' 476 CALL SLARHS( PATH, XTYPE, UPLO, ' ', N, N, KD, 477 $ KD, NRHS, A, LDAB, XACT, LDA, B, 478 $ LDA, ISEED, INFO ) 479 CALL SLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) 480* 481 SRNAMT = 'SPBTRS' 482 CALL SPBTRS( UPLO, N, KD, NRHS, AFAC, LDAB, X, 483 $ LDA, INFO ) 484* 485* Check error code from SPBTRS. 486* 487 IF( INFO.NE.0 ) 488 $ CALL ALAERH( PATH, 'SPBTRS', INFO, 0, UPLO, 489 $ N, N, KD, KD, NRHS, IMAT, NFAIL, 490 $ NERRS, NOUT ) 491* 492 CALL SLACPY( 'Full', N, NRHS, B, LDA, WORK, 493 $ LDA ) 494 CALL SPBT02( UPLO, N, KD, NRHS, A, LDAB, X, LDA, 495 $ WORK, LDA, RWORK, RESULT( 2 ) ) 496* 497*+ TEST 3 498* Check solution from generated exact solution. 499* 500 CALL SGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 501 $ RESULT( 3 ) ) 502* 503*+ TESTS 4, 5, and 6 504* Use iterative refinement to improve the solution. 505* 506 SRNAMT = 'SPBRFS' 507 CALL SPBRFS( UPLO, N, KD, NRHS, A, LDAB, AFAC, 508 $ LDAB, B, LDA, X, LDA, RWORK, 509 $ RWORK( NRHS+1 ), WORK, IWORK, 510 $ INFO ) 511* 512* Check error code from SPBRFS. 513* 514 IF( INFO.NE.0 ) 515 $ CALL ALAERH( PATH, 'SPBRFS', INFO, 0, UPLO, 516 $ N, N, KD, KD, NRHS, IMAT, NFAIL, 517 $ NERRS, NOUT ) 518* 519 CALL SGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 520 $ RESULT( 4 ) ) 521 CALL SPBT05( UPLO, N, KD, NRHS, A, LDAB, B, LDA, 522 $ X, LDA, XACT, LDA, RWORK, 523 $ RWORK( NRHS+1 ), RESULT( 5 ) ) 524* 525* Print information about the tests that did not 526* pass the threshold. 527* 528 DO 30 K = 2, 6 529 IF( RESULT( K ).GE.THRESH ) THEN 530 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 531 $ CALL ALAHD( NOUT, PATH ) 532 WRITE( NOUT, FMT = 9998 )UPLO, N, KD, 533 $ NRHS, IMAT, K, RESULT( K ) 534 NFAIL = NFAIL + 1 535 END IF 536 30 CONTINUE 537 NRUN = NRUN + 5 538 40 CONTINUE 539* 540*+ TEST 7 541* Get an estimate of RCOND = 1/CNDNUM. 542* 543 SRNAMT = 'SPBCON' 544 CALL SPBCON( UPLO, N, KD, AFAC, LDAB, ANORM, RCOND, 545 $ WORK, IWORK, INFO ) 546* 547* Check error code from SPBCON. 548* 549 IF( INFO.NE.0 ) 550 $ CALL ALAERH( PATH, 'SPBCON', INFO, 0, UPLO, N, 551 $ N, KD, KD, -1, IMAT, NFAIL, NERRS, 552 $ NOUT ) 553* 554 RESULT( 7 ) = SGET06( RCOND, RCONDC ) 555* 556* Print the test ratio if it is .GE. THRESH. 557* 558 IF( RESULT( 7 ).GE.THRESH ) THEN 559 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 560 $ CALL ALAHD( NOUT, PATH ) 561 WRITE( NOUT, FMT = 9997 )UPLO, N, KD, IMAT, 7, 562 $ RESULT( 7 ) 563 NFAIL = NFAIL + 1 564 END IF 565 NRUN = NRUN + 1 566 50 CONTINUE 567 60 CONTINUE 568 70 CONTINUE 569 80 CONTINUE 570 90 CONTINUE 571* 572* Print a summary of the results. 573* 574 CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS ) 575* 576 9999 FORMAT( ' UPLO=''', A1, ''', N=', I5, ', KD=', I5, ', NB=', I4, 577 $ ', type ', I2, ', test ', I2, ', ratio= ', G12.5 ) 578 9998 FORMAT( ' UPLO=''', A1, ''', N=', I5, ', KD=', I5, ', NRHS=', I3, 579 $ ', type ', I2, ', test(', I2, ') = ', G12.5 ) 580 9997 FORMAT( ' UPLO=''', A1, ''', N=', I5, ', KD=', I5, ',', 10X, 581 $ ' type ', I2, ', test(', I2, ') = ', G12.5 ) 582 RETURN 583* 584* End of SCHKPB 585* 586 END 587