LAPACK: TESTING/LIN/zchktb.f Source File

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