LAPACK: TESTING/LIN/dchkgt.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DCHKGT( DOTYPE, NN, NVAL, NNS, NSVAL, THRESH, TSTERR,
00002      $                   A, AF, B, X, XACT, WORK, RWORK, IWORK, NOUT )
00003 *
00004 *  -- LAPACK test routine (version 3.1) --
00005 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       LOGICAL            TSTERR
00010       INTEGER            NN, NNS, NOUT
00011       DOUBLE PRECISION   THRESH
00012 *     ..
00013 *     .. Array Arguments ..
00014       LOGICAL            DOTYPE( * )
00015       INTEGER            IWORK( * ), NSVAL( * ), NVAL( * )
00016       DOUBLE PRECISION   A( * ), AF( * ), B( * ), RWORK( * ), WORK( * ),
00017      $                   X( * ), XACT( * )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  DCHKGT tests DGTTRF, -TRS, -RFS, and -CON
00024 *
00025 *  Arguments
00026 *  =========
00027 *
00028 *  DOTYPE  (input) LOGICAL array, dimension (NTYPES)
00029 *          The matrix types to be used for testing.  Matrices of type j
00030 *          (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
00031 *          .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
00032 *
00033 *  NN      (input) INTEGER
00034 *          The number of values of N contained in the vector NVAL.
00035 *
00036 *  NVAL    (input) INTEGER array, dimension (NN)
00037 *          The values of the matrix dimension N.
00038 *
00039 *  NNS     (input) INTEGER
00040 *          The number of values of NRHS contained in the vector NSVAL.
00041 *
00042 *  NSVAL   (input) INTEGER array, dimension (NNS)
00043 *          The values of the number of right hand sides NRHS.
00044 *
00045 *  THRESH  (input) DOUBLE PRECISION
00046 *          The threshold value for the test ratios.  A result is
00047 *          included in the output file if RESULT >= THRESH.  To have
00048 *          every test ratio printed, use THRESH = 0.
00049 *
00050 *  TSTERR  (input) LOGICAL
00051 *          Flag that indicates whether error exits are to be tested.
00052 *
00053 *  A       (workspace) DOUBLE PRECISION array, dimension (NMAX*4)
00054 *
00055 *  AF      (workspace) DOUBLE PRECISION array, dimension (NMAX*4)
00056 *
00057 *  B       (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
00058 *          where NSMAX is the largest entry in NSVAL.
00059 *
00060 *  X       (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
00061 *
00062 *  XACT    (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
00063 *
00064 *  WORK    (workspace) DOUBLE PRECISION array, dimension
00065 *                      (NMAX*max(3,NSMAX))
00066 *
00067 *  RWORK   (workspace) DOUBLE PRECISION array, dimension
00068 *                      (max(NMAX,2*NSMAX))
00069 *
00070 *  IWORK   (workspace) INTEGER array, dimension (2*NMAX)
00071 *
00072 *  NOUT    (input) INTEGER
00073 *          The unit number for output.
00074 *
00075 *  =====================================================================
00076 *
00077 *     .. Parameters ..
00078       DOUBLE PRECISION   ONE, ZERO
00079       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00080       INTEGER            NTYPES
00081       PARAMETER          ( NTYPES = 12 )
00082       INTEGER            NTESTS
00083       PARAMETER          ( NTESTS = 7 )
00084 *     ..
00085 *     .. Local Scalars ..
00086       LOGICAL            TRFCON, ZEROT
00087       CHARACTER          DIST, NORM, TRANS, TYPE
00088       CHARACTER*3        PATH
00089       INTEGER            I, IMAT, IN, INFO, IRHS, ITRAN, IX, IZERO, J,
00090      $                   K, KL, KOFF, KU, LDA, M, MODE, N, NERRS, NFAIL,
00091      $                   NIMAT, NRHS, NRUN
00092       DOUBLE PRECISION   AINVNM, ANORM, COND, RCOND, RCONDC, RCONDI,
00093      $                   RCONDO
00094 *     ..
00095 *     .. Local Arrays ..
00096       CHARACTER          TRANSS( 3 )
00097       INTEGER            ISEED( 4 ), ISEEDY( 4 )
00098       DOUBLE PRECISION   RESULT( NTESTS ), Z( 3 )
00099 *     ..
00100 *     .. External Functions ..
00101       DOUBLE PRECISION   DASUM, DGET06, DLANGT
00102       EXTERNAL           DASUM, DGET06, DLANGT
00103 *     ..
00104 *     .. External Subroutines ..
00105       EXTERNAL           ALAERH, ALAHD, ALASUM, DCOPY, DERRGE, DGET04,
00106      $                   DGTCON, DGTRFS, DGTT01, DGTT02, DGTT05, DGTTRF,
00107      $                   DGTTRS, DLACPY, DLAGTM, DLARNV, DLATB4, DLATMS,
00108      $                   DSCAL
00109 *     ..
00110 *     .. Intrinsic Functions ..
00111       INTRINSIC          MAX
00112 *     ..
00113 *     .. Scalars in Common ..
00114       LOGICAL            LERR, OK
00115       CHARACTER*32       SRNAMT
00116       INTEGER            INFOT, NUNIT
00117 *     ..
00118 *     .. Common blocks ..
00119       COMMON             / INFOC / INFOT, NUNIT, OK, LERR
00120       COMMON             / SRNAMC / SRNAMT
00121 *     ..
00122 *     .. Data statements ..
00123       DATA               ISEEDY / 0, 0, 0, 1 / , TRANSS / 'N', 'T',
00124      $                   'C' /
00125 *     ..
00126 *     .. Executable Statements ..
00127 *
00128       PATH( 1: 1 ) = 'Double precision'
00129       PATH( 2: 3 ) = 'GT'
00130       NRUN = 0
00131       NFAIL = 0
00132       NERRS = 0
00133       DO 10 I = 1, 4
00134          ISEED( I ) = ISEEDY( I )
00135    10 CONTINUE
00136 *
00137 *     Test the error exits
00138 *
00139       IF( TSTERR )
00140      $   CALL DERRGE( PATH, NOUT )
00141       INFOT = 0
00142 *
00143       DO 110 IN = 1, NN
00144 *
00145 *        Do for each value of N in NVAL.
00146 *
00147          N = NVAL( IN )
00148          M = MAX( N-1, 0 )
00149          LDA = MAX( 1, N )
00150          NIMAT = NTYPES
00151          IF( N.LE.0 )
00152      $      NIMAT = 1
00153 *
00154          DO 100 IMAT = 1, NIMAT
00155 *
00156 *           Do the tests only if DOTYPE( IMAT ) is true.
00157 *
00158             IF( .NOT.DOTYPE( IMAT ) )
00159      $         GO TO 100
00160 *
00161 *           Set up parameters with DLATB4.
00162 *
00163             CALL DLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
00164      $                   COND, DIST )
00165 *
00166             ZEROT = IMAT.GE.8 .AND. IMAT.LE.10
00167             IF( IMAT.LE.6 ) THEN
00168 *
00169 *              Types 1-6:  generate matrices of known condition number.
00170 *
00171                KOFF = MAX( 2-KU, 3-MAX( 1, N ) )
00172                SRNAMT = 'DLATMS'
00173                CALL DLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, COND,
00174      $                      ANORM, KL, KU, 'Z', AF( KOFF ), 3, WORK,
00175      $                      INFO )
00176 *
00177 *              Check the error code from DLATMS.
00178 *
00179                IF( INFO.NE.0 ) THEN
00180                   CALL ALAERH( PATH, 'DLATMS', INFO, 0, ' ', N, N, KL,
00181      $                         KU, -1, IMAT, NFAIL, NERRS, NOUT )
00182                   GO TO 100
00183                END IF
00184                IZERO = 0
00185 *
00186                IF( N.GT.1 ) THEN
00187                   CALL DCOPY( N-1, AF( 4 ), 3, A, 1 )
00188                   CALL DCOPY( N-1, AF( 3 ), 3, A( N+M+1 ), 1 )
00189                END IF
00190                CALL DCOPY( N, AF( 2 ), 3, A( M+1 ), 1 )
00191             ELSE
00192 *
00193 *              Types 7-12:  generate tridiagonal matrices with
00194 *              unknown condition numbers.
00195 *
00196                IF( .NOT.ZEROT .OR. .NOT.DOTYPE( 7 ) ) THEN
00197 *
00198 *                 Generate a matrix with elements from [-1,1].
00199 *
00200                   CALL DLARNV( 2, ISEED, N+2*M, A )
00201                   IF( ANORM.NE.ONE )
00202      $               CALL DSCAL( N+2*M, ANORM, A, 1 )
00203                ELSE IF( IZERO.GT.0 ) THEN
00204 *
00205 *                 Reuse the last matrix by copying back the zeroed out
00206 *                 elements.
00207 *
00208                   IF( IZERO.EQ.1 ) THEN
00209                      A( N ) = Z( 2 )
00210                      IF( N.GT.1 )
00211      $                  A( 1 ) = Z( 3 )
00212                   ELSE IF( IZERO.EQ.N ) THEN
00213                      A( 3*N-2 ) = Z( 1 )
00214                      A( 2*N-1 ) = Z( 2 )
00215                   ELSE
00216                      A( 2*N-2+IZERO ) = Z( 1 )
00217                      A( N-1+IZERO ) = Z( 2 )
00218                      A( IZERO ) = Z( 3 )
00219                   END IF
00220                END IF
00221 *
00222 *              If IMAT > 7, set one column of the matrix to 0.
00223 *
00224                IF( .NOT.ZEROT ) THEN
00225                   IZERO = 0
00226                ELSE IF( IMAT.EQ.8 ) THEN
00227                   IZERO = 1
00228                   Z( 2 ) = A( N )
00229                   A( N ) = ZERO
00230                   IF( N.GT.1 ) THEN
00231                      Z( 3 ) = A( 1 )
00232                      A( 1 ) = ZERO
00233                   END IF
00234                ELSE IF( IMAT.EQ.9 ) THEN
00235                   IZERO = N
00236                   Z( 1 ) = A( 3*N-2 )
00237                   Z( 2 ) = A( 2*N-1 )
00238                   A( 3*N-2 ) = ZERO
00239                   A( 2*N-1 ) = ZERO
00240                ELSE
00241                   IZERO = ( N+1 ) / 2
00242                   DO 20 I = IZERO, N - 1
00243                      A( 2*N-2+I ) = ZERO
00244                      A( N-1+I ) = ZERO
00245                      A( I ) = ZERO
00246    20             CONTINUE
00247                   A( 3*N-2 ) = ZERO
00248                   A( 2*N-1 ) = ZERO
00249                END IF
00250             END IF
00251 *
00252 *+    TEST 1
00253 *           Factor A as L*U and compute the ratio
00254 *              norm(L*U - A) / (n * norm(A) * EPS )
00255 *
00256             CALL DCOPY( N+2*M, A, 1, AF, 1 )
00257             SRNAMT = 'DGTTRF'
00258             CALL DGTTRF( N, AF, AF( M+1 ), AF( N+M+1 ), AF( N+2*M+1 ),
00259      $                   IWORK, INFO )
00260 *
00261 *           Check error code from DGTTRF.
00262 *
00263             IF( INFO.NE.IZERO )
00264      $         CALL ALAERH( PATH, 'DGTTRF', INFO, IZERO, ' ', N, N, 1,
00265      $                      1, -1, IMAT, NFAIL, NERRS, NOUT )
00266             TRFCON = INFO.NE.0
00267 *
00268             CALL DGTT01( N, A, A( M+1 ), A( N+M+1 ), AF, AF( M+1 ),
00269      $                   AF( N+M+1 ), AF( N+2*M+1 ), IWORK, WORK, LDA,
00270      $                   RWORK, RESULT( 1 ) )
00271 *
00272 *           Print the test ratio if it is .GE. THRESH.
00273 *
00274             IF( RESULT( 1 ).GE.THRESH ) THEN
00275                IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00276      $            CALL ALAHD( NOUT, PATH )
00277                WRITE( NOUT, FMT = 9999 )N, IMAT, 1, RESULT( 1 )
00278                NFAIL = NFAIL + 1
00279             END IF
00280             NRUN = NRUN + 1
00281 *
00282             DO 50 ITRAN = 1, 2
00283                TRANS = TRANSS( ITRAN )
00284                IF( ITRAN.EQ.1 ) THEN
00285                   NORM = 'O'
00286                ELSE
00287                   NORM = 'I'
00288                END IF
00289                ANORM = DLANGT( NORM, N, A, A( M+1 ), A( N+M+1 ) )
00290 *
00291                IF( .NOT.TRFCON ) THEN
00292 *
00293 *                 Use DGTTRS to solve for one column at a time of inv(A)
00294 *                 or inv(A^T), computing the maximum column sum as we
00295 *                 go.
00296 *
00297                   AINVNM = ZERO
00298                   DO 40 I = 1, N
00299                      DO 30 J = 1, N
00300                         X( J ) = ZERO
00301    30                CONTINUE
00302                      X( I ) = ONE
00303                      CALL DGTTRS( TRANS, N, 1, AF, AF( M+1 ),
00304      $                            AF( N+M+1 ), AF( N+2*M+1 ), IWORK, X,
00305      $                            LDA, INFO )
00306                      AINVNM = MAX( AINVNM, DASUM( N, X, 1 ) )
00307    40             CONTINUE
00308 *
00309 *                 Compute RCONDC = 1 / (norm(A) * norm(inv(A))
00310 *
00311                   IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
00312                      RCONDC = ONE
00313                   ELSE
00314                      RCONDC = ( ONE / ANORM ) / AINVNM
00315                   END IF
00316                   IF( ITRAN.EQ.1 ) THEN
00317                      RCONDO = RCONDC
00318                   ELSE
00319                      RCONDI = RCONDC
00320                   END IF
00321                ELSE
00322                   RCONDC = ZERO
00323                END IF
00324 *
00325 *+    TEST 7
00326 *              Estimate the reciprocal of the condition number of the
00327 *              matrix.
00328 *
00329                SRNAMT = 'DGTCON'
00330                CALL DGTCON( NORM, N, AF, AF( M+1 ), AF( N+M+1 ),
00331      $                      AF( N+2*M+1 ), IWORK, ANORM, RCOND, WORK,
00332      $                      IWORK( N+1 ), INFO )
00333 *
00334 *              Check error code from DGTCON.
00335 *
00336                IF( INFO.NE.0 )
00337      $            CALL ALAERH( PATH, 'DGTCON', INFO, 0, NORM, N, N, -1,
00338      $                         -1, -1, IMAT, NFAIL, NERRS, NOUT )
00339 *
00340                RESULT( 7 ) = DGET06( RCOND, RCONDC )
00341 *
00342 *              Print the test ratio if it is .GE. THRESH.
00343 *
00344                IF( RESULT( 7 ).GE.THRESH ) THEN
00345                   IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00346      $               CALL ALAHD( NOUT, PATH )
00347                   WRITE( NOUT, FMT = 9997 )NORM, N, IMAT, 7,
00348      $               RESULT( 7 )
00349                   NFAIL = NFAIL + 1
00350                END IF
00351                NRUN = NRUN + 1
00352    50       CONTINUE
00353 *
00354 *           Skip the remaining tests if the matrix is singular.
00355 *
00356             IF( TRFCON )
00357      $         GO TO 100
00358 *
00359             DO 90 IRHS = 1, NNS
00360                NRHS = NSVAL( IRHS )
00361 *
00362 *              Generate NRHS random solution vectors.
00363 *
00364                IX = 1
00365                DO 60 J = 1, NRHS
00366                   CALL DLARNV( 2, ISEED, N, XACT( IX ) )
00367                   IX = IX + LDA
00368    60          CONTINUE
00369 *
00370                DO 80 ITRAN = 1, 3
00371                   TRANS = TRANSS( ITRAN )
00372                   IF( ITRAN.EQ.1 ) THEN
00373                      RCONDC = RCONDO
00374                   ELSE
00375                      RCONDC = RCONDI
00376                   END IF
00377 *
00378 *                 Set the right hand side.
00379 *
00380                   CALL DLAGTM( TRANS, N, NRHS, ONE, A, A( M+1 ),
00381      $                         A( N+M+1 ), XACT, LDA, ZERO, B, LDA )
00382 *
00383 *+    TEST 2
00384 *                 Solve op(A) * X = B and compute the residual.
00385 *
00386                   CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
00387                   SRNAMT = 'DGTTRS'
00388                   CALL DGTTRS( TRANS, N, NRHS, AF, AF( M+1 ),
00389      $                         AF( N+M+1 ), AF( N+2*M+1 ), IWORK, X,
00390      $                         LDA, INFO )
00391 *
00392 *                 Check error code from DGTTRS.
00393 *
00394                   IF( INFO.NE.0 )
00395      $               CALL ALAERH( PATH, 'DGTTRS', INFO, 0, TRANS, N, N,
00396      $                            -1, -1, NRHS, IMAT, NFAIL, NERRS,
00397      $                            NOUT )
00398 *
00399                   CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
00400                   CALL DGTT02( TRANS, N, NRHS, A, A( M+1 ), A( N+M+1 ),
00401      $                         X, LDA, WORK, LDA, RWORK, RESULT( 2 ) )
00402 *
00403 *+    TEST 3
00404 *                 Check solution from generated exact solution.
00405 *
00406                   CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
00407      $                         RESULT( 3 ) )
00408 *
00409 *+    TESTS 4, 5, and 6
00410 *                 Use iterative refinement to improve the solution.
00411 *
00412                   SRNAMT = 'DGTRFS'
00413                   CALL DGTRFS( TRANS, N, NRHS, A, A( M+1 ), A( N+M+1 ),
00414      $                         AF, AF( M+1 ), AF( N+M+1 ),
00415      $                         AF( N+2*M+1 ), IWORK, B, LDA, X, LDA,
00416      $                         RWORK, RWORK( NRHS+1 ), WORK,
00417      $                         IWORK( N+1 ), INFO )
00418 *
00419 *                 Check error code from DGTRFS.
00420 *
00421                   IF( INFO.NE.0 )
00422      $               CALL ALAERH( PATH, 'DGTRFS', INFO, 0, TRANS, N, N,
00423      $                            -1, -1, NRHS, IMAT, NFAIL, NERRS,
00424      $                            NOUT )
00425 *
00426                   CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
00427      $                         RESULT( 4 ) )
00428                   CALL DGTT05( TRANS, N, NRHS, A, A( M+1 ), A( N+M+1 ),
00429      $                         B, LDA, X, LDA, XACT, LDA, RWORK,
00430      $                         RWORK( NRHS+1 ), RESULT( 5 ) )
00431 *
00432 *                 Print information about the tests that did not pass
00433 *                 the threshold.
00434 *
00435                   DO 70 K = 2, 6
00436                      IF( RESULT( K ).GE.THRESH ) THEN
00437                         IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00438      $                     CALL ALAHD( NOUT, PATH )
00439                         WRITE( NOUT, FMT = 9998 )TRANS, N, NRHS, IMAT,
00440      $                     K, RESULT( K )
00441                         NFAIL = NFAIL + 1
00442                      END IF
00443    70             CONTINUE
00444                   NRUN = NRUN + 5
00445    80          CONTINUE
00446    90       CONTINUE
00447 *
00448   100    CONTINUE
00449   110 CONTINUE
00450 *
00451 *     Print a summary of the results.
00452 *
00453       CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
00454 *
00455  9999 FORMAT( 12X, 'N =', I5, ',', 10X, ' type ', I2, ', test(', I2,
00456      $      ') = ', G12.5 )
00457  9998 FORMAT( ' TRANS=''', A1, ''', N =', I5, ', NRHS=', I3, ', type ',
00458      $      I2, ', test(', I2, ') = ', G12.5 )
00459  9997 FORMAT( ' NORM =''', A1, ''', N =', I5, ',', 10X, ' type ', I2,
00460      $      ', test(', I2, ') = ', G12.5 )
00461       RETURN
00462 *
00463 *     End of DCHKGT
00464 *
00465       END