LAPACK: TESTING/LIN/dchkpp.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DCHKPP( DOTYPE, NN, NVAL, NNS, NSVAL, THRESH, TSTERR,
00002      $                   NMAX, A, AFAC, AINV, B, X, XACT, WORK, RWORK,
00003      $                   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, NNS, NOUT
00012       DOUBLE PRECISION   THRESH
00013 *     ..
00014 *     .. Array Arguments ..
00015       LOGICAL            DOTYPE( * )
00016       INTEGER            IWORK( * ), NSVAL( * ), NVAL( * )
00017       DOUBLE PRECISION   A( * ), AFAC( * ), AINV( * ), B( * ),
00018      $                   RWORK( * ), WORK( * ), X( * ), XACT( * )
00019 *     ..
00020 *
00021 *  Purpose
00022 *  =======
00023 *
00024 *  DCHKPP tests DPPTRF, -TRI, -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 *  NNS     (input) INTEGER
00041 *          The number of values of NRHS contained in the vector NSVAL.
00042 *
00043 *  NSVAL   (input) INTEGER array, dimension (NNS)
00044 *          The values of the number of right hand sides NRHS.
00045 *
00046 *  THRESH  (input) DOUBLE PRECISION
00047 *          The threshold value for the test ratios.  A result is
00048 *          included in the output file if RESULT >= THRESH.  To have
00049 *          every test ratio printed, use THRESH = 0.
00050 *
00051 *  TSTERR  (input) LOGICAL
00052 *          Flag that indicates whether error exits are to be tested.
00053 *
00054 *  NMAX    (input) INTEGER
00055 *          The maximum value permitted for N, used in dimensioning the
00056 *          work arrays.
00057 *
00058 *  A       (workspace) DOUBLE PRECISION array, dimension
00059 *                      (NMAX*(NMAX+1)/2)
00060 *
00061 *  AFAC    (workspace) DOUBLE PRECISION array, dimension
00062 *                      (NMAX*(NMAX+1)/2)
00063 *
00064 *  AINV    (workspace) DOUBLE PRECISION array, dimension
00065 *                      (NMAX*(NMAX+1)/2)
00066 *
00067 *  B       (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
00068 *          where NSMAX is the largest entry in NSVAL.
00069 *
00070 *  X       (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
00071 *
00072 *  XACT    (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
00073 *
00074 *  WORK    (workspace) DOUBLE PRECISION array, dimension
00075 *                      (NMAX*max(3,NSMAX))
00076 *
00077 *  RWORK   (workspace) DOUBLE PRECISION array, dimension
00078 *                      (max(NMAX,2*NSMAX))
00079 *
00080 *  IWORK   (workspace) INTEGER array, dimension (NMAX)
00081 *
00082 *  NOUT    (input) INTEGER
00083 *          The unit number for output.
00084 *
00085 *  =====================================================================
00086 *
00087 *     .. Parameters ..
00088       DOUBLE PRECISION   ZERO
00089       PARAMETER          ( ZERO = 0.0D+0 )
00090       INTEGER            NTYPES
00091       PARAMETER          ( NTYPES = 9 )
00092       INTEGER            NTESTS
00093       PARAMETER          ( NTESTS = 8 )
00094 *     ..
00095 *     .. Local Scalars ..
00096       LOGICAL            ZEROT
00097       CHARACTER          DIST, PACKIT, TYPE, UPLO, XTYPE
00098       CHARACTER*3        PATH
00099       INTEGER            I, IMAT, IN, INFO, IOFF, IRHS, IUPLO, IZERO, K,
00100      $                   KL, KU, LDA, MODE, N, NERRS, NFAIL, NIMAT, NPP,
00101      $                   NRHS, NRUN
00102       DOUBLE PRECISION   ANORM, CNDNUM, RCOND, RCONDC
00103 *     ..
00104 *     .. Local Arrays ..
00105       CHARACTER          PACKS( 2 ), UPLOS( 2 )
00106       INTEGER            ISEED( 4 ), ISEEDY( 4 )
00107       DOUBLE PRECISION   RESULT( NTESTS )
00108 *     ..
00109 *     .. External Functions ..
00110       DOUBLE PRECISION   DGET06, DLANSP
00111       EXTERNAL           DGET06, DLANSP
00112 *     ..
00113 *     .. External Subroutines ..
00114       EXTERNAL           ALAERH, ALAHD, ALASUM, DCOPY, DERRPO, DGET04,
00115      $                   DLACPY, DLARHS, DLATB4, DLATMS, DPPCON, DPPRFS,
00116      $                   DPPT01, DPPT02, DPPT03, DPPT05, DPPTRF, DPPTRI,
00117      $                   DPPTRS
00118 *     ..
00119 *     .. Scalars in Common ..
00120       LOGICAL            LERR, OK
00121       CHARACTER*32       SRNAMT
00122       INTEGER            INFOT, NUNIT
00123 *     ..
00124 *     .. Common blocks ..
00125       COMMON             / INFOC / INFOT, NUNIT, OK, LERR
00126       COMMON             / SRNAMC / SRNAMT
00127 *     ..
00128 *     .. Intrinsic Functions ..
00129       INTRINSIC          MAX
00130 *     ..
00131 *     .. Data statements ..
00132       DATA               ISEEDY / 1988, 1989, 1990, 1991 /
00133       DATA               UPLOS / 'U', 'L' / , PACKS / 'C', 'R' /
00134 *     ..
00135 *     .. Executable Statements ..
00136 *
00137 *     Initialize constants and the random number seed.
00138 *
00139       PATH( 1: 1 ) = 'Double precision'
00140       PATH( 2: 3 ) = 'PP'
00141       NRUN = 0
00142       NFAIL = 0
00143       NERRS = 0
00144       DO 10 I = 1, 4
00145          ISEED( I ) = ISEEDY( I )
00146    10 CONTINUE
00147 *
00148 *     Test the error exits
00149 *
00150       IF( TSTERR )
00151      $   CALL DERRPO( PATH, NOUT )
00152       INFOT = 0
00153 *
00154 *     Do for each value of N in NVAL
00155 *
00156       DO 110 IN = 1, NN
00157          N = NVAL( IN )
00158          LDA = MAX( N, 1 )
00159          XTYPE = 'N'
00160          NIMAT = NTYPES
00161          IF( N.LE.0 )
00162      $      NIMAT = 1
00163 *
00164          DO 100 IMAT = 1, NIMAT
00165 *
00166 *           Do the tests only if DOTYPE( IMAT ) is true.
00167 *
00168             IF( .NOT.DOTYPE( IMAT ) )
00169      $         GO TO 100
00170 *
00171 *           Skip types 3, 4, or 5 if the matrix size is too small.
00172 *
00173             ZEROT = IMAT.GE.3 .AND. IMAT.LE.5
00174             IF( ZEROT .AND. N.LT.IMAT-2 )
00175      $         GO TO 100
00176 *
00177 *           Do first for UPLO = 'U', then for UPLO = 'L'
00178 *
00179             DO 90 IUPLO = 1, 2
00180                UPLO = UPLOS( IUPLO )
00181                PACKIT = PACKS( IUPLO )
00182 *
00183 *              Set up parameters with DLATB4 and generate a test matrix
00184 *              with DLATMS.
00185 *
00186                CALL DLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
00187      $                      CNDNUM, DIST )
00188 *
00189                SRNAMT = 'DLATMS'
00190                CALL DLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
00191      $                      CNDNUM, ANORM, KL, KU, PACKIT, A, LDA, WORK,
00192      $                      INFO )
00193 *
00194 *              Check error code from DLATMS.
00195 *
00196                IF( INFO.NE.0 ) THEN
00197                   CALL ALAERH( PATH, 'DLATMS', INFO, 0, UPLO, N, N, -1,
00198      $                         -1, -1, IMAT, NFAIL, NERRS, NOUT )
00199                   GO TO 90
00200                END IF
00201 *
00202 *              For types 3-5, zero one row and column of the matrix to
00203 *              test that INFO is returned correctly.
00204 *
00205                IF( ZEROT ) THEN
00206                   IF( IMAT.EQ.3 ) THEN
00207                      IZERO = 1
00208                   ELSE IF( IMAT.EQ.4 ) THEN
00209                      IZERO = N
00210                   ELSE
00211                      IZERO = N / 2 + 1
00212                   END IF
00213 *
00214 *                 Set row and column IZERO of A to 0.
00215 *
00216                   IF( IUPLO.EQ.1 ) THEN
00217                      IOFF = ( IZERO-1 )*IZERO / 2
00218                      DO 20 I = 1, IZERO - 1
00219                         A( IOFF+I ) = ZERO
00220    20                CONTINUE
00221                      IOFF = IOFF + IZERO
00222                      DO 30 I = IZERO, N
00223                         A( IOFF ) = ZERO
00224                         IOFF = IOFF + I
00225    30                CONTINUE
00226                   ELSE
00227                      IOFF = IZERO
00228                      DO 40 I = 1, IZERO - 1
00229                         A( IOFF ) = ZERO
00230                         IOFF = IOFF + N - I
00231    40                CONTINUE
00232                      IOFF = IOFF - IZERO
00233                      DO 50 I = IZERO, N
00234                         A( IOFF+I ) = ZERO
00235    50                CONTINUE
00236                   END IF
00237                ELSE
00238                   IZERO = 0
00239                END IF
00240 *
00241 *              Compute the L*L' or U'*U factorization of the matrix.
00242 *
00243                NPP = N*( N+1 ) / 2
00244                CALL DCOPY( NPP, A, 1, AFAC, 1 )
00245                SRNAMT = 'DPPTRF'
00246                CALL DPPTRF( UPLO, N, AFAC, INFO )
00247 *
00248 *              Check error code from DPPTRF.
00249 *
00250                IF( INFO.NE.IZERO ) THEN
00251                   CALL ALAERH( PATH, 'DPPTRF', INFO, IZERO, UPLO, N, N,
00252      $                         -1, -1, -1, IMAT, NFAIL, NERRS, NOUT )
00253                   GO TO 90
00254                END IF
00255 *
00256 *              Skip the tests if INFO is not 0.
00257 *
00258                IF( INFO.NE.0 )
00259      $            GO TO 90
00260 *
00261 *+    TEST 1
00262 *              Reconstruct matrix from factors and compute residual.
00263 *
00264                CALL DCOPY( NPP, AFAC, 1, AINV, 1 )
00265                CALL DPPT01( UPLO, N, A, AINV, RWORK, RESULT( 1 ) )
00266 *
00267 *+    TEST 2
00268 *              Form the inverse and compute the residual.
00269 *
00270                CALL DCOPY( NPP, AFAC, 1, AINV, 1 )
00271                SRNAMT = 'DPPTRI'
00272                CALL DPPTRI( UPLO, N, AINV, INFO )
00273 *
00274 *              Check error code from DPPTRI.
00275 *
00276                IF( INFO.NE.0 )
00277      $            CALL ALAERH( PATH, 'DPPTRI', INFO, 0, UPLO, N, N, -1,
00278      $                         -1, -1, IMAT, NFAIL, NERRS, NOUT )
00279 *
00280                CALL DPPT03( UPLO, N, A, AINV, WORK, LDA, RWORK, RCONDC,
00281      $                      RESULT( 2 ) )
00282 *
00283 *              Print information about the tests that did not pass
00284 *              the threshold.
00285 *
00286                DO 60 K = 1, 2
00287                   IF( RESULT( K ).GE.THRESH ) THEN
00288                      IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00289      $                  CALL ALAHD( NOUT, PATH )
00290                      WRITE( NOUT, FMT = 9999 )UPLO, N, IMAT, K,
00291      $                  RESULT( K )
00292                      NFAIL = NFAIL + 1
00293                   END IF
00294    60          CONTINUE
00295                NRUN = NRUN + 2
00296 *
00297                DO 80 IRHS = 1, NNS
00298                   NRHS = NSVAL( IRHS )
00299 *
00300 *+    TEST 3
00301 *              Solve and compute residual for  A * X = B.
00302 *
00303                   SRNAMT = 'DLARHS'
00304                   CALL DLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU,
00305      $                         NRHS, A, LDA, XACT, LDA, B, LDA, ISEED,
00306      $                         INFO )
00307                   CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
00308 *
00309                   SRNAMT = 'DPPTRS'
00310                   CALL DPPTRS( UPLO, N, NRHS, AFAC, X, LDA, INFO )
00311 *
00312 *              Check error code from DPPTRS.
00313 *
00314                   IF( INFO.NE.0 )
00315      $               CALL ALAERH( PATH, 'DPPTRS', INFO, 0, UPLO, N, N,
00316      $                            -1, -1, NRHS, IMAT, NFAIL, NERRS,
00317      $                            NOUT )
00318 *
00319                   CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
00320                   CALL DPPT02( UPLO, N, NRHS, A, X, LDA, WORK, LDA,
00321      $                         RWORK, RESULT( 3 ) )
00322 *
00323 *+    TEST 4
00324 *              Check solution from generated exact solution.
00325 *
00326                   CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
00327      $                         RESULT( 4 ) )
00328 *
00329 *+    TESTS 5, 6, and 7
00330 *              Use iterative refinement to improve the solution.
00331 *
00332                   SRNAMT = 'DPPRFS'
00333                   CALL DPPRFS( UPLO, N, NRHS, A, AFAC, B, LDA, X, LDA,
00334      $                         RWORK, RWORK( NRHS+1 ), WORK, IWORK,
00335      $                         INFO )
00336 *
00337 *              Check error code from DPPRFS.
00338 *
00339                   IF( INFO.NE.0 )
00340      $               CALL ALAERH( PATH, 'DPPRFS', INFO, 0, UPLO, N, N,
00341      $                            -1, -1, NRHS, IMAT, NFAIL, NERRS,
00342      $                            NOUT )
00343 *
00344                   CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
00345      $                         RESULT( 5 ) )
00346                   CALL DPPT05( UPLO, N, NRHS, A, B, LDA, X, LDA, XACT,
00347      $                         LDA, RWORK, RWORK( NRHS+1 ),
00348      $                         RESULT( 6 ) )
00349 *
00350 *                 Print information about the tests that did not pass
00351 *                 the threshold.
00352 *
00353                   DO 70 K = 3, 7
00354                      IF( RESULT( K ).GE.THRESH ) THEN
00355                         IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00356      $                     CALL ALAHD( NOUT, PATH )
00357                         WRITE( NOUT, FMT = 9998 )UPLO, N, NRHS, IMAT,
00358      $                     K, RESULT( K )
00359                         NFAIL = NFAIL + 1
00360                      END IF
00361    70             CONTINUE
00362                   NRUN = NRUN + 5
00363    80          CONTINUE
00364 *
00365 *+    TEST 8
00366 *              Get an estimate of RCOND = 1/CNDNUM.
00367 *
00368                ANORM = DLANSP( '1', UPLO, N, A, RWORK )
00369                SRNAMT = 'DPPCON'
00370                CALL DPPCON( UPLO, N, AFAC, ANORM, RCOND, WORK, IWORK,
00371      $                      INFO )
00372 *
00373 *              Check error code from DPPCON.
00374 *
00375                IF( INFO.NE.0 )
00376      $            CALL ALAERH( PATH, 'DPPCON', INFO, 0, UPLO, N, N, -1,
00377      $                         -1, -1, IMAT, NFAIL, NERRS, NOUT )
00378 *
00379                RESULT( 8 ) = DGET06( RCOND, RCONDC )
00380 *
00381 *              Print the test ratio if greater than or equal to THRESH.
00382 *
00383                IF( RESULT( 8 ).GE.THRESH ) THEN
00384                   IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00385      $               CALL ALAHD( NOUT, PATH )
00386                   WRITE( NOUT, FMT = 9999 )UPLO, N, IMAT, 8,
00387      $               RESULT( 8 )
00388                   NFAIL = NFAIL + 1
00389                END IF
00390                NRUN = NRUN + 1
00391    90       CONTINUE
00392   100    CONTINUE
00393   110 CONTINUE
00394 *
00395 *     Print a summary of the results.
00396 *
00397       CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
00398 *
00399  9999 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', type ', I2, ', test ',
00400      $      I2, ', ratio =', G12.5 )
00401  9998 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NRHS=', I3, ', type ',
00402      $      I2, ', test(', I2, ') =', G12.5 )
00403       RETURN
00404 *
00405 *     End of DCHKPP
00406 *
00407       END