/* pde46h_nu.c  homogeneous Biharmonic PDE in six dimensions
 * solve Uvvvv(v,w,x,y,z,t) + Uwwww(v,w,x,y,z,t) + Uxxxx(v,w,x,y,z,t) + 
 *       Uyyyy(v,w,x,y,z,t) + Uzzzz(v,w,x,y,z,t) + Utttt(v,w,x,y,z,t) + 
 *       2 Uvv(v,w,x,y,z,t) + 2 Uww(v,w,x,y,z,t) + 2 Uxx(v,w,x,y,z,t) +
 *       2 Uyy(v,w,x,y,z,t) + 2 Uzz(v,w,x,y,z,t) + 2 Utt(v,w,x,y,z,t) + 
 *         6 U(v,w,x,y,z,t) = 0
 *
 * boundary conditions computed using uana(v,w,x,y,z,t)
 * analytic solution is U(v,w,x,y,z,t) = sin(v+w+x+y+z+t)
 *
 * subscripts in code use v[i], w[ii], x[iii], y[iiii], z[iiiii], t[iiiiii]
 */

#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 nv  6
#define nw  6
#define nx  6
#define ny  6
#define nz  6
#define nt  6
#define nvwxyzt (nv-2)*(nw-2)*(nx-2)*(ny-2)*(nz-2)*(nt-2)

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

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

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

} /* end sk */

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

static double A[nvwxyzt*(nvwxyzt+1)]; /* last column is RHS */
static double U[nvwxyzt];
static double vg[nv];
static double wg[nw];
static double xg[nx];
static double yg[ny];
static double zg[nz];
static double tg[nt];
static double Ua[nvwxyzt];
static double vmin = 0.0; /* problem parameters */
static double vmax = 1.5708;
static double wmin = 0.0;
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 hv, hw, hx, hy, hz, ht;


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

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

  nuderiv(vord, nv, i, vg, cv);
  v = vg[i];
  nuderiv(word, nw, ii, wg, cw);
  w = wg[ii];
  nuderiv(xord, nx, iii, xg, cx);
  x = xg[iii];
  nuderiv(yord, ny, iiii, yg, cy);
  y = yg[iiii];
  nuderiv(zord, nz, iiiii, zg, cz);
  z = zg[iiiii];
  nuderiv(tord, nt, iiiiii, tg, ct);
  t = tg[iiiiii];
  discrete = 0.0;


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

