LAPACK: TESTING/LIN/dchksp.f Source File

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