/* dspdchk.f -- translated by f2c (version of 20 August 1993 13:15:44).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
#include "f2c.h"
/* Common Block Declarations */
struct {
doublereal a[40000], m[200];
} system_;
#define system_1 system_
struct {
integer n, lda;
} matdim_;
#define matdim_1 matdim_
struct {
char curpform[5];
} forms_;
#define forms_1 forms_
/* Table of constant values */
static integer c__3 = 3;
static integer c__1 = 1;
static integer c__9 = 9;
static integer c__4 = 4;
/* Subroutine */ int dspdchk_(x, ldx, b, x0, work, ldw, pform, matvec,
matvectrans, psolve, psolvetrans, psolveq, psolvetransq, backsolve,
tol, scaledtol, ltest, maxws, spdres, numtests, numsusp, criterr,
pform_len)
doublereal *x;
integer *ldx;
doublereal *b, *x0, *work;
integer *ldw;
char *pform;
/* Subroutine */ int (*matvec) (), (*matvectrans) (), (*psolve) (), (*
psolvetrans) (), (*psolveq) (), (*psolvetransq) (), (*backsolve) ();
doublereal *tol, *scaledtol;
logical *ltest;
integer *maxws;
logical *spdres;
integer *numtests, *numsusp, *criterr;
ftnlen pform_len;
{
/* Format strings */
static char fmt_81[] = "(\002WARNING: COULD NOT FORM 2-D POISSON (N=\002\
,i4,\002),INFO=\002,i2)";
static char fmt_82[] = "(\002WARNING: COULD NOT FORM 3-D POISSON (N=\002\
,i4,\002),INFO=\002,i2)";
static char fmt_83[] = "(\002WARNING: COULD NOT FORM ORDER\002,i4,\002 W\
ATHEN MATRIX\002)";
static char fmt_92[] = "(\002ERROR: RHS FOR TEST MATRIX\002,a4,\002 NOT \
FORMED\002)";
static char fmt_93[] = "(\002ERROR: INITIAL GUESS FOR TEST MATRIX\002,\
a4,\002 NOT FORMED\002)";
static char fmt_94[] = "(\002ERROR: PRECONDITIONER\002,a5,\002 NOT FOR\
MED\002)";
/* System generated locals */
integer x_dim1, x_offset, i__1, i__2;
/* Builtin functions */
integer s_rsle(), do_lio(), e_rsle(), s_wsle(), e_wsle(), s_wsfe(),
do_fio(), e_wsfe();
/* Subroutine */ int s_copy();
/* Local variables */
extern /* Subroutine */ int bicg_();
static logical flag_;
static integer ival, info[9];
extern /* Subroutine */ int bicgstab_();
static integer iter[9], nsys;
static logical tstbicgs, tstcheby, tstgmres;
static char initialguess[4];
static integer i, j, k;
extern /* Subroutine */ int cheby_();
static char aform[4];
static doublereal resid[9];
static integer iindx;
static doublereal anorm;
extern /* Subroutine */ int gmres_(), dcopy_();
static integer maxit;
static logical tstcg, tstjacobi;
extern doublereal negonefun_();
extern /* Subroutine */ int cg_(), comp2dense_(), jacobi_();
static integer nx, ny, nz;
extern /* Subroutine */ int vecgen_();
extern logical lsamen_();
static integer needws;
extern /* Subroutine */ int wathen_(), precon_(), gen57pt_();
static logical tstcgs;
extern /* Subroutine */ int result_();
static integer restrt;
static logical tstqmr, tstsor;
extern /* Subroutine */ int cgs_();
static char rhs[4];
extern /* Subroutine */ int qmr_(), sor_();
static logical tstbicg;
static integer matform;
extern doublereal matnorm_();
static integer npforms, ipointr;
extern doublereal zerofun_();
/* Fortran I/O blocks */
static cilist io___2 = { 0, 9, 0, 0, 0 };
static cilist io___4 = { 0, 6, 0, 0, 0 };
static cilist io___17 = { 0, 9, 0, 0, 0 };
static cilist io___28 = { 0, 6, 0, fmt_81, 0 };
static cilist io___29 = { 0, 10, 0, fmt_81, 0 };
static cilist io___30 = { 0, 6, 0, fmt_82, 0 };
static cilist io___31 = { 0, 10, 0, fmt_82, 0 };
static cilist io___33 = { 0, 6, 0, fmt_83, 0 };
static cilist io___34 = { 0, 10, 0, fmt_83, 0 };
static cilist io___35 = { 0, 6, 0, 0, 0 };
static cilist io___38 = { 0, 6, 0, fmt_92, 0 };
static cilist io___39 = { 0, 10, 0, fmt_92, 0 };
static cilist io___40 = { 0, 6, 0, fmt_93, 0 };
static cilist io___41 = { 0, 10, 0, fmt_93, 0 };
static cilist io___45 = { 0, 6, 0, fmt_94, 0 };
static cilist io___46 = { 0, 10, 0, fmt_94, 0 };
static cilist io___49 = { 0, 6, 0, 0, 0 };
static cilist io___50 = { 0, 10, 0, 0, 0 };
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. Common Blocks .. */
/* MAXDIM2 = MAXDIM*MAXDIM. */
/* Purpose */
/* ======= */
/* Subroutine to test the performance of the template kernels */
/* on symmetric positivie definite matrices. */
/* Generates, solves, and checks accuracy of linear systems. */
/* Algorithms tested: */
/* 1. CG */
/* 2. Cheby */
/* 3. SOR */
/* 4. BiCG */
/* 5. CGS */
/* 6. BiCGSTAB */
/* 7. GMRES */
/* 8. QMR */
/* 9. Jacobi ( for diagonally dominant Poisson matrices only ) */
/* Various systems are generated. Each method attempts to solve the */
/* system to the input TOL in MAXIT iterations. Each method iterates */
/* using various preconditioners. */
/* The result is compared with the normalized scaled residual */
/* || b - A*x || / ( ||A||||x||*N*TOL ). */
/* In order to do this, the solution vectors are stored in matrix */
/* X( LDX,* ). Column j contains the solution vector for algorithm j, */
/* j as defined above. */
/* ================================================================= */
/* .. Parameters .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. Local Arrays .. */
/* .. */
/* .. External Functions .. */
/* PDE Coefficient functions. */
/* .. */
/* .. Executable Statements .. */
/* Parameter adjustments */
--ltest;
pform -= 5;
--work;
--x0;
--b;
x_dim1 = *ldx;
x_offset = x_dim1 + 1;
x -= x_offset;
/* Function Body */
npforms = 2;
*numtests = 0;
*numsusp = 0;
*criterr = 0;
s_rsle(&io___2);
do_lio(&c__3, &c__1, (char *)&nsys, (ftnlen)sizeof(integer));
e_rsle();
/* Check for quick return. */
if (nsys < 0) {
s_wsle(&io___4);
do_lio(&c__9, &c__1, "ERROR IN NONSYMMETRIC TESTER: NUMBER OF SYSTEM\
S TO BE GENERATED IS LESS THAN 0", 78L);
e_wsle();
return 0;
}
flag_ = FALSE_;
for (i = 1; i <= 9; ++i) {
if (ltest[i]) {
flag_ = TRUE_;
}
/* L5: */
}
if (! flag_) {
return 0;
}
tstcg = ltest[1];
tstcheby = ltest[2];
tstsor = ltest[3];
tstbicg = ltest[4];
tstcgs = ltest[5];
tstbicgs = ltest[6];
tstgmres = ltest[7];
tstqmr = ltest[8];
tstjacobi = ltest[9];
/* L10: */
i__1 = nsys;
for (matform = 1; matform <= i__1; ++matform) {
s_rsle(&io___17);
do_lio(&c__9, &c__1, aform, 4L);
do_lio(&c__3, &c__1, (char *)&nx, (ftnlen)sizeof(integer));
do_lio(&c__3, &c__1, (char *)&ny, (ftnlen)sizeof(integer));
do_lio(&c__3, &c__1, (char *)&nz, (ftnlen)sizeof(integer));
do_lio(&c__9, &c__1, rhs, 4L);
do_lio(&c__9, &c__1, initialguess, 4L);
e_rsle();
/* The following two matrices are generated using a 5- or 7-poi
nt */
/* stencil using centered differences on a 1d, 2d, or 3d grid,
*/
/* with Dirichlet boundary conditions. */
/* The last 7 arguments to this routine are the coefficient */
/* functions for the PDE: */
/* delx( a delx u ) + dely ( b dely u) + delz ( c delz u ) +
*/
/* delx ( d u ) + dely (e u) + delz( f u ) + g u */
if (lsamen_(&c__4, aform, "F2SH", 4L, 4L)) {
/* -u_xx - u_yy */
matdim_1.n = nx * ny * nz;
ipointr = 1;
iindx = ipointr + matdim_1.n + 1;
ival = iindx + matdim_1.n * 5;
gen57pt_(&nx, &ny, &nz, &work[ival], &work[iindx], &work[ipointr],
negonefun_, negonefun_, zerofun_, zerofun_, zerofun_,
zerofun_, zerofun_);
comp2dense_(&work[ival], &work[ipointr], &work[iindx], &
matdim_1.n, system_1.a, &matdim_1.lda, "ROW", info, 3L);
if (info[0] != 0) {
s_wsfe(&io___28);
do_fio(&c__1, (char *)&matdim_1.n, (ftnlen)sizeof(integer));
do_fio(&c__1, (char *)&info[0], (ftnlen)sizeof(integer));
e_wsfe();
s_wsfe(&io___29);
do_fio(&c__1, (char *)&matdim_1.n, (ftnlen)sizeof(integer));
do_fio(&c__1, (char *)&info[0], (ftnlen)sizeof(integer));
e_wsfe();
goto L60;
}
} else if (lsamen_(&c__4, aform, "F3SH", 4L, 4L)) {
/* -u_xx - u_yy - u_zz */
matdim_1.n = nx * ny * nz;
ipointr = 1;
iindx = ipointr + matdim_1.n + 1;
ival = iindx + matdim_1.n * 10;
gen57pt_(&nx, &ny, &nz, &work[ival], &work[iindx], &work[ipointr],
negonefun_, negonefun_, negonefun_, zerofun_, zerofun_,
zerofun_, zerofun_);
comp2dense_(&work[ival], &work[ipointr], &work[iindx], &
matdim_1.n, system_1.a, &matdim_1.lda, "ROW", info, 3L);
if (info[0] != 0) {
s_wsfe(&io___30);
do_fio(&c__1, (char *)&matdim_1.n, (ftnlen)sizeof(integer));
do_fio(&c__1, (char *)&info[0], (ftnlen)sizeof(integer));
e_wsfe();
s_wsfe(&io___31);
do_fio(&c__1, (char *)&matdim_1.n, (ftnlen)sizeof(integer));
do_fio(&c__1, (char *)&info[0], (ftnlen)sizeof(integer));
e_wsfe();
goto L60;
}
} else if (lsamen_(&c__4, aform, "WATHEN", 4L, 6L)) {
/* Form Wathen matrix. */
/* ( Matrix order will be 3*NX*NY + 2*NX + 2*NY + 1 ) */
k = 0;
wathen_(&nx, &ny, &k, &matdim_1.n, system_1.a, &matdim_1.lda, &
work[1], ldw, info);
if (info[0] != 0) {
s_wsfe(&io___33);
i__2 = nx * 3 * ny + (nx << 1) + (ny << 1) + 1;
do_fio(&c__1, (char *)&i__2, (ftnlen)sizeof(integer));
do_fio(&c__1, (char *)&info[0], (ftnlen)sizeof(integer));
e_wsfe();
s_wsfe(&io___34);
i__2 = nx * 3 * ny + (nx << 1) + (ny << 1) + 1;
do_fio(&c__1, (char *)&i__2, (ftnlen)sizeof(integer));
do_fio(&c__1, (char *)&info[0], (ftnlen)sizeof(integer));
e_wsfe();
goto L60;
}
} else {
s_wsle(&io___35);
do_lio(&c__9, &c__1, aform, 4L);
do_lio(&c__9, &c__1, "IS AN UNKNOWM MATRIX TYPE", 25L);
e_wsle();
goto L60;
}
maxit = matdim_1.n << 2;
anorm = matnorm_(&matdim_1.n, system_1.a, &matdim_1.lda);
/* Form RHS and initial guess vectors. */
vecgen_(rhs, &matdim_1.n, system_1.a, &matdim_1.lda, &b[1], &info[1],
4L);
vecgen_(initialguess, &matdim_1.n, system_1.a, &matdim_1.lda, &x0[1],
&info[2], 4L);
if (info[1] != 0) {
s_wsfe(&io___38);
do_fio(&c__1, rhs, 4L);
e_wsfe();
s_wsfe(&io___39);
do_fio(&c__1, rhs, 4L);
e_wsfe();
goto L60;
} else if (info[2] != 0) {
s_wsfe(&io___40);
do_fio(&c__1, initialguess, 4L);
e_wsfe();
s_wsfe(&io___41);
do_fio(&c__1, initialguess, 4L);
e_wsfe();
goto L60;
}
/* L20: */
/* Solve system using the various algorithms, using no */
/* preconditioning, then diagonal preconditioning. */
i__2 = npforms;
for (i = 1; i <= i__2; ++i) {
for (j = 1; j <= 9; ++j) {
info[j - 1] = 0;
iter[j - 1] = maxit;
resid[j - 1] = *tol;
dcopy_(&matdim_1.n, &x0[1], &c__1, &x[j * x_dim1 + 1], &c__1);
/* L30: */
}
precon_(&matdim_1.n, system_1.a, &matdim_1.lda, pform + i * 5,
system_1.m, info, 5L);
if (info[0] != 0) {
s_wsfe(&io___45);
do_fio(&c__1, pform + i * 5, 5L);
e_wsfe();
s_wsfe(&io___46);
do_fio(&c__1, pform + i * 5, 5L);
e_wsfe();
goto L50;
}
s_copy(forms_1.curpform, pform + i * 5, 5L, 5L);
if (tstcg) {
cg_(&matdim_1.n, &b[1], &x[x_dim1 + 1], &work[1], ldw, iter,
resid, matvec, psolve, info);
}
if (tstcheby) {
cheby_(&matdim_1.n, &b[1], &x[(x_dim1 << 1) + 1], &work[1],
ldw, &iter[1], &resid[1], matvec, psolve, &info[1]);
}
if (tstsor) {
/* Set OMEGA */
if (matform == 1) {
/* Guass-Seidel */
work[1] = 1.;
} else {
work[1] = 1.2;
}
sor_(&matdim_1.n, &b[1], &x[x_dim1 * 3 + 1], &work[1], ldw, &
iter[2], &resid[2], matvec, backsolve, &info[2]);
}
if (tstbicg) {
bicg_(&matdim_1.n, &b[1], &x[(x_dim1 << 2) + 1], &work[1],
ldw, &iter[3], &resid[3], matvec, matvectrans, psolve,
psolvetrans, &info[3]);
}
if (tstcgs) {
cgs_(&matdim_1.n, &b[1], &x[x_dim1 * 5 + 1], &work[1], ldw, &
iter[4], &resid[4], matvec, psolve, &info[4]);
}
if (tstbicgs) {
bicgstab_(&matdim_1.n, &b[1], &x[x_dim1 * 6 + 1], &work[1],
ldw, &iter[5], &resid[5], matvec, psolve, &info[5]);
}
if (tstgmres) {
/* For the symmetric case, restarts = N. */
restrt = matdim_1.n;
needws = matdim_1.n * (restrt + 4) + (restrt + 1) * (restrt +
2);
if (needws > *maxws) {
info[6] = 100;
s_wsle(&io___49);
do_lio(&c__9, &c__1, "WARNING: NOT ENOUGH WORKSPACE FOR \
GMRES, (TEST MATRIX", 53L);
do_lio(&c__3, &c__1, (char *)&matform, (ftnlen)sizeof(
integer));
do_lio(&c__9, &c__1, ")", 1L);
e_wsle();
s_wsle(&io___50);
do_lio(&c__9, &c__1, "NOT ENOUGH WORKSPACE FOR GMRES, (T\
EST MATRIX", 44L);
do_lio(&c__3, &c__1, (char *)&matform, (ftnlen)sizeof(
integer));
do_lio(&c__9, &c__1, " REQUIRES MAXLEN=", 17L);
do_lio(&c__3, &c__1, (char *)&needws, (ftnlen)sizeof(
integer));
do_lio(&c__9, &c__1, ")", 1L);
e_wsle();
} else {
gmres_(&matdim_1.n, &b[1], &x[x_dim1 * 7 + 1], &restrt, &
work[1], ldw, &work[matdim_1.n * (restrt + 5) + 1]
, ldw, &iter[6], &resid[6], matvec, psolve, &info[
6]);
}
}
if (tstqmr) {
qmr_(&matdim_1.n, &b[1], &x[(x_dim1 << 3) + 1], &work[1], ldw,
&iter[7], &resid[7], matvec, matvectrans, psolveq,
psolvetransq, &info[7]);
}
if (tstjacobi) {
if (! lsamen_(&c__3, aform, "WATH", 4L, 4L) && i == 1) {
/* Since preconditioning does not apply to
Jacobi, it is */
/* only called once per test matrix. (Wath
en matrix is */
/* not diagonally dominant.) */
jacobi_(&matdim_1.n, &b[1], &x[x_dim1 * 9 + 1], &work[1],
ldw, &iter[8], &resid[8], matvec, &info[8]);
if (info[8] != 0) {
++(*numsusp);
}
} else {
/* Flag not to check accuracy. */
info[8] = 100;
}
}
/* Check for convergence. */
for (j = 1; j <= 8; ++j) {
if (info[j - 1] != 0) {
++(*numsusp);
}
/* L40: */
}
/* Write results to file. */
result_(&matdim_1.n, system_1.a, &matdim_1.lda, &x[x_offset], ldx,
&b[1], &work[1], "SPD", pform + i * 5, iter, resid, tol,
info, aform, &anorm, <est[1], scaledtol, spdres,
criterr, 3L, 5L, 4L);
++(*numtests);
L50:
;
}
L60:
;
}
/* L200: */
return 0;
/* -- End of DSPDCHK */
} /* dspdchk_ */