/* pde64sb_eq.c  discretize, build linear equations, solve biharmonic
 * nabla^4 P(x,y,z,t,u,v) + 2 nabla^2 P(x,y,z,t,u,v) = F(x,y,z,t,u,v) 
 *
 * solve Pxxxx(x,y,z,t,u,v) + Pyyyy(x,y,z,t,u,v) + Pzzzz(x,y,z,t,u,v) +
 *       Ptttt(x,y,z,t,u,v) + Puuuu(x,y,z,t,u,v) + Pvvvv(x,y,z,t,u,v) + 
 *       2*(Pxx(x,y,z,t,u,v) + Pyy(x,y,z,t,u,v) + Pzz(x,y,z,t,u,v) +
 *       Ptt(x,y,z,t,u,v) + Puu(x,y,z,t,u,v) + Pvv(x,y,z,t,u,v)) = 
 *       F(x,y,z,t,u,v)
 * F(x,y,z,t,u,v) = see d6p4sb_rhs.jpg 
 *              
 * boundary conditions computed using P(x,y,z,t,u,v)
 * analytic solution is P(x,y,z,t,u,v) = sin(x*u+y*v+z+t)
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.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
#define ny  5
#define nz  5
#define nt  5
#define nu  5
#define nv  5
#define nxyztuv (nx-2)*(ny-2)*(nz-2)*(nt-2)*(nu-2)*(nv-2)

static double A[nxyztuv*(nxyztuv+1)]; /* last column is RHS */
static double U[nxyztuv];  /* computed solution */
static double Ua[nxyztuv]; /* known solution when testing */
static double xg[nx];      /* grid coordinates */
static double yg[ny];
static double zg[nz];
static double tg[nt];
static double ug[nt];
static double vg[nt];
static double xmin = -0.5; /* problem parameters */
static double xmax =  0.5; /* for nx = 5 */
static double ymin = -0.5;
static double ymax =  0.5;
static double zmin = -0.5;
static double zmax =  0.5;
static double tmin = -0.5;
static double tmax =  0.5;
static double umin = -0.5;
static double umax =  0.5;
static double vmin = -0.5;
static double vmax =  0.5;

/* derivative point arrays, static in this version, [point][term] */
static double cxx[nx][nx], cyy[ny][ny], czz[nz][nz], 
              ctt[nt][nt], cuu[nu][nu], cvv[nv][nv]; 
static double cxxxx[nx][nx], cyyyy[ny][ny], czzzz[nz][nz], 
              ctttt[nt][nt], cuuuu[nu][nu], cvvvv[nv][nv];
static double hx, hy, hz, ht, hu, hv;

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

int sk(int k, int i, int ii, int iii, int iiii, int iiiii, int iiiiii) 
       /* for x,y,z,t,u,v */
{
  return k*(nxyztuv+1) + 
         (i-1)*(ny-2)*(nz-2)*(nt-2)*(nu-2)*(nv-2) + 
         (ii-1)*(nz-2)*(nt-2)*(nu-2)*(nv-2) +
         (iii-1)*(nt-2)*(nu-2)*(nv-2) + 
         (iiii-1)*(nu-2)*(nv-2) + 
         (iiiii-1)*(nv-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)*(nxyztuv+1) + 
         s(j,jj,jjj,jjjj,jjjjj,jjjjjj);
} /* end ss */


