/* pde_bihar44d_eq.c  discretization solution biharmonic 4th order 4D 
 * solve  uxxxx(x,y,z,t) + uyyyy(x,y,z,t) + uzzzz(x,y,z,t) + utttt(x,y,z,t) + 
 *        2*uxxyy(x,y,z,t) + 2*uxxzz(x,y,z,t)  + 2*uxxtt(x,y,z,t)  + 
 *        2*uyyzz(x,y,z,t) + 2*uyytt(x,y,z,t)  + 2*uzztt(x,y,z,t)  + 
 *        2*u(x,y,z,t) = f(x,y,z,t)
 *
 *  f(x,y,z,t) := (y^4*z^4*t^4       +
 *                 x^4*z^4*t^4       +
 *                 x^4*y^4*t^4       +
 *                 x^4*y^4*z^4       +
 *                 4*z^2*t^2         +
 *                 8*y*z^3*t^3*x     +
 *                 2*y^2*z^4*t^4*x^2 +
 *                 4*y^2*t^2         +
 *                 8*y^3*z*t^3*x     +
 *                 2*y^4*z^2*t^4*x^2 +
 *                 4*y^2*z^2         +
 *                 8*y^3*z^3*t*x     +
 *                 2*y^4*z^4*t^2*x^2 +
 *                 4*x^2*t^2         +
 *                 8*x^3*z*t^3*y     +
 *                 2*x^4*z^2*t^4*y^2 +
 *                 4*x^2*z^2         +
 *                 8*x^3*z^3*t*y     +
 *                 2*x^4*z^4*t^2*y^2 +
 *                 4*x^2*y^2         +
 *                 8*x^3*y^3*t*z     +
 *                 2*x^4*y^4*t^2*z^2 +
 *                 2                )*exp(x*y*z*t)
 *
 *  in the hyper cube  xmin< x <xmax and  ymin< y <ymax  ,
 *                     zmin< z <zmax and  tmin< t <tmax 
 *
 *     subject to boundary conditions of the first
 *     kind, Dirichlet
 *
 *  Dirichlet boundary conditions are computed using:
 *  ub(x,y,z,t) := exp(x*y*z*t)
 *
 *  uses   simeq.c  nderiv.c 
 */

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include "simeq.h"
#include "deriv.h"

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


#define nx  7
#define ny  7
#define nz  7
#define nt  7
/* DOF not including boundary */
#define nxyzt (nx-2)*(ny-2)*(nz-2)*(nt-2)

static int debug = 0;
static double xmin = 0.0; /* problem parameters */
static double xmax = 1.1; 
static double ymin = 0.0;
static double ymax = 1.2;
static double zmin = 0.0;
static double zmax = 1.3;
static double tmin = 0.0;
static double tmax = 1.4;

/* global nx, hx, xg, ny, hy, yg, nz, hz, zg, nt, ht, tg  needed by functions*/
static double xg[nx]; /* X grig */
static double yg[ny]; /* Y grid */
static double zg[ny]; /* z grid */
static double tg[ny]; /* t grid */
static double hx, hy, hz, ht;
static double ug[nxyzt];          /* solution, does not include boundary */
static double kg[nxyzt][nxyzt+1]; /* simultaneous equations, incl RHS */ 
static double Ua[nxyzt];          /* analytic solution for error check */
static double maxerr;

static double cxxxx[nx][nx];
static double cyyyy[ny][ny];
static double czzzz[nz][nz];
static double ctttt[nt][nt];
static double cxx[nx][nx];
static double cyy[ny][ny];
static double czz[nz][nz];
static double ctt[nt][nt];
/* cxxyy etc. built dynamically */

/* function prototypes */
static double F(double x, double y, double z, double t); /* RHS */
static double ub(double x, double y, double z, double t);
static void discretize(int i, int ii, int iii, int iiii);
static void dsetup();
static void write_soln();
static void check_soln();
static int s(int i, int ii, int iii, int iiii) /* index 1..nx-1  etc */
{
  return (i-1) + (ii-1)*(nx-2) + (iii-1)*(nx-2)*(ny-2) +
         (iiii-1)*(nx-2)*(ny-2)*(nz-2);
} /* end s */

