/* pde44_sparse.c  discretize, build sparse linear equations, solve
 * solve uxxxx(x,y,z,t) + 2 uyyyy(x,y,z,t) + 3 uzzzz(x,y,z,t) +
 *     4 utttt(x,y,z,t) + 5  uxyt(x,y,z,t) + 6 uyyzz(x,y,z,t) +
 *     7  uztt(x,y,z,t) + 8  uxxx(x,y,z,t) + 9  uyyt(x,y,z,t) +
 *     10  uzt(x,y,z,t) + 11   ut(x,y,z,t) + 12    u(x,y,z,t) =
 *     F(x,y,z,t)
 *     F(x,y,z,t) = 706 + 120 t^2 + 140 y^2z + 768 x +
 *                  180 z^2*t + 200 y^2z*t + 44 t^3 +
 *                  110 y^2z^2t + 66 xy + 12t^4 +
 *                  24 z^4 + 36 y^4 + 48 x^4 +
 *                  60 y^2z^2t^2 + 72 xyt + 84 xz + 96 t 
 *
 * boundary conditions computed using u(x,y,z,t)
 * analytic solution is u(x,y,z,t) = t^4 + 2 z^4 + 3 y^4 + 4 x^4 +
 *                                   5 y^2 z^2 t^2 + 6 x y t + 
 *                                   7 x z + 8 t + 9
 *
 * replace continuous derivatives with discrete derivatives
 * solve for Uijkl numerical approximation U(x_i,y_j,z_k,t_l)
 * A * U = F, solve simultaneous equations for U
 */

/* fourth order numerical difference order=4, npoints=5, term=0
 * d^4u/dx^4 = (1/(b4[0]*hx^4))*(a4[0][0]*u(x,y,z,t) + a4[0][1]*u(x+hx,y,z,t) +
 * a4[0][2]*u(x+2hx,y,z,t) + a4[0][3]*u(x+3hx,y,z,t) +a4[0][4]*u(x+4hx,y,z,t))
 *
 * d^4u/dx^4 using subscripts i,j,k,l for i a minimum subscript of x
 * d^4u/dx^4 =(1/(b4[0]*hx^4))*(a4[0][0]*u(i,j,k,l) + a4[0][1]*u(i+1,j,k,l) +
 * a4[0][2]*u(i+2,j,k,l) + a4[0][3]*u(i+3,j,k,l) +a4[0][4]*u(i+4,j,k,l))
 *
 * d^4u/dy^4 using subscripts i,j,k,l for j a minimum subscript of y, etc 
 * d^4u/dy^4 =(1/(b4[0]*hy^4))*(a4[0][0]*u(i,j,k,l) + a4[0][1]*u(i,j+1,k,l) +
 * a4[0][2]*u(i,j+2,k,l) + a4[0][3]*u(i,j+3,k,l) +a4[0][4]*u(i,j+4,k,l))
 *
 * d^4u/dx^4 using subscripts i,j,k,l for i-2 and i+2 allowed subscripts of x
 * d^4u/dx^4 =(1/(b4[2]*hx^4))*(a4[2][0]*u(i-2,j,k,l) + a4[2][1]*u(i-1,j,k,l) +
 * a4[2][2]*u(i,j,k,l) + a4[2][3]*u(i+1,j,k,l) +a4[2][4]*u(i+2,j,k,l))
 *
 * d^4u/dx^4 using subscripts i,j,k,l for i a maximum subscript of x
 * d^4u/dx^4 =(1/(b4[4]*hx^4))*(a4[4][0]*u(i-4,j,k,l) + a4[4][1]*u(i-3,j,k,l) +
 * a4[4][2]*u(i-2,j,k,l) + a4[4][3]*u(i-1,j,k,l) +a4[4][4]*u(i,j,k,l))
 *
 * subscripts in code use x[i], y[ii], z[iii], t[iiii]
 * derivatives are computed with b,hx,a included as in
 * cx[j] for first derivative of x at point j
 * cxxxx[j] for fourth derivative of x at point j
 * cyy[jj]*czz[jjj] for d^4u/dy^2dz^2 at jj, jjj
 */