double f(double x, double y, double z, double t, double u, double v)
{
  double x2, x4, y2, y4, z2, z4, t2, t4, u2, u4, v2, v4;
  double cscc, csss, ccsc, cccs, scss, sssc, sscs, sccc;

  x2 = x*x;
  x4 = x2*x2;
  y2 = y*y;
  y4 = y2*y2;
  z2 = z*z;
  z4 = z2*z2;
  t2 = t*t;
  t4 = t2*t2;
  u2 = u*u;
  u4 = u2*u2;
  v2 = v*v;
  v4 = v2*v2;

  cscc = cos(x*u)*sin(y*v)*cos(z)*cos(t);
  csss = cos(x*u)*sin(y*v)*sin(z)*sin(t);
  ccsc = cos(x*u)*cos(y*v)*sin(z)*cos(t);
  cccs = cos(x*u)*cos(y*v)*cos(z)*sin(t);
  scss = sin(x*u)*cos(y*v)*sin(z)*sin(t);
  sssc = sin(x*u)*sin(y*v)*sin(z)*cos(t);
  sscs = sin(x*u)*sin(y*v)*cos(z)*sin(t);
  sccc = sin(x*u)*cos(y*v)*cos(z)*cos(t);

  return
     2.0*sssc + 2.0*scss - 2.0*sccc + 2.0*sscs
   - 2.0*cscc - 2.0*ccsc + 2.0*csss - 2.0*cccs
   - x4*sssc - x4*sscs + x4*cscc - x4*csss
   + x4*ccsc + x4*cccs - 2.0*u2*sccc  + 2.0*u2*scss
   + 2.0*u2*sssc  + 2.0*u2*sscs  - 2.0*u2*cscc  + 2.0*u2*csss
   - 2.0*u2*ccsc  - 2.0*u2*cccs  - 2.0*y2*cccs  - 2.0*v2*sccc
   + 2.0*v2*scss  + 2.0*v2*sssc  + 2.0*v2*sscs  - 2.0*v2*cscc
   + 2.0*v2*csss  - 2.0*v2*ccsc  - 2.0*v2*cccs  - 2.0*x2*sccc
   + 2.0*x2*scss  + 2.0*x2*sssc  + 2.0*x2*sscs  - 2.0*x2*cscc
   + 2.0*x2*csss  - 2.0*x2*ccsc  - 2.0*x2*cccs  - 2.0*y2*sccc
   + 2.0*y2*scss  + 2.0*y2*sssc  + 2.0*y2*sscs  - 2.0*y2*cscc
   + 2.0*y2*csss
   - 2.0*y2*ccsc + y4*sccc  - y4*scss - y4*sssc
   - y4*sscs + y4*cscc  - y4*csss + y4*ccsc
   + y4*cccs + u4*sccc  - u4*scss - u4*sssc
   - u4*sscs + u4*cscc  - u4*csss + u4*ccsc
   + u4*cccs + v4*sccc  - v4*scss - v4*sssc
   - v4*sscs + v4*cscc  - v4*csss + v4*ccsc
   + v4*cccs + x4*sccc  - x4*scss;
}

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

void printA()
{
  int i, ii, iii, iiii, iiiii, iiiiii, j, jj, jjj, jjjj, jjjjj, jjjjjj, cs, k;
  cs = (nx-2)*(ny-2)*(nz-2)*(nt-2)*(nu-2)*(nv-2); /* constant RHS column */
  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++)
            {
              k = s(i,ii,iii,iiii,iiiii,iiiiii); /* row index */
              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++)
                        {
                          printf(
  "A[row(%1d,%1d,%1d,%1d,%1d,%1d,%1d),col(%1d,%1d,%1d,%1d,%1d,%1d)] =%f\n",
                             i,ii,iii,iiii,iiiii,iiiiii,
                             j,jj,jjj,jjjj,jjjjj,jjjjjj,
                             A[ss(i,ii,iii,iiii,iiiii,iiiiii,
                             j,jj,jjj,jjjj,jjjjj,jjjjjj)]);
	                } /* jjjjjj */
	              } /* jjjjj */
	            } /* jjjj */
	          } /* jjj */
	        } /* jj */
              } /* j */
              printf("RHS A[%d,%d]=%f\n\n", k, cs, A[k*(nxyztuv+1)+cs]);
              /* equivalent A[sk(k,i,ii,iii,iiii,iiiii,iiiiii+1)] */
            } /* iiiiii */
          } /* iiiii */
        } /* iiii */
      } /* iii */
    } /* ii */
  } /* i */
  printf("\n");
} /* end printA */