int main(int argc, char * argv[])
{
  double x, y, z, t;
  int i, ii, iii, iiii, j, jj, jjj, jjjj;
  int i4; /* s(i,ii,iii,iiii) */
  int j4; /* s(j,jj,jjj,jjjj) */
  double val, err, avgerr;
  int stat;

  printf("pde_bihar44d_eq.c running.\n");
  printf("solve  uxxxx(x,y,z,t) + uyyyy(x,y,z,t) +\n");
  printf("       uzzzz(x,y,z,t) + utttt(x,y,z,t) +\n"); 
  printf("     2*uxxyy(x,y,z,t) + 2*uxxzz(x,y,z,t)  + 2*uxxtt(x,y,z,t)  +\n"); 
  printf("     2*uyyzz(x,y,z,t) + 2*uyytt(x,y,z,t)  + 2*uzztt(x,y,z,t)  +\n"); 
  printf("     2*u(x,y,z,t) = f(x,y,z,t) \n");
  printf(" \n");
  printf("f(x,y,z,t) := (y^4*z^4*t^4 + x^4*z^4*t^4 + x^4*y^4*t^4 +\n");
  printf("               x^4*y^4*z^4 + 4*z^2*t^2 + 8*y*z^3*t^3*x +\n");
  printf("               2*y^2*z^4*t^4*x^2 + 4*y^2*t^2 + 8*y^3*z*t^3*x +\n");
  printf("               2*y^4*z^2*t^4*x^2 + 4*y^2*z^2 + 8*y^3*z^3*t*x +\n");
  printf("               2*y^4*z^4*t^2*x^2 + 4*x^2*t^2 + 8*x^3*z*t^3*y +\n");
  printf("               2*x^4*z^2*t^4*y^2 + 4*x^2*z^2 + 8*x^3*z^3*t*y +\n");
  printf("               2*x^4*z^4*t^2*y^2 + 4*x^2*y^2 + 8*x^3*y^3*t*z +\n");
  printf("               2*x^4*y^4*t^2*z^2 + 2) * exp(x*y*z*t) \n");       
  printf(" \n");
  printf("in the hyper cube   xmin< x <xmax  and  ymin< y <ymax,\n");
  printf("                    zmin< z <zmax  and  tmin< t <tmax,\n");
  printf(" \n");
  printf("   subject to boundary conditions of the first kind, Dirichlet\n");
  printf("   ub(x,y,z,t) := exp(x*y*z*t) \n");
  printf("xmin=%7.3f, xmax=%7.3f, ymin=%7.3f, ymax=%7.3f\n",
         xmin, xmax, ymin, ymax);
  printf("zmin=%7.3f, zmax=%7.3f, tmin=%7.3f, tmax=%7.3f\n",
         zmin, zmax, tmin, tmax);
  printf("nx=%d, ny=%d, nz=%d, nt=%d \n",nx, ny, nz, nt);
  hx = (xmax-xmin)/(double)(nx-1);
  printf("x grid  hx=%f \n", hx);
  for(i=0; i<nx; i++) /* does not need to be uniform */
  {
    xg[i] = xmin+i*hx;
    printf("xg[%d]=%f\n", i, xg[i]);
  }  
  printf(" \n");
  hy = (ymax-ymin)/(double)(ny-1);
  printf("y grid hy=%f \n", hy);
  for(ii=0; ii<ny; ii++)
  {
    yg[ii] = ymin+ii*hy;
    printf("yg[%d]=%f\n", ii, yg[ii]);
  }  
  printf(" \n");
  hz = (zmax-zmin)/(double)(nz-1);
  printf("z grid hz=%f \n", hz);
  for(iii=0; iii<nz; iii++)
  {
    zg[iii] = zmin+iii*hz;
    printf("zg[%d]=%f\n", iii, zg[iii]);
  }  
  printf(" \n");
  ht = (tmax-tmin)/(double)(nt-1);
  printf("t grid  ht=%f \n", ht);
  for(iiii=0; iiii<nt; iiii++)
  {
    tg[iiii] = tmin+iiii*ht;
    printf("tg[%d]=%f\n", iiii, tg[iiii]);
  }  
  printf(" \n");

  /* put solution in Ua vector, not usually known */ 
  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];
          i4 = s(i,ii,iii,iiii);
          Ua[i4] = ub(x,y,z,t);
        }
      }
    }
  }
  printf(" \n");

  /* initialize kg  */ 
  for(i=1; i<nx-1; i++)
  {
    for(j=1; j<nx-1; j++)
    {
      for(ii=1; ii<ny-1; ii++)
      {
        for(jj=1; jj<ny-1; jj++)
        {
          for(iii=1; iii<nz-1; iii++)
          {
            for(jjj=1; jjj<nz-1; jjj++)
            {
              for(iiii=1; iiii<nt-1; iiii++)
              {
                for(jjjj=1; jjjj<nt-1; jjjj++)
                {
                  kg[s(i,ii,iii,iiii)][s(j,jj,jjj,jjj)] = 0.0;
		}
	      }
	    }
	  }
        }
      }
    }
  }
  
  dsetup(); /* setup discrete coefficients */
 
  /* compute entries in stiffness, equation matrix */ 
  printf("compute stiffness, equation matrix \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++)
        {
          /* compute discrete derivatives for s(i,ii,iii,iiii) equation */  
          discretize(i, ii, iii, iiii);
                        /* sets kg[s(i,ii,iii,iiii)][s(j,jj,jjj,jjjj)] */
        } /* end iiii loop */
      } /* end iii loop */
    } /* end ii loop */
  } /* end i loop */

  printf("k computed discrete matrix \n");
  printf("f computed forcing function, RHS \n");

  simeqb(nxyzt, (double*)kg, ug); /* B is kg fg */

  write_soln(); /* to file x y z t Ua soln */

  /* check_soln(); needs fixing */
  /* independent of nowing analytic solution */

  avgerr = 0.0;
  maxerr = 0.0;
  printf("ug computed by discretization, 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++)
        {
          i4 = s(i,ii,iii,iiii);
          err = abs(ug[i4]-Ua[i4]);
          if(err>maxerr) maxerr = err;
          avgerr = avgerr + err;
          printf("ug[%d,%d,%d,%d]=%f, Ua=%f, err=%g \n",
                 i, ii, iii, iiii, ug[i4], Ua[i4], ug[i4]-Ua[i4]);
        }
      }
    }
  }
  printf("xmax=%f, ymax=%f, xmax=%f, tmax=%f \n", xmax, ymax, zmax, tmax);
  printf("nx=%d, ny=%d, nz=%d, nt=%d \n", nx, ny, nz, nt);
  printf("                  maxerr=%g, avgerr=%g \n",
         maxerr,avgerr/(double)(nxyzt));
  printf("\n");

  fflush(stdout);
  printf("use plot4d_gl to see solution \n");
  
  return 0;
} /* end pde_bihar44d_eq main */

