/* pde45h_eq.c  homogeneous Biharmonic PDE in five dimensions
 * solve uwwww(w,x,y,z,t) + uxxxx(w,x,y,z,t) + uyyyy(w,x,y,z,t) + 
 *       uzzzz(w,x,y,z,t) + utttt(w,x,y,z,t) + 2 uww(w,x,y,z,t) +
 *       2 uxx(w,x,y,z,t) + 2 uyy(w,x,y,z,t) + 2 uzz(w,x,y,z,t) +
 *       2 utt(w,x,y,z,t) + 5 u(w,x,y,z,t) = 0
 *
 * boundary conditions computed using u(w,x,y,z,t)
 * analytic solution is u(w,x,y,z,t) = sin(w+x+y+z+t)
 *
 * subscripts in code use w[i], x[ii], y[iii], z[iiii], t[iiiii]
 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>
#include "deriv.h"
#include "simeq.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 nw  6
#define nx  6
#define ny  6
#define nz  6
#define nt  6
#define nwxyzt (nw-2)*(nx-2)*(ny-2)*(nz-2)*(nt-2)

int s(int i, int ii, int iii, int iiii, int iiiii) /* for w,x,y,z,t */
{
  int k;

  k = (i-1)*(nx-2)*(ny-2)*(nz-2)*(nt-2) +
      (ii-1)*(ny-2)*(nz-2)*(nt-2) + (iii-1)*(nz-2)*(nt-2) +
      (iiii-1)*(nt-2) + (iiiii-1);
  if(k<0 || k>=nwxyzt) printf("bad s(%d %d %d %d %d\n",i,ii,iii,iiii,iiiii);  
  return k;
} /* end s */

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

} /* end sk */

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

static double A[nwxyzt*(nwxyzt+1)]; /* last column is RHS */
static double U[nwxyzt];
static double wg[nw];
static double xg[nx];
static double yg[ny];
static double zg[nz];
static double tg[nt];
static double Ua[nwxyzt];
static double wmin = 0.0; /* problem parameters */
static double wmax = 1.5708;
static double xmin = 0.0;
static double xmax = 1.5708;
static double ymin = 0.0;
static double ymax = 1.5708;
static double zmin = 0.0;
static double zmax = 1.5708;
static double tmin = 0.0;
static double tmax = 1.5708;
static double hw, hx, hy, hz, ht;


static double uana(double w, double x, double y, double z, double t)
                                        /*analytic solution and boundaries*/
{
  return sin(w+x+y+z+t);
}

double eval_derivative(int word, int xord, int yord, int zord, int tord,
                       int i, int ii, int iii, int iiii, int iiiii, double u[])
{
  int j, jj, jjj, jjjj, jjjjj;
  double cw[nw];
  double cx[nx];
  double cy[ny];
  double cz[nz];
  double ct[nt];
  double discrete;
  double w, x, y, z, t;

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


  /* five cases of one partial w, x, y, z, t to any order */
  if(word!=0 && xord==0 && yord==0 && zord==0 && tord==0)
    for(j=0; j<nw; j++)
      if(j==0 || j==nw-1) discrete += cw[j]*uana(wg[j],x,y,z,t);
      else discrete += cw[j]*u[s(j,ii,iii,iiii,iiiii)];
  else if(word==0 && xord!=0 && yord==0 && zord==0 && tord==0)
    for(jj=0; jj<nx; jj++)
      if(jj==0 || jj==nx-1) discrete += cx[jj]*uana(w,xg[jj],y,z,t);
      else discrete += cx[jj]*u[s(i,jj,iii,iiii,iiiii)];
  else if(word==0 && xord==0 && yord!=0 && zord==0 && tord==0)
    for(jjj=0; jjj<ny; jjj++)
      if(jjj==0 || jjj==ny-1) discrete += cy[jjj]*uana(w,x,yg[jjj],z,t);
      else discrete += cy[jjj]*u[s(i,ii,jjj,iiii,iiiii)];
  else if(word==0 && xord==0 && yord==0 && zord!=0 && tord==0)
    for(jjjj=0; jjjj<nz; jjjj++)
      if(jjjj==0 || jjjj==nz-1) discrete += cz[jjjj]*uana(w,x,y,zg[jjjj],t);
      else discrete += cz[jjjj]*u[s(i,ii,iii,jjjj,iiiii)];
  else if(word==0 && xord==0 && yord==0 && zord==0 && tord!=0)
    for(jjjjj=0; jjjjj<nt; jjjjj++)
      if(jjjjj==0 || jjjjj==nt-1) discrete += ct[jjjjj]*uana(w,x,y,z,tg[jjjjj]);
      else discrete += ct[jjjjj]*u[s(i,ii,iii,iiii,jjjjj)];
  return discrete;
} /* end eval_deriv */

