LAPACK: TESTING/LIN/dchkpb.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DCHKPB( 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 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       LOGICAL            TSTERR
00011       INTEGER            NMAX, NN, NNB, NNS, NOUT
00012       DOUBLE PRECISION   THRESH
00013 *     ..
00014 *     .. Array Arguments ..
00015       LOGICAL            DOTYPE( * )
00016       INTEGER            IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * )
00017       DOUBLE PRECISION   A( * ), AFAC( * ), AINV( * ), B( * ),
00018      $                   RWORK( * ), WORK( * ), X( * ), XACT( * )
00019 *     ..
00020 *
00021 *  Purpose
00022 *  =======
00023 *
00024 *  DCHKPB tests DPBTRF, -TRS, -RFS, and -CON.
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 *  NNB     (input) INTEGER
00041 *          The number of values of NB contained in the vector NBVAL.
00042 *
00043 *  NBVAL   (input) INTEGER array, dimension (NBVAL)
00044 *          The values of the blocksize NB.
00045 *
00046 *  NNS     (input) INTEGER
00047 *          The number of values of NRHS contained in the vector NSVAL.
00048 *
00049 *  NSVAL   (input) INTEGER array, dimension (NNS)
00050 *          The values of the number of right hand sides NRHS.
00051 *
00052 *  THRESH  (input) DOUBLE PRECISION
00053 *          The threshold value for the test ratios.  A result is
00054 *          included in the output file if RESULT >= THRESH.  To have
00055 *          every test ratio printed, use THRESH = 0.
00056 *
00057 *  TSTERR  (input) LOGICAL
00058 *          Flag that indicates whether error exits are to be tested.
00059 *
00060 *  NMAX    (input) INTEGER
00061 *          The maximum value permitted for N, used in dimensioning the
00062 *          work arrays.
00063 *
00064 *  A       (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
00065 *
00066 *  AFAC    (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
00067 *
00068 *  AINV    (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
00069 *
00070 *  B       (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
00071 *          where NSMAX is the largest entry in NSVAL.
00072 *
00073 *  X       (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
00074 *
00075 *  XACT    (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
00076 *
00077 *  WORK    (workspace) DOUBLE PRECISION array, dimension
00078 *                      (NMAX*max(3,NSMAX))
00079 *
00080 *  RWORK   (workspace) DOUBLE PRECISION array, dimension
00081 *                      (max(NMAX,2*NSMAX))
00082 *
00083 *  IWORK   (workspace) INTEGER array, dimension (NMAX)
00084 *
00085 *  NOUT    (input) INTEGER
00086 *          The unit number for output.
00087 *
00088 *  =====================================================================
00089 *
00090 *     .. Parameters ..
00091       DOUBLE PRECISION   ONE, ZERO
00092       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00093       INTEGER            NTYPES, NTESTS
00094       PARAMETER          ( NTYPES = 8, NTESTS = 7 )
00095       INTEGER            NBW
00096       PARAMETER          ( NBW = 4 )
00097 *     ..
00098 *     .. Local Scalars ..
00099       LOGICAL            ZEROT
00100       CHARACTER          DIST, PACKIT, TYPE, UPLO, XTYPE
00101       CHARACTER*3        PATH
00102       INTEGER            I, I1, I2, IKD, IMAT, IN, INB, INFO, IOFF,
00103      $                   IRHS, IUPLO, IW, IZERO, K, KD, KL, KOFF, KU,
00104      $                   LDA, LDAB, MODE, N, NB, NERRS, NFAIL, NIMAT,
00105      $                   NKD, NRHS, NRUN
00106       DOUBLE PRECISION   AINVNM, ANORM, CNDNUM, RCOND, RCONDC
00107 *     ..
00108 *     .. Local Arrays ..
00109       INTEGER            ISEED( 4 ), ISEEDY( 4 ), KDVAL( NBW )
00110       DOUBLE PRECISION   RESULT( NTESTS )
00111 *     ..
00112 *     .. External Functions ..
00113       DOUBLE PRECISION   DGET06, DLANGE, DLANSB
00114       EXTERNAL           DGET06, DLANGE, DLANSB
00115 *     ..
00116 *     .. External Subroutines ..
00117       EXTERNAL           ALAERH, ALAHD, ALASUM, DCOPY, DERRPO, DGET04,
00118      $                   DLACPY, DLARHS, DLASET, DLATB4, DLATMS, DPBCON,
00119      $                   DPBRFS, DPBT01, DPBT02, DPBT05, DPBTRF, DPBTRS,
00120      $                   DSWAP, XLAENV
00121 *     ..
00122 *     .. Intrinsic Functions ..
00123       INTRINSIC          MAX, MIN
00124 *     ..
00125 *     .. Scalars in Common ..
00126       LOGICAL            LERR, OK
00127       CHARACTER*32       SRNAMT
00128       INTEGER            INFOT, NUNIT
00129 *     ..
00130 *     .. Common blocks ..
00131       COMMON             / INFOC / INFOT, NUNIT, OK, LERR
00132       COMMON             / SRNAMC / SRNAMT
00133 *     ..
00134 *     .. Data statements ..
00135       DATA               ISEEDY / 1988, 1989, 1990, 1991 /
00136 *     ..
00137 *     .. Executable Statements ..
00138 *
00139 *     Initialize constants and the random number seed.
00140 *
00141       PATH( 1: 1 ) = 'Double precision'
00142       PATH( 2: 3 ) = 'PB'
00143       NRUN = 0
00144       NFAIL = 0
00145       NERRS = 0
00146       DO 10 I = 1, 4
00147          ISEED( I ) = ISEEDY( I )
00148    10 CONTINUE
00149 *
00150 *     Test the error exits
00151 *
00152       IF( TSTERR )
00153      $   CALL DERRPO( PATH, NOUT )
00154       INFOT = 0
00155       CALL XLAENV( 2, 2 )
00156       KDVAL( 1 ) = 0
00157 *
00158 *     Do for each value of N in NVAL
00159 *
00160       DO 90 IN = 1, NN
00161          N = NVAL( IN )
00162          LDA = MAX( N, 1 )
00163          XTYPE = 'N'
00164 *
00165 *        Set limits on the number of loop iterations.
00166 *
00167          NKD = MAX( 1, MIN( N, 4 ) )
00168          NIMAT = NTYPES
00169          IF( N.EQ.0 )
00170      $      NIMAT = 1
00171 *
00172          KDVAL( 2 ) = N + ( N+1 ) / 4
00173          KDVAL( 3 ) = ( 3*N-1 ) / 4
00174          KDVAL( 4 ) = ( N+1 ) / 4
00175 *
00176          DO 80 IKD = 1, NKD
00177 *
00178 *           Do for KD = 0, (5*N+1)/4, (3N-1)/4, and (N+1)/4. This order
00179 *           makes it easier to skip redundant values for small values
00180 *           of N.
00181 *
00182             KD = KDVAL( IKD )
00183             LDAB = KD + 1
00184 *
00185 *           Do first for UPLO = 'U', then for UPLO = 'L'
00186 *
00187             DO 70 IUPLO = 1, 2
00188                KOFF = 1
00189                IF( IUPLO.EQ.1 ) THEN
00190                   UPLO = 'U'
00191                   KOFF = MAX( 1, KD+2-N )
00192                   PACKIT = 'Q'
00193                ELSE
00194                   UPLO = 'L'
00195                   PACKIT = 'B'
00196                END IF
00197 *
00198                DO 60 IMAT = 1, NIMAT
00199 *
00200 *                 Do the tests only if DOTYPE( IMAT ) is true.
00201 *
00202                   IF( .NOT.DOTYPE( IMAT ) )
00203      $               GO TO 60
00204 *
00205 *                 Skip types 2, 3, or 4 if the matrix size is too small.
00206 *
00207                   ZEROT = IMAT.GE.2 .AND. IMAT.LE.4
00208                   IF( ZEROT .AND. N.LT.IMAT-1 )
00209      $               GO TO 60
00210 *
00211                   IF( .NOT.ZEROT .OR. .NOT.DOTYPE( 1 ) ) THEN
00212 *
00213 *                    Set up parameters with DLATB4 and generate a test
00214 *                    matrix with DLATMS.
00215 *
00216                      CALL DLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM,
00217      $                            MODE, CNDNUM, DIST )
00218 *
00219                      SRNAMT = 'DLATMS'
00220                      CALL DLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
00221      $                            CNDNUM, ANORM, KD, KD, PACKIT,
00222      $                            A( KOFF ), LDAB, WORK, INFO )
00223 *
00224 *                    Check error code from DLATMS.
00225 *
00226                      IF( INFO.NE.0 ) THEN
00227                         CALL ALAERH( PATH, 'DLATMS', INFO, 0, UPLO, N,
00228      $                               N, KD, KD, -1, IMAT, NFAIL, NERRS,
00229      $                               NOUT )
00230                         GO TO 60
00231                      END IF
00232                   ELSE IF( IZERO.GT.0 ) THEN
00233 *
00234 *                    Use the same matrix for types 3 and 4 as for type
00235 *                    2 by copying back the zeroed out column,
00236 *
00237                      IW = 2*LDA + 1
00238                      IF( IUPLO.EQ.1 ) THEN
00239                         IOFF = ( IZERO-1 )*LDAB + KD + 1
00240                         CALL DCOPY( IZERO-I1, WORK( IW ), 1,
00241      $                              A( IOFF-IZERO+I1 ), 1 )
00242                         IW = IW + IZERO - I1
00243                         CALL DCOPY( I2-IZERO+1, WORK( IW ), 1,
00244      $                              A( IOFF ), MAX( LDAB-1, 1 ) )
00245                      ELSE
00246                         IOFF = ( I1-1 )*LDAB + 1
00247                         CALL DCOPY( IZERO-I1, WORK( IW ), 1,
00248      $                              A( IOFF+IZERO-I1 ),
00249      $                              MAX( LDAB-1, 1 ) )
00250                         IOFF = ( IZERO-1 )*LDAB + 1
00251                         IW = IW + IZERO - I1
00252                         CALL DCOPY( I2-IZERO+1, WORK( IW ), 1,
00253      $                              A( IOFF ), 1 )
00254                      END IF
00255                   END IF
00256 *
00257 *                 For types 2-4, zero one row and column of the matrix
00258 *                 to test that INFO is returned correctly.
00259 *
00260                   IZERO = 0
00261                   IF( ZEROT ) THEN
00262                      IF( IMAT.EQ.2 ) THEN
00263                         IZERO = 1
00264                      ELSE IF( IMAT.EQ.3 ) THEN
00265                         IZERO = N
00266                      ELSE
00267                         IZERO = N / 2 + 1
00268                      END IF
00269 *
00270 *                    Save the zeroed out row and column in WORK(*,3)
00271 *
00272                      IW = 2*LDA
00273                      DO 20 I = 1, MIN( 2*KD+1, N )
00274                         WORK( IW+I ) = ZERO
00275    20                CONTINUE
00276                      IW = IW + 1
00277                      I1 = MAX( IZERO-KD, 1 )
00278                      I2 = MIN( IZERO+KD, N )
00279 *
00280                      IF( IUPLO.EQ.1 ) THEN
00281                         IOFF = ( IZERO-1 )*LDAB + KD + 1
00282                         CALL DSWAP( IZERO-I1, A( IOFF-IZERO+I1 ), 1,
00283      $                              WORK( IW ), 1 )
00284                         IW = IW + IZERO - I1
00285                         CALL DSWAP( I2-IZERO+1, A( IOFF ),
00286      $                              MAX( LDAB-1, 1 ), WORK( IW ), 1 )
00287                      ELSE
00288                         IOFF = ( I1-1 )*LDAB + 1
00289                         CALL DSWAP( IZERO-I1, A( IOFF+IZERO-I1 ),
00290      $                              MAX( LDAB-1, 1 ), WORK( IW ), 1 )
00291                         IOFF = ( IZERO-1 )*LDAB + 1
00292                         IW = IW + IZERO - I1
00293                         CALL DSWAP( I2-IZERO+1, A( IOFF ), 1,
00294      $                              WORK( IW ), 1 )
00295                      END IF
00296                   END IF
00297 *
00298 *                 Do for each value of NB in NBVAL
00299 *
00300                   DO 50 INB = 1, NNB
00301                      NB = NBVAL( INB )
00302                      CALL XLAENV( 1, NB )
00303 *
00304 *                    Compute the L*L' or U'*U factorization of the band
00305 *                    matrix.
00306 *
00307                      CALL DLACPY( 'Full', KD+1, N, A, LDAB, AFAC, LDAB )
00308                      SRNAMT = 'DPBTRF'
00309                      CALL DPBTRF( UPLO, N, KD, AFAC, LDAB, INFO )
00310 *
00311 *                    Check error code from DPBTRF.
00312 *
00313                      IF( INFO.NE.IZERO ) THEN
00314                         CALL ALAERH( PATH, 'DPBTRF', INFO, IZERO, UPLO,
00315      $                               N, N, KD, KD, NB, IMAT, NFAIL,
00316      $                               NERRS, NOUT )
00317                         GO TO 50
00318                      END IF
00319 *
00320 *                    Skip the tests if INFO is not 0.
00321 *
00322                      IF( INFO.NE.0 )
00323      $                  GO TO 50
00324 *
00325 *+    TEST 1
00326 *                    Reconstruct matrix from factors and compute
00327 *                    residual.
00328 *
00329                      CALL DLACPY( 'Full', KD+1, N, AFAC, LDAB, AINV,
00330      $                            LDAB )
00331                      CALL DPBT01( UPLO, N, KD, A, LDAB, AINV, LDAB,
00332      $                            RWORK, RESULT( 1 ) )
00333 *
00334 *                    Print the test ratio if it is .GE. THRESH.
00335 *
00336                      IF( RESULT( 1 ).GE.THRESH ) THEN
00337                         IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00338      $                     CALL ALAHD( NOUT, PATH )
00339                         WRITE( NOUT, FMT = 9999 )UPLO, N, KD, NB, IMAT,
00340      $                     1, RESULT( 1 )
00341                         NFAIL = NFAIL + 1
00342                      END IF
00343                      NRUN = NRUN + 1
00344 *
00345 *                    Only do other tests if this is the first blocksize.
00346 *
00347                      IF( INB.GT.1 )
00348      $                  GO TO 50
00349 *
00350 *                    Form the inverse of A so we can get a good estimate
00351 *                    of RCONDC = 1/(norm(A) * norm(inv(A))).
00352 *
00353                      CALL DLASET( 'Full', N, N, ZERO, ONE, AINV, LDA )
00354                      SRNAMT = 'DPBTRS'
00355                      CALL DPBTRS( UPLO, N, KD, N, AFAC, LDAB, AINV, LDA,
00356      $                            INFO )
00357 *
00358 *                    Compute RCONDC = 1/(norm(A) * norm(inv(A))).
00359 *
00360                      ANORM = DLANSB( '1', UPLO, N, KD, A, LDAB, RWORK )
00361                      AINVNM = DLANGE( '1', N, N, AINV, LDA, RWORK )
00362                      IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
00363                         RCONDC = ONE
00364                      ELSE
00365                         RCONDC = ( ONE / ANORM ) / AINVNM
00366                      END IF
00367 *
00368                      DO 40 IRHS = 1, NNS
00369                         NRHS = NSVAL( IRHS )
00370 *
00371 *+    TEST 2
00372 *                    Solve and compute residual for A * X = B.
00373 *
00374                         SRNAMT = 'DLARHS'
00375                         CALL DLARHS( PATH, XTYPE, UPLO, ' ', N, N, KD,
00376      $                               KD, NRHS, A, LDAB, XACT, LDA, B,
00377      $                               LDA, ISEED, INFO )
00378                         CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
00379 *
00380                         SRNAMT = 'DPBTRS'
00381                         CALL DPBTRS( UPLO, N, KD, NRHS, AFAC, LDAB, X,
00382      $                               LDA, INFO )
00383 *
00384 *                    Check error code from DPBTRS.
00385 *
00386                         IF( INFO.NE.0 )
00387      $                     CALL ALAERH( PATH, 'DPBTRS', INFO, 0, UPLO,
00388      $                                  N, N, KD, KD, NRHS, IMAT, NFAIL,
00389      $                                  NERRS, NOUT )
00390 *
00391                         CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK,
00392      $                               LDA )
00393                         CALL DPBT02( UPLO, N, KD, NRHS, A, LDAB, X, LDA,
00394      $                               WORK, LDA, RWORK, RESULT( 2 ) )
00395 *
00396 *+    TEST 3
00397 *                    Check solution from generated exact solution.
00398 *
00399                         CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
00400      $                               RESULT( 3 ) )
00401 *
00402 *+    TESTS 4, 5, and 6
00403 *                    Use iterative refinement to improve the solution.
00404 *
00405                         SRNAMT = 'DPBRFS'
00406                         CALL DPBRFS( UPLO, N, KD, NRHS, A, LDAB, AFAC,
00407      $                               LDAB, B, LDA, X, LDA, RWORK,
00408      $                               RWORK( NRHS+1 ), WORK, IWORK,
00409      $                               INFO )
00410 *
00411 *                    Check error code from DPBRFS.
00412 *
00413                         IF( INFO.NE.0 )
00414      $                     CALL ALAERH( PATH, 'DPBRFS', INFO, 0, UPLO,
00415      $                                  N, N, KD, KD, NRHS, IMAT, NFAIL,
00416      $                                  NERRS, NOUT )
00417 *
00418                         CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
00419      $                               RESULT( 4 ) )
00420                         CALL DPBT05( UPLO, N, KD, NRHS, A, LDAB, B, LDA,
00421      $                               X, LDA, XACT, LDA, RWORK,
00422      $                               RWORK( NRHS+1 ), RESULT( 5 ) )
00423 *
00424 *                       Print information about the tests that did not
00425 *                       pass the threshold.
00426 *
00427                         DO 30 K = 2, 6
00428                            IF( RESULT( K ).GE.THRESH ) THEN
00429                               IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00430      $                           CALL ALAHD( NOUT, PATH )
00431                               WRITE( NOUT, FMT = 9998 )UPLO, N, KD,
00432      $                           NRHS, IMAT, K, RESULT( K )
00433                               NFAIL = NFAIL + 1
00434                            END IF
00435    30                   CONTINUE
00436                         NRUN = NRUN + 5
00437    40                CONTINUE
00438 *
00439 *+    TEST 7
00440 *                    Get an estimate of RCOND = 1/CNDNUM.
00441 *
00442                      SRNAMT = 'DPBCON'
00443                      CALL DPBCON( UPLO, N, KD, AFAC, LDAB, ANORM, RCOND,
00444      $                            WORK, IWORK, INFO )
00445 *
00446 *                    Check error code from DPBCON.
00447 *
00448                      IF( INFO.NE.0 )
00449      $                  CALL ALAERH( PATH, 'DPBCON', INFO, 0, UPLO, N,
00450      $                               N, KD, KD, -1, IMAT, NFAIL, NERRS,
00451      $                               NOUT )
00452 *
00453                      RESULT( 7 ) = DGET06( RCOND, RCONDC )
00454 *
00455 *                    Print the test ratio if it is .GE. THRESH.
00456 *
00457                      IF( RESULT( 7 ).GE.THRESH ) THEN
00458                         IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00459      $                     CALL ALAHD( NOUT, PATH )
00460                         WRITE( NOUT, FMT = 9997 )UPLO, N, KD, IMAT, 7,
00461      $                     RESULT( 7 )
00462                         NFAIL = NFAIL + 1
00463                      END IF
00464                      NRUN = NRUN + 1
00465    50             CONTINUE
00466    60          CONTINUE
00467    70       CONTINUE
00468    80    CONTINUE
00469    90 CONTINUE
00470 *
00471 *     Print a summary of the results.
00472 *
00473       CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
00474 *
00475  9999 FORMAT( ' UPLO=''', A1, ''', N=', I5, ', KD=', I5, ', NB=', I4,
00476      $      ', type ', I2, ', test ', I2, ', ratio= ', G12.5 )
00477  9998 FORMAT( ' UPLO=''', A1, ''', N=', I5, ', KD=', I5, ', NRHS=', I3,
00478      $      ', type ', I2, ', test(', I2, ') = ', G12.5 )
00479  9997 FORMAT( ' UPLO=''', A1, ''', N=', I5, ', KD=', I5, ',', 10X,
00480      $      ' type ', I2, ', test(', I2, ') = ', G12.5 )
00481       RETURN
00482 *
00483 *     End of DCHKPB
00484 *
00485       END