void check_soln(double u[])
{
  /* check Uvvvv(v,w,x,y,z,t) + Uwwww(v,w,x,y,z,t) + Uxxxx(v,w,x,y,z,t) +
           Uyyyy(v,w,x,y,z,t) + Uzzzz(v,w,x,y,z,t) + Utttt(v,w,x,y,z,t) +
           2 Uvv(v,w,x,y,z,t) + 2 Uww(v,w,x,y,z,t) + 2 Uxx(v,w,x,y,z,t) +
           2 Uyy(v,w,x,y,z,t) + 2 Uzz(v,w,x,y,z,t) + 2 Utt(v,w,x,y,z,t) +
           6 u(v,w,x,y,z,t) = 0
  */
  double val, err, smaxerr;
  int i, ii, iii, iiii, iiiii, iiiiii;
  int ib=0, iib=0, iiib=0, iiiib=0, iiiiib=0, iiiiiib=0;

  smaxerr = 0.0;
  for(i=1; i<nv-1; i++)
  {
    for(ii=1; ii<nw-1; ii++)
    {
      for(iii=1; iii<nx-1; iii++)
      {
        for(iiii=1; iiii<ny-1; iiii++)
        {
          for(iiiii=1; iiiii<nz-1; iiiii++)
          {
            for(iiiiii=1; iiiiii<nt-1; iiiiii++)
            {
      val = 0.0;
      /* Uvvvv(v,w,x,y,z,t) */
      val += eval_derivative(4,0,0,0,0,0,i,ii,iii,iiii,iiiii,iiiiii,u);
      /* + Uwwww(v,w,x,y,z,t) */
      val += eval_derivative(0,4,0,0,0,0,i,ii,iii,iiii,iiiii,iiiiii,u);
      /* Uxxxx(v,w,x,y,z,t) */
      val += eval_derivative(0,0,4,0,0,0,i,ii,iii,iiii,iiiii,iiiiii,u);
      /* + Uyyyy(v,w,x,y,z,t) */
      val += eval_derivative(0,0,0,4,0,0,i,ii,iii,iiii,iiiii,iiiiii,u);
      /* + Uzzzz(v,w,x,y,z,t) */
      val += eval_derivative(0,0,0,0,4,0,i,ii,iii,iiii,iiiii,iiiiii,u);
      /* + Utttt(v,w,x,y,z,t) */
      val += eval_derivative(0,0,0,0,0,4,i,ii,iii,iiii,iiiii,iiiiii,u);
      /* + 2 Uvv(v,w,x,y,z,t) */
      val += 2.0*eval_derivative(2,0,0,0,0,0,i,ii,iii,iiii,iiiii,iiiiii,u);
      /* + 2 Uww(v,w,x,y,z,t) */
      val += 2.0*eval_derivative(0,2,0,0,0,0,i,ii,iii,iiii,iiiii,iiiiii,u);
      /* + 2 Uxx(v,w,x,y,z,t) */
      val += 2.0*eval_derivative(0,0,2,0,0,0,i,ii,iii,iiii,iiiii,iiiiii,u);
      /* + 2 Uyy(v,w,x,y,z,t) */
      val += 2.0*eval_derivative(0,0,0,2,0,0,i,ii,iii,iiii,iiiii,iiiiii,u);
      /* + 2 Uzz(v,w,x,y,z,t) */
      val += 2.0*eval_derivative(0,0,0,0,2,0,i,ii,iii,iiii,iiiii,iiiiii,u);
      /* + 2 Utt(v,w,x,y,z,t) */
      val += 2.0*eval_derivative(0,0,0,0,0,2,i,ii,iii,iiii,iiiii,iiiiii,u);
      /* + 6 u(v,w,x,y,z,t) */
      val += 6.0*uana(vg[i],wg[ii],xg[iii],yg[iiii],zg[iiiii],tg[iiiiii]);
      /* -       F(v,w,x,y,z,t)  zero */
      err = abs(val);
      if(err>smaxerr)
      {
        smaxerr = err;
        ib = i;
        iib = ii;
        iiib = iii;
        iiiib = iiii;
        iiiiib = iiiii;
        iiiiiib = iiiiii;
      }
           } /* end iiiii */
          } /* end iiiii */
        } /* end iiii */
      } /* end iii */
    } /* end ii */
  } /* end i */
  printf("check_soln maxerr=%g \n",smaxerr);
  printf("at(%1d,%1d,%1d,%1d,%1d,%1d) \n",ib,iib, iiib, iiiib, iiiiib,iiiiiib);
} /* end check_soln */

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

  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);
    return;
    exit(1);
  }
  fprintf(outp, "6 %2d %2d %2d %2d %2d %2d\n", nv, nw, nx, ny, nz, nt);
  for(i=0; i<nv; i++)
  {
    for(ii=0; ii<nw; ii++)
    {
      for(iii=0; iii<nx; iii++)
      {
        for(iiii=0; iiii<ny; iiii++)
        {
          for(iiiii=0; iiiii<nz; iiiii++)
          {
            for(iiiiii=0; iiiiii<nt; iiiiii++)
            {
              if(i==0 || i==nv-1 || ii==0 || ii==nw-1 ||
               iii==0 || iii==nx-1 || iiii==0 || iiii==ny-1 ||
		 iiiii==0 || iiiii==nz-1 || iiiiii==0 || iiiiii==nt-1)
                fprintf(outp, "%f %f %f %f %f %f %f \n",
		        vg[i],wg[ii],xg[iii],yg[iiii],zg[iiiii],tg[iiiiii],
                        uana(vg[i],wg[ii],xg[iii],yg[iiii],zg[iiiii],tg[iiiiii]));
              else
                fprintf(outp, "%f %f %f %f %f %f %f \n",
		        vg[i],wg[ii],xg[iii],yg[iiii],zg[iiiii],tg[iiiiii],
                        U[s(i,ii,iii,iiii,iiiii,iiiiii)]);
	    }
          }  
        }
      }
    }
  }
  fclose(outp);
  printf("solution file %s written \n", name);
} /* end write_soln */
 