/* PDE problem definition functions  f and ub */

static double F(double x, double y, double z, double t)
{
  return (y*y*y*y*z*z*z*z*t*t*t*t + x*x*x*x*z*z*z*z*t*t*t*t +
          x*x*x*x*y*y*y*y*t*t*t*t + x*x*x*x*y*y*y*y*z*z*z*z +
          4*z*z*t*t + 8*y*z*z*z*t*t*t*x + 2*y*y*z*z*z*z*t*t*t*t*x*x +
          4*y*y*t*t + 8*y*y*y*z*t*t*t*x + 2*y*y*y*y*z*z*t*t*t*t*x*x +
          4*y*y*z*z + 8*y*y*y*z*z*z*t*x + 2*y*y*y*y*z*z*z*z*t*t*x*x +
          4*x*x*t*t + 8*x*x*x*z*t*t*t*y + 2*x*x*x*x*z*z*t*t*t*t*y*y +
          4*x*x*z*z + 8*x*x*x*z*z*z*t*y + 2*x*x*x*x*z*z*z*z*t*t*y*y +
          4*x*x*y*y + 8*x*x*x*y*y*y*t*z + 2*x*x*x*x*y*y*y*y*t*t*z*z +
          2       ) * exp(x*y*z*t);

}

static double ub(double x, double y, double z, double t)
{
  return exp(x*y*z*t);
}

