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

#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 nx  5 /* 6 makes 4^7 = 2^14 = 16384 eqn in 16384 unk */
#define ny  5
#define nz  5
#define nt  5
#define nu  5
#define nv  5
#define nw  5
#define nxyztuvw (nx-2)*(ny-2)*(nz-2)*(nt-2)*(nu-2)*(nv-2)*(nw-2)

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

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

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

} /* end sk */

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

static double A[nxyztuvw*(nxyztuvw+1)]; /* last column is RHS */
static double U[nxyztuvw];
static double xg[nx];
static double yg[ny];
static double zg[nz];
static double tg[nt];
static double ug[nu];
static double vg[nv];
static double wg[nw];
static double Ua[nxyztuvw];
static double xmin = 0.0; /* problem parameters */
static double xmax = 0.5;
static double ymin = 0.0;
static double ymax = 0.5;
static double zmin = 0.0;
static double zmax = 0.5;
static double tmin = 0.0;
static double tmax = 0.5;
static double umin = 0.0;
static double umax = 0.5;
static double vmin = 0.0;
static double vmax = 0.5;
static double wmin = 0.0;
static double wmax = 0.5;
static double hx, hy, hz, ht, hu, hv, hw;

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

static double f(double x, double y, double z, double t, double u,
                double v, double w)
                              /* RHS */
{
  return sin(x+y+z+t+u+v+w);
}

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

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

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

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

  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++)
        {
          for(iiiii=1; iiiii<nu-1; iiiii++)
          {
            for(iiiiii=1; iiiiii<nv-1; iiiiii++)
            {
              for(iiiiiii=1; iiiiiii<nw-1; iiiiiii++)
              {
    val = 0.0;
    /* Uxxxx(x,y,z,t,u,v,w) */
    val += eval_derivative(4,0,0,0,0,0,0,i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,u);
    /* + Uyyyy(x,y,z,t,u,v,w) */
    val += eval_derivative(0,4,0,0,0,0,0,i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,u);
    /* + Uzzzz(x,y,z,t,u,v,w) */
    val += eval_derivative(0,0,4,0,0,0,0,i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,u);
    /* + Utttt(x,y,z,t,u,v,w) */
    val += eval_derivative(0,0,0,4,0,0,0,i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,u);
    /* + Uuuuu(x,y,z,t,u,v,w) */
    val += eval_derivative(0,0,0,0,4,0,0,i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,u);
    /* + Uvvvv(x,y,z,t,u,v,w) */
    val += eval_derivative(0,0,0,0,0,4,0,i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,u);
    /* + Uwwww(x,y,z,t,u,v,w) */
    val += eval_derivative(0,0,0,0,0,0,4,i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,u);
    /* + 2 Uxx(x,y,z,t,u,v,w) */
    val += 2.0*eval_derivative(2,0,0,0,0,0,0,i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,u);
    /* + 2 Uyy(x,y,z,t,u,v,w) */
    val += 2.0*eval_derivative(0,2,0,0,0,0,0,i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,u);
    /* + 2 Uzz(x,y,z,t,u,v,w) */
    val += 2.0*eval_derivative(0,0,2,0,0,0,0,i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,u);
    /* + 2 Utt(x,y,z,t,u,v,w) */
    val += 2.0*eval_derivative(0,0,0,2,0,0,0,i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,u);
    /* + 2 Uuu(x,y,z,t,u,v,w) */
    val += 2.0*eval_derivative(0,0,0,0,2,0,0,i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,u);
    /* + 2 Uvv(x,y,z,t,u,v,w) */
    val += 2.0*eval_derivative(0,0,0,0,0,2,0,i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,u);
    /* + 2 Uww(x,y,z,t,u,v,w) */
    val += 2.0*eval_derivative(0,0,0,0,0,0,2,i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,u);
    /* + 8 u(x,y,z,t,u,v,w) */
    val += 8.0*Ub(xg[i],yg[ii],zg[iii],tg[iiii],ug[iiiii],
                  vg[iiiiii],wg[iiiiiii]);
    /* - f(x,y,z,t,u,v,w) */
    val -= f(xg[i],yg[ii],zg[iii],tg[iiii],ug[iiiii],
                  vg[iiiiii],wg[iiiiiii]);
    err = abs(val);
    if(err>smaxerr)
    {
      smaxerr = err;
      ib = i;
      iib = ii;
      iiib = iii;
      iiiib = iiii;
      iiiiib = iiiii;
      iiiiiib = iiiiii;
      iiiiiiib = iiiiiii;
    }
             } /* end iiiiiii */
           } /* end iiiiii */
          } /* 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] = "pde47h_nu_c";
  char extn[] = ".dat";
  char * name;
  int i, ii, iii, iiii, iiiii, iiiiii, iiiiiii;

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

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

  printf("pde47h_nu.c running \n");
  printf("Given Uxxxx(x,y,z,t,u,v,w) + Uyyyy(x,y,z,t,u,v,w) + Uzzzz(x,y,z,t,u,v,w) +\n");
  printf("      Utttt(x,y,z,t,u,v,w) + Uuuuu(x,y,z,t,u,v,w) + Uvvvv(x,y,z,t,u,v,w) +\n");
  printf("      Uwwww(x,y,z,t,u,v,w) +\n");
  printf("      2 Uxx(x,y,z,t,u,v,w) + 2 Uyy(x,y,z,t,u,v,w) + 2 Uzz(x,y,z,t,u,v,w) +\n");
  printf("      2 Utt(x,y,z,t,u,v,w) + 2 Uuu(x,y,z,t,u,v,w) + 2 Uvv(x,y,z,t,u,v,w) +\n");
  printf("      2 Uww(x,y,z,t,u,v,w) +\n");
  printf("      8   u(x,y,z,t,u,v,w) = f(x,y,z,t,u,v,w) \n");
  printf("\n"); 
  printf("xmin<=x<=xmax ymin<=y<=ymax zmin<=z<=zmax \n");
  printf("tmin<=t<=tmax umin<=u<=umax vmin<=v<=vmax \n");
  printf("wmin<=w<=wmax  Boundaries \n");
  printf("Analytic solution U(x,y,z,t,u,v,w) =  sin(x+y+z+t+u+v+w) \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("umin=%g, umax=%g, vmin=%g, vmax=%g \n", umin, umax, vmin, vmax);
  printf("wmin=%g, wmax=%g \n", wmin, wmax);
  printf("nx=%d, ny=%d, nz=%d, nt=%d, nu=%d, nv=%d, nw=%d \n",
         nx, ny, nz, nt, nu, nv, nw);
  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 minimums \n");
  hx = (double)(nx*(nx-1))/2.0;
  i = 0;
  xg[i] = xmin;
  printf("i=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
         i, xg[i], ymin, zmin, tmin, umin, vmin, wmin, Ua[i]);
  for(i=1; i<nx; i++)
  {
    xg[i] = xg[i-1]+(xmax-xmin)*(double)(nx-i)/hx;
    Ua[i] = Ub(xg[i],ymin,zmin,tmin,umin,vmin,wmin); /* temp */
    printf("i=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
	   i, xg[i], ymin, zmin, tmin, umin, vmin, wmin, Ua[i]);
  } /* i  x */  
  printf("\n");

  printf("y grid and analytic solution at minimums \n");
  hy = (double)(ny*(ny-1))/2.0;
  ii = 0;
  yg[ii] = ymin;
  printf("ii=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
         ii, xmin, yg[ii], zmin, tmin, umin, vmin, wmin, Ua[ii]);
  for(ii=1; ii<ny; ii++)
  {
    yg[ii] = yg[ii-1]+(ymax-ymin)*(double)(ny-ii)/hy;
    Ua[ii] = Ub(xmin, yg[ii], zmin, tmin, umin, vmin, wmin); /* temp */
    printf("ii=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
	   ii, xmin, yg[ii], zmin, tmin, umin, vmin, wmin, Ua[ii]);
  } /* ii  y */
  printf("\n");

  printf("z grid and analytic solution at minimums \n");
  hz = (double)(nz*(nz-1))/2.0;
  iii = 0;
  zg[iii] = zmin;
  for(iii=1; iii<nz; iii++)
  {
    zg[iii] = zg[iii-1]+(zmax-zmin)*(double)(nz-iii)/hz;
    Ua[iii] = Ub(xmin, ymin, zg[iii], tmin, umin, vmin, wmin); /* temp */
    printf("iii=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
	   iii, xmin, ymin, zg[iii], tmin, umin, vmin, wmin, Ua[iii]);
  } /* iii   z */
  printf("\n");

  printf("t grid and analytic solution at minimums \n");
  ht = (double)(nt*(nt-1))/2.0;
  iiii = 0;
  tg[iiii] = tmin;
  printf("iiii=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
         iiii, xmin, ymin, zmin, tg[iiii], umin, vmin, wmin, Ua[iiii]);
  for(iiii=1; iiii<nt; iiii++)
  {
    tg[iiii] = tg[iiii-1]+(tmax-tmin)*(double)(nt-iiii)/ht;
    Ua[iiii] = Ub(xmin, ymin, zmin, tg[iiii], umin, vmin, wmin); /* temp */
    printf("iiii=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
	   iiii, xmin, ymin, zmin, tg[iiii], umin, vmin, wmin, Ua[iiii]);
  } /* iiii  t */
  printf("\n");

  printf("u grid and analytic solution at minimums \n");
  hu = (double)(nu*(nu-1))/2.0;
  iiiii = 0;
  ug[iiiii] = umin;
  printf("iiiii=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
         iiiii, xmin, ymin, zmin, tmin, ug[iiiii], vmin, wmin, Ua[iiiii]);
  for(iiiii=1; iiiii<nu; iiiii++)
  {
    ug[iiiii] = ug[iiiii-1]+(umax-umin)*(double)(nu-iiiii)/hu;
    Ua[iiiii] = Ub(xmin,ymin,zmin,tmin,ug[iiiii],vmin,wmin); /* temp */
    printf("iiiii=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
	   iiiii, xmin, ymin, zmin, tmin, ug[iiiii], vmin, wmin, Ua[iiiii]);
  } /* iiiii  u */  
  printf("\n");

  printf("v grid and analytic solution at minimums \n");
  hv = (double)(nv*(nv-1))/2.0;
  iiiiii = 0;
  vg[iiiiii] = vmin;
  printf("iiiiii=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
         iiiiii, xmin, ymin, zmin, tmin, umin, vg[iiiiii], wmin, Ua[iiiiii]);
  for(iiiiii=1; iiiiii<nv; iiiiii++)
  {
    vg[iiiiii] = vg[iiiiii-1]+(vmax-vmin)*(double)(nv-iiiiii)/hv;
    Ua[iiiiii] = Ub(xmin,ymin,zmin,tmin,umin,vg[iiiiii],wmin); /* temp */
    printf("iiiiii=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
	   iiiiii, xmin, ymin, zmin, tmin, umin, vg[iiiiii], wmin, Ua[iiiiii]);
  } /* iiiiii  v */  
  printf("\n");

  printf("w grid and analytic solution at minimums \n");
  hw = (double)(nw*(nw-1))/2.0;
  iiiiiii = 0;
  wg[iiiiiii] = wmin;
  printf("iiiiiii=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
         iiiiiii, vmin, wg[ii], xmin, ymin, zmin, tmin, Ua[ii]);
  for(iiiiiii=1; iiiiiii<nw; iiiiiii++)
  {
    wg[iiiiiii] = wg[iiiiiii-1]+(wmax-wmin)*(double)(nw-iiiiiii)/hw;
    Ua[iiiiiii] = Ub(xmin,ymin,zmin,tmin,umin,vmin,wg[iiiiiii]); /* temp */
    printf("iiiiiii=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
	   iiiiiii, xmin, ymin, zmin, tmin, umin, vmin, wg[iiiiiii], Ua[iiiiiii]);
  } /* iiiiiii  w */  
  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<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];
          for(iiiii=1; iiiii<nu-1; iiiii++)
          {
            u = ug[iiiii];
            for(iiiiii=1; iiiiii<nv-1; iiiiii++)
            {
              v = vg[iiiiii];
              for(iiiiiii=1; iiiiiii<nw-1; iiiiiii++)
              {
                w = wg[iiiiiii];
                Ua[s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)] = Ub(x,y,z,t,u,v,w);
              } /* iiiiiii */
            } /* 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");

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

                /* for each term in first equation */

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

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

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

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

	        /* for d^4U/du^4 */
                nuderiv(4,nu,iiiii,ug,cuuuu);
                for(jjjjj=0; jjjjj<nu; jjjjj++)
                {
                  val = cuuuu[jjjjj];
                  if(jjjjj==0 || jjjjj==nu-1)
	            A[k*(nxyztuvw+1)+cs] -= val*
		      Ub(xg[i],yg[ii],zg[iii],tg[iiii],ug[jjjjj],
		         vg[iiiiii],wg[iiiiiii]);
                  else A[sk(k,i,ii,iii,iiii,jjjjj,iiiiii,iiiiiii)] += val;
	        }

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

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

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

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

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

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

	        /* for 2 d^2U/du^2 */
                nuderiv(2,nu,iiiii,ug,cuu);
                for(jjjjj=0; jjjjj<nu; jjjjj++)
                { 
                  val = 2.0*cuu[jjjjj];
                  if(jjjjj==0 || jjjjj==nu-1)
	            A[k*(nxyztuvw+1)+cs] -= val*
		      Ub(xg[i],yg[ii],zg[iii],tg[iiii],ug[jjjjj],
		         vg[iiiiii],wg[iiiiiii]);
                  else A[sk(k,i,ii,iii,iiii,jjjjj,iiiiii,iiiiiii)] += val;
	        }

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

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

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

                /* for f(x,y,z,t,u,v,w) RHS */
	        A[k*(nxyztuvw+1)+cs] += f(x,y,z,t,u,v,w);
 
                /* end terms */

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

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

/* end pde47h_nu.c */