void build_deriv()
{
  int i, ii, iii, iiii, iiiii, iiiiii;

  printf("build_deriv \n");
  for(i=0; i<nx; i++)
  {
    rderiv(4,nx,i,hx,cxxxx[i]);
    rderiv(2,nx,i,hx,cxx[i]);
  }
  for(ii=0; ii<ny; ii++)
  {
    rderiv(4,ny,ii,hy,cyyyy[ii]);
    rderiv(2,ny,ii,hy,cyy[ii]);
  }
  for(iii=0; iii<nz; iii++)
  {
    rderiv(4,nz,iii,hz,czzzz[iii]);
    rderiv(2,nz,iii,hz,czz[iii]);
  }
  for(iiii=0; iiii<nt; iiii++)
  {
    rderiv(4,nt,iiii,ht,ctttt[iiii]);
    rderiv(2,nt,iiii,ht,ctt[iiii]);
  }
  for(iiiii=0; iiiii<nu; iiiii++)
  {
    rderiv(4,nu,iiiii,hu,cuuuu[iiiii]);
    rderiv(2,nu,iiiii,hu,cuu[iiiii]);
  }
  for(iiiiii=0; iiiiii<nv; iiiiii++)
  {
    rderiv(4,nv,iiiiii,hv,cuuuu[iiiiii]);
    rderiv(2,nv,iiiiii,hv,cuu[iiiiii]);
  }
} /* end build_deriv */

void init_matrix()
{
  int i, ii, iii, iiii, iiiii, iiiiii, j, jj, jjj, jjjj, jjjjj, jjjjjj;

  printf("initialize matrix A \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(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++)
                        {
                          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 */
} /* end init matrix */

void save_known_soln()
{
  int i, ii, iii, iiii, iiiii, iiiiii;
  double x, y, z, t, u, v;

  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];
              Ua[s(i,ii,iii,iiii,iiiii,iiiiii)] = pana(x,y,z,t,u,v);
            } /* iiiiii */
          } /* iiiii */
        } /* iiii */
      } /* iii */
    } /* ii */
  } /* i */
} /* end save_known_soln */