static void discretize(int i, int ii, int iii, int iiii)
                   /* sets kg[s(i,ii,iii,iiii)][s(j,jj,jjj,jjjj)...] */
                   /* this code unique to this specific PDE */
{
  int cs = nxyzt;        /* RHS index */
  int i4;
  int j, jj, jjj, jjjj;
  double val;

  i4 = s(i,ii,iii,iiii); /* equation index */
  printf("discretize equation i4=%d \n", i4);
  /* for each term in each equation */

  /* for uxxxx  d^4u/dx^4 */
  for(j=0; j<nx; j++)
  {
    val = cxxxx[i][j];
    if(j==0 || j==nx-1)
      kg[i4][cs] -= val*ub(xg[j],yg[ii],zg[iii],tg[iiii]);
    else kg[i4][s(j,ii,iii,iiii)] += val;
  }

  /* for uyyyy  d^4u/dy^4 */
  for(jj=0; jj<ny; jj++)
  {
    val = cyyyy[ii][jj];
    if(jj==0 || jj==ny-1)
      kg[i4][cs] -=  val*ub(xg[i],yg[jj],zg[iii],tg[iiii]);
    else kg[i4][s(i,jj,iii,iiii)] += val;
  }

  /* for uzzzz  d^4u/dz^4 */
  for(jjj=0; jjj<nz; jjj++)
  {
    val = cyyyy[iii][jjj];
    if(jjj==0 || jjj==nz-1)
      kg[i4][cs] -=  val*ub(xg[i],yg[ii],zg[jjj],tg[iiii]);
    else kg[i4][s(i,ii,jjj,iiii)] += val;
  }

  /* for utttt  d^4u/dt^4 */
  for(jjjj=0; jjjj<nt; jjjj++)
  {
    val = cyyyy[iiii][jjjj];
    if(jjjj==0 || jjjj==nt-1)
      kg[i4][cs] -=  val*ub(xg[i],yg[ii],zg[iii],tg[jjjj]);
    else kg[i4][s(i,ii,iii,jjjj)] += val;
  }

  /* for uxxyy  d^4u/dx^2dy^2 */
  for(j=0; j<nx; j++)
  {
    for(jj=0; jj<ny; jj++)
    {
      val = 2.0*cxx[i][j]*cyy[ii][jj];
      if(j==0 || j==nx-1 || jj==0 || jj==ny-1)
	kg[i4][cs] -= val*ub(xg[j],yg[jj],zg[iii],tg[iiii]);
      else kg[i4][s(j,jj,iii,iiii)] += val;
    } /* end jj */
  } /* end j */

  /* for uxxzz  d^4u/dx^2dz^2 */
  for(j=0; j<nx; j++)
  {
    for(jjj=0; jjj<nz; jjj++)
    {
      val = 2.0*cxx[i][j]*czz[iii][jjj];
      if(j==0 || j==nx-1 || jjj==0 || jjj==nz-1)
	kg[i4][cs] -= val*ub(xg[j],yg[ii],zg[jjj],tg[iiii]);
      else kg[i4][s(j,ii,jjj,iiii)] += val;
    } /* end jjj */
  } /* end j */

  /* for uxxtt  d^4u/dx^2dt^2 */
  for(j=0; j<nx; j++)
  {
    for(jjjj=0; jjjj<nt; jjjj++)
    {
      val = 2.0*cxx[i][j]*ctt[iiii][jjjj];
      if(j==0 || j==nx-1 || jjjj==0 || jjjj==nt-1)
	kg[i4][cs] -= val*ub(xg[j],yg[ii],zg[iii],tg[jjjj]);
      else kg[i4][s(j,ii,iii,jjjj)] += val;
    } /* end jjjj */
  } /* end j */

  /* for uyyzz  d^4u/dy^2dz^2 */
  for(jj=0; jj<ny; jj++)
  {
    for(jjj=0; jjj<nz; jjj++)
    {
      val = 2.0*cyy[ii][jj]*czz[iii][jjj];
      if(jj==0 || jj==ny-1 || jjj==0 || jjj==nz-1)
	kg[i4][cs] -= val*ub(xg[i],yg[jj],zg[jjj],tg[iiii]);
      else kg[i4][s(i,jj,jjj,iiii)] += val;
    } /* end jjj */
  } /* end jj */

  /* for uyytt  d^4u/dy^2dt^2 */
  for(jj=0; jj<ny; jj++)
  {
    for(jjjj=0; jjjj<nt; jjjj++)
    {
      val = 2.0*cyy[ii][jj]*ctt[iiii][jjjj];
      if(jj==0 || jj==ny-1 || jjjj==0 || jjjj==nt-1)
	kg[i4][cs] -= val*ub(xg[i],yg[jj],zg[iii],tg[jjjj]);
      else kg[i4][s(i,jj,iii,jjjj)] += val;
    } /* end jjjj */
  } /* end jj */

  /* for uzztt  d^4u/dz^2dt^2 */
  for(jjj=0; jjj<nz; jjj++)
  {
    for(jjjj=0; jjjj<nt; jjjj++)
    {
      val = 2.0*czz[iii][jjj]*ctt[iiii][jjjj];
      if(jjj==0 || jjj==nz-1 || jjjj==0 || jjjj==nt-1)
	kg[i4][cs] -= val*ub(xg[i],yg[ii],zg[jjj],tg[jjjj]);
      else kg[i4][s(i,ii,jjj,jjjj)] += val;
    } /* end jjjj */
  } /* end jjj */

  kg[i4][i4] += 2.0; /* U */ 
} /* end discretize */

