LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dchksp.f
Go to the documentation of this file.
00001 *> \brief \b DCHKSP
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *  Definition:
00009 *  ===========
00010 *
00011 *       SUBROUTINE DCHKSP( DOTYPE, NN, NVAL, NNS, NSVAL, THRESH, TSTERR,
00012 *                          NMAX, A, AFAC, AINV, B, X, XACT, WORK, RWORK,
00013 *                          IWORK, NOUT )
00014 * 
00015 *       .. Scalar Arguments ..
00016 *       LOGICAL            TSTERR
00017 *       INTEGER            NMAX, NN, NNS, NOUT
00018 *       DOUBLE PRECISION   THRESH
00019 *       ..
00020 *       .. Array Arguments ..
00021 *       LOGICAL            DOTYPE( * )
00022 *       INTEGER            IWORK( * ), NSVAL( * ), NVAL( * )
00023 *       DOUBLE PRECISION   A( * ), AFAC( * ), AINV( * ), B( * ),
00024 *      $                   RWORK( * ), WORK( * ), X( * ), XACT( * )
00025 *       ..
00026 *  
00027 *
00028 *> \par Purpose:
00029 *  =============
00030 *>
00031 *> \verbatim
00032 *>
00033 *> DCHKSP tests DSPTRF, -TRI, -TRS, -RFS, and -CON
00034 *> \endverbatim
00035 *
00036 *  Arguments:
00037 *  ==========
00038 *
00039 *> \param[in] DOTYPE
00040 *> \verbatim
00041 *>          DOTYPE is LOGICAL array, dimension (NTYPES)
00042 *>          The matrix types to be used for testing.  Matrices of type j
00043 *>          (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
00044 *>          .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
00045 *> \endverbatim
00046 *>
00047 *> \param[in] NN
00048 *> \verbatim
00049 *>          NN is INTEGER
00050 *>          The number of values of N contained in the vector NVAL.
00051 *> \endverbatim
00052 *>
00053 *> \param[in] NVAL
00054 *> \verbatim
00055 *>          NVAL is INTEGER array, dimension (NN)
00056 *>          The values of the matrix dimension N.
00057 *> \endverbatim
00058 *>
00059 *> \param[in] NNS
00060 *> \verbatim
00061 *>          NNS is INTEGER
00062 *>          The number of values of NRHS contained in the vector NSVAL.
00063 *> \endverbatim
00064 *>
00065 *> \param[in] NSVAL
00066 *> \verbatim
00067 *>          NSVAL is INTEGER array, dimension (NNS)
00068 *>          The values of the number of right hand sides NRHS.
00069 *> \endverbatim
00070 *>
00071 *> \param[in] THRESH
00072 *> \verbatim
00073 *>          THRESH is DOUBLE PRECISION
00074 *>          The threshold value for the test ratios.  A result is
00075 *>          included in the output file if RESULT >= THRESH.  To have
00076 *>          every test ratio printed, use THRESH = 0.
00077 *> \endverbatim
00078 *>
00079 *> \param[in] TSTERR
00080 *> \verbatim
00081 *>          TSTERR is LOGICAL
00082 *>          Flag that indicates whether error exits are to be tested.
00083 *> \endverbatim
00084 *>
00085 *> \param[in] NMAX
00086 *> \verbatim
00087 *>          NMAX is INTEGER
00088 *>          The maximum value permitted for N, used in dimensioning the
00089 *>          work arrays.
00090 *> \endverbatim
00091 *>
00092 *> \param[out] A
00093 *> \verbatim
00094 *>          A is DOUBLE PRECISION array, dimension
00095 *>                      (NMAX*(NMAX+1)/2)
00096 *> \endverbatim
00097 *>
00098 *> \param[out] AFAC
00099 *> \verbatim
00100 *>          AFAC is DOUBLE PRECISION array, dimension
00101 *>                      (NMAX*(NMAX+1)/2)
00102 *> \endverbatim
00103 *>
00104 *> \param[out] AINV
00105 *> \verbatim
00106 *>          AINV is DOUBLE PRECISION array, dimension
00107 *>                      (NMAX*(NMAX+1)/2)
00108 *> \endverbatim
00109 *>
00110 *> \param[out] B
00111 *> \verbatim
00112 *>          B is DOUBLE PRECISION array, dimension (NMAX*NSMAX)
00113 *>          where NSMAX is the largest entry in NSVAL.
00114 *> \endverbatim
00115 *>
00116 *> \param[out] X
00117 *> \verbatim
00118 *>          X is DOUBLE PRECISION array, dimension (NMAX*NSMAX)
00119 *> \endverbatim
00120 *>
00121 *> \param[out] XACT
00122 *> \verbatim
00123 *>          XACT is DOUBLE PRECISION array, dimension (NMAX*NSMAX)
00124 *> \endverbatim
00125 *>
00126 *> \param[out] WORK
00127 *> \verbatim
00128 *>          WORK is DOUBLE PRECISION array, dimension
00129 *>                      (NMAX*max(2,NSMAX))
00130 *> \endverbatim
00131 *>
00132 *> \param[out] RWORK
00133 *> \verbatim
00134 *>          RWORK is DOUBLE PRECISION array,
00135 *>                                 dimension (NMAX+2*NSMAX)
00136 *> \endverbatim
00137 *>
00138 *> \param[out] IWORK
00139 *> \verbatim
00140 *>          IWORK is INTEGER array, dimension (2*NMAX)
00141 *> \endverbatim
00142 *>
00143 *> \param[in] NOUT
00144 *> \verbatim
00145 *>          NOUT is INTEGER
00146 *>          The unit number for output.
00147 *> \endverbatim
00148 *
00149 *  Authors:
00150 *  ========
00151 *
00152 *> \author Univ. of Tennessee 
00153 *> \author Univ. of California Berkeley 
00154 *> \author Univ. of Colorado Denver 
00155 *> \author NAG Ltd. 
00156 *
00157 *> \date November 2011
00158 *
00159 *> \ingroup double_lin
00160 *
00161 *  =====================================================================
00162       SUBROUTINE DCHKSP( DOTYPE, NN, NVAL, NNS, NSVAL, THRESH, TSTERR,
00163      $                   NMAX, A, AFAC, AINV, B, X, XACT, WORK, RWORK,
00164      $                   IWORK, NOUT )
00165 *
00166 *  -- LAPACK test routine (version 3.4.0) --
00167 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00168 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00169 *     November 2011
00170 *
00171 *     .. Scalar Arguments ..
00172       LOGICAL            TSTERR
00173       INTEGER            NMAX, NN, NNS, NOUT
00174       DOUBLE PRECISION   THRESH
00175 *     ..
00176 *     .. Array Arguments ..
00177       LOGICAL            DOTYPE( * )
00178       INTEGER            IWORK( * ), NSVAL( * ), NVAL( * )
00179       DOUBLE PRECISION   A( * ), AFAC( * ), AINV( * ), B( * ),
00180      $                   RWORK( * ), WORK( * ), X( * ), XACT( * )
00181 *     ..
00182 *
00183 *  =====================================================================
00184 *
00185 *     .. Parameters ..
00186       DOUBLE PRECISION   ZERO
00187       PARAMETER          ( ZERO = 0.0D+0 )
00188       INTEGER            NTYPES
00189       PARAMETER          ( NTYPES = 10 )
00190       INTEGER            NTESTS
00191       PARAMETER          ( NTESTS = 8 )
00192 *     ..
00193 *     .. Local Scalars ..
00194       LOGICAL            TRFCON, ZEROT
00195       CHARACTER          DIST, PACKIT, TYPE, UPLO, XTYPE
00196       CHARACTER*3        PATH
00197       INTEGER            I, I1, I2, IMAT, IN, INFO, IOFF, IRHS, IUPLO,
00198      $                   IZERO, J, K, KL, KU, LDA, MODE, N, NERRS,
00199      $                   NFAIL, NIMAT, NPP, NRHS, NRUN, NT
00200       DOUBLE PRECISION   ANORM, CNDNUM, RCOND, RCONDC
00201 *     ..
00202 *     .. Local Arrays ..
00203       CHARACTER          UPLOS( 2 )
00204       INTEGER            ISEED( 4 ), ISEEDY( 4 )
00205       DOUBLE PRECISION   RESULT( NTESTS )
00206 *     ..
00207 *     .. External Functions ..
00208       LOGICAL            LSAME
00209       DOUBLE PRECISION   DGET06, DLANSP
00210       EXTERNAL           LSAME, DGET06, DLANSP
00211 *     ..
00212 *     .. External Subroutines ..
00213       EXTERNAL           ALAERH, ALAHD, ALASUM, DCOPY, DERRSY, DGET04,
00214      $                   DLACPY, DLARHS, DLATB4, DLATMS, DPPT02, DPPT03,
00215      $                   DPPT05, DSPCON, DSPRFS, DSPT01, DSPTRF, DSPTRI,
00216      $                   DSPTRS
00217 *     ..
00218 *     .. Intrinsic Functions ..
00219       INTRINSIC          MAX, MIN
00220 *     ..
00221 *     .. Scalars in Common ..
00222       LOGICAL            LERR, OK
00223       CHARACTER*32       SRNAMT
00224       INTEGER            INFOT, NUNIT
00225 *     ..
00226 *     .. Common blocks ..
00227       COMMON             / INFOC / INFOT, NUNIT, OK, LERR
00228       COMMON             / SRNAMC / SRNAMT
00229 *     ..
00230 *     .. Data statements ..
00231       DATA               ISEEDY / 1988, 1989, 1990, 1991 /
00232       DATA               UPLOS / 'U', 'L' /
00233 *     ..
00234 *     .. Executable Statements ..
00235 *
00236 *     Initialize constants and the random number seed.
00237 *
00238       PATH( 1: 1 ) = 'Double precision'
00239       PATH( 2: 3 ) = 'SP'
00240       NRUN = 0
00241       NFAIL = 0
00242       NERRS = 0
00243       DO 10 I = 1, 4
00244          ISEED( I ) = ISEEDY( I )
00245    10 CONTINUE
00246 *
00247 *     Test the error exits
00248 *
00249       IF( TSTERR )
00250      $   CALL DERRSY( PATH, NOUT )
00251       INFOT = 0
00252 *
00253 *     Do for each value of N in NVAL
00254 *
00255       DO 170 IN = 1, NN
00256          N = NVAL( IN )
00257          LDA = MAX( N, 1 )
00258          XTYPE = 'N'
00259          NIMAT = NTYPES
00260          IF( N.LE.0 )
00261      $      NIMAT = 1
00262 *
00263          IZERO = 0
00264          DO 160 IMAT = 1, NIMAT
00265 *
00266 *           Do the tests only if DOTYPE( IMAT ) is true.
00267 *
00268             IF( .NOT.DOTYPE( IMAT ) )
00269      $         GO TO 160
00270 *
00271 *           Skip types 3, 4, 5, or 6 if the matrix size is too small.
00272 *
00273             ZEROT = IMAT.GE.3 .AND. IMAT.LE.6
00274             IF( ZEROT .AND. N.LT.IMAT-2 )
00275      $         GO TO 160
00276 *
00277 *           Do first for UPLO = 'U', then for UPLO = 'L'
00278 *
00279             DO 150 IUPLO = 1, 2
00280                UPLO = UPLOS( IUPLO )
00281                IF( LSAME( UPLO, 'U' ) ) THEN
00282                   PACKIT = 'C'
00283                ELSE
00284                   PACKIT = 'R'
00285                END IF
00286 *
00287 *              Set up parameters with DLATB4 and generate a test matrix
00288 *              with DLATMS.
00289 *
00290                CALL DLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
00291      $                      CNDNUM, DIST )
00292 *
00293                SRNAMT = 'DLATMS'
00294                CALL DLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
00295      $                      CNDNUM, ANORM, KL, KU, PACKIT, A, LDA, WORK,
00296      $                      INFO )
00297 *
00298 *              Check error code from DLATMS.
00299 *
00300                IF( INFO.NE.0 ) THEN
00301                   CALL ALAERH( PATH, 'DLATMS', INFO, 0, UPLO, N, N, -1,
00302      $                         -1, -1, IMAT, NFAIL, NERRS, NOUT )
00303                   GO TO 150
00304                END IF
00305 *
00306 *              For types 3-6, zero one or more rows and columns of
00307 *              the matrix to test that INFO is returned correctly.
00308 *
00309                IF( ZEROT ) THEN
00310                   IF( IMAT.EQ.3 ) THEN
00311                      IZERO = 1
00312                   ELSE IF( IMAT.EQ.4 ) THEN
00313                      IZERO = N
00314                   ELSE
00315                      IZERO = N / 2 + 1
00316                   END IF
00317 *
00318                   IF( IMAT.LT.6 ) THEN
00319 *
00320 *                    Set row and column IZERO to zero.
00321 *
00322                      IF( IUPLO.EQ.1 ) THEN
00323                         IOFF = ( IZERO-1 )*IZERO / 2
00324                         DO 20 I = 1, IZERO - 1
00325                            A( IOFF+I ) = ZERO
00326    20                   CONTINUE
00327                         IOFF = IOFF + IZERO
00328                         DO 30 I = IZERO, N
00329                            A( IOFF ) = ZERO
00330                            IOFF = IOFF + I
00331    30                   CONTINUE
00332                      ELSE
00333                         IOFF = IZERO
00334                         DO 40 I = 1, IZERO - 1
00335                            A( IOFF ) = ZERO
00336                            IOFF = IOFF + N - I
00337    40                   CONTINUE
00338                         IOFF = IOFF - IZERO
00339                         DO 50 I = IZERO, N
00340                            A( IOFF+I ) = ZERO
00341    50                   CONTINUE
00342                      END IF
00343                   ELSE
00344                      IOFF = 0
00345                      IF( IUPLO.EQ.1 ) THEN
00346 *
00347 *                       Set the first IZERO rows and columns to zero.
00348 *
00349                         DO 70 J = 1, N
00350                            I2 = MIN( J, IZERO )
00351                            DO 60 I = 1, I2
00352                               A( IOFF+I ) = ZERO
00353    60                      CONTINUE
00354                            IOFF = IOFF + J
00355    70                   CONTINUE
00356                      ELSE
00357 *
00358 *                       Set the last IZERO rows and columns to zero.
00359 *
00360                         DO 90 J = 1, N
00361                            I1 = MAX( J, IZERO )
00362                            DO 80 I = I1, N
00363                               A( IOFF+I ) = ZERO
00364    80                      CONTINUE
00365                            IOFF = IOFF + N - J
00366    90                   CONTINUE
00367                      END IF
00368                   END IF
00369                ELSE
00370                   IZERO = 0
00371                END IF
00372 *
00373 *              Compute the L*D*L' or U*D*U' factorization of the matrix.
00374 *
00375                NPP = N*( N+1 ) / 2
00376                CALL DCOPY( NPP, A, 1, AFAC, 1 )
00377                SRNAMT = 'DSPTRF'
00378                CALL DSPTRF( UPLO, N, AFAC, IWORK, INFO )
00379 *
00380 *              Adjust the expected value of INFO to account for
00381 *              pivoting.
00382 *
00383                K = IZERO
00384                IF( K.GT.0 ) THEN
00385   100             CONTINUE
00386                   IF( IWORK( K ).LT.0 ) THEN
00387                      IF( IWORK( K ).NE.-K ) THEN
00388                         K = -IWORK( K )
00389                         GO TO 100
00390                      END IF
00391                   ELSE IF( IWORK( K ).NE.K ) THEN
00392                      K = IWORK( K )
00393                      GO TO 100
00394                   END IF
00395                END IF
00396 *
00397 *              Check error code from DSPTRF.
00398 *
00399                IF( INFO.NE.K )
00400      $            CALL ALAERH( PATH, 'DSPTRF', INFO, K, UPLO, N, N, -1,
00401      $                         -1, -1, IMAT, NFAIL, NERRS, NOUT )
00402                IF( INFO.NE.0 ) THEN
00403                   TRFCON = .TRUE.
00404                ELSE
00405                   TRFCON = .FALSE.
00406                END IF
00407 *
00408 *+    TEST 1
00409 *              Reconstruct matrix from factors and compute residual.
00410 *
00411                CALL DSPT01( UPLO, N, A, AFAC, IWORK, AINV, LDA, RWORK,
00412      $                      RESULT( 1 ) )
00413                NT = 1
00414 *
00415 *+    TEST 2
00416 *              Form the inverse and compute the residual.
00417 *
00418                IF( .NOT.TRFCON ) THEN
00419                   CALL DCOPY( NPP, AFAC, 1, AINV, 1 )
00420                   SRNAMT = 'DSPTRI'
00421                   CALL DSPTRI( UPLO, N, AINV, IWORK, WORK, INFO )
00422 *
00423 *              Check error code from DSPTRI.
00424 *
00425                   IF( INFO.NE.0 )
00426      $               CALL ALAERH( PATH, 'DSPTRI', INFO, 0, UPLO, N, N,
00427      $                            -1, -1, -1, IMAT, NFAIL, NERRS, NOUT )
00428 *
00429                   CALL DPPT03( UPLO, N, A, AINV, WORK, LDA, RWORK,
00430      $                         RCONDC, RESULT( 2 ) )
00431                   NT = 2
00432                END IF
00433 *
00434 *              Print information about the tests that did not pass
00435 *              the threshold.
00436 *
00437                DO 110 K = 1, NT
00438                   IF( RESULT( K ).GE.THRESH ) THEN
00439                      IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00440      $                  CALL ALAHD( NOUT, PATH )
00441                      WRITE( NOUT, FMT = 9999 )UPLO, N, IMAT, K,
00442      $                  RESULT( K )
00443                      NFAIL = NFAIL + 1
00444                   END IF
00445   110          CONTINUE
00446                NRUN = NRUN + NT
00447 *
00448 *              Do only the condition estimate if INFO is not 0.
00449 *
00450                IF( TRFCON ) THEN
00451                   RCONDC = ZERO
00452                   GO TO 140
00453                END IF
00454 *
00455                DO 130 IRHS = 1, NNS
00456                   NRHS = NSVAL( IRHS )
00457 *
00458 *+    TEST 3
00459 *              Solve and compute residual for  A * X = B.
00460 *
00461                   SRNAMT = 'DLARHS'
00462                   CALL DLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU,
00463      $                         NRHS, A, LDA, XACT, LDA, B, LDA, ISEED,
00464      $                         INFO )
00465                   CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
00466 *
00467                   SRNAMT = 'DSPTRS'
00468                   CALL DSPTRS( UPLO, N, NRHS, AFAC, IWORK, X, LDA,
00469      $                         INFO )
00470 *
00471 *              Check error code from DSPTRS.
00472 *
00473                   IF( INFO.NE.0 )
00474      $               CALL ALAERH( PATH, 'DSPTRS', INFO, 0, UPLO, N, N,
00475      $                            -1, -1, NRHS, IMAT, NFAIL, NERRS,
00476      $                            NOUT )
00477 *
00478                   CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
00479                   CALL DPPT02( UPLO, N, NRHS, A, X, LDA, WORK, LDA,
00480      $                         RWORK, RESULT( 3 ) )
00481 *
00482 *+    TEST 4
00483 *              Check solution from generated exact solution.
00484 *
00485                   CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
00486      $                         RESULT( 4 ) )
00487 *
00488 *+    TESTS 5, 6, and 7
00489 *              Use iterative refinement to improve the solution.
00490 *
00491                   SRNAMT = 'DSPRFS'
00492                   CALL DSPRFS( UPLO, N, NRHS, A, AFAC, IWORK, B, LDA, X,
00493      $                         LDA, RWORK, RWORK( NRHS+1 ), WORK,
00494      $                         IWORK( N+1 ), INFO )
00495 *
00496 *              Check error code from DSPRFS.
00497 *
00498                   IF( INFO.NE.0 )
00499      $               CALL ALAERH( PATH, 'DSPRFS', INFO, 0, UPLO, N, N,
00500      $                            -1, -1, NRHS, IMAT, NFAIL, NERRS,
00501      $                            NOUT )
00502 *
00503                   CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
00504      $                         RESULT( 5 ) )
00505                   CALL DPPT05( UPLO, N, NRHS, A, B, LDA, X, LDA, XACT,
00506      $                         LDA, RWORK, RWORK( NRHS+1 ),
00507      $                         RESULT( 6 ) )
00508 *
00509 *                 Print information about the tests that did not pass
00510 *                 the threshold.
00511 *
00512                   DO 120 K = 3, 7
00513                      IF( RESULT( K ).GE.THRESH ) THEN
00514                         IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00515      $                     CALL ALAHD( NOUT, PATH )
00516                         WRITE( NOUT, FMT = 9998 )UPLO, N, NRHS, IMAT,
00517      $                     K, RESULT( K )
00518                         NFAIL = NFAIL + 1
00519                      END IF
00520   120             CONTINUE
00521                   NRUN = NRUN + 5
00522   130          CONTINUE
00523 *
00524 *+    TEST 8
00525 *              Get an estimate of RCOND = 1/CNDNUM.
00526 *
00527   140          CONTINUE
00528                ANORM = DLANSP( '1', UPLO, N, A, RWORK )
00529                SRNAMT = 'DSPCON'
00530                CALL DSPCON( UPLO, N, AFAC, IWORK, ANORM, RCOND, WORK,
00531      $                      IWORK( N+1 ), INFO )
00532 *
00533 *              Check error code from DSPCON.
00534 *
00535                IF( INFO.NE.0 )
00536      $            CALL ALAERH( PATH, 'DSPCON', INFO, 0, UPLO, N, N, -1,
00537      $                         -1, -1, IMAT, NFAIL, NERRS, NOUT )
00538 *
00539                RESULT( 8 ) = DGET06( RCOND, RCONDC )
00540 *
00541 *              Print the test ratio if it is .GE. THRESH.
00542 *
00543                IF( RESULT( 8 ).GE.THRESH ) THEN
00544                   IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00545      $               CALL ALAHD( NOUT, PATH )
00546                   WRITE( NOUT, FMT = 9999 )UPLO, N, IMAT, 8,
00547      $               RESULT( 8 )
00548                   NFAIL = NFAIL + 1
00549                END IF
00550                NRUN = NRUN + 1
00551   150       CONTINUE
00552   160    CONTINUE
00553   170 CONTINUE
00554 *
00555 *     Print a summary of the results.
00556 *
00557       CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
00558 *
00559  9999 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', type ', I2, ', test ',
00560      $      I2, ', ratio =', G12.5 )
00561  9998 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NRHS=', I3, ', type ',
00562      $      I2, ', test(', I2, ') =', G12.5 )
00563       RETURN
00564 *
00565 *     End of DCHKSP
00566 *
00567       END
 All Files Functions