int main(int argc, char *argv[])
{
  double v, w, x, y, z, t;
  int i,ii,iii,iiii,iiiii,iiiiii,j,jj,jjj,jjjj,jjjjj,jjjjjj;
  int k, cs;
  /* coefficients of terms in first equation */
  /* derivative point arrays, dynamic in this version */
  double cvv[nv], cww[nw], cxx[nx], cyy[ny], czz[nz],  ctt[nt]; 
  double cvvvv[nv], 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("pde46h_nu.c running \n");
  printf("Given Uvvvv(v,w,x,y,z,t) + Uwwww(v,w,x,y,z,t) + Uxxxx(v,w,x,y,z,t) +\n");
  printf("      Uyyyy(v,w,x,y,z,t) + Uzzzz(v,w,x,y,z,t) + Utttt(v,w,x,y,z,t) +\n");
  printf("      2 Uvv(v,w,x,y,z,t) + 2 Uww(v,w,x,y,z,t) + 2 Uxx(v,w,x,y,z,t) +\n");
  printf("      2 Uyy(v,w,x,y,z,t) + 2 Uzz(v,w,x,y,z,t) + 2 Utt(v,w,x,y,z,t) +\n");
  printf("      6   u(v,w,x,y,z,t) = 0 \n");
  printf("\n"); 
  printf("vmin<=v<=vmax wmin<=w<=wmax xmin<=x<=xmax \n");
  printf("ymin<=y<=ymax zmin<=z<=zmax tmin<=t<=tmax Boundaries \n");
  printf("Analytic solution U(v,w,x,y,z,t) =  sin(v+w+x+y+z+t) \n");
  printf("\n");

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

  printf("check nuderiv(4,nv,0,vg,cvvvv \n");
  nuderiv(4,nv,0,vg,cvvvv);
  for(i=0; i<nv; i++) printf("cvvvv[%d]=%f \n", i, cvvvv[i]);
  printf("check nuderiv(4,nv,nv-1,vg,cvvvv \n");
  nuderiv(4,nv,nv-1,vg,cvvvv);
  for(i=0; i<nv; i++) printf("cvvvv[%d]=%f \n", i, cvvvv[i]);
  printf("\n");

  printf("put solution in Ua vector \n");
  /* put solution in Ua vector */
  for(i=1; i<nv-1; i++)
  {
    v = vg[i];
    for(ii=1; ii<nw-1; ii++)
    {
      w = wg[ii];
      for(iii=1; iii<nx-1; iii++)
      {
        x = xg[iii];
        for(iiii=1; iiii<ny-1; iiii++)
        {
          y = yg[iiii];
          for(iiiii=1; iiiii<nz-1; iiiii++)
          {
            z = zg[iiiii];
            for(iiiiii=1; iiiiii<nt-1; iiiiii++)
            {
              t = tg[iiiiii];
              Ua[s(i,ii,iii,iiii,iiiii,iiiiii)] = uana(v,w,x,y,z,t);
            } /* 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<nv-1; i++)
  {
    for(ii=1; ii<nw-1; ii++)
    {
      for(iii=1; iii<nx-1; iii++)
      {
        for(iiii=1; iiii<ny-1; iiii++)
        {
          for(iiiii=1; iiiii<nz-1; iiiii++)
          {
            for(iiiiii=1; iiiiii<nt-1; iiiiii++)
            {
              for(j=1; j<nv-1; j++)
              {
                for(jj=1; jj<nw-1; jj++)
                {
                  for(jjj=1; jjj<nx-1; jjj++)
                  {
                    for(jjjj=1; jjjj<ny-1; jjjj++)
                    {
                      for(jjjjj=1; jjjjj<nz-1; jjjjj++)
                      {
                        for(jjjjjj=1; jjjjjj<nt-1; jjjjjj++)
                        {
                          A[ss(i,ii,iii,iiii,iiiii,iiiiii,
                               j,jj,jjj,jjjj,jjjjj,jjjjjj)] = 0.0;
	                } /* jjjjjj */
	              } /* jjjjj */
	            } /* jjjj */
	          } /* jjj */
	        } /* jj */
              } /* j */
            } /* iiiiii */
          } /* 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 = nvwxyzt; /* constant RHS column */
  val = 0.0;
  for(i=1; i<nv-1; i++)
    {
    for(ii=1; ii<nw-1; ii++)
    {
      for(iii=1; iii<nx-1; iii++)
      {
        for(iiii=1; iiii<ny-1; iiii++)
        {
          for(iiiii=1; iiiii<nz-1; iiiii++)
          {
            for(iiiiii=1; iiiiii<nt-1; iiiiii++)
            { 
              k = s(i,ii,iii,iiii,iiiii,iiiiii); /* row index */

              /* for each term in first equation */

	      /* for d^4u/dv^4 */
              nuderiv(4,nv,i,vg,cvvvv);
              for(j=0; j<nv; j++)
              {
                val = cvvvv[j];
                if(j==0 || j==nv-1)
	          A[k*(nvwxyzt+1)+cs] -= val*
		    uana(vg[j],wg[ii],xg[iii],yg[iiii],zg[iiiii],tg[iiiiii]);
                else A[sk(k,j,ii,iii,iiii,iiiii,iiiiii)] += val;
	      }

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

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

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

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

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

	      /* for 2 d^2u/dv^2 */
              nuderiv(2,nv,i,vg,cvv);
              for(j=0; j<nv; j++)
              {
                val = 2.0*cvv[j];
                if(j==0 || j==nv-1)
	          A[k*(nvwxyzt+1)+cs] -= val*
		    uana(vg[i],wg[jj],xg[iii],yg[iiii],zg[iiiii],tg[iiiiii]);
                else A[sk(k,j,ii,iii,iiii,iiiii,iiiiii)] += val;
	      }

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

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

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

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

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

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

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

              /* end terms */

	    } /* end iiiiii loop */
	  } /* 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", nvwxyzt, nvwxyzt);

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

/* end pde46h_nu.c */