LAPACK: TESTING/LIN/dchktr.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DCHKTR( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL,
00002      $                   THRESH, TSTERR, NMAX, A, AINV, B, X, XACT,
00003      $                   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( * ), AINV( * ), B( * ), RWORK( * ),
00018      $                   WORK( * ), X( * ), XACT( * )
00019 *     ..
00020 *
00021 *  Purpose
00022 *  =======
00023 *
00024 *  DCHKTR tests DTRTRI, -TRS, -RFS, and -CON, and DLATRS
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 column 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 (NNB)
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 leading dimension of the work arrays.
00062 *          NMAX >= the maximum value of N in NVAL.
00063 *
00064 *  A       (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
00065 *
00066 *  AINV    (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
00067 *
00068 *  B       (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
00069 *          where NSMAX is the largest entry in NSVAL.
00070 *
00071 *  X       (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
00072 *
00073 *  XACT    (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
00074 *
00075 *  WORK    (workspace) DOUBLE PRECISION array, dimension
00076 *                      (NMAX*max(3,NSMAX))
00077 *
00078 *  RWORK   (workspace) DOUBLE PRECISION array, dimension
00079 *                      (max(NMAX,2*NSMAX))
00080 *
00081 *  IWORK   (workspace) INTEGER array, dimension (NMAX)
00082 *
00083 *  NOUT    (input) INTEGER
00084 *          The unit number for output.
00085 *
00086 *  =====================================================================
00087 *
00088 *     .. Parameters ..
00089       INTEGER            NTYPE1, NTYPES
00090       PARAMETER          ( NTYPE1 = 10, NTYPES = 18 )
00091       INTEGER            NTESTS
00092       PARAMETER          ( NTESTS = 9 )
00093       INTEGER            NTRAN
00094       PARAMETER          ( NTRAN = 3 )
00095       DOUBLE PRECISION   ONE, ZERO
00096       PARAMETER          ( ONE = 1.0D0, ZERO = 0.0D0 )
00097 *     ..
00098 *     .. Local Scalars ..
00099       CHARACTER          DIAG, NORM, TRANS, UPLO, XTYPE
00100       CHARACTER*3        PATH
00101       INTEGER            I, IDIAG, IMAT, IN, INB, INFO, IRHS, ITRAN,
00102      $                   IUPLO, K, LDA, N, NB, NERRS, NFAIL, NRHS, NRUN
00103       DOUBLE PRECISION   AINVNM, ANORM, DUMMY, RCOND, RCONDC, RCONDI,
00104      $                   RCONDO, SCALE
00105 *     ..
00106 *     .. Local Arrays ..
00107       CHARACTER          TRANSS( NTRAN ), UPLOS( 2 )
00108       INTEGER            ISEED( 4 ), ISEEDY( 4 )
00109       DOUBLE PRECISION   RESULT( NTESTS )
00110 *     ..
00111 *     .. External Functions ..
00112       LOGICAL            LSAME
00113       DOUBLE PRECISION   DLANTR
00114       EXTERNAL           LSAME, DLANTR
00115 *     ..
00116 *     .. External Subroutines ..
00117       EXTERNAL           ALAERH, ALAHD, ALASUM, DCOPY, DERRTR, DGET04,
00118      $                   DLACPY, DLARHS, DLATRS, DLATTR, DTRCON, DTRRFS,
00119      $                   DTRT01, DTRT02, DTRT03, DTRT05, DTRT06, DTRTRI,
00120      $                   DTRTRS, XLAENV
00121 *     ..
00122 *     .. Scalars in Common ..
00123       LOGICAL            LERR, OK
00124       CHARACTER*32       SRNAMT
00125       INTEGER            INFOT, IOUNIT
00126 *     ..
00127 *     .. Common blocks ..
00128       COMMON             / INFOC / INFOT, IOUNIT, 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' / , TRANSS / 'N', 'T', 'C' /
00137 *     ..
00138 *     .. Executable Statements ..
00139 *
00140 *     Initialize constants and the random number seed.
00141 *
00142       PATH( 1: 1 ) = 'Double precision'
00143       PATH( 2: 3 ) = 'TR'
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 DERRTR( PATH, NOUT )
00155       INFOT = 0
00156       CALL XLAENV( 2, 2 )
00157 *
00158       DO 120 IN = 1, NN
00159 *
00160 *        Do for each value of N in NVAL
00161 *
00162          N = NVAL( IN )
00163          LDA = MAX( 1, N )
00164          XTYPE = 'N'
00165 *
00166          DO 80 IMAT = 1, NTYPE1
00167 *
00168 *           Do the tests only if DOTYPE( IMAT ) is true.
00169 *
00170             IF( .NOT.DOTYPE( IMAT ) )
00171      $         GO TO 80
00172 *
00173             DO 70 IUPLO = 1, 2
00174 *
00175 *              Do first for UPLO = 'U', then for UPLO = 'L'
00176 *
00177                UPLO = UPLOS( IUPLO )
00178 *
00179 *              Call DLATTR to generate a triangular test matrix.
00180 *
00181                SRNAMT = 'DLATTR'
00182                CALL DLATTR( IMAT, UPLO, 'No transpose', DIAG, ISEED, N,
00183      $                      A, LDA, X, WORK, INFO )
00184 *
00185 *              Set IDIAG = 1 for non-unit matrices, 2 for unit.
00186 *
00187                IF( LSAME( DIAG, 'N' ) ) THEN
00188                   IDIAG = 1
00189                ELSE
00190                   IDIAG = 2
00191                END IF
00192 *
00193                DO 60 INB = 1, NNB
00194 *
00195 *                 Do for each blocksize in NBVAL
00196 *
00197                   NB = NBVAL( INB )
00198                   CALL XLAENV( 1, NB )
00199 *
00200 *+    TEST 1
00201 *                 Form the inverse of A.
00202 *
00203                   CALL DLACPY( UPLO, N, N, A, LDA, AINV, LDA )
00204                   SRNAMT = 'DTRTRI'
00205                   CALL DTRTRI( UPLO, DIAG, N, AINV, LDA, INFO )
00206 *
00207 *                 Check error code from DTRTRI.
00208 *
00209                   IF( INFO.NE.0 )
00210      $               CALL ALAERH( PATH, 'DTRTRI', INFO, 0, UPLO // DIAG,
00211      $                            N, N, -1, -1, NB, IMAT, NFAIL, NERRS,
00212      $                            NOUT )
00213 *
00214 *                 Compute the infinity-norm condition number of A.
00215 *
00216                   ANORM = DLANTR( 'I', UPLO, DIAG, N, N, A, LDA, RWORK )
00217                   AINVNM = DLANTR( 'I', UPLO, DIAG, N, N, AINV, LDA,
00218      $                     RWORK )
00219                   IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
00220                      RCONDI = ONE
00221                   ELSE
00222                      RCONDI = ( ONE / ANORM ) / AINVNM
00223                   END IF
00224 *
00225 *                 Compute the residual for the triangular matrix times
00226 *                 its inverse.  Also compute the 1-norm condition number
00227 *                 of A.
00228 *
00229                   CALL DTRT01( UPLO, DIAG, N, A, LDA, AINV, LDA, RCONDO,
00230      $                         RWORK, RESULT( 1 ) )
00231 *
00232 *                 Print the test ratio if it is .GE. THRESH.
00233 *
00234                   IF( RESULT( 1 ).GE.THRESH ) THEN
00235                      IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00236      $                  CALL ALAHD( NOUT, PATH )
00237                      WRITE( NOUT, FMT = 9999 )UPLO, DIAG, N, NB, IMAT,
00238      $                  1, RESULT( 1 )
00239                      NFAIL = NFAIL + 1
00240                   END IF
00241                   NRUN = NRUN + 1
00242 *
00243 *                 Skip remaining tests if not the first block size.
00244 *
00245                   IF( INB.NE.1 )
00246      $               GO TO 60
00247 *
00248                   DO 40 IRHS = 1, NNS
00249                      NRHS = NSVAL( IRHS )
00250                      XTYPE = 'N'
00251 *
00252                      DO 30 ITRAN = 1, NTRAN
00253 *
00254 *                    Do for op(A) = A, A**T, or A**H.
00255 *
00256                         TRANS = TRANSS( ITRAN )
00257                         IF( ITRAN.EQ.1 ) THEN
00258                            NORM = 'O'
00259                            RCONDC = RCONDO
00260                         ELSE
00261                            NORM = 'I'
00262                            RCONDC = RCONDI
00263                         END IF
00264 *
00265 *+    TEST 2
00266 *                       Solve and compute residual for op(A)*x = b.
00267 *
00268                         SRNAMT = 'DLARHS'
00269                         CALL DLARHS( PATH, XTYPE, UPLO, TRANS, N, N, 0,
00270      $                               IDIAG, NRHS, A, LDA, XACT, LDA, B,
00271      $                               LDA, ISEED, INFO )
00272                         XTYPE = 'C'
00273                         CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
00274 *
00275                         SRNAMT = 'DTRTRS'
00276                         CALL DTRTRS( UPLO, TRANS, DIAG, N, NRHS, A, LDA,
00277      $                               X, LDA, INFO )
00278 *
00279 *                       Check error code from DTRTRS.
00280 *
00281                         IF( INFO.NE.0 )
00282      $                     CALL ALAERH( PATH, 'DTRTRS', INFO, 0,
00283      $                                  UPLO // TRANS // DIAG, N, N, -1,
00284      $                                  -1, NRHS, IMAT, NFAIL, NERRS,
00285      $                                  NOUT )
00286 *
00287 *                       This line is needed on a Sun SPARCstation.
00288 *
00289                         IF( N.GT.0 )
00290      $                     DUMMY = A( 1 )
00291 *
00292                         CALL DTRT02( UPLO, TRANS, DIAG, N, NRHS, A, LDA,
00293      $                               X, LDA, B, LDA, WORK, RESULT( 2 ) )
00294 *
00295 *+    TEST 3
00296 *                       Check solution from generated exact solution.
00297 *
00298                         CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
00299      $                               RESULT( 3 ) )
00300 *
00301 *+    TESTS 4, 5, and 6
00302 *                       Use iterative refinement to improve the solution
00303 *                       and compute error bounds.
00304 *
00305                         SRNAMT = 'DTRRFS'
00306                         CALL DTRRFS( UPLO, TRANS, DIAG, N, NRHS, A, LDA,
00307      $                               B, LDA, X, LDA, RWORK,
00308      $                               RWORK( NRHS+1 ), WORK, IWORK,
00309      $                               INFO )
00310 *
00311 *                       Check error code from DTRRFS.
00312 *
00313                         IF( INFO.NE.0 )
00314      $                     CALL ALAERH( PATH, 'DTRRFS', INFO, 0,
00315      $                                  UPLO // TRANS // DIAG, N, N, -1,
00316      $                                  -1, NRHS, IMAT, NFAIL, NERRS,
00317      $                                  NOUT )
00318 *
00319                         CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
00320      $                               RESULT( 4 ) )
00321                         CALL DTRT05( UPLO, TRANS, DIAG, N, NRHS, A, LDA,
00322      $                               B, LDA, X, LDA, XACT, LDA, RWORK,
00323      $                               RWORK( NRHS+1 ), RESULT( 5 ) )
00324 *
00325 *                       Print information about the tests that did not
00326 *                       pass the threshold.
00327 *
00328                         DO 20 K = 2, 6
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 = 9998 )UPLO, TRANS,
00333      $                           DIAG, N, NRHS, IMAT, K, RESULT( K )
00334                               NFAIL = NFAIL + 1
00335                            END IF
00336    20                   CONTINUE
00337                         NRUN = NRUN + 5
00338    30                CONTINUE
00339    40             CONTINUE
00340 *
00341 *+    TEST 7
00342 *                       Get an estimate of RCOND = 1/CNDNUM.
00343 *
00344                   DO 50 ITRAN = 1, 2
00345                      IF( ITRAN.EQ.1 ) THEN
00346                         NORM = 'O'
00347                         RCONDC = RCONDO
00348                      ELSE
00349                         NORM = 'I'
00350                         RCONDC = RCONDI
00351                      END IF
00352                      SRNAMT = 'DTRCON'
00353                      CALL DTRCON( NORM, UPLO, DIAG, N, A, LDA, RCOND,
00354      $                            WORK, IWORK, INFO )
00355 *
00356 *                       Check error code from DTRCON.
00357 *
00358                      IF( INFO.NE.0 )
00359      $                  CALL ALAERH( PATH, 'DTRCON', INFO, 0,
00360      $                               NORM // UPLO // DIAG, N, N, -1, -1,
00361      $                               -1, IMAT, NFAIL, NERRS, NOUT )
00362 *
00363                      CALL DTRT06( RCOND, RCONDC, UPLO, DIAG, N, A, LDA,
00364      $                            RWORK, RESULT( 7 ) )
00365 *
00366 *                    Print the test ratio if it is .GE. THRESH.
00367 *
00368                      IF( RESULT( 7 ).GE.THRESH ) THEN
00369                         IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00370      $                     CALL ALAHD( NOUT, PATH )
00371                         WRITE( NOUT, FMT = 9997 )NORM, UPLO, N, IMAT,
00372      $                     7, RESULT( 7 )
00373                         NFAIL = NFAIL + 1
00374                      END IF
00375                      NRUN = NRUN + 1
00376    50             CONTINUE
00377    60          CONTINUE
00378    70       CONTINUE
00379    80    CONTINUE
00380 *
00381 *        Use pathological test matrices to test DLATRS.
00382 *
00383          DO 110 IMAT = NTYPE1 + 1, NTYPES
00384 *
00385 *           Do the tests only if DOTYPE( IMAT ) is true.
00386 *
00387             IF( .NOT.DOTYPE( IMAT ) )
00388      $         GO TO 110
00389 *
00390             DO 100 IUPLO = 1, 2
00391 *
00392 *              Do first for UPLO = 'U', then for UPLO = 'L'
00393 *
00394                UPLO = UPLOS( IUPLO )
00395                DO 90 ITRAN = 1, NTRAN
00396 *
00397 *                 Do for op(A) = A, A**T, and A**H.
00398 *
00399                   TRANS = TRANSS( ITRAN )
00400 *
00401 *                 Call DLATTR to generate a triangular test matrix.
00402 *
00403                   SRNAMT = 'DLATTR'
00404                   CALL DLATTR( IMAT, UPLO, TRANS, DIAG, ISEED, N, A,
00405      $                         LDA, X, WORK, INFO )
00406 *
00407 *+    TEST 8
00408 *                 Solve the system op(A)*x = b.
00409 *
00410                   SRNAMT = 'DLATRS'
00411                   CALL DCOPY( N, X, 1, B, 1 )
00412                   CALL DLATRS( UPLO, TRANS, DIAG, 'N', N, A, LDA, B,
00413      $                         SCALE, RWORK, INFO )
00414 *
00415 *                 Check error code from DLATRS.
00416 *
00417                   IF( INFO.NE.0 )
00418      $               CALL ALAERH( PATH, 'DLATRS', INFO, 0,
00419      $                            UPLO // TRANS // DIAG // 'N', N, N,
00420      $                            -1, -1, -1, IMAT, NFAIL, NERRS, NOUT )
00421 *
00422                   CALL DTRT03( UPLO, TRANS, DIAG, N, 1, A, LDA, SCALE,
00423      $                         RWORK, ONE, B, LDA, X, LDA, WORK,
00424      $                         RESULT( 8 ) )
00425 *
00426 *+    TEST 9
00427 *                 Solve op(A)*X = b again with NORMIN = 'Y'.
00428 *
00429                   CALL DCOPY( N, X, 1, B( N+1 ), 1 )
00430                   CALL DLATRS( UPLO, TRANS, DIAG, 'Y', N, A, LDA,
00431      $                         B( N+1 ), SCALE, RWORK, INFO )
00432 *
00433 *                 Check error code from DLATRS.
00434 *
00435                   IF( INFO.NE.0 )
00436      $               CALL ALAERH( PATH, 'DLATRS', INFO, 0,
00437      $                            UPLO // TRANS // DIAG // 'Y', N, N,
00438      $                            -1, -1, -1, IMAT, NFAIL, NERRS, NOUT )
00439 *
00440                   CALL DTRT03( UPLO, TRANS, DIAG, N, 1, A, LDA, SCALE,
00441      $                         RWORK, ONE, B( N+1 ), LDA, X, LDA, WORK,
00442      $                         RESULT( 9 ) )
00443 *
00444 *                 Print information about the tests that did not pass
00445 *                 the threshold.
00446 *
00447                   IF( RESULT( 8 ).GE.THRESH ) THEN
00448                      IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00449      $                  CALL ALAHD( NOUT, PATH )
00450                      WRITE( NOUT, FMT = 9996 )'DLATRS', UPLO, TRANS,
00451      $                  DIAG, 'N', N, IMAT, 8, RESULT( 8 )
00452                      NFAIL = NFAIL + 1
00453                   END IF
00454                   IF( RESULT( 9 ).GE.THRESH ) THEN
00455                      IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00456      $                  CALL ALAHD( NOUT, PATH )
00457                      WRITE( NOUT, FMT = 9996 )'DLATRS', UPLO, TRANS,
00458      $                  DIAG, 'Y', N, IMAT, 9, RESULT( 9 )
00459                      NFAIL = NFAIL + 1
00460                   END IF
00461                   NRUN = NRUN + 2
00462    90          CONTINUE
00463   100       CONTINUE
00464   110    CONTINUE
00465   120 CONTINUE
00466 *
00467 *     Print a summary of the results.
00468 *
00469       CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
00470 *
00471  9999 FORMAT( ' UPLO=''', A1, ''', DIAG=''', A1, ''', N=', I5, ', NB=',
00472      $      I4, ', type ', I2, ', test(', I2, ')= ', G12.5 )
00473  9998 FORMAT( ' UPLO=''', A1, ''', TRANS=''', A1, ''', DIAG=''', A1,
00474      $      ''', N=', I5, ', NB=', I4, ', type ', I2, 
00475 ',     $      test(', I2, ')= ', G12.5 )
00476  9997 FORMAT( ' NORM=''', A1, ''', UPLO =''', A1, ''', N=', I5, ',',
00477      $      11X, ' type ', I2, ', test(', I2, ')=', G12.5 )
00478  9996 FORMAT( 1X, A, '( ''', A1, ''', ''', A1, ''', ''', A1, ''', ''',
00479      $      A1, ''',', I5, ', ... ), type ', I2, ', test(', I2, ')=',
00480      $      G12.5 )
00481       RETURN
00482 *
00483 *     End of DCHKTR
00484 *
00485       END