void build_equations() /* specific to this PDE */
{
  double val;
  int k, cs;
  int i, ii, iii, iiii, iiiii, iiiiii, j, jj, jjj, jjjj, jjjjj, jjjjjj;

  /* compute entries in A matrix, inside boundary */
  printf("build equations, A matrix \n");
  cs = (nx-2)*(ny-2)*(nz-2)*(nt-2)*(nu-2)*(nv-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++)
        {
          for(iiiii=1; iiiii<nu-1; iiiii++)
          {
            for(iiiiii=1; iiiiii<nv-1; iiiiii++)
            {
              k = s(i,ii,iii,iiii,iiiii,iiiiii); /* row index */

              /* for each term in first equation */

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

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

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

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

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

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

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

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

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

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

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

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

              /* for f(x,y,z,t,u,v) RHS constant */
	      A[k*(nxyztuv+1)+cs] += 
              f(xg[i],yg[ii],zg[iii],tg[iiii],ug[iiiii],vg[iiiiii]);

              /* end terms */

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

void check_soln(double UU[]) /* specific to this PDE */
{
  double val, smaxerr, err, est;
  int k;
  int i, ii, iii, iiii, iiiii, iiiiii, j, jj, jjj, jjjj, jjjjj, jjjjjj;

  /* compute entries in A matrix, inside boundary */
  printf("check UU against PDE \n");
  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++)
            {
              k = s(i,ii,iii,iiii,iiiii,iiiiii); /* row index */
              est = 0.0;
              /* for each term in equation */

	      /* for d^4P/dx^4 */
              for(j=0; j<nx; j++)
              {
                val = cxxxx[i][j];
                if(j==0 || j==nx-1)
	          est += 
		  val*pana(xg[j],yg[ii],zg[iii],tg[iiii],ug[iiiii],vg[iiiiii]);
		else est += val*UU[s(j,ii,iii,iiii,iiiii,iiiiii)];
	      }

	      /* for d^4P/dy^4 */
              for(jj=0; jj<ny; jj++)
              {
                val = cyyyy[ii][jj];
                if(jj==0 || jj==ny-1)
	          est +=  
		  val*pana(xg[i],yg[jj],zg[iii],tg[iiii],ug[iiiii],vg[iiiiii]);
		else est += val*UU[s(i,jj,iii,iiii,iiiii,iiiiii)];
	      }

	      /* for d^4P/dz^4 */
              for(jjj=0; jjj<nz; jjj++)
              {
                val = czzzz[iii][jjj];
                if(jjj==0 || jjj==nz-1)
	          est += 
		  val*pana(xg[i],yg[ii],zg[jjj],tg[iiii],ug[iiiii],vg[iiiiii]);
                else est += val*UU[s(i,ii,jjj,iiii,iiiii,iiiiii)];
	      }

	      /* for d^4P/dt^4 */
              for(jjjj=0; jjjj<nt; jjjj++)
              {
                val = ctttt[iiii][jjjj];
                if(jjjj==0 || jjjj==nt-1)
	          est += 
		  val*pana(xg[i],yg[ii],zg[iii],tg[jjjj],ug[iiiii],vg[iiiiii]);
                else est += val*UU[s(i,ii,iii,jjjj,iiiii,iiiiii)];
	      }

	      /* for d^4P/du^4 */
              for(jjjjj=0; jjjjj<nu; jjjjj++)
              {
                val = cuuuu[iiiii][jjjjj];
                if(jjjjj==0 || jjjjj==nu-1)
	          est += 
		  val*pana(xg[i],yg[ii],zg[iii],tg[iiii],ug[jjjjj],vg[iiiiii]);
                else est += val*UU[s(i,ii,iii,iiii,jjjjj,iiiiii)];
	      }

	      /* for d^4P/dv^4 */
              for(jjjjjj=0; jjjjjj<nv; jjjjjj++)
              {
                val = cvvvv[iiiiii][jjjjjj];
                if(jjjjjj==0 || jjjjjj==nv-1)
	          est += 
		  val*pana(xg[i],yg[ii],zg[iii],tg[iiii],ug[iiiii],vg[jjjjjj]);
                else est += val*UU[s(i,ii,iii,iiii,iiiii,jjjjjj)];
	      }

	      /* for d^2P/dx^4 */
              for(j=0; j<nx; j++)
              {
                val = 2.0*cxx[i][j];
                if(j==0 || j==nx-1)
	          est += 
		  val*pana(xg[j],yg[ii],zg[iii],tg[iiii],ug[iiiii],vg[iiiiii]);
		else est += val*UU[s(j,ii,iii,iiii,iiiii,iiiiii)];
	      }

	      /* for d^2P/dy^4 */
              for(jj=0; jj<ny; jj++)
              {
                val = 2.0*cyy[ii][jj];
                if(jj==0 || jj==ny-1)
	          est +=  
		  val*pana(xg[i],yg[jj],zg[iii],tg[iiii],ug[iiiii],vg[iiiiii]);
		else est += val*UU[s(i,jj,iii,iiii,iiiii,iiiiii)];
	      }

	      /* for d^2P/dz^4 */
              for(jjj=0; jjj<nz; jjj++)
              {
                val = 2.0*czz[iii][jjj];
                if(jjj==0 || jjj==nz-1)
	          est += 
		  val*pana(xg[i],yg[ii],zg[jjj],tg[iiii],ug[iiiii],vg[iiiiii]);
                else est += val*UU[s(i,ii,jjj,iiii,iiiii,iiiiii)];
	      }

	      /* for d^2P/dt^4 */
              for(jjjj=0; jjjj<nt; jjjj++)
              {
                val = 2.0*ctt[iiii][jjjj];
                if(jjjj==0 || jjjj==nt-1)
	          est += 
		  val*pana(xg[i],yg[ii],zg[iii],tg[jjjj],ug[iiiii],vg[iiiiii]);
                else est += val*UU[s(i,ii,iii,jjjj,iiiii,iiiiii)];
	      }

	      /* for d^2P/du^4 */
              for(jjjjj=0; jjjjj<nu; jjjjj++)
              {
                val = 2.0*cuu[iiiii][jjjjj];
                if(jjjjj==0 || jjjjj==nu-1)
	          est += 
		  val*pana(xg[i],yg[ii],zg[iii],tg[iiii],ug[jjjjj],vg[iiiiii]);
                else est += val*UU[s(i,ii,iii,iiii,jjjjj,iiiiii)];
	      }

	      /* for d^2P/dv^4 */
              for(jjjjjj=0; jjjjjj<nv; jjjjjj++)
              {
                val = 2.0*cvv[iiiiii][jjjjjj];
                if(jjjjjj==0 || jjjjjj==nv-1)
	          est += 
		  val*pana(xg[i],yg[ii],zg[iii],tg[iiii],ug[iiiii],vg[jjjjjj]);
                else est += val*UU[s(i,ii,iii,iiii,iiiii,jjjjjj)];
	      }

              /* for f(x,y,z,t,u,v) RHS constant */
	      err = est - 
              f(xg[i],yg[ii],zg[iii],tg[iiii],ug[iiiii],vg[iiiiii]);
              err = abs(err);
              smaxerr = max(smaxerr,err);
              /* if(err>0.001)printf("%1d %1d %1d %1d %1d %1d est=%e err=%e \n",
		 i,ii,iii,iiii,iiiii,iiiiii,est,err); */
              /* end terms */

	    } /* end iiiiii loop */
	  } /* end iiiii loop */
	} /* end iiii loop */
      } /* end iii loop */
    } /* end ii loop */
  } /* end i loop */
  printf("check against PDE, smaxerr=%e \n", smaxerr);
} /* end check_soln */

