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 (NBVAL) 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*> \date December 2016 167* 168*> \ingroup single_lin 169* 170* ===================================================================== 171 SUBROUTINE SCHKPB( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, 172 $ THRESH, TSTERR, NMAX, A, AFAC, AINV, B, X, 173 $ XACT, WORK, RWORK, IWORK, NOUT ) 174* 175* -- LAPACK test routine (version 3.7.0) -- 176* -- LAPACK is a software package provided by Univ. of Tennessee, -- 177* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 178* December 2016 179* 180* .. Scalar Arguments .. 181 LOGICAL TSTERR 182 INTEGER NMAX, NN, NNB, NNS, NOUT 183 REAL THRESH 184* .. 185* .. Array Arguments .. 186 LOGICAL DOTYPE( * ) 187 INTEGER IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * ) 188 REAL A( * ), AFAC( * ), AINV( * ), B( * ), 189 $ RWORK( * ), WORK( * ), X( * ), XACT( * ) 190* .. 191* 192* ===================================================================== 193* 194* .. Parameters .. 195 REAL ONE, ZERO 196 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 197 INTEGER NTYPES, NTESTS 198 PARAMETER ( NTYPES = 8, NTESTS = 7 ) 199 INTEGER NBW 200 PARAMETER ( NBW = 4 ) 201* .. 202* .. Local Scalars .. 203 LOGICAL ZEROT 204 CHARACTER DIST, PACKIT, TYPE, UPLO, XTYPE 205 CHARACTER*3 PATH 206 INTEGER I, I1, I2, IKD, IMAT, IN, INB, INFO, IOFF, 207 $ IRHS, IUPLO, IW, IZERO, K, KD, KL, KOFF, KU, 208 $ LDA, LDAB, MODE, N, NB, NERRS, NFAIL, NIMAT, 209 $ NKD, NRHS, NRUN 210 REAL AINVNM, ANORM, CNDNUM, RCOND, RCONDC 211* .. 212* .. Local Arrays .. 213 INTEGER ISEED( 4 ), ISEEDY( 4 ), KDVAL( NBW ) 214 REAL RESULT( NTESTS ) 215* .. 216* .. External Functions .. 217 REAL SGET06, SLANGE, SLANSB 218 EXTERNAL SGET06, SLANGE, SLANSB 219* .. 220* .. External Subroutines .. 221 EXTERNAL ALAERH, ALAHD, ALASUM, SCOPY, SERRPO, SGET04, 222 $ SLACPY, SLARHS, SLASET, SLATB4, SLATMS, SPBCON, 223 $ SPBRFS, SPBT01, SPBT02, SPBT05, SPBTRF, SPBTRS, 224 $ SSWAP, XLAENV 225* .. 226* .. Intrinsic Functions .. 227 INTRINSIC MAX, MIN 228* .. 229* .. Scalars in Common .. 230 LOGICAL LERR, OK 231 CHARACTER*32 SRNAMT 232 INTEGER INFOT, NUNIT 233* .. 234* .. Common blocks .. 235 COMMON / INFOC / INFOT, NUNIT, OK, LERR 236 COMMON / SRNAMC / SRNAMT 237* .. 238* .. Data statements .. 239 DATA ISEEDY / 1988, 1989, 1990, 1991 / 240* .. 241* .. Executable Statements .. 242* 243* Initialize constants and the random number seed. 244* 245 PATH( 1: 1 ) = 'Single precision' 246 PATH( 2: 3 ) = 'PB' 247 NRUN = 0 248 NFAIL = 0 249 NERRS = 0 250 DO 10 I = 1, 4 251 ISEED( I ) = ISEEDY( I ) 252 10 CONTINUE 253* 254* Test the error exits 255* 256 IF( TSTERR ) 257 $ CALL SERRPO( PATH, NOUT ) 258 INFOT = 0 259 CALL XLAENV( 2, 2 ) 260 KDVAL( 1 ) = 0 261* 262* Do for each value of N in NVAL 263* 264 DO 90 IN = 1, NN 265 N = NVAL( IN ) 266 LDA = MAX( N, 1 ) 267 XTYPE = 'N' 268* 269* Set limits on the number of loop iterations. 270* 271 NKD = MAX( 1, MIN( N, 4 ) ) 272 NIMAT = NTYPES 273 IF( N.EQ.0 ) 274 $ NIMAT = 1 275* 276 KDVAL( 2 ) = N + ( N+1 ) / 4 277 KDVAL( 3 ) = ( 3*N-1 ) / 4 278 KDVAL( 4 ) = ( N+1 ) / 4 279* 280 DO 80 IKD = 1, NKD 281* 282* Do for KD = 0, (5*N+1)/4, (3N-1)/4, and (N+1)/4. This order 283* makes it easier to skip redundant values for small values 284* of N. 285* 286 KD = KDVAL( IKD ) 287 LDAB = KD + 1 288* 289* Do first for UPLO = 'U', then for UPLO = 'L' 290* 291 DO 70 IUPLO = 1, 2 292 KOFF = 1 293 IF( IUPLO.EQ.1 ) THEN 294 UPLO = 'U' 295 KOFF = MAX( 1, KD+2-N ) 296 PACKIT = 'Q' 297 ELSE 298 UPLO = 'L' 299 PACKIT = 'B' 300 END IF 301* 302 DO 60 IMAT = 1, NIMAT 303* 304* Do the tests only if DOTYPE( IMAT ) is true. 305* 306 IF( .NOT.DOTYPE( IMAT ) ) 307 $ GO TO 60 308* 309* Skip types 2, 3, or 4 if the matrix size is too small. 310* 311 ZEROT = IMAT.GE.2 .AND. IMAT.LE.4 312 IF( ZEROT .AND. N.LT.IMAT-1 ) 313 $ GO TO 60 314* 315 IF( .NOT.ZEROT .OR. .NOT.DOTYPE( 1 ) ) THEN 316* 317* Set up parameters with SLATB4 and generate a test 318* matrix with SLATMS. 319* 320 CALL SLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, 321 $ MODE, CNDNUM, DIST ) 322* 323 SRNAMT = 'SLATMS' 324 CALL SLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, 325 $ CNDNUM, ANORM, KD, KD, PACKIT, 326 $ A( KOFF ), LDAB, WORK, INFO ) 327* 328* Check error code from SLATMS. 329* 330 IF( INFO.NE.0 ) THEN 331 CALL ALAERH( PATH, 'SLATMS', INFO, 0, UPLO, N, 332 $ N, KD, KD, -1, IMAT, NFAIL, NERRS, 333 $ NOUT ) 334 GO TO 60 335 END IF 336 ELSE IF( IZERO.GT.0 ) THEN 337* 338* Use the same matrix for types 3 and 4 as for type 339* 2 by copying back the zeroed out column, 340* 341 IW = 2*LDA + 1 342 IF( IUPLO.EQ.1 ) THEN 343 IOFF = ( IZERO-1 )*LDAB + KD + 1 344 CALL SCOPY( IZERO-I1, WORK( IW ), 1, 345 $ A( IOFF-IZERO+I1 ), 1 ) 346 IW = IW + IZERO - I1 347 CALL SCOPY( I2-IZERO+1, WORK( IW ), 1, 348 $ A( IOFF ), MAX( LDAB-1, 1 ) ) 349 ELSE 350 IOFF = ( I1-1 )*LDAB + 1 351 CALL SCOPY( IZERO-I1, WORK( IW ), 1, 352 $ A( IOFF+IZERO-I1 ), 353 $ MAX( LDAB-1, 1 ) ) 354 IOFF = ( IZERO-1 )*LDAB + 1 355 IW = IW + IZERO - I1 356 CALL SCOPY( I2-IZERO+1, WORK( IW ), 1, 357 $ A( IOFF ), 1 ) 358 END IF 359 END IF 360* 361* For types 2-4, zero one row and column of the matrix 362* to test that INFO is returned correctly. 363* 364 IZERO = 0 365 IF( ZEROT ) THEN 366 IF( IMAT.EQ.2 ) THEN 367 IZERO = 1 368 ELSE IF( IMAT.EQ.3 ) THEN 369 IZERO = N 370 ELSE 371 IZERO = N / 2 + 1 372 END IF 373* 374* Save the zeroed out row and column in WORK(*,3) 375* 376 IW = 2*LDA 377 DO 20 I = 1, MIN( 2*KD+1, N ) 378 WORK( IW+I ) = ZERO 379 20 CONTINUE 380 IW = IW + 1 381 I1 = MAX( IZERO-KD, 1 ) 382 I2 = MIN( IZERO+KD, N ) 383* 384 IF( IUPLO.EQ.1 ) THEN 385 IOFF = ( IZERO-1 )*LDAB + KD + 1 386 CALL SSWAP( IZERO-I1, A( IOFF-IZERO+I1 ), 1, 387 $ WORK( IW ), 1 ) 388 IW = IW + IZERO - I1 389 CALL SSWAP( I2-IZERO+1, A( IOFF ), 390 $ MAX( LDAB-1, 1 ), WORK( IW ), 1 ) 391 ELSE 392 IOFF = ( I1-1 )*LDAB + 1 393 CALL SSWAP( IZERO-I1, A( IOFF+IZERO-I1 ), 394 $ MAX( LDAB-1, 1 ), WORK( IW ), 1 ) 395 IOFF = ( IZERO-1 )*LDAB + 1 396 IW = IW + IZERO - I1 397 CALL SSWAP( I2-IZERO+1, A( IOFF ), 1, 398 $ WORK( IW ), 1 ) 399 END IF 400 END IF 401* 402* Do for each value of NB in NBVAL 403* 404 DO 50 INB = 1, NNB 405 NB = NBVAL( INB ) 406 CALL XLAENV( 1, NB ) 407* 408* Compute the L*L' or U'*U factorization of the band 409* matrix. 410* 411 CALL SLACPY( 'Full', KD+1, N, A, LDAB, AFAC, LDAB ) 412 SRNAMT = 'SPBTRF' 413 CALL SPBTRF( UPLO, N, KD, AFAC, LDAB, INFO ) 414* 415* Check error code from SPBTRF. 416* 417 IF( INFO.NE.IZERO ) THEN 418 CALL ALAERH( PATH, 'SPBTRF', INFO, IZERO, UPLO, 419 $ N, N, KD, KD, NB, IMAT, NFAIL, 420 $ NERRS, NOUT ) 421 GO TO 50 422 END IF 423* 424* Skip the tests if INFO is not 0. 425* 426 IF( INFO.NE.0 ) 427 $ GO TO 50 428* 429*+ TEST 1 430* Reconstruct matrix from factors and compute 431* residual. 432* 433 CALL SLACPY( 'Full', KD+1, N, AFAC, LDAB, AINV, 434 $ LDAB ) 435 CALL SPBT01( UPLO, N, KD, A, LDAB, AINV, LDAB, 436 $ RWORK, RESULT( 1 ) ) 437* 438* Print the test ratio if it is .GE. THRESH. 439* 440 IF( RESULT( 1 ).GE.THRESH ) THEN 441 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 442 $ CALL ALAHD( NOUT, PATH ) 443 WRITE( NOUT, FMT = 9999 )UPLO, N, KD, NB, IMAT, 444 $ 1, RESULT( 1 ) 445 NFAIL = NFAIL + 1 446 END IF 447 NRUN = NRUN + 1 448* 449* Only do other tests if this is the first blocksize. 450* 451 IF( INB.GT.1 ) 452 $ GO TO 50 453* 454* Form the inverse of A so we can get a good estimate 455* of RCONDC = 1/(norm(A) * norm(inv(A))). 456* 457 CALL SLASET( 'Full', N, N, ZERO, ONE, AINV, LDA ) 458 SRNAMT = 'SPBTRS' 459 CALL SPBTRS( UPLO, N, KD, N, AFAC, LDAB, AINV, LDA, 460 $ INFO ) 461* 462* Compute RCONDC = 1/(norm(A) * norm(inv(A))). 463* 464 ANORM = SLANSB( '1', UPLO, N, KD, A, LDAB, RWORK ) 465 AINVNM = SLANGE( '1', N, N, AINV, LDA, RWORK ) 466 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 467 RCONDC = ONE 468 ELSE 469 RCONDC = ( ONE / ANORM ) / AINVNM 470 END IF 471* 472 DO 40 IRHS = 1, NNS 473 NRHS = NSVAL( IRHS ) 474* 475*+ TEST 2 476* Solve and compute residual for A * X = B. 477* 478 SRNAMT = 'SLARHS' 479 CALL SLARHS( PATH, XTYPE, UPLO, ' ', N, N, KD, 480 $ KD, NRHS, A, LDAB, XACT, LDA, B, 481 $ LDA, ISEED, INFO ) 482 CALL SLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) 483* 484 SRNAMT = 'SPBTRS' 485 CALL SPBTRS( UPLO, N, KD, NRHS, AFAC, LDAB, X, 486 $ LDA, INFO ) 487* 488* Check error code from SPBTRS. 489* 490 IF( INFO.NE.0 ) 491 $ CALL ALAERH( PATH, 'SPBTRS', INFO, 0, UPLO, 492 $ N, N, KD, KD, NRHS, IMAT, NFAIL, 493 $ NERRS, NOUT ) 494* 495 CALL SLACPY( 'Full', N, NRHS, B, LDA, WORK, 496 $ LDA ) 497 CALL SPBT02( UPLO, N, KD, NRHS, A, LDAB, X, LDA, 498 $ WORK, LDA, RWORK, RESULT( 2 ) ) 499* 500*+ TEST 3 501* Check solution from generated exact solution. 502* 503 CALL SGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 504 $ RESULT( 3 ) ) 505* 506*+ TESTS 4, 5, and 6 507* Use iterative refinement to improve the solution. 508* 509 SRNAMT = 'SPBRFS' 510 CALL SPBRFS( UPLO, N, KD, NRHS, A, LDAB, AFAC, 511 $ LDAB, B, LDA, X, LDA, RWORK, 512 $ RWORK( NRHS+1 ), WORK, IWORK, 513 $ INFO ) 514* 515* Check error code from SPBRFS. 516* 517 IF( INFO.NE.0 ) 518 $ CALL ALAERH( PATH, 'SPBRFS', INFO, 0, UPLO, 519 $ N, N, KD, KD, NRHS, IMAT, NFAIL, 520 $ NERRS, NOUT ) 521* 522 CALL SGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 523 $ RESULT( 4 ) ) 524 CALL SPBT05( UPLO, N, KD, NRHS, A, LDAB, B, LDA, 525 $ X, LDA, XACT, LDA, RWORK, 526 $ RWORK( NRHS+1 ), RESULT( 5 ) ) 527* 528* Print information about the tests that did not 529* pass the threshold. 530* 531 DO 30 K = 2, 6 532 IF( RESULT( K ).GE.THRESH ) THEN 533 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 534 $ CALL ALAHD( NOUT, PATH ) 535 WRITE( NOUT, FMT = 9998 )UPLO, N, KD, 536 $ NRHS, IMAT, K, RESULT( K ) 537 NFAIL = NFAIL + 1 538 END IF 539 30 CONTINUE 540 NRUN = NRUN + 5 541 40 CONTINUE 542* 543*+ TEST 7 544* Get an estimate of RCOND = 1/CNDNUM. 545* 546 SRNAMT = 'SPBCON' 547 CALL SPBCON( UPLO, N, KD, AFAC, LDAB, ANORM, RCOND, 548 $ WORK, IWORK, INFO ) 549* 550* Check error code from SPBCON. 551* 552 IF( INFO.NE.0 ) 553 $ CALL ALAERH( PATH, 'SPBCON', INFO, 0, UPLO, N, 554 $ N, KD, KD, -1, IMAT, NFAIL, NERRS, 555 $ NOUT ) 556* 557 RESULT( 7 ) = SGET06( RCOND, RCONDC ) 558* 559* Print the test ratio if it is .GE. THRESH. 560* 561 IF( RESULT( 7 ).GE.THRESH ) THEN 562 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 563 $ CALL ALAHD( NOUT, PATH ) 564 WRITE( NOUT, FMT = 9997 )UPLO, N, KD, IMAT, 7, 565 $ RESULT( 7 ) 566 NFAIL = NFAIL + 1 567 END IF 568 NRUN = NRUN + 1 569 50 CONTINUE 570 60 CONTINUE 571 70 CONTINUE 572 80 CONTINUE 573 90 CONTINUE 574* 575* Print a summary of the results. 576* 577 CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS ) 578* 579 9999 FORMAT( ' UPLO=''', A1, ''', N=', I5, ', KD=', I5, ', NB=', I4, 580 $ ', type ', I2, ', test ', I2, ', ratio= ', G12.5 ) 581 9998 FORMAT( ' UPLO=''', A1, ''', N=', I5, ', KD=', I5, ', NRHS=', I3, 582 $ ', type ', I2, ', test(', I2, ') = ', G12.5 ) 583 9997 FORMAT( ' UPLO=''', A1, ''', N=', I5, ', KD=', I5, ',', 10X, 584 $ ' type ', I2, ', test(', I2, ') = ', G12.5 ) 585 RETURN 586* 587* End of SCHKPB 588* 589 END 590