LAPACK: TESTING/LIN/schkge.f Source File

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