void write_soln(char tag[], double u[]) /* for plot6d_gl */
{
  FILE* outp;
  char name1[32] = "pde64sb_eq_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", nx, ny, nz, nt, nu, nv);
  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++)
            {
              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)
                fprintf(outp, "%f %f %f %f %f %f %f \n",
		        xg[i],yg[ii],zg[iii],tg[iiii],ug[iiiii],vg[iiiiii],
                        pana(xg[i],yg[ii],zg[iii],tg[iiii],ug[iiiii],vg[iiiiii]));
              else
                fprintf(outp, "%f %f %f %f %f %f %f \n",
		        xg[i],yg[ii],zg[iii],tg[iiii],ug[iiiii],vg[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 x, y, z, t, u, v;
  int i, ii, iii, iiii, iiiii, iiiiii, j, jj, jjj, jjjj, jjjjj, jjjjjj;
  int k, cs;
  double val, err, avgerr, maxerr;
  double time_start;
  double time_now;
  double time_last;

  printf("pde64sb_eq.c running \n");
  printf("Given nabla^4 P(x,y,z,t,u,v) + 2 nabla^2 P(x,y,z,t,u,v) = F(x,y,z,t,u,v) \n");
  printf("  F(x,y,z,t,u,v) = exp(x+y+z*t+u+v)*(12.0+t*t*t*t+z*z*z*z+2.0*t*t+2.0*z*z); \n");
  printf(" solve Pxxxx(x,y,z,t,u,v) + Pyyyy(x,y,z,t,u,v) + Pzzzz(x,y,z,t,u,v)\n");
  printf("       Ptttt(x,y,z,t,u,v) + Puuuu(x,y,z,t,u,v) + Pvvvv(x,y,z,t,u,v) +\n");
  printf("       2*(Pxx(x,y,z,t,u,v) + Pyy(x,y,z,t,u,v) + Pzz(x,y,z,t,u,v) +\n");
  printf("       Ptt(x,y,z,t,u,v) + Puu(x,y,z,t,u,v) + Pvv(x,y,z,t,u,v)) =\n");  printf("       F(x,y,z,t,u,v) \n");
  printf(" \n");
  printf("F(x,y,z,t,u,v) = see d6p4sb_rhs.jpg \n");
  printf(" \n");
  printf("xmin<=x<=xmax ymin<=y<=ymax zmin<=z<=zmax Boundaries \n");
  printf("tmin<=t<=tmax umin<=u<=umax vmin<=v<=vmax Boundaries \n");
  printf("Analytic solution P(x,y,z,t,u,v) =  sin(x*u+y*v+z+t) \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("nx=%1d, ny=%1d, nz=%1d, nt=%1d, nu=%1d, nv=%1d \n", 
	 nx, ny, nz, nt, nu, nv);
  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,umin,vmin \n");
  hx = (xmax-xmin)/(double)(nx-1);
  for(i=0; i<nx; i++)
  {
    xg[i] = xmin+(double)i*hx;
    Ua[i] = pana(xg[i],ymin,zmin,tmin,umin,vmin); /* temp */
    printf("i=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
	   i, xg[i], ymin, zmin, tmin, umin, vmin, Ua[i]);
  } /* i */  
  printf("\n");
  printf("y grid and analytic solution at xmin,zmin,tmin,umin,vmin \n");
  hy = (ymax-ymin)/(double)(ny-1);
  for(ii=0; ii<ny; ii++)
  {
    yg[ii] = ymin+(double)ii*hy;
    Ua[ii] = pana(xmin, yg[ii],zmin,tmin,umin,vmin); /* temp */
    printf("ii=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
	   ii, xmin, yg[ii], zmin, tmin, umin, vmin, Ua[ii]);
  } /* ii */
  printf("\n");
  printf("z grid and analytic solution at xmin,ymin,tmin,umin,vmin \n");
  hz = (zmax-zmin)/(double)(nz-1);
  for(iii=0; iii<nz; iii++)
  {
    zg[iii] = zmin+(double)iii*hz;
    Ua[iii] = pana(xmin, ymin, zg[iii],tmin,umin,vmin); /* temp */
    printf("iii=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
	   iii, xmin, ymin, zg[iii], tmin, umin, vmin, Ua[iii]);
  } /* iii */
  printf("\n");
  printf("t grid and analytic solution at xmin,ymin,zmin,umin,vmin \n");
  ht = (tmax-tmin)/(double)(nt-1);
  for(iiii=0; iiii<nt; iiii++)
  {
    tg[iiii] = tmin+(double)iiii*ht;
    Ua[iiii] = pana(xmin, ymin, zmin, tg[iiii], umin, vmin); /* temp */
    printf("iiii=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
	   iiii, xmin, ymin, zmin, tg[iiii], umin, vmin, Ua[iiii]);
  } /* iiii */
  printf("\n");
  printf("u grid and analytic solution at xmin,ymin,zmin,tmin,vmin \n");
  hu = (umax-umin)/(double)(nu-1);
  for(iiiii=0; iiiii<nu; iiiii++)
  {
    ug[iiiii] = umin+(double)iiiii*hu;
    Ua[iiiii] = pana(xmin, ymin, zmin, tmin, ug[iiiii], vmin); /* temp */
    printf("iiiii=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
	   iiiii, xmin, ymin, zmin, tmin, ug[iiiii], vmin, Ua[iiiii]);
  } /* iiiii */
  printf("\n");
  printf("v grid and analytic solution at xmin,ymin,zmin,tmin,umin \n");
  hv = (vmax-vmin)/(double)(nv-1);
  for(iiiiii=0; iiiiii<nv; iiiiii++)
  {
    vg[iiiiii] = vmin+(double)iiiiii*hv;
    Ua[iiiiii] = pana(xmin, ymin, zmin, tmin, umin, vg[iiiiii]); /* temp */
    printf("iiiiii=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f \n",
	   iiiiii, xmin, ymin, zmin, tmin, umin, vg[iiiiii], Ua[iiiiii]);
  } /* iiiiii */
  printf("\n");

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

  save_known_soln();

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

  init_matrix(); /* A */

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

  build_equations(); /* specific to this PDE */

  /* printA(); lots of output */

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

  simeqb(nxyztuv, 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");


  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++)
            {
              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[%1d,%1d,%1d,%1d,%1d,%1d]=%9.4f, Ua=%9.4f, err=%g \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 */

  write_soln("", U);
  printf("                                        maxerr=%g, avgerr=%g \n",
         maxerr, avgerr/(double)(nx*ny*nz*nt));
  printf("\n");
  printf("check known against PDE \n");
  check_soln(Ua);
  printf("check solution against PDE \n");
  check_soln(U);

  time_now = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("total CPU time  = %f seconds \n", time_now-time_start);
  printf("\n");
  printf("end pde64sb_eq.c \n");
  return 0;
} /* end main */
/* end pde64sb_eq.c */