static void dsetup()
{
  int i, ii, iii, iiii, j;

  /* each such as cxxxx[i][j] is j_th coefficient at DOF i */
  for(i=0; i<nx; i++) rderiv(4, nx, i, hx, cxxxx[i]);
  if(debug)
  {
    printf("debug print of cxxxx[i][j] \n");
    for(i=0; i<nx; i++)
    {
      printf("i=%d ",i);
      for(j=0; j<nx; j++) printf("%f ", cxxxx[i][j]);
      printf(" \n");
    }
  }
  for(i=0; i<nx; i++) rderiv(2, nx, i, hx, cxx[i]);
  if(debug)
  {
    printf("debug print of cxx[i][j] \n");
    for(i=0; i<nx; i++)
    {
      printf("i=%d ",i);
      for(j=0; j<nx; j++) printf("%f ", cxx[i][j]);
      printf(" \n");
    }
  }
  for(ii=0; ii<ny; ii++) rderiv(4, ny, ii, hy, cyyyy[ii]);
  for(ii=0; ii<ny; ii++) rderiv(2, ny, ii, hy, cyy[ii]);
  for(iii=0; iii<nz; iii++) rderiv(4, nz, iii, hz, czzzz[iii]);
  for(iii=0; iii<nz; iii++) rderiv(2, nz, iii, hz, czz[iii]);
  for(iiii=0; iiii<nt; iiii++) rderiv(4, nt, iiii, ht, ctttt[iiii]);
  for(iiii=0; iiii<nt; iiii++) rderiv(2, nt, iiii, ht, ctt[iiii]);
} /* end dsetup */

static void write_soln()
{
  FILE * fout;
  int i4;
  int i, ii, iii, iiii;

  fout = fopen("pde_bihar44d_eq_c.dat", "w");
  if(fout==NULL)
  {
    printf("file pde_bihar44d_eq_c.dat can not be opened for writing\n");
    return;
  }
  printf("writing pde_bihar44d_eq_c.dat \n"); /* including boundary */
  fprintf(fout,"4 %d %d %d %d\n", nx, ny, nz, nt);
  for(i=0; i<nx; i++)
  {
    for(ii=0; ii<ny; ii++)
    {
      for(iii=0; iii<nz; iii++)
      {
        for(iiii=0; iiii<nt; iiii++)
        {
          if(i==0 || i==nx-1 || ii==0 || ii==ny-1 ||
             iii==0 || iii==nz-1 || iiii==0 || iiii==nt-1)
	  {
            fprintf(fout," %f %f %f %f %f \n",xg[i],yg[ii],zg[iii],tg[iiii],
                    ub(xg[i],yg[ii],zg[iii],tg[iiii]));
	  }
          else
	  {
            i4 = s(i,ii,iii,iiii);
            fprintf(fout," %f %f %f %f %f \n",xg[i],yg[ii],zg[iii],tg[iiii],
                    ug[i4]);
	  }
	} /* end t */
      } /* end z */
    } /* end y */
  } /* end x */
  fflush(fout);
  fclose(fout);
  printf("finished writing pde_bihar44d_eq_c.dat \n");
} /* end write_soln */

