LAPACK: TESTING/LIN/dchkpo.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DCHKPO( 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 *  DCHKPO tests DPOTRF, -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 *  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   ZERO
00092       PARAMETER          ( ZERO = 0.0D+0 )
00093       INTEGER            NTYPES
00094       PARAMETER          ( NTYPES = 9 )
00095       INTEGER            NTESTS
00096       PARAMETER          ( NTESTS = 8 )
00097 *     ..
00098 *     .. Local Scalars ..
00099       LOGICAL            ZEROT
00100       CHARACTER          DIST, TYPE, UPLO, XTYPE
00101       CHARACTER*3        PATH
00102       INTEGER            I, IMAT, IN, INB, INFO, IOFF, IRHS, IUPLO,
00103      $                   IZERO, K, KL, KU, LDA, MODE, N, NB, NERRS,
00104      $                   NFAIL, NIMAT, NRHS, NRUN
00105       DOUBLE PRECISION   ANORM, CNDNUM, RCOND, RCONDC
00106 *     ..
00107 *     .. Local Arrays ..
00108       CHARACTER          UPLOS( 2 )
00109       INTEGER            ISEED( 4 ), ISEEDY( 4 )
00110       DOUBLE PRECISION   RESULT( NTESTS )
00111 *     ..
00112 *     .. External Functions ..
00113       DOUBLE PRECISION   DGET06, DLANSY
00114       EXTERNAL           DGET06, DLANSY
00115 *     ..
00116 *     .. External Subroutines ..
00117       EXTERNAL           ALAERH, ALAHD, ALASUM, DERRPO, DGET04, DLACPY,
00118      $                   DLARHS, DLATB4, DLATMS, DPOCON, DPORFS, DPOT01,
00119      $                   DPOT02, DPOT03, DPOT05, DPOTRF, DPOTRI, DPOTRS,
00120      $                   XLAENV
00121 *     ..
00122 *     .. Scalars in Common ..
00123       LOGICAL            LERR, OK
00124       CHARACTER*32       SRNAMT
00125       INTEGER            INFOT, NUNIT
00126 *     ..
00127 *     .. Common blocks ..
00128       COMMON             / INFOC / INFOT, NUNIT, OK, LERR
00129       COMMON             / SRNAMC / SRNAMT
00130 *     ..
00131 *     .. Intrinsic Functions ..
00132       INTRINSIC          MAX
00133 *     ..
00134 *     .. Data statements ..
00135       DATA               ISEEDY / 1988, 1989, 1990, 1991 /
00136       DATA               UPLOS / 'U', 'L' /
00137 *     ..
00138 *     .. Executable Statements ..
00139 *
00140 *     Initialize constants and the random number seed.
00141 *
00142       PATH( 1: 1 ) = 'Double precision'
00143       PATH( 2: 3 ) = 'PO'
00144       NRUN = 0
00145       NFAIL = 0
00146       NERRS = 0
00147       DO 10 I = 1, 4
00148          ISEED( I ) = ISEEDY( I )
00149    10 CONTINUE
00150 *
00151 *     Test the error exits
00152 *
00153       IF( TSTERR )
00154      $   CALL DERRPO( PATH, NOUT )
00155       INFOT = 0
00156       CALL XLAENV( 2, 2 )
00157 *
00158 *     Do for each value of N in NVAL
00159 *
00160       DO 120 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 110 IMAT = 1, NIMAT
00170 *
00171 *           Do the tests only if DOTYPE( IMAT ) is true.
00172 *
00173             IF( .NOT.DOTYPE( IMAT ) )
00174      $         GO TO 110
00175 *
00176 *           Skip types 3, 4, or 5 if the matrix size is too small.
00177 *
00178             ZEROT = IMAT.GE.3 .AND. IMAT.LE.5
00179             IF( ZEROT .AND. N.LT.IMAT-2 )
00180      $         GO TO 110
00181 *
00182 *           Do first for UPLO = 'U', then for UPLO = 'L'
00183 *
00184             DO 100 IUPLO = 1, 2
00185                UPLO = UPLOS( IUPLO )
00186 *
00187 *              Set up parameters with DLATB4 and generate a test matrix
00188 *              with DLATMS.
00189 *
00190                CALL DLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
00191      $                      CNDNUM, DIST )
00192 *
00193                SRNAMT = 'DLATMS'
00194                CALL DLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
00195      $                      CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK,
00196      $                      INFO )
00197 *
00198 *              Check error code from DLATMS.
00199 *
00200                IF( INFO.NE.0 ) THEN
00201                   CALL ALAERH( PATH, 'DLATMS', INFO, 0, UPLO, N, N, -1,
00202      $                         -1, -1, IMAT, NFAIL, NERRS, NOUT )
00203                   GO TO 100
00204                END IF
00205 *
00206 *              For types 3-5, zero one row and column of the matrix to
00207 *              test that INFO is returned correctly.
00208 *
00209                IF( ZEROT ) THEN
00210                   IF( IMAT.EQ.3 ) THEN
00211                      IZERO = 1
00212                   ELSE IF( IMAT.EQ.4 ) THEN
00213                      IZERO = N
00214                   ELSE
00215                      IZERO = N / 2 + 1
00216                   END IF
00217                   IOFF = ( IZERO-1 )*LDA
00218 *
00219 *                 Set row and column IZERO of A to 0.
00220 *
00221                   IF( IUPLO.EQ.1 ) THEN
00222                      DO 20 I = 1, IZERO - 1
00223                         A( IOFF+I ) = ZERO
00224    20                CONTINUE
00225                      IOFF = IOFF + IZERO
00226                      DO 30 I = IZERO, N
00227                         A( IOFF ) = ZERO
00228                         IOFF = IOFF + LDA
00229    30                CONTINUE
00230                   ELSE
00231                      IOFF = IZERO
00232                      DO 40 I = 1, IZERO - 1
00233                         A( IOFF ) = ZERO
00234                         IOFF = IOFF + LDA
00235    40                CONTINUE
00236                      IOFF = IOFF - IZERO
00237                      DO 50 I = IZERO, N
00238                         A( IOFF+I ) = ZERO
00239    50                CONTINUE
00240                   END IF
00241                ELSE
00242                   IZERO = 0
00243                END IF
00244 *
00245 *              Do for each value of NB in NBVAL
00246 *
00247                DO 90 INB = 1, NNB
00248                   NB = NBVAL( INB )
00249                   CALL XLAENV( 1, NB )
00250 *
00251 *                 Compute the L*L' or U'*U factorization of the matrix.
00252 *
00253                   CALL DLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
00254                   SRNAMT = 'DPOTRF'
00255                   CALL DPOTRF( UPLO, N, AFAC, LDA, INFO )
00256 *
00257 *                 Check error code from DPOTRF.
00258 *
00259                   IF( INFO.NE.IZERO ) THEN
00260                      CALL ALAERH( PATH, 'DPOTRF', INFO, IZERO, UPLO, N,
00261      $                            N, -1, -1, NB, IMAT, NFAIL, NERRS,
00262      $                            NOUT )
00263                      GO TO 90
00264                   END IF
00265 *
00266 *                 Skip the tests if INFO is not 0.
00267 *
00268                   IF( INFO.NE.0 )
00269      $               GO TO 90
00270 *
00271 *+    TEST 1
00272 *                 Reconstruct matrix from factors and compute residual.
00273 *
00274                   CALL DLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
00275                   CALL DPOT01( UPLO, N, A, LDA, AINV, LDA, RWORK,
00276      $                         RESULT( 1 ) )
00277 *
00278 *+    TEST 2
00279 *                 Form the inverse and compute the residual.
00280 *
00281                   CALL DLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
00282                   SRNAMT = 'DPOTRI'
00283                   CALL DPOTRI( UPLO, N, AINV, LDA, INFO )
00284 *
00285 *                 Check error code from DPOTRI.
00286 *
00287                   IF( INFO.NE.0 )
00288      $               CALL ALAERH( PATH, 'DPOTRI', INFO, 0, UPLO, N, N,
00289      $                            -1, -1, -1, IMAT, NFAIL, NERRS, NOUT )
00290 *
00291                   CALL DPOT03( UPLO, N, A, LDA, AINV, LDA, WORK, LDA,
00292      $                         RWORK, RCONDC, RESULT( 2 ) )
00293 *
00294 *                 Print information about the tests that did not pass
00295 *                 the threshold.
00296 *
00297                   DO 60 K = 1, 2
00298                      IF( RESULT( K ).GE.THRESH ) THEN
00299                         IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00300      $                     CALL ALAHD( NOUT, PATH )
00301                         WRITE( NOUT, FMT = 9999 )UPLO, N, NB, IMAT, K,
00302      $                     RESULT( K )
00303                         NFAIL = NFAIL + 1
00304                      END IF
00305    60             CONTINUE
00306                   NRUN = NRUN + 2
00307 *
00308 *                 Skip the rest of the tests unless this is the first
00309 *                 blocksize.
00310 *
00311                   IF( INB.NE.1 )
00312      $               GO TO 90
00313 *
00314                   DO 80 IRHS = 1, NNS
00315                      NRHS = NSVAL( IRHS )
00316 *
00317 *+    TEST 3
00318 *                 Solve and compute residual for A * X = B .
00319 *
00320                      SRNAMT = 'DLARHS'
00321                      CALL DLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU,
00322      $                            NRHS, A, LDA, XACT, LDA, B, LDA,
00323      $                            ISEED, INFO )
00324                      CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
00325 *
00326                      SRNAMT = 'DPOTRS'
00327                      CALL DPOTRS( UPLO, N, NRHS, AFAC, LDA, X, LDA,
00328      $                            INFO )
00329 *
00330 *                 Check error code from DPOTRS.
00331 *
00332                      IF( INFO.NE.0 )
00333      $                  CALL ALAERH( PATH, 'DPOTRS', INFO, 0, UPLO, N,
00334      $                               N, -1, -1, NRHS, IMAT, NFAIL,
00335      $                               NERRS, NOUT )
00336 *
00337                      CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
00338                      CALL DPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK,
00339      $                            LDA, RWORK, RESULT( 3 ) )
00340 *
00341 *+    TEST 4
00342 *                 Check solution from generated exact solution.
00343 *
00344                      CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
00345      $                            RESULT( 4 ) )
00346 *
00347 *+    TESTS 5, 6, and 7
00348 *                 Use iterative refinement to improve the solution.
00349 *
00350                      SRNAMT = 'DPORFS'
00351                      CALL DPORFS( UPLO, N, NRHS, A, LDA, AFAC, LDA, B,
00352      $                            LDA, X, LDA, RWORK, RWORK( NRHS+1 ),
00353      $                            WORK, IWORK, INFO )
00354 *
00355 *                 Check error code from DPORFS.
00356 *
00357                      IF( INFO.NE.0 )
00358      $                  CALL ALAERH( PATH, 'DPORFS', INFO, 0, UPLO, N,
00359      $                               N, -1, -1, NRHS, IMAT, NFAIL,
00360      $                               NERRS, NOUT )
00361 *
00362                      CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
00363      $                            RESULT( 5 ) )
00364                      CALL DPOT05( UPLO, N, NRHS, A, LDA, B, LDA, X, LDA,
00365      $                            XACT, LDA, RWORK, RWORK( NRHS+1 ),
00366      $                            RESULT( 6 ) )
00367 *
00368 *                    Print information about the tests that did not pass
00369 *                    the threshold.
00370 *
00371                      DO 70 K = 3, 7
00372                         IF( RESULT( K ).GE.THRESH ) THEN
00373                            IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00374      $                        CALL ALAHD( NOUT, PATH )
00375                            WRITE( NOUT, FMT = 9998 )UPLO, N, NRHS,
00376      $                        IMAT, K, RESULT( K )
00377                            NFAIL = NFAIL + 1
00378                         END IF
00379    70                CONTINUE
00380                      NRUN = NRUN + 5
00381    80             CONTINUE
00382 *
00383 *+    TEST 8
00384 *                 Get an estimate of RCOND = 1/CNDNUM.
00385 *
00386                   ANORM = DLANSY( '1', UPLO, N, A, LDA, RWORK )
00387                   SRNAMT = 'DPOCON'
00388                   CALL DPOCON( UPLO, N, AFAC, LDA, ANORM, RCOND, WORK,
00389      $                         IWORK, INFO )
00390 *
00391 *                 Check error code from DPOCON.
00392 *
00393                   IF( INFO.NE.0 )
00394      $               CALL ALAERH( PATH, 'DPOCON', INFO, 0, UPLO, N, N,
00395      $                            -1, -1, -1, IMAT, NFAIL, NERRS, NOUT )
00396 *
00397                   RESULT( 8 ) = DGET06( RCOND, RCONDC )
00398 *
00399 *                 Print the test ratio if it is .GE. THRESH.
00400 *
00401                   IF( RESULT( 8 ).GE.THRESH ) THEN
00402                      IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00403      $                  CALL ALAHD( NOUT, PATH )
00404                      WRITE( NOUT, FMT = 9997 )UPLO, N, IMAT, 8,
00405      $                  RESULT( 8 )
00406                      NFAIL = NFAIL + 1
00407                   END IF
00408                   NRUN = NRUN + 1
00409    90          CONTINUE
00410   100       CONTINUE
00411   110    CONTINUE
00412   120 CONTINUE
00413 *
00414 *     Print a summary of the results.
00415 *
00416       CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
00417 *
00418  9999 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NB =', I4, ', type ',
00419      $      I2, ', test ', I2, ', ratio =', G12.5 )
00420  9998 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NRHS=', I3, ', type ',
00421      $      I2, ', test(', I2, ') =', G12.5 )
00422  9997 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ',', 10X, ' type ', I2,
00423      $      ', test(', I2, ') =', G12.5 )
00424       RETURN
00425 *
00426 *     End of DCHKPO
00427 *
00428       END