LAPACK: TESTING/LIN/zchkge.f Source File

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