/* The subscript versions of the derivatives are substituted into the
 * equation to be solved. The resulting equation is analytically solved
 * for u(i,j,k,l) = other u terms with coefficients and + f(xi,yj,zk,tl)
 * The boundaries xmin..xmax, ymin..ymax, zmin..zmax, tmin..tmax are known
 * The number of grid points in each dimension is known nx, ny, nz, nt
 * hx=(xmax-xmin)/(nx-1), hy=(ymax-ymin)/(ny-1),
 * hz=(zmax-zmin)/(nz-1), ht=(tmax-tmin)/(tx-1), 
 * thus xi=xmin+i*hx, yj=ymin+j*hy, zk=zmin+k*hz, tl=tmin+l*ht
 * then the linear system of equations is built,
 * nxyzt=(nx-2)*(ny-2)*(nz-2)*(nt-2) equations because the boundary value
 *                                   is known on each end.
 * The matrix A is nxyzt rows by nxyzt columns (nxyzt+1) columns including F
 * The forcing function vector F is nxyzt elements adjunct to A
 * The solution vector U is nxyzt elements.
 * The building of the A matix requires 4 nested loops for rows
 * i  in 1..nx-1, j  in 1..ny-1,  k  in 1..nz-1, l  in 1..nt-1 for rows of A
 * then each term in the differential equation fills in that row
 * and cs = (nx-2)*(ny-2)*(nz-2)*(nt-2) the RHS, fills the last column of A
 * subscripting functions are used to compute linear subscripts of A, U
 * row=(i-1)*(ny-2)*(nz-2)*(nt-2)+(ii-1)*(nz-2)*(nt-2)+(iii-1)*(nt-2)+(iiii-1)
 * row*(nxyzt+1)+col for A
 *
 * care must be taken to move the negative of boundary elements to
 * the right hand size of the equation to be solved. These are subtracted
 * from the computed value of the forcing function, cs, for that row.
 * tests are j==0 || j==nx-1   ...  jjjj==0 || jjjj==nt-1 for boundaries
 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include "deriv.h"
#include "sparse.h"

#undef  abs
#define abs(x) ((x)<0.0?(-(x)):(x))
#undef  min
#define min(a,b) ((a)<(b)?(a):(b))
#undef  max
#define max(a,b) ((a)>(b)?(a):(b))

#define nx  6
#define ny  6
#define nz  6
#define nt  6
#define nxyzt (nx-2)*(ny-2)*(nz-2)*(nt-2)
#define nn max(nx, max(ny, max(nz, nt)))

int s(int i, int ii, int iii, int iiii) /* for x,y,z,t */
{
  return (i-1)*(ny-2)*(nz-2)*(nt-2) + (ii-1)*(nz-2)*(nt-2) +
         (iii-1)*(nt-2) + (iiii-1);
} /* end s */

int sk(int k, int i, int ii, int iii, int iiii) /* for x,y,z,t */
{
  return k*(nxyzt+1) + (i-1)*(ny-2)*(nz-2)*(nt-2) + (ii-1)*(nz-2)*(nt-2) +
         (iii-1)*(nt-2) + (iiii-1);

} /* end sk */

int ss(int i, int ii, int iii, int iiii, int j, int jj, int jjj, int jjjj)
{
  return s(i,ii,iii,iiii)*(nxyzt+1) + s(j,jj,jjj,jjjj);
} /* end ss */

/* was: static double A[nxyzt*(nxyzt+1)]; */
static struct sparse A; /* linear equations including RHS */
static double U[nxyzt];
static double xg[nx];
static double yg[ny];
static double zg[nz];
static double tg[nt];
static double Ua[nxyzt];
static double xmin = 0.0; /* problem parameters */
static double xmax = 1.0;
static double ymin = 0.0;
static double ymax = 1.0;
static double zmin = 0.0;
static double zmax = 1.0;
static double tmin = 0.0;
static double tmax = 1.0;
static double hx, hy, hz, ht;


static double f(double x, double y, double z, double t)
{
  return 706.0 + 120.0*t*t + 140.0*y*y*z + 768.0*x + 180.0*z*z*t +
         200.0*y*y*z*t + 44.0*t*t*t + 110.0*y*y*z*z*t + 66.0*x*y +
         12.0*t*t*t*t + 24.0*z*z*z*z + 36.0*y*y*y*y + 48.0*x*x*x*x +
         60.0*y*y*z*z*t*t + 72.0*x*y*t + 84.0*x*z + 96.0*t; 
}

static double uana(double x, double y, double z, double t)
                                        /*analytic solution and boundaries*/
{
  return t*t*t*t + 2.0*z*z*z*z + 3.0*y*y*y*y + 4.0*x*x*x*x +
         5.0*y*y*z*z*t*t + 6.0*x*y*t + 7.0*x*z + 8.0*t + 9.0;
}