static void check_soln()
{
               /* for all interior points */
  double x, y, z, t, err, maxerr, avgerr, rmserr;
  double sumerr, sumsqerr, eval_deriv;
  int cnt = 0;
  int i4;
  double tcxxyy[ny];
  int check=0;
  int i, ii, iii, iiii, j, jj, jjj, jjjj;

  printf("check_soln values \n");
  sumerr = 0.0;
  sumsqerr = 0.0;
  maxerr = 0.0;
  for(i=1; i<nx-1; i++) /* full equation */
  {
    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];
          i4 = s(i,ii,iii,iiii);

          if(i==nx-2 && ii==ny-2 && iii==nz-2 && iiii==nt-2) check = 1;
          if(check)
	  {
            printf("\nchecking on check against PDE i4=%d \n", i4);
            printf("xg[i]=%f, yg[ii]=%f, zg[iii]=%f, tg[iiii]=%f \n",
                   xg[i], yg[ii], zg[iii], tg[iiii]);
	  }
          /* PDE 2 u(x,y,z,t) term */
          eval_deriv = 2.0*ug[i4];
          if(check) printf("2.0*ug = %f \n", 2.0*ug[i4]);

          /* PDE uxxxx term */          
          for(j=0; j<nx; j++)
          {
	    eval_deriv += cxxxx[i][j]*ug[s(j,ii,iii,iiii)];
            if(check) printf("cxxxx[i][%d]*ug[s(j,ii,iii,iiii)] = %f \n",
			     j, cxxxx[i][j]*ug[s(j,ii,iii,iiii)]);
          } 
          /* PDE uyyyy term */
          for(jj=0; jj<ny; jj++)
          {
	    eval_deriv += cyyyy[ii][jj]*ug[s(i,jj,iii,iiii)];
            if(check) printf("cyyyy[ii][%d]*ug[s(i,jj,iii,iiii)] = %f \n",
			     jj, cxxxx[ii][jj]*ug[s(i,jj,iii,iiii)]);
          }
          /* PDE uzzzz term */
          for(jjj=0; jjj<nz; jjj++)
          {
	    eval_deriv += czzzz[iii][jjj]*ug[s(i,ii,jjj,iiii)];
            if(check) printf("czzzz[iii][%d]*ug[s(i,ii,jjj,iiii)] = %f \n",
			     jjj, czzzz[iii][jjj]*ug[s(i,ii,jjj,iiii)]);
          } 
          /* PDE utttt term */
          for(jjjj=0; jjjj<nt; jjjj++)
          {
	    eval_deriv += ctttt[iiii][jjjj]*ug[s(i,ii,iii,jjjj)];
            if(check) printf("ctttt[iiii][%d]*ug[s(i,ii,iii,jjjj)] = %f \n",
			     jjjj, ctttt[iiii][jjjj]*ug[s(i,ii,iii,jjjj)]);
          }
          /* PDE  2*uxxyy term */
          for(j=0; j<nx; j++)
          {
            for(jj=0; jj<ny; jj++)
            {
              eval_deriv += 2.0*cxx[i][j]*cyy[ii][jj];
            } /* end jj */
          } /* end j */
          /* PDE  2*uxxzz term */
          for(j=0; j<nx; j++)
          {
            for(jjj=0; jjj<nz; jjj++)
            {
              eval_deriv += 2.0*cxx[i][j]*czz[iii][jjj];
            } /* end jjj */
          } /* end j */
          /* PDE  2*uxxtt term */
          for(j=0; j<nx; j++)
          {
            for(jjjj=0; jjjj<nt; jjjj++)
            {
              eval_deriv += 2.0*cxx[i][j]*ctt[iiii][jjjj];
            } /* end jjjj */
          } /* end j */
          /* PDE  2*uyyzz term */
          for(jj=0; jj<ny; jj++)
          {
            for(jjj=0; jjj<nz; jjj++)
            {
              eval_deriv += 2.0*cyy[ii][jj]*czz[iii][jjj];
            } /* end jjj */
          } /* end jj */
          /* PDE  2*uyytt term */
          for(jj=0; jj<ny; jj++)
          {
            for(jjjj=0; jjjj<nt; jjjj++)
            {
              eval_deriv += 2.0*cyy[ii][jj]*ctt[iiii][jjjj];
            } /* end jjjj */
          } /* end jj */
          /* PDE  2*uzztt term */
          for(jjj=0; jjj<nz; jjj++)
          {
            for(jjjj=0; jjjj<nt; jjjj++)
            {
              eval_deriv += 2.0*czz[iii][jjj]*ctt[iiii][jjjj];
            } /* end jjjj */
          } /* end jjj */


          err = eval_deriv-F(x,y,z,t);
          err = abs(err);
          sumerr += err;
          sumsqerr += err*err;
          maxerr = max(err,maxerr);
          cnt++;
	} /* end t */
      } /* end z */ 
    } /* end y */
  } /* end x */
  rmserr = sqrt(sumsqerr/(double)cnt);
  avgerr = sumerr/(double)cnt;
  printf("maxerr=%g, rmserr=%g, avgerr=%g \n",maxerr,rmserr,avgerr);
} /* end check_soln */ 
/* end pde_bihar44d_eq.c */