schkpp.c
Go to the documentation of this file.
00001 /* schkpp.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Common Block Declarations */
00017 
00018 struct {
00019     integer infot, nunit;
00020     logical ok, lerr;
00021 } infoc_;
00022 
00023 #define infoc_1 infoc_
00024 
00025 struct {
00026     char srnamt[32];
00027 } srnamc_;
00028 
00029 #define srnamc_1 srnamc_
00030 
00031 /* Table of constant values */
00032 
00033 static integer c__0 = 0;
00034 static integer c_n1 = -1;
00035 static integer c__1 = 1;
00036 static integer c__8 = 8;
00037 
00038 /* Subroutine */ int schkpp_(logical *dotype, integer *nn, integer *nval, 
00039         integer *nns, integer *nsval, real *thresh, logical *tsterr, integer *
00040         nmax, real *a, real *afac, real *ainv, real *b, real *x, real *xact, 
00041         real *work, real *rwork, integer *iwork, integer *nout)
00042 {
00043     /* Initialized data */
00044 
00045     static integer iseedy[4] = { 1988,1989,1990,1991 };
00046     static char uplos[1*2] = "U" "L";
00047     static char packs[1*2] = "C" "R";
00048 
00049     /* Format strings */
00050     static char fmt_9999[] = "(\002 UPLO = '\002,a1,\002', N =\002,i5,\002, "
00051             "type \002,i2,\002, test \002,i2,\002, ratio =\002,g12.5)";
00052     static char fmt_9998[] = "(\002 UPLO = '\002,a1,\002', N =\002,i5,\002, "
00053             "NRHS=\002,i3,\002, type \002,i2,\002, test(\002,i2,\002) =\002,g"
00054             "12.5)";
00055 
00056     /* System generated locals */
00057     integer i__1, i__2, i__3;
00058 
00059     /* Builtin functions */
00060     /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
00061     integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
00062 
00063     /* Local variables */
00064     integer i__, k, n, in, kl, ku, lda, npp, ioff, mode, imat, info;
00065     char path[3], dist[1];
00066     integer irhs, nrhs;
00067     char uplo[1], type__[1];
00068     integer nrun;
00069     extern /* Subroutine */ int alahd_(integer *, char *);
00070     integer nfail, iseed[4];
00071     real rcond;
00072     extern /* Subroutine */ int sget04_(integer *, integer *, real *, integer 
00073             *, real *, integer *, real *, real *);
00074     integer nimat;
00075     extern doublereal sget06_(real *, real *);
00076     real anorm;
00077     extern /* Subroutine */ int sppt01_(char *, integer *, real *, real *, 
00078             real *, real *);
00079     integer iuplo, izero, nerrs;
00080     extern /* Subroutine */ int sppt02_(char *, integer *, integer *, real *, 
00081             real *, integer *, real *, integer *, real *, real *), 
00082             scopy_(integer *, real *, integer *, real *, integer *), sppt03_(
00083             char *, integer *, real *, real *, real *, integer *, real *, 
00084             real *, real *), sppt05_(char *, integer *, integer *, 
00085             real *, real *, integer *, real *, integer *, real *, integer *, 
00086             real *, real *, real *);
00087     logical zerot;
00088     char xtype[1];
00089     extern /* Subroutine */ int slatb4_(char *, integer *, integer *, integer 
00090             *, char *, integer *, integer *, real *, integer *, real *, char *
00091 ), alaerh_(char *, char *, integer *, 
00092             integer *, char *, integer *, integer *, integer *, integer *, 
00093             integer *, integer *, integer *, integer *, integer *);
00094     real rcondc;
00095     char packit[1];
00096     extern /* Subroutine */ int alasum_(char *, integer *, integer *, integer 
00097             *, integer *);
00098     real cndnum;
00099     extern /* Subroutine */ int slacpy_(char *, integer *, integer *, real *, 
00100             integer *, real *, integer *), slarhs_(char *, char *, 
00101             char *, char *, integer *, integer *, integer *, integer *, 
00102             integer *, real *, integer *, real *, integer *, real *, integer *
00103 , integer *, integer *);
00104     extern doublereal slansp_(char *, char *, integer *, real *, real *);
00105     extern /* Subroutine */ int sppcon_(char *, integer *, real *, real *, 
00106             real *, real *, integer *, integer *), slatms_(integer *, 
00107             integer *, char *, integer *, char *, real *, integer *, real *, 
00108             real *, integer *, integer *, char *, real *, integer *, real *, 
00109             integer *), serrpo_(char *, integer *), spprfs_(char *, integer *, integer *, real *, real *, 
00110             real *, integer *, real *, integer *, real *, real *, real *, 
00111             integer *, integer *);
00112     real result[8];
00113     extern /* Subroutine */ int spptrf_(char *, integer *, real *, integer *), spptri_(char *, integer *, real *, integer *), 
00114             spptrs_(char *, integer *, integer *, real *, real *, integer *, 
00115             integer *);
00116 
00117     /* Fortran I/O blocks */
00118     static cilist io___34 = { 0, 0, 0, fmt_9999, 0 };
00119     static cilist io___37 = { 0, 0, 0, fmt_9998, 0 };
00120     static cilist io___39 = { 0, 0, 0, fmt_9999, 0 };
00121 
00122 
00123 
00124 /*  -- LAPACK test routine (version 3.1) -- */
00125 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00126 /*     November 2006 */
00127 
00128 /*     .. Scalar Arguments .. */
00129 /*     .. */
00130 /*     .. Array Arguments .. */
00131 /*     .. */
00132 
00133 /*  Purpose */
00134 /*  ======= */
00135 
00136 /*  SCHKPP tests SPPTRF, -TRI, -TRS, -RFS, and -CON */
00137 
00138 /*  Arguments */
00139 /*  ========= */
00140 
00141 /*  DOTYPE  (input) LOGICAL array, dimension (NTYPES) */
00142 /*          The matrix types to be used for testing.  Matrices of type j */
00143 /*          (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = */
00144 /*          .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. */
00145 
00146 /*  NN      (input) INTEGER */
00147 /*          The number of values of N contained in the vector NVAL. */
00148 
00149 /*  NVAL    (input) INTEGER array, dimension (NN) */
00150 /*          The values of the matrix dimension N. */
00151 
00152 /*  NNS     (input) INTEGER */
00153 /*          The number of values of NRHS contained in the vector NSVAL. */
00154 
00155 /*  NSVAL   (input) INTEGER array, dimension (NNS) */
00156 /*          The values of the number of right hand sides NRHS. */
00157 
00158 /*  THRESH  (input) REAL */
00159 /*          The threshold value for the test ratios.  A result is */
00160 /*          included in the output file if RESULT >= THRESH.  To have */
00161 /*          every test ratio printed, use THRESH = 0. */
00162 
00163 /*  TSTERR  (input) LOGICAL */
00164 /*          Flag that indicates whether error exits are to be tested. */
00165 
00166 /*  NMAX    (input) INTEGER */
00167 /*          The maximum value permitted for N, used in dimensioning the */
00168 /*          work arrays. */
00169 
00170 /*  A       (workspace) REAL array, dimension */
00171 /*                      (NMAX*(NMAX+1)/2) */
00172 
00173 /*  AFAC    (workspace) REAL array, dimension */
00174 /*                      (NMAX*(NMAX+1)/2) */
00175 
00176 /*  AINV    (workspace) REAL array, dimension */
00177 /*                      (NMAX*(NMAX+1)/2) */
00178 
00179 /*  B       (workspace) REAL array, dimension (NMAX*NSMAX) */
00180 /*          where NSMAX is the largest entry in NSVAL. */
00181 
00182 /*  X       (workspace) REAL array, dimension (NMAX*NSMAX) */
00183 
00184 /*  XACT    (workspace) REAL array, dimension (NMAX*NSMAX) */
00185 
00186 /*  WORK    (workspace) REAL array, dimension */
00187 /*                      (NMAX*max(3,NSMAX)) */
00188 
00189 /*  RWORK   (workspace) REAL array, dimension */
00190 /*                      (max(NMAX,2*NSMAX)) */
00191 
00192 /*  IWORK   (workspace) INTEGER array, dimension (NMAX) */
00193 
00194 /*  NOUT    (input) INTEGER */
00195 /*          The unit number for output. */
00196 
00197 /*  ===================================================================== */
00198 
00199 /*     .. Parameters .. */
00200 /*     .. */
00201 /*     .. Local Scalars .. */
00202 /*     .. */
00203 /*     .. Local Arrays .. */
00204 /*     .. */
00205 /*     .. External Functions .. */
00206 /*     .. */
00207 /*     .. External Subroutines .. */
00208 /*     .. */
00209 /*     .. Scalars in Common .. */
00210 /*     .. */
00211 /*     .. Common blocks .. */
00212 /*     .. */
00213 /*     .. Intrinsic Functions .. */
00214 /*     .. */
00215 /*     .. Data statements .. */
00216     /* Parameter adjustments */
00217     --iwork;
00218     --rwork;
00219     --work;
00220     --xact;
00221     --x;
00222     --b;
00223     --ainv;
00224     --afac;
00225     --a;
00226     --nsval;
00227     --nval;
00228     --dotype;
00229 
00230     /* Function Body */
00231 /*     .. */
00232 /*     .. Executable Statements .. */
00233 
00234 /*     Initialize constants and the random number seed. */
00235 
00236     s_copy(path, "Single precision", (ftnlen)1, (ftnlen)16);
00237     s_copy(path + 1, "PP", (ftnlen)2, (ftnlen)2);
00238     nrun = 0;
00239     nfail = 0;
00240     nerrs = 0;
00241     for (i__ = 1; i__ <= 4; ++i__) {
00242         iseed[i__ - 1] = iseedy[i__ - 1];
00243 /* L10: */
00244     }
00245 
00246 /*     Test the error exits */
00247 
00248     if (*tsterr) {
00249         serrpo_(path, nout);
00250     }
00251     infoc_1.infot = 0;
00252 
00253 /*     Do for each value of N in NVAL */
00254 
00255     i__1 = *nn;
00256     for (in = 1; in <= i__1; ++in) {
00257         n = nval[in];
00258         lda = max(n,1);
00259         *(unsigned char *)xtype = 'N';
00260         nimat = 9;
00261         if (n <= 0) {
00262             nimat = 1;
00263         }
00264 
00265         i__2 = nimat;
00266         for (imat = 1; imat <= i__2; ++imat) {
00267 
00268 /*           Do the tests only if DOTYPE( IMAT ) is true. */
00269 
00270             if (! dotype[imat]) {
00271                 goto L100;
00272             }
00273 
00274 /*           Skip types 3, 4, or 5 if the matrix size is too small. */
00275 
00276             zerot = imat >= 3 && imat <= 5;
00277             if (zerot && n < imat - 2) {
00278                 goto L100;
00279             }
00280 
00281 /*           Do first for UPLO = 'U', then for UPLO = 'L' */
00282 
00283             for (iuplo = 1; iuplo <= 2; ++iuplo) {
00284                 *(unsigned char *)uplo = *(unsigned char *)&uplos[iuplo - 1];
00285                 *(unsigned char *)packit = *(unsigned char *)&packs[iuplo - 1]
00286                         ;
00287 
00288 /*              Set up parameters with SLATB4 and generate a test matrix */
00289 /*              with SLATMS. */
00290 
00291                 slatb4_(path, &imat, &n, &n, type__, &kl, &ku, &anorm, &mode, 
00292                         &cndnum, dist);
00293 
00294                 s_copy(srnamc_1.srnamt, "SLATMS", (ftnlen)32, (ftnlen)6);
00295                 slatms_(&n, &n, dist, iseed, type__, &rwork[1], &mode, &
00296                         cndnum, &anorm, &kl, &ku, packit, &a[1], &lda, &work[
00297                         1], &info);
00298 
00299 /*              Check error code from SLATMS. */
00300 
00301                 if (info != 0) {
00302                     alaerh_(path, "SLATMS", &info, &c__0, uplo, &n, &n, &c_n1, 
00303                              &c_n1, &c_n1, &imat, &nfail, &nerrs, nout);
00304                     goto L90;
00305                 }
00306 
00307 /*              For types 3-5, zero one row and column of the matrix to */
00308 /*              test that INFO is returned correctly. */
00309 
00310                 if (zerot) {
00311                     if (imat == 3) {
00312                         izero = 1;
00313                     } else if (imat == 4) {
00314                         izero = n;
00315                     } else {
00316                         izero = n / 2 + 1;
00317                     }
00318 
00319 /*                 Set row and column IZERO of A to 0. */
00320 
00321                     if (iuplo == 1) {
00322                         ioff = (izero - 1) * izero / 2;
00323                         i__3 = izero - 1;
00324                         for (i__ = 1; i__ <= i__3; ++i__) {
00325                             a[ioff + i__] = 0.f;
00326 /* L20: */
00327                         }
00328                         ioff += izero;
00329                         i__3 = n;
00330                         for (i__ = izero; i__ <= i__3; ++i__) {
00331                             a[ioff] = 0.f;
00332                             ioff += i__;
00333 /* L30: */
00334                         }
00335                     } else {
00336                         ioff = izero;
00337                         i__3 = izero - 1;
00338                         for (i__ = 1; i__ <= i__3; ++i__) {
00339                             a[ioff] = 0.f;
00340                             ioff = ioff + n - i__;
00341 /* L40: */
00342                         }
00343                         ioff -= izero;
00344                         i__3 = n;
00345                         for (i__ = izero; i__ <= i__3; ++i__) {
00346                             a[ioff + i__] = 0.f;
00347 /* L50: */
00348                         }
00349                     }
00350                 } else {
00351                     izero = 0;
00352                 }
00353 
00354 /*              Compute the L*L' or U'*U factorization of the matrix. */
00355 
00356                 npp = n * (n + 1) / 2;
00357                 scopy_(&npp, &a[1], &c__1, &afac[1], &c__1);
00358                 s_copy(srnamc_1.srnamt, "SPPTRF", (ftnlen)32, (ftnlen)6);
00359                 spptrf_(uplo, &n, &afac[1], &info);
00360 
00361 /*              Check error code from SPPTRF. */
00362 
00363                 if (info != izero) {
00364                     alaerh_(path, "SPPTRF", &info, &izero, uplo, &n, &n, &
00365                             c_n1, &c_n1, &c_n1, &imat, &nfail, &nerrs, nout);
00366                     goto L90;
00367                 }
00368 
00369 /*              Skip the tests if INFO is not 0. */
00370 
00371                 if (info != 0) {
00372                     goto L90;
00373                 }
00374 
00375 /* +    TEST 1 */
00376 /*              Reconstruct matrix from factors and compute residual. */
00377 
00378                 scopy_(&npp, &afac[1], &c__1, &ainv[1], &c__1);
00379                 sppt01_(uplo, &n, &a[1], &ainv[1], &rwork[1], result);
00380 
00381 /* +    TEST 2 */
00382 /*              Form the inverse and compute the residual. */
00383 
00384                 scopy_(&npp, &afac[1], &c__1, &ainv[1], &c__1);
00385                 s_copy(srnamc_1.srnamt, "SPPTRI", (ftnlen)32, (ftnlen)6);
00386                 spptri_(uplo, &n, &ainv[1], &info);
00387 
00388 /*              Check error code from SPPTRI. */
00389 
00390                 if (info != 0) {
00391                     alaerh_(path, "SPPTRI", &info, &c__0, uplo, &n, &n, &c_n1, 
00392                              &c_n1, &c_n1, &imat, &nfail, &nerrs, nout);
00393                 }
00394 
00395                 sppt03_(uplo, &n, &a[1], &ainv[1], &work[1], &lda, &rwork[1], 
00396                         &rcondc, &result[1]);
00397 
00398 /*              Print information about the tests that did not pass */
00399 /*              the threshold. */
00400 
00401                 for (k = 1; k <= 2; ++k) {
00402                     if (result[k - 1] >= *thresh) {
00403                         if (nfail == 0 && nerrs == 0) {
00404                             alahd_(nout, path);
00405                         }
00406                         io___34.ciunit = *nout;
00407                         s_wsfe(&io___34);
00408                         do_fio(&c__1, uplo, (ftnlen)1);
00409                         do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00410                         do_fio(&c__1, (char *)&imat, (ftnlen)sizeof(integer));
00411                         do_fio(&c__1, (char *)&k, (ftnlen)sizeof(integer));
00412                         do_fio(&c__1, (char *)&result[k - 1], (ftnlen)sizeof(
00413                                 real));
00414                         e_wsfe();
00415                         ++nfail;
00416                     }
00417 /* L60: */
00418                 }
00419                 nrun += 2;
00420 
00421                 i__3 = *nns;
00422                 for (irhs = 1; irhs <= i__3; ++irhs) {
00423                     nrhs = nsval[irhs];
00424 
00425 /* +    TEST 3 */
00426 /*              Solve and compute residual for  A * X = B. */
00427 
00428                     s_copy(srnamc_1.srnamt, "SLARHS", (ftnlen)32, (ftnlen)6);
00429                     slarhs_(path, xtype, uplo, " ", &n, &n, &kl, &ku, &nrhs, &
00430                             a[1], &lda, &xact[1], &lda, &b[1], &lda, iseed, &
00431                             info);
00432                     slacpy_("Full", &n, &nrhs, &b[1], &lda, &x[1], &lda);
00433 
00434                     s_copy(srnamc_1.srnamt, "SPPTRS", (ftnlen)32, (ftnlen)6);
00435                     spptrs_(uplo, &n, &nrhs, &afac[1], &x[1], &lda, &info);
00436 
00437 /*              Check error code from SPPTRS. */
00438 
00439                     if (info != 0) {
00440                         alaerh_(path, "SPPTRS", &info, &c__0, uplo, &n, &n, &
00441                                 c_n1, &c_n1, &nrhs, &imat, &nfail, &nerrs, 
00442                                 nout);
00443                     }
00444 
00445                     slacpy_("Full", &n, &nrhs, &b[1], &lda, &work[1], &lda);
00446                     sppt02_(uplo, &n, &nrhs, &a[1], &x[1], &lda, &work[1], &
00447                             lda, &rwork[1], &result[2]);
00448 
00449 /* +    TEST 4 */
00450 /*              Check solution from generated exact solution. */
00451 
00452                     sget04_(&n, &nrhs, &x[1], &lda, &xact[1], &lda, &rcondc, &
00453                             result[3]);
00454 
00455 /* +    TESTS 5, 6, and 7 */
00456 /*              Use iterative refinement to improve the solution. */
00457 
00458                     s_copy(srnamc_1.srnamt, "SPPRFS", (ftnlen)32, (ftnlen)6);
00459                     spprfs_(uplo, &n, &nrhs, &a[1], &afac[1], &b[1], &lda, &x[
00460                             1], &lda, &rwork[1], &rwork[nrhs + 1], &work[1], &
00461                             iwork[1], &info);
00462 
00463 /*              Check error code from SPPRFS. */
00464 
00465                     if (info != 0) {
00466                         alaerh_(path, "SPPRFS", &info, &c__0, uplo, &n, &n, &
00467                                 c_n1, &c_n1, &nrhs, &imat, &nfail, &nerrs, 
00468                                 nout);
00469                     }
00470 
00471                     sget04_(&n, &nrhs, &x[1], &lda, &xact[1], &lda, &rcondc, &
00472                             result[4]);
00473                     sppt05_(uplo, &n, &nrhs, &a[1], &b[1], &lda, &x[1], &lda, 
00474                             &xact[1], &lda, &rwork[1], &rwork[nrhs + 1], &
00475                             result[5]);
00476 
00477 /*                 Print information about the tests that did not pass */
00478 /*                 the threshold. */
00479 
00480                     for (k = 3; k <= 7; ++k) {
00481                         if (result[k - 1] >= *thresh) {
00482                             if (nfail == 0 && nerrs == 0) {
00483                                 alahd_(nout, path);
00484                             }
00485                             io___37.ciunit = *nout;
00486                             s_wsfe(&io___37);
00487                             do_fio(&c__1, uplo, (ftnlen)1);
00488                             do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer))
00489                                     ;
00490                             do_fio(&c__1, (char *)&nrhs, (ftnlen)sizeof(
00491                                     integer));
00492                             do_fio(&c__1, (char *)&imat, (ftnlen)sizeof(
00493                                     integer));
00494                             do_fio(&c__1, (char *)&k, (ftnlen)sizeof(integer))
00495                                     ;
00496                             do_fio(&c__1, (char *)&result[k - 1], (ftnlen)
00497                                     sizeof(real));
00498                             e_wsfe();
00499                             ++nfail;
00500                         }
00501 /* L70: */
00502                     }
00503                     nrun += 5;
00504 /* L80: */
00505                 }
00506 
00507 /* +    TEST 8 */
00508 /*              Get an estimate of RCOND = 1/CNDNUM. */
00509 
00510                 anorm = slansp_("1", uplo, &n, &a[1], &rwork[1]);
00511                 s_copy(srnamc_1.srnamt, "SPPCON", (ftnlen)32, (ftnlen)6);
00512                 sppcon_(uplo, &n, &afac[1], &anorm, &rcond, &work[1], &iwork[
00513                         1], &info);
00514 
00515 /*              Check error code from SPPCON. */
00516 
00517                 if (info != 0) {
00518                     alaerh_(path, "SPPCON", &info, &c__0, uplo, &n, &n, &c_n1, 
00519                              &c_n1, &c_n1, &imat, &nfail, &nerrs, nout);
00520                 }
00521 
00522                 result[7] = sget06_(&rcond, &rcondc);
00523 
00524 /*              Print the test ratio if greater than or equal to THRESH. */
00525 
00526                 if (result[7] >= *thresh) {
00527                     if (nfail == 0 && nerrs == 0) {
00528                         alahd_(nout, path);
00529                     }
00530                     io___39.ciunit = *nout;
00531                     s_wsfe(&io___39);
00532                     do_fio(&c__1, uplo, (ftnlen)1);
00533                     do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00534                     do_fio(&c__1, (char *)&imat, (ftnlen)sizeof(integer));
00535                     do_fio(&c__1, (char *)&c__8, (ftnlen)sizeof(integer));
00536                     do_fio(&c__1, (char *)&result[7], (ftnlen)sizeof(real));
00537                     e_wsfe();
00538                     ++nfail;
00539                 }
00540                 ++nrun;
00541 L90:
00542                 ;
00543             }
00544 L100:
00545             ;
00546         }
00547 /* L110: */
00548     }
00549 
00550 /*     Print a summary of the results. */
00551 
00552     alasum_(path, nout, &nfail, &nrun, &nerrs);
00553 
00554     return 0;
00555 
00556 /*     End of SCHKPP */
00557 
00558 } /* schkpp_ */

swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:59