double eval_derivative(int xord, int yord, int zord, int tord,
                       int i, int ii, int iii, int iiii, double u[])
{
  int j, jj, jjj, jjjj;
  double cx[nx];
  double cy[ny];
  double cz[nz];
  double ct[nt];
  double p1[nn];
  double p2[nn];
  double p3[nn];
  double discrete;
  double x, y, z, t;

  rderiv(xord, nx, i, hx, cx);
  x = xg[i];
  rderiv(yord, ny, ii, hy, cy);
  y = yg[ii];
  rderiv(zord, nz, iii, hz, cz);
  z = zg[iii];
  rderiv(tord, nt, iiii, ht, ct);
  t = tg[iiii];
  discrete = 0.0;

  /* four cases of one partial x, y, z, t to any order */
  if(xord!=0 && yord==0 && zord==0 && tord==0)
    for(j=0; j<nx; j++)
      if(j==0 || j==nx-1) discrete += cx[j]*uana(xg[j],y,z,t);
      else discrete += cx[j]*u[s(j,ii,iii,iiii)];
  else if(xord==0 && yord!=0 && zord==0 && tord==0)
    for(jj=0; jj<ny; jj++)
      if(jj==0 || jj==ny-1) discrete += cy[jj]*uana(x,yg[jj],z,t);
      else discrete += cy[jj]*u[s(i,jj,iii,iiii)];
  else if(xord==0 && yord==0 && zord!=0 && tord==0)
    for(jjj=0; jjj<nz; jjj++)
      if(jjj==0 || jjj==nz-1) discrete += cz[jjj]*uana(x,y,zg[jjj],t);
      else discrete += cz[jjj]*u[s(i,ii,jjj,iiii)];
  else if(xord==0 && yord==0 && zord==0 && tord!=0)
    for(jjjj=0; jjjj<nt; jjjj++)
      if(jjjj==0 || jjjj==nt-1) discrete += ct[jjjj]*uana(x,y,z,tg[jjjj]);
      else discrete += ct[jjjj]*u[s(i,ii,iii,jjjj)];
  /* six cases of two partials xy, xz, xt, yz, yt, zt to any order */
  else if(xord!=0 && yord!=0 && zord==0 && tord==0)
  {
    for(j=0; j<nx; j++)
    {
      p1[j] = 0.0;
      for(jj=0; jj<ny; jj++)
        if(j==0 || j==nx-1 || jj==0 || jj==ny-1)
             p1[j] += cy[jj]*uana(xg[j],yg[jj],zg[iii],tg[iiii]);
        else p1[j] += cy[jj]*u[s(j,jj,iii,iiii)];
      discrete += cx[j]*p1[j];
    }
  }
  else if(xord!=0 && yord==0 && zord!=0 && tord==0)
  {
    for(j=0; j<nx; j++)
    {
      p1[j] = 0.0;
      for(jjj=0; jjj<nz; jjj++)
        if(j==0 || j==nx-1 || jjj==0 || jjj==nz-1)
	     p1[j] += ct[jjj]*uana(xg[j],yg[ii],zg[jjj],tg[iiii]);
        else p1[j] += cz[jjj]*u[s(j,ii,jjj,iiii)];
      discrete += cx[j]*p1[j];
    }
  }
  else if(xord!=0 && yord==0 && zord==0 && tord!=0)
  {
    for(j=0; j<nx; j++)
    {
      p1[j] = 0.0;
      for(jjjj=0; jjjj<nt; jjjj++)
        if(j==0 || j==nx-1 || jjjj==0 || jjjj==nt-1)
             p1[j] += ct[jjjj]*uana(xg[j],yg[ii],zg[iii],tg[jjjj]);
        else p1[j] += ct[jjjj]*u[s(j,ii,iii,jjjj)];
      discrete += cx[j]*p1[j];
    }
  }
  else if(xord==0 && yord!=0 && zord!=0 && tord==0)
  {
    for(jj=0; jj<ny; jj++)
    {
      p1[jj] = 0.0;
      for(jjj=0; jjj<nz; jjj++)
        if(jj==0 || jj==ny-1 || jjj==0 || jjj==nz-1)
             p1[jj] += cz[jjj]*uana(xg[i],yg[jj],zg[jjj],tg[iiii]);
        else p1[jj] += cz[jjj]*u[s(i,jj,jjj,iiii)];
      discrete += cy[jj]*p1[jj];
    }
  }
  else if(xord==0 && yord!=0 && zord==0 && tord!=0)
  {
    for(jj=0; jj<ny; jj++)
    {
      p1[jj] = 0.0;
      for(jjjj=0; jjjj<nt; jjjj++)
        if(jj==0 || jj==ny-1 || jjjj==0 || jjjj==nt-1)
             p1[jj] += ct[jjjj]*uana(xg[i],yg[jj],zg[iii],tg[jjjj]);
        else p1[jj] += ct[jjjj]*u[s(i,jj,iii,jjjj)];
      discrete += cy[jj]*p1[jj];
    }
  }
  else if(xord==0 && yord==0 && zord!=0 && tord!=0)
  {
    for(jjj=0; jjj<nz; jjj++)
    {
      p1[jjj] = 0.0;
      for(jjjj=0; jjjj<nt; jjjj++)
        if(jjj==0 || jjj==nz-1 || jjjj==0 || jjjj==nt-1)
             p1[jjj] += ct[jjjj]*uana(xg[i],yg[ii],zg[jjj],tg[jjjj]);
        else p1[jjj] += ct[jjjj]*u[s(i,ii,jjj,jjjj)];
      discrete += cz[jjj]*p1[jjj];
    }
  }
  /* four cases of three partials xyz, xyt, xzt, yzt to any order */
  else if(xord!=0 && yord!=0 && zord!=0 && tord==0)
  {
    for(j=0; j<nx; j++)
    {
      p1[j] = 0.0;
      for(jj=0; jj<ny; jj++)
      {
        p2[jj] = 0.0;
        for(jjj=0; jjj<nz; jjj++)
          if(j==0 || j==nx-1 ||
             jj==0 || jj==ny-1 ||
             jjj==0 || jjj==nz-1)
               p2[jj] += cz[jjj]*uana(xg[j],yg[jj],zg[jjj],tg[iiii]);
          else p2[jj] += cz[jjj]*u[s(j,jj,jjj,iiii)];
        p1[j] += cy[jj]*p2[jj];
      }
      discrete += cx[j]*p1[j];
    }
  }
  else if(xord!=0 && yord!=0 && zord==0 && tord!=0)
  {
    for(j=0; j<nx; j++)
    {
      p1[j] = 0.0;
      for(jj=0; jj<ny; jj++)
      {
        p2[jj] = 0.0;
        for(jjjj=0; jjjj<nt; jjjj++)
          if(j==0 || j==nx-1 ||
             jj==0 || jj==ny-1 ||
             jjjj==0 || jjjj==nt-1)
               p2[jj] += ct[jjjj]*uana(xg[j],yg[jj],zg[iii],tg[jjjj]);
          else p2[jj] += ct[jjjj]*u[s(j,jj,iii,jjjj)];
        p1[j] += cy[jj]*p2[jj];
      }
      discrete += cx[j]*p1[j];
    }
  }
  else if(xord!=0 && yord==0 && zord!=0 && tord!=0)
  {
    for(j=0; j<nx; j++)
    {
      p1[j] = 0.0;
      for(jjj=0; jjj<nz; jjj++)
      {
        p2[jjj] = 0.0;
        for(jjjj=0; jjjj<nt; jjjj++)
          if(j==0 || j==nx-1 ||
             jjj==0 || jjj==nz-1 ||
             jjjj==0 || jjjj==nt-1)
               p2[jjj] += ct[jjjj]*uana(xg[j],yg[ii],zg[jjj],tg[jjjj]);
          else p2[jjj] += ct[jjjj]*u[s(j,ii,jjj,jjjj)];
        p1[j] += cz[jjj]*p2[jjj];
      }
      discrete += cx[j]*p1[j];
    }
  }
  else if(xord==0 && yord!=0 && zord!=0 && tord!=0)
  {
    for(jj=0; jj<ny; jj++)
    {
      p1[jj] = 0.0;
      for(jjj=0; jjj<nz; jjj++)
      {
        p2[jjj] = 0.0;
        for(jjjj=0; jjjj<nt; jjjj++)
          if(jj==0 || jj==ny-1 ||
             jjj==0 || jjj==nz-1 ||
             jjjj==0 || jjjj==nt-1)
               p2[jjj] += ct[jjjj]*uana(xg[i],yg[jj],zg[jjj],tg[jjjj]);
          else p2[jjj] += ct[jjjj]*u[s(i,jj,jjj,jjjj)];
        p1[jj] += cz[jjj]*p2[jjj];
      }
      discrete += cy[jj]*p1[jj];
    }
  }
  /* one case of four partials */
  else
  {
    for(j=0; j<nx; j++)
    {
      p1[j] = 0.0;
      for(jj=0; jj<ny; jj++)
      {
        p2[jj] = 0.0;
        for(jjj=0; jjj<nz; jjj++)
        {
          p3[jjj] = 0.0;
          for(jjjj=0; jjjj<nt; jjjj++)
	    if(j==0 || j==nx-1 ||
               jj==0 || jj==ny-1 ||
               jjj==0 || jjj==nz-1 ||
               jjjj==0 || jjjj==nt-1)
               p3[jjj] += ct[jjjj]*uana(xg[j],yg[jj],zg[jjj],tg[jjjj]);
	    else p3[jjj] += ct[jjjj]*u[s(j,jj,jjj,jjjj)];
          p2[jj] += cz[jjj]*p3[jjj];
        } /* end jjj */
        p1[j] += cy[jj]*p2[jj];
      } /* end jj */
      discrete += cx[j]*p1[j];
    } /* end j */
  } /* end if */
  return discrete;
} /* end eval_deriv */

