LAPACK: TESTING/LIN/dchkge.f Source File

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