void check_soln(double u[])
{
  /* check uwwww(w,x,y,z,t) + uxxxx(w,x,y,z,t) + uyyyy(w,x,y,z,t) +
           uzzzz(w,x,y,z,t) + utttt(w,x,y,z,t) + 2 uww(w,x,y,z,t) + 
           2 uxx(w,x,y,z,t) + 2 uyy(w,x,y,z,t) + 2 uzz(w,x,y,z,t) + 
           2 utt(w,x,y,z,t) + 5 u(w,x,y,z,t) = 0
  */
  double val, err, smaxerr;
  int i, ii, iii, iiii, iiiii;

  smaxerr = 0.0;
  for(i=1; i<nw-1; i++)
  {
    for(ii=1; ii<nx-1; ii++)
    {
      for(iii=1; iii<ny-1; iii++)
      {
        for(iiii=1; iiii<nz-1; iiii++)
        {
          for(iiiii=1; iiiii<nt-1; iiiii++)
          {
            val = 0.0;
            /* uwwww(w,x,y,z,t) */
            val += eval_derivative(4,0,0,0,0,i,ii,iii,iiii,iiiii,u);
            /* uxxxx(w,x,y,z,t) */
            val += eval_derivative(0,4,0,0,0,i,ii,iii,iiii,iiiii,u);
	    /* + uyyyy(w,x,y,z,t) */
            val += eval_derivative(0,0,4,0,0,i,ii,iii,iiii,iiiii,u);
	    /* + uzzzz(w,x,y,z,t) */
            val += eval_derivative(0,0,0,4,0,i,ii,iii,iiii,iiiii,u);
            /* + utttt(w,x,y,z,t) */
            val += eval_derivative(0,0,0,0,4,i,ii,iii,iiii,iiiii,u);
            /* + 2 uww(w,x,y,z,t) */
            val += 2.0*eval_derivative(2,0,0,0,0,i,ii,iii,iiii,iiiii,u);
            /* + 2 uxx(w,x,y,z,t) */
            val += 2.0*eval_derivative(0,2,0,0,0,i,ii,iii,iiii,iiiii,u);
	    /* + 2 uyy(w,x,y,z,t) */
            val += 2.0*eval_derivative(0,0,2,0,0,i,ii,iii,iiii,iiiii,u);
	    /* + 2 uzz(w,x,y,z,t) */
            val += 2.0*eval_derivative(0,0,0,2,0,i,ii,iii,iiii,iiiii,u);
	    /* + 2 utt(w,x,y,z,t) */
            val += 2.0*eval_derivative(0,0,0,0,2,i,ii,iii,iiii,iiiii,u);
	    /* + 5 u(w,x,y,z,t) */
            val += 5.0*uana(wg[i],xg[ii],yg[iii],zg[iiii],tg[iiiii]);
            /* -       F(w,x,y,z,t)  zero */
	    err = abs(val);
            if(err>smaxerr) smaxerr = err;
          } /* end iiiii */
        } /* end iiii */
      } /* end iii */
    } /* end ii */
  } /* end i */
  printf("check_soln maxerr=%g \n",smaxerr);
} /* end check_soln */

void write_soln(char tag[], double u[]) /* for plot5d_gl */
{
  FILE* outp;
  char name1[32] = "pde45h_eq_c";
  char extn[] = ".dat";
  char * name;
  int i, ii, iii, iiii, iiiii;

  strcat(name1,tag);
  name = strcat(name1,extn);
  printf("writing solution to %s \n", name);
  outp = fopen(name,"w");
  if(outp==NULL)
  {
    printf("ERROR unable to open %s for writing \n", name);
    exit(1);
  }
  fprintf(outp, "5 %2d %2d %2d %2d %2d \n", nw, nx, ny, nz, nt);
  for(i=1; i<nw-1; i++)
  {
    for(ii=1; ii<nx-1; ii++)
    {
      for(iii=1; iii<ny-1; iii++)
      {
        for(iiii=1; iiii<nz-1; iiii++)
        {
          for(iiiii=1; iiiii<nt-1; iiiii++)
          {
            fprintf(outp, "%f %f %f %f %f %f\n",
	    wg[i],xg[ii],yg[iii],zg[iiii],tg[iiiii],U[s(i,ii,iii,iiii,iiiii)]);
          }  
        }
      }
    }
  }
  fclose(outp);
  printf("solution file %s written \n", name);
} /* end write_soln */
 