void check_soln(double u[])
{
  /* check uxxxx(x,y,z,t) + 2 uyyyy(x,y,z,t) + 3 uzzzz(x,y,z,t) +
         4 utttt(x,y,z,t) + 5  uxyt(x,y,z,t) + 6 uyyzz(x,y,z,t) +
         7  uztt(x,y,z,t) + 8  uxxx(x,y,z,t) + 9  uyyt(x,y,z,t) +
         10  uzt(x,y,z,t) + 11   ut(x,y,z,t) + 12    u(x,y,z,t) = F(x,y,z,t)
  */
  double val, err, smaxerr;
  int i, ii, iii, iiii;

  smaxerr = 0.0;
  for(i=1; i<nx-1; i++)
  {
    for(ii=1; ii<ny-1; ii++)
    {
      for(iii=1; iii<nz-1; iii++)
      {
        for(iiii=1; iiii<nt-1; iiii++)
        {
          val = 0.0;
          /* uxxxx(x,y,z,t) */
          val += eval_derivative(4, 0, 0, 0, i, ii, iii, iiii, u);
	  /* + 2 uyyyy(x,y,z,t) */
          val += 2.0*eval_derivative(0, 4, 0, 0, i, ii, iii, iiii, u);
	  /* + 3 uzzzz(x,y,z,t) */
          val += 3.0*eval_derivative(0, 0, 4, 0, i, ii, iii, iiii, u);
          /* + 4 utttt(x,y,z,t) */
          val += 4.0*eval_derivative(0, 0, 0, 4, i, ii, iii, iiii, u);
	  /* + 5  uxyt(x,y,z,t) */
          val += 5.0*eval_derivative(1, 1, 0, 1, i, ii, iii, iiii, u);
	  /* + 6 uyyzz(x,y,z,t) */
          val += 6.0*eval_derivative(0, 2, 2, 0, i, ii, iii, iiii, u);
          /* + 7  uztt(x,y,z,t) */
          val += 7.0*eval_derivative(0, 0, 1, 2, i, ii, iii, iiii, u);
	  /* + 8  uxxx(x,y,z,t) */
          val += 8.0*eval_derivative(3, 0, 0, 0, i, ii, iii, iiii, u);
	  /* + 9  uyyt(x,y,z,t) */
          val += 9.0*eval_derivative(0, 2, 0, 1, i, ii, iii, iiii, u);
          /* +10   uzt(x,y,z,t) */
          val += 10.0*eval_derivative(0, 0, 1, 1, i, ii, iii, iiii, u);
	  /* +11    ut(x,y,z,t) */
          val += 11.0*eval_derivative(0, 0, 0, 1, i, ii, iii, iiii, u);
	  /* +12     u(x,y,z,t) */
          val += 12.0*uana(xg[i],yg[ii],zg[iii],tg[iiii]);
          /* -       F(x,y,z,t) should be zero */
	  err = abs(val - f(xg[i],yg[ii],zg[iii],tg[iiii]));
          if(err>smaxerr) smaxerr = err;
        } /* end iiii */
      } /* end iii */
    } /* end ii */
  } /* end i */
  printf("check_soln maxerr=%g \n",smaxerr);
} /* end check_soln */
 

