LAPACK: TESTING/LIN/zchksp.f Source File

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