int main(int argc, char *argv[])
{
  double w, x, y, z, t;
  int i,ii,iii,iiii,iiiii, j, jj, jjj, jjjj, jjjjj;
  int k, cs;
  /* coefficients of terms in first equation */
  /* derivative point arrays, dynamic in this version */
  double cww[nw], cxx[nx], cyy[ny], czz[nz],  ctt[nt]; 
  double cwwww[nw], cxxxx[nx], cyyyy[ny], czzzz[nz], ctttt[nt];
  double val, err, avgerr, maxerr;
  double time_start;
  double time_now;
  double time_last;

  printf("pde45h_eq.c running \n");
  printf("Given uwwww(w,x,y,z,t) + uxxxx(w,x,y,z,t) + uyyyy(w,x,y,z,t) +\n");
  printf("      uzzzz(w,x,y,z,t) + utttt(w,x,y,z,t) + 2 uww(w,x,y,z,t) +\n");
  printf("      2 uxx(w,x,y,z,t) + 2 uyy(w,x,y,z,t) + 2 uzz(w,x,y,z,t) +\n");
  printf("      2 utt(w,x,y,z,t) + 5 u(w,x,y,z,t) = 0 \n");
  printf("\n"); 
  printf("wmin<=w<=wmax xmin<=x<=xmax ymin<=y<=ymax zmin<=z<=zmax \n");
  printf("tmin<=t<=tmax Boundaries \n");
  printf("Analytic solution u(w,x,y,z,t) =  sin(w+x+y+z+t) \n");
  printf("\n");

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

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

  printf("put solution in Ua vector \n");
  /* put solution in Ua vector */
  for(i=1; i<nw-1; i++)
  {
    w = wg[i];
    for(ii=1; ii<nx-1; ii++)
    {
      x = xg[ii];
      for(iii=1; iii<ny-1; iii++)
      {
        y = yg[iii];
        for(iiii=1; iiii<nz-1; iiii++)
        {
          z = zg[iiii];
          for(iiiii=1; iiiii<nt-1; iiiii++)
          {
            t = tg[iiiii];
            Ua[s(i,ii,iii,iiii,iiiii)] = uana(w,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, w, ii, x, iii, y, iiii, z, iiiii, t, Ua[s(i,ii,iii,iiii,iiiii)]); */
          } /* iiiii */
        } /* 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");

  printf("initialize A \n");
  /* initialize A */
  for(i=1; i<nw-1; i++)
  {
    for(ii=1; ii<nx-1; ii++)
    {
      for(iii=1; iii<ny-1; iii++)
      {
        for(iiii=1; iiii<nz-1; iiii++)
        {
          for(iiiii=1; iiiii<nt-1; iiiii++)
          {
            for(j=1; j<nw-1; j++)
            {
              for(jj=1; jj<nx-1; jj++)
              {
                for(jjj=1; jjj<ny-1; jjj++)
                {
                  for(jjjj=1; jjjj<nz-1; jjjj++)
                  {
                    for(jjjjj=1; jjjjj<nt-1; jjjjj++)
                    {
                      A[ss(i,ii,iii,iiii,iiiii,j,jj,jjj,jjjj,jjjjj)] = 0.0;
	            } /* jjjjj */
	          } /* jjjj */
	        } /* jjj */
	      } /* jj */
            } /* j */
          } /* iiiii */
        } /* 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");

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

            /* for each term in first equation */

	    /* for d^4u/dw^4 */
            rderiv(4,nw,i,hw,cwwww);
            for(j=0; j<nw; j++)
            {
              val = cwwww[j];
              if(j==0 || j==nw-1)
	        A[k*(nwxyzt+1)+cs] -= val*
                              uana(wg[j],xg[ii],yg[iii],zg[iiii],tg[iiiii]);
              else A[sk(k,j,ii,iii,iiii,iiiii)] += val;
	    }

	    /* for d^4u/dx^4 */
            rderiv(4,nx,ii,hx,cxxxx);
            for(jj=0; jj<nx; jj++)
            {
              val = cxxxx[jj];
              if(jj==0 || jj==nx-1)
		A[k*(nwxyzt+1)+cs] -= val*
                              uana(wg[i],xg[jj],yg[iii],zg[iiii],tg[iiiii]);
              else A[sk(k,i,jj,iii,iiii,iiiii)] += val;
	    }

	    /* for d^4u/dy^4 */
            rderiv(4,ny,iii,hy,cyyyy);
            for(jjj=0; jjj<ny; jjj++)
            {
              val = cyyyy[jjj];
              if(jjj==0 || jjj==ny-1)
	        A[k*(nwxyzt+1)+cs] -=  val*
                              uana(wg[i],xg[ii],yg[jjj],zg[iiii],tg[iiiii]);
              else A[sk(k,i,ii,jjj,iiii,iiiii)] += val;
	    }

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

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

	    /* for 2 d^2u/dw^2 */
            rderiv(2,nw,i,hw,cww);
            for(j=0; j<nw; j++)
            {
              val = 2.0*cww[j];
              if(j==0 || j==nw-1)
	        A[k*(nwxyzt+1)+cs] -= val*
		              uana(wg[j],xg[ii],yg[iii],zg[iiii],tg[iiiii]);
              else A[sk(k,j,ii,iii,iiii,iiiii)] += val;
	    }

	    /* for 2 d^2u/dx^2 */
            rderiv(2,nx,ii,hx,cxx);
            for(jj=0; jj<nx; jj++)
            {
              val = 2.0*cxx[jj];
              if(jj==0 || jj==nx-1)
	        A[k*(nwxyzt+1)+cs] -= val*
	                      uana(wg[i],xg[jj],yg[iii],zg[iiii],tg[iiiii]);
              else A[sk(k,i,jj,iii,iiii,iiiii)] += val;
	    }

	    /* for 2 d^2u/dy^2 */
            rderiv(2,ny,iii,hy,cyy);
            for(jjj=0; jjj<ny; jjj++)
            {
              val = 2.0*cyy[jjj];
              if(jjj==0 || jjj==ny-1)
	        A[k*(nwxyzt+1)+cs] -=  val*
                              uana(wg[i],xg[ii],yg[jjj],zg[iiii],tg[iiiii]);
              else A[sk(k,i,ii,jjj,iiii,iiiii)] += val;
	    }

	    /* for 2 d^2u/dz^2 */
            rderiv(2,nz,iiii,hz,czz);
            for(jjjj=0; jjjj<nz; jjjj++)
            { 
              val = 2.0*czz[jjjj];
              if(jjjj==0 || jjjj==nz-1)
	        A[k*(nwxyzt+1)+cs] -= val*
                              uana(wg[i],xg[ii],yg[iii],zg[jjjj],tg[iiiii]);
              else A[sk(k,i,ii,iii,jjjj,iiiii)] += val;
	    }

	    /* for 2 d^2u/dt^2 */
            rderiv(2,nt,iiiii,ht,ctt);
            for(jjjjj=0; jjjjj<nt; jjjjj++)
            {
              val = 2.0*ctt[jjjjj];
              if(jjjjj==0 || jjjjj==nt-1)
	        A[k*(nwxyzt+1)+cs] -= val*
                              uana(wg[i],xg[ii],yg[iii],zg[iiii],tg[jjjjj]);
              else A[sk(k,i,ii,iii,iiii,jjjjj)] += val;
	    }

	    /* for 5 u(w,x,y,z,t) */
            A[sk(k,i,ii,iii,iiii,iiiii)] += 5.0;

            /* for f(w,x,y,z,t) RHS zero */

            /* end terms */

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

  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", nwxyzt, nwxyzt);

  simeqb(nwxyzt, A, U); /* A destroyed */

  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);
  write_soln("1",Ua);

  printf("check PDE against computed solution\n");
  check_soln(U);
  write_soln("2",U);
  printf("\n");

  avgerr = 0.0;
  maxerr = 0.0;
  printf("          U computed, Ua analytic, error \n");
  for(i=1; i<nw-1; i++)
  {
    for(ii=1; ii<nx-1; ii++)
    {
      for(iii=1; iii<ny-1; iii++)
      {
        for(iiii=1; iiii<nz-1; iiii++)
        {
          for(iiiii=1; iiiii<nt-1; iiiii++)
          {
            err = abs(U[s(i,ii,iii,iiii,iiiii)]-Ua[s(i,ii,iii,iiii,iiiii)]);
            if(err>maxerr) maxerr = err;
            avgerr = avgerr + err;
            printf("ug[%d,%d,%d,%d,%d]=%9.4f, Ua=%9.4f, err=%g \n",
                   i,ii,iii,iiii,iiiii,U[s(i,ii,iii,iiii,iiiii)],
		   Ua[s(i,ii,iii,iiii,iiiii)], err);
          } /* iiiii */
        } /* 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=%e, avgerr=%e \n",
         maxerr, avgerr/(double)(nwxyzt));
  printf("\n");
  printf("end pde45h_eq.c \n");
  return 0;
} /* end main */

/* end pde45h_eq.c */