int main(int argc, char *argv[])
{
  double x, y, z, t;
  int i, ii, iii, iiii, j, jj, jjj, jjjj;
  int k, cs;
  /* coeficients of terms in first equation */
  double coefdxxxx = 1.0; /* computed below if not constant */
  double coefdyyyy = 2.0;
  double coefdzzzz = 3.0;
  double coefdtttt = 4.0;
  double coefdxyt  = 5.0;
  double coefdyyzz = 6.0;
  double coefdztt  = 7.0;
  double coefdxxx  = 8.0;
  double coefdyyt  = 9.0;
  double coefdzt   = 10.0;
  double coefdt    = 11.0;
  double coefu     = 12.0;
  /* derivative point arrays, dynamic in this version */
  double cx[nx], cy[ny], cz[nz], ct[nt]; 
  double cxx[nx], cyy[ny], czz[nz], ctt[nt]; 
  double cxxx[nx], cyyy[ny], czzz[nz], cttt[nt]; 
  double cxxxx[nx], cyyyy[ny], czzzz[nz], ctttt[nt];
  double val, err, avgerr, maxerr;
  double time_start;
  double time_now;
  double time_last;

  printf("pde44_sparse.c running \n");
  printf("Given uxxxx(x,y,z,t) + 2 uyyyy(x,y,z,t) + 3 uzzzz(x,y,z,t) + \n");
  printf("    4 utttt(x,y,z,t) + 5  uxyt(x,y,z,t) + 6 uyyzz(x,y,z,t) + \n");
  printf("    7  uztt(x,y,z,t) + 8  uxxx(x,y,z,t) + 9  uyyt(x,y,z,t) + \n");
  printf("    10  uzt(x,y,z,t) + 11   ut(x,y,z,t) + 12    u(x,y,z,t) = \n");
  printf("  F(x,y,z,t) \n");
  printf("  F(x,y,z,t) = 706 + 120 t^2 + 140 y^2z + 768 x +   \n");
  printf("               180 z^2*t + 200 y^2z*t + 44 t^3 +    \n");
  printf("               110 y^2z^2t + 66 xy + 12t^4 +        \n");
  printf("               24 z^4 + 36 y^4 + 48 x^4 +           \n");
  printf("               60 y^2z^2t^2 + 72 xyt + 84 xz + 96 t \n"); 
  printf("xmin<=x<=xmax ymin<=y<=ymax zmin<=z<=zmax Boundaries \n");
  printf("Analytic solution u(x,y,z,t) =  t^4 + 2 z^4 + 3 y^4 + \n");
  printf("     4 x^4 + 5 y^2 z^2 t^2 + 6 x y t + 7 x z + 8 t + 9 \n");
  printf("\n");

  printf("xmin=%g, xmax=%g, ymin=%g, ymax=%g \n", xmin, xmax, ymin, ymax);
  printf("zmin=%g, zmax=%g, tmin=%g, tmax=%g \n", zmin, zmax, tmin, tmax);
  printf("nx=%d, ny=%d, nz=%d, nt=%d \n", nx, ny, nz, nt);
  time_start = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("time start = %f seconds \n", time_start);
  time_last = time_start;

  printf("x grid and analytic solution at ymin,zmin,tmin \n");
  hx = (xmax-xmin)/(double)(nx-1);
  for(i=0; i<nx; i++)
  {
    xg[i] = xmin+(double)i*hx;
    Ua[i] = uana(xg[i],ymin,zmin,tmin); /* temp */
    printf("i=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
            i, xg[i], ymin, zmin, tmin, Ua[i]);
  } /* i */  
  printf("\n");
  printf("y grid and analytic solution at xmin,zmin,tmin \n");
  hy = (ymax-ymin)/(double)(ny-1);
  for(ii=0; ii<ny; ii++)
  {
    yg[ii] = ymin+(double)ii*hy;
    Ua[ii] = uana(xmin, yg[ii],zmin,tmin); /* temp */
    printf("ii=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
            ii, xmin, yg[ii], zmin, tmin, Ua[ii]);
  } /* ii */
  printf("\n");
  printf("z grid and analytic solution at xmin,ymin,tmin \n");
  hz = (zmax-zmin)/(double)(nz-1);
  for(iii=0; iii<nz; iii++)
  {
    zg[iii] = zmin+(double)iii*hz;
    Ua[iii] = uana(xmin, ymin, zg[iii],tmin); /* temp */
    printf("iii=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
            iii, xmin, ymin, zg[iii], tmin, Ua[iii]);
  } /* iii */
  printf("\n");
  printf("t grid and analytic solution at xmin,ymin,zmin \n");
  ht = (tmax-tmin)/(double)(nt-1);
  for(iiii=0; iiii<nt; iiii++)
  {
    tg[iiii] = tmin+(double)iiii*ht;
    Ua[iiii] = uana(xmin, ymin, zmin, tg[iiii]); /* temp */
    printf("iiii=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
            iiii, xmin, ymin, zmin, tg[iiii], Ua[iiii]);
  } /* iiii */
  printf("\n");

  printf("check rderiv(4,nx,0,hx,cxxxx \n");
  rderiv(4,nx,0,hx,cxxxx);
  for(i=0; i<nx; i++) printf("cxxxx[%d]=%f \n", i, cxxxx[i]);
  printf("check rderiv(4,nx,0,hx,cxxxx \n");
  rderiv(4,nx,nx-1,hx,cxxxx);
  for(i=0; i<nx; i++) printf("cxxxx[%d]=%f \n", i, cxxxx[i]);
  printf("\n");

  printf("put solution in Ua vector \n");
  /* put solution in Ua vector */
  for(i=1; i<nx-1; i++)
  {
    x = xg[i];
    for(ii=1; ii<ny-1; ii++)
    {
      y = yg[ii];
      for(iii=1; iii<nz-1; iii++)
      {
        z = zg[iii];
        for(iiii=1; iiii<nt-1; iiii++)
        {
          t = tg[iiii];
          Ua[s(i,ii,iii,iiii)] = uana(x,y,z,t);
          /* printf("solution at i=%d,x=%g, ii=%d,y=%g, iii=%d,z=%g, iiii=%d,t=%g is %g \n",
	     i, x, ii, y, iii, z, iiii, t, Ua[s(i,ii,iii,iiii)]); */
        } /* iiii */
      } /* iii */
    } /* ii */
  } /* i */
  time_now = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("time now = %f seconds, time for previous section = %f seconds \n",
         time_now, time_now-time_last);
  time_last = time_now;
  printf("\n");

  /* initialize A */

  printf("initialize A \n");
  sparseInit(&A, nxyzt);

  time_now = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("time now = %f seconds, time for previous section = %f seconds \n",
         time_now, time_now-time_last);
  time_last = time_now;
  printf("\n");

  /* compute entries in A matrix, inside boundary */
  printf("compute A matrix \n");
  cs = (nx-2)*(ny-2)*(nz-2)*(nt-2); /* constant RHS column */
  val = 0.0;
  for(i=1; i<nx-1; i++)
  {
    for(ii=1; ii<ny-1; ii++)
    {
      for(iii=1; iii<nz-1; iii++)
      {
        for(iiii=1; iiii<nt-1; iiii++)
        {
          k = s(i,ii,iii,iiii); /* row index */

          /* for each term in first equation */

	  /* for d^4u/dx^4 */
          coefdxxxx = 1.0;
          rderiv(4,nx,i,hx,cxxxx);
          for(j=0; j<nx; j++)
          {
            val = coefdxxxx*cxxxx[j];
            if(j==0 || j==nx-1)
	     /* A[k*(nxyzt+1)+cs] -= val*uana(xg[j],yg[ii],zg[iii],tg[iiii]);*/
	      sparseAdd(A,k,cs,-val*uana(xg[j],yg[ii],zg[iii],tg[iiii]));
	    else /* A[sk(k,j,ii,iii,iiii)] += val; */
	      sparseAdd(A,k,s(j,ii,iii,iiii),val);
	  }
	  /* for d^4u/dy^4 */
          coefdyyyy = 2.0;
          rderiv(4,ny,ii,hy,cyyyy);
          for(jj=0; jj<ny; jj++)
          {
            val = coefdyyyy*cyyyy[jj];
            if(jj==0 || jj==ny-1)
	      sparseAdd(A,k,cs,-val*uana(xg[i],yg[jj],zg[iii],tg[iiii]));
            else sparseAdd(A,k,s(i,jj,iii,iiii),val);
	  }

	  /* for d^4u/dz^4 */
          coefdzzzz = 3.0;
          rderiv(4,nz,iii,hz,czzzz);
          for(jjj=0; jjj<nz; jjj++)
          {
            val = coefdzzzz*czzzz[jjj];
            if(jjj==0 || jjj==nz-1)
	      sparseAdd(A,k,cs,-val*uana(xg[i],yg[ii],zg[jjj],tg[iiii]));
            else sparseAdd(A,k,s(i,ii,jjj,iiii),val);
	  }

	  /* for d^4u/dt^4 */
          coefdtttt = 4.0;
          rderiv(4,nt,iiii,ht,ctttt);
          for(jjjj=0; jjjj<nt; jjjj++)
          {
            val = coefdtttt*ctttt[jjjj];
            if(jjjj==0 || jjjj==nt-1)
	      sparseAdd(A,k,cs,-val*uana(xg[i],yg[ii],zg[iii],tg[jjjj]));
            else sparseAdd(A,k,s(i,ii,iii,jjjj),val);
	  }

	  /* for d^3u/dxdydt */
          coefdxyt  = 5.0;
          rderiv(1,nx,i,hx,cx);
          rderiv(1,ny,ii,hy,cy);
          rderiv(1,nt,iiii,ht,ct);
          for(j=0; j<nx; j++)
	  {
            for(jj=0; jj<ny; jj++)
	    {
              for(jjjj=0; jjjj<nt; jjjj++)
              {
                val = coefdxyt*cx[j]*cy[jj]*ct[jjjj];
                if(j==0 || j==nx-1 || jj==0 || jj==ny-1 ||
                   jjjj==0 || jjjj==nt-1)
		  sparseAdd(A,k,cs,-val*uana(xg[j],yg[jj],zg[iii],tg[jjjj]));
                else sparseAdd(A,k,s(j,jj,iii,jjjj),val);
	      } /* end jjjj */
	    } /* end jj */
	  } /* end j */

	  /* for d^4u/dy^2dz^2 */
          coefdyyzz = 6.0;
          rderiv(2,ny,ii,hy,cyy);
          rderiv(2,nz,iii,hz,czz);
          for(jj=0; jj<ny; jj++)
	  {
            for(jjj=0; jjj<nz; jjj++)
            {
              val = coefdyyzz*cyy[jj]*czz[jjj];
              if(jj==0 || jj==ny-1 || jjj==0 || jjj==nz-1)
		sparseAdd(A,k,cs,-val*uana(xg[i],yg[jj],zg[jjj],tg[iiii]));
              else sparseAdd(A,k,s(i,jj,jjj,iiii),val);
	    } /* end jjj */
	  } /* end jj */

	  /* for d^3u/dzdt^2 */
          coefdztt  = 7.0;
          rderiv(1,nz,iii,hz,cz);
          rderiv(2,nt,iiii,ht,ctt);
          for(jjj=0; jjj<nz; jjj++)
	  {
            for(jjjj=0; jjjj<nt; jjjj++)
            {
              val = coefdztt*cz[jjj]*ctt[jjjj];
              if(jjj==0 || jjj==nz-1 || jjjj==0 || jjjj==nt-1)
		sparseAdd(A,k,cs,-val*uana(xg[i],yg[ii],zg[jjj],tg[jjjj]));
              else sparseAdd(A,k,s(i,ii,jjj,jjjj),val);
	    } /* end jjjj */
	  } /* end jjj */

	  /* for d^3u/dx^3 */
          coefdxxx  = 8.0;
          rderiv(3,nx,i,hx,cxxx);
          for(j=0; j<nx; j++)
          {
            val = coefdxxx*cxxx[j];
            if(j==0 || j==nx-1)
	      sparseAdd(A,k,cs,-val*uana(xg[j],yg[ii],zg[iii],tg[iiii]));
            else sparseAdd(A,k,s(j,ii,iii,iiii),val);
	  }

	  /* for d^3u/dy^2dt */
          coefdyyt  = 9.0;
          rderiv(2,ny,ii,hy,cyy);
          rderiv(1,nt,iiii,ht,ct);
          for(jj=0; jj<ny; jj++)
	  {
            for(jjjj=0; jjjj<nt; jjjj++)
            {
              val = coefdyyt*cyy[jj]*ct[jjjj];
              if(jj==0 || jj==ny-1 || jjjj==0 || jjjj==nt-1)
		sparseAdd(A,k,cs,-val*uana(xg[i],yg[jj],zg[iii],tg[jjjj]));
              else sparseAdd(A,k,s(i,jj,iii,jjjj),val);
	    } /* end jjjj */
	  } /* end jj */

	  /* for d^2u/dzdt */
          coefdzt   = 10.0;
          rderiv(1,nz,iii,hz,cz);
          rderiv(1,nt,iiii,ht,ct);
          for(jjj=0; jjj<nz; jjj++)
	  {
            for(jjjj=0; jjjj<nt; jjjj++)
            {
              val = coefdzt*cz[jjj]*ct[jjjj];
              if(jjj==0 || jjj==nz-1 || jjjj==0 || jjjj==nt-1)
	        sparseAdd(A,k,cs,-val*uana(xg[i],yg[ii],zg[jjj],tg[jjjj]));
              else sparseAdd(A,k,s(i,ii,jjj,jjjj),val);
	    } /* end jjjj */
	  } /* end jjj */

	  /* for du/dt */
          coefdt    = 11.0;
          rderiv(1,nt,iiii,ht,ct);
          for(jjjj=0; jjjj<nt; jjjj++)
          {
            val = coefdt*ct[jjjj];
            if(jjjj==0 || jjjj==nt-1)
	      sparseAdd(A,k,cs,-val*uana(xg[i],yg[ii],zg[iii],tg[jjjj]));
            else sparseAdd(A,k,s(i,ii,iii,jjjj),val);
	  }

	  /* for u(x,y,z,t) */
          coefu     = 12.0;
          sparseAdd(A,k,s(i,ii,iii,iiii),coefu);

          /* for f(x,y,z,t) RHS constant */
          sparseAdd(A,k,cs,f(xg[i],yg[ii],zg[iii],tg[iiii]));

          /* end terms */

	} /* end iiii loop */
      } /* end iii loop */
    } /* end ii loop */
  } /* end i loop */

  /* printA(); lots of output */
  /* sparseWrite_All(A); */ /* optional */

  printf("A computed stiffness matrix \n");
  time_now = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("time now = %f seconds, time for previous section = %f seconds \n",
         time_now, time_now-time_last);
  time_last = time_now;
  printf("\n");

  printf("solve %d equations in %d unknowns \n", nxyzt, nxyzt);

  /* simeqb(nxyzt, A, U); */
  sparseSimeq(A,U);

  printf("equations solved \n");
  time_now = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("time now = %f seconds, time for previous section = %f seconds \n",
         time_now, time_now-time_last);
  time_last = time_now;
  printf("\n");
  
  printf("check PDE against known solution\n");
  check_soln(Ua);
  printf("check PDE against computer solution\n");
  check_soln(U);
  printf("\n");

  avgerr = 0.0;
  maxerr = 0.0;
  printf("          U computed, Ua analytic, error \n");
  for(i=1; i<nx-1; i++)
  {
    for(ii=1; ii<ny-1; ii++)
    {
      for(iii=1; iii<nz-1; iii++)
      {
        for(iiii=1; iiii<nt-1; iiii++)
        {
          err = abs(U[s(i,ii,iii,iiii)]-Ua[s(i,ii,iii,iiii)]);
          if(err>maxerr) maxerr = err;
          avgerr = avgerr + err;
          printf("ug[%d,%d,%d,%d]=%9.4f, Ua=%9.4f, err=%g \n",
                  i, ii, iii, iiii, U[s(i,ii,iii,iiii)],
                  Ua[s(i,ii,iii,iiii)], err);
        } /* iiii */
      } /* iii */
    } /* ii */
  } /* i */
  time_now = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("total CPU time  = %f seconds \n", time_now-time_start);
  printf("\n");


  printf("                                        maxerr=%g, avgerr=%g \n",
         maxerr, avgerr/(double)(nx*ny*nz*nt));
  printf("\n");
  printf("end pde44_sparse.c \n");
  return 0;
} /* end main */

/* end pde44_sparse.c */