// pde_bihar44e_eq.java  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*uxx(x,y,z,t) + 2*uyy(x,y,z,t) + 2*uzz(x,y,z,t) + 
//        2*utt(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+
//                2*y^2*z^2*t^2+
//                2*x^2*z^2*t^2+
//                2*x^2*y^2*t^2+
//                2*x^2*y^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 
 
import java.io.*;
//                uses rderiv.java and simeq.java

class pde_bihar44e_eq
{
  //  nx, ny, ny, nz, xg, yg, zg, tg needed by many, available at call 
  int nx, ny, nz, nt;
  int i4;
  double xg[];   // X grig
  double yg[];   // Y grid
  double zg[];   // Z grig
  double tg[];   // T grid
  double ug[];   // solution
  double kg[][]; // matrix
  double fg[];   // RHS
  double maxerr;
  double cxxxx[][];
  double cyyyy[][];
  double czzzz[][];
  double ctttt[][];
  double cxx[][];
  double cyy[][];
  double czz[][];
  double ctt[][];
  boolean debug = true;
  double hx, hy, hz, ht;
  int nxyzt;

  pde_bihar44e_eq(){}; // construct the class once, run many times

  void run(double xmax, double ymax, double zmax, double tmax,
           int anx, int any, int anz, int ant)
  {
    nx = anx; // needed global
    ny = any;
    nz = anz;
    nt = ant;
    nxyzt = (nx-2)*(ny-2)*(nz-2)*(nt-2);
    xg = new double[nx];
    yg = new double[ny];
    zg = new double[nz];
    tg = new double[nt];
    double x, y, z, t;
    double val, err, avgerr;

    double xmin = 0.0; // problem parameters 
    double ymin = 0.0;
    double zmin = 0.0;
    double tmin = 0.0;

    kg = new double[nxyzt][nxyzt];  
    fg = new double[nxyzt];
    ug = new double[nxyzt];
    double Ua[] = new double[nxyzt];

    cxxxx = new double[nx][nx];
    cyyyy = new double[ny][ny];
    czzzz = new double[nz][nz];
    ctttt = new double[nt][nt];
    cxx = new double[nx][nx];
    cyy = new double[ny][ny];
    czz = new double[nz][nz];
    ctt = new double[nt][nt];

    System.out.println("pde_bihar44e_eq.java running.");
    System.out.println("solve  uxxxx(x,y,z,t) + uyyyy(x,y,z,t) + uzzzz(x,y,z,t) + utttt(x,y,z,t) + ");
    System.out.println("       2*uxx(x,y,z,t) + 2*uyy(x,y,z,t)  + 2*uzz(x,y,z,t)  + ");
    System.out.println("       2*utt(x,y,z,t) + 2*u(x,y,z,t) = f(x,y,z,t)");
    System.out.println(" ");
    System.out.println("f(x,y,z,t) := (y^4*z^4*t^4 + x^4*z^4*t^4 +");
    System.out.println("               x^4*y^4*t^4 + x^4*y^4*z^4 +" );
    System.out.println("               2*y^2*z^2*t^2 + 2*x^2*z^2*t^2 +");
    System.out.println("               2*x^2*y^2*t^2 + 2*x^2*y^2*z^2 +");
    System.out.println("               2)                            *exp(x*y*z*t)");
    System.out.println("in the hyper cube   xmin<x<xmax  and  ymin<y<ymax,");
    System.out.println("                    zmin<z<zmax  and  tmin<t<tmax,");
    System.out.println(" ");
    System.out.println("   subject to boundary conditions of the first");
    System.out.println("   kind, Dirichlet.");
    System.out.println("   ub(x,y,z,t) := exp(x*y*z*t)");
    System.out.println(" ");
    System.out.println("xmin="+xmin+", xmax="+xmax+
                       ", ymin="+ymin+", ymax="+ymax);
    System.out.println("zmin="+zmin+", zmax="+zmax+
                       ", tmin="+tmin+", tmax="+tmax);
    System.out.println("nx="+nx+", ny="+ny+", nz="+nz+", nt="+nt);

    hx = (xmax-xmin)/(double)(nx-1);
    if(debug) System.out.println("x grid  hx="+hx);
    for(int i=0; i<nx; i++) // does not need to be uniform 
    {
      xg[i] = xmin+(double)i*hx;
      if(debug) System.out.println("xg["+i+"]="+xg[i]);
    }  
    if(debug) System.out.println(" ");
    hy = (ymax-ymin)/(double)(ny-1);
    if(debug) System.out.println("y grid hy="+hy);
    for(int ii=0; ii<ny; ii++)
    {
      yg[ii] = ymin+(double)ii*hy;
      if(debug) System.out.println("yg["+ii+"]="+yg[ii]);
    }  
    if(debug) System.out.println(" ");
    hz = (zmax-zmin)/(double)(nz-1);
    if(debug) System.out.println("z grid hz="+hz);
    for(int iii=0; iii<nz; iii++)
    {
      zg[iii] = zmin+(double)iii*hz;
      if(debug) System.out.println("zg["+iii+"]="+zg[iii]);
    }  
    if(debug) System.out.println(" ");
    ht = (tmax-tmin)/(double)(nt-1);
    if(debug) System.out.println("t grid  ht="+ht);
    for(int iiii=0; iiii<nt; iiii++)
    {
      tg[iiii] = tmin+(double)iiii*ht;
      if(debug) System.out.println("tg["+iiii+"]="+tg[iiii]);
    }  
    if(debug) System.out.println(" ");

    // put solution in Ua vector 
    for(int i=1; i<nx-1; i++)
    {
      x = xg[i];
      for(int ii=1; ii<ny-1; ii++)
      {
	y = yg[ii];
        for(int iii=1; iii<nz-1; iii++)
        {
          z = zg[iii];
          for(int iiii=1; iiii<nt-1; iiii++)
          {
	    t = tg[iiii];
            Ua[s(i,ii,iii,iiii)] = ub(x,y,z,t);
	  }
	}
      }
    }
    if(debug) System.out.println(" ");

    // initialize k and f 
    for(int i=1; i<nx-1; i++)
    {
      for(int j=1; j<nx-1; j++)
      {
        for(int ii=1; ii<ny-1; ii++)
        {
          for(int jj=1; jj<ny-1; jj++)
          {
            for(int iii=1; iii<nz-1; iii++)
            {
              for(int jjj=1; jjj<nz-1; jjj++)
              {
                for(int iiii=1; iiii<nt-1; iiii++)
                {
                  for(int jjjj=1; jjjj<nt-1; jjjj++)
                  {
	            kg[s(i,ii,iii,iiii)][s(j,jj,jjj,jjjj)] = 0.0;
		  }
		}
	      }
	    }
	  }
        }
      }
    }
    for(int i=1; i<nx-1; i++)
    {
      x = xg[i];
      for(int ii=1; ii<ny-1; ii++)
      {
        y = yg[ii];
        for(int iii=1; iii<nz-1; iii++)
        {
          z = zg[iii];
          for(int iiii=1; iiii<nt-1; iiii++)
          {
            t = tg[iiii];
	    fg[s(i,ii,iii,iiii)] = F(x,y,z,t);
	  }
	}
      }
    }

    dsetup(); // setup discrete coefficients 

    // compute entries in kg (stiffness matrix) discretization matrix 
    System.out.println("compute discretization matrix  ");
    for(int i=1; i<nx-1; i++) // not boundary rows
    {
      for(int ii=1; ii<ny-1; ii++) // not boundary rows
      {
        for(int iii=1; iii<nz-1; iii++) // not boundary rows
        {
          for(int iiii=1; iiii<nt-1; iiii++) // not boundary rows
          {
            // compute discrete derivatives around (i,ii) DOF  
	    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 

    System.out.println("k computed discrete matrix");
    System.out.println("f computed RHS");

    // write_eqn(); // all
    write_eqnij(); // i,j,val non zero

    new simeq(kg, fg, ug);

    write_soln(); // to file x y soln

    check_soln(); // independent of nowing analytic solution

    avgerr = 0.0;
    maxerr = 0.0;
    System.out.println("ug computed, Ua analytic, error  ");
    for(int i=1; i<nx-1; i++)
    {
      for(int ii=1; ii<ny-1; ii++)
      {
        for(int iii=1; iii<nz-1; iii++)
         {
           for(int iiii=1; iiii<nt-1; iiii++)
           {
	     i4 = s(i,ii,iii,iiii);
             err = Math.abs(ug[i4]-Ua[i4]);
             if(err>maxerr) maxerr = err;
             avgerr = avgerr + err;
             System.out.println("ug["+i+","+ii+","+iii+","+iiii+"]="+
                               ug[i4]+", Ua="+Ua[i4]+", err="+(ug[i4]-Ua[i4]));
	   }
	 }
      }
    }
    System.out.println("xmax="+xmax+", ymax="+ymax+", nx="+nx+", ny="+ny);
    System.out.println("zmax="+zmax+", tmax="+tmax+", nz="+nz+", nt="+nt);
    System.out.println("                  maxerr="+maxerr+", avgerr="+
                       avgerr/(double)(nxyzt));
    System.out.println("");
  } // end pde_bihar44e_eq constructor

  int s(int i, int ii, int iii, int iiii) // index 1..nx-1  etc
  {
      return (i-1)*(ny-2)*(nz-2)*(nt-2) + (ii-1)*(nz-2)*(nt-2) +
             (iii-1)*(nt-2) + (iiii-1);
  } // end s

  // PDE problem definition functions  f and ub

  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 +
            2.0*y*y*z*z*t*t + 2.0*x*x*z*z*t*t +
            2.0*x*x*y*y*t*t + 2.0*x*x*y*y*z*z +
            2.0       ) * Math.exp(x*y*z*t);
  }

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

  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 i4;
    double val;

    i4 = s(i,ii,iii,iiii); // equation index 
    if(debug) System.out.println("discretize equation i4="+i4);
    // for each term in each equation 

    // for uxxxx  d^4u/dx^4 
    for(int j=0; j<nx; j++)
    {
      val = cxxxx[i][j];
      if(j==0 || j==nx-1)
        fg[i4] -= 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(int jj=0; jj<ny; jj++)
    {
      val = cyyyy[ii][jj];
      if(jj==0 || jj==ny-1)
        fg[i4] -=  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(int jjj=0; jjj<nz; jjj++)
    {
      val = czzzz[iii][jjj];
      if(jjj==0 || jjj==nz-1)
        fg[i4] -=  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(int jjjj=0; jjjj<nt; jjjj++)
    {
      val = ctttt[iiii][jjjj];
      if(jjjj==0 || jjjj==nt-1)
        fg[i4] -=  val*ub(xg[i],yg[ii],zg[iii],tg[jjjj]);
      else kg[i4][s(i,ii,iii,jjjj)] += val;
    }

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

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

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

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

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

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

  void write_soln()
  {
    double x, y, z, t;
    int dim = 4;
    try
    {
      FileWriter fout = new FileWriter("pde_bihar44e_eq.dat");
      BufferedWriter fileout = new BufferedWriter(fout);
      System.out.println("writing pde_bihar44e_eq.dat");
      fileout.write(dim+" "+nx+" "+ny+" "+nz+" "+nt+" \n");

      for(int i=0; i<nx; i++)
      {
	x = xg[i];
	for(int ii=0; ii<ny; ii++)
	{
	  y = yg[ii];
	  for(int iii=0; iii<nz; iii++)
	  {
	    z = zg[iii];
	    for(int iiii=0; iiii<nt; iiii++)
	    {
	      t = tg[iiii];
	      if(i==0 || i==nx-1 || ii==0 || ii==ny-1 ||
		 iii==0 || iii==nz-1 || iiii==0 || iiii==nt-1)
		   fileout.write(" "+x+" "+y+" "+z+" "+t+" "+
				 ub(x,y,z,t)+" \n");
	      else
		fileout.write(" "+x+" "+y+" "+z+" "+t+" "+
			      ug[s(i,ii,iii,iiii)]+" \n");
	    } // end t
 	  } // end z
	} // end y
      } // end x
      fileout.close();
      System.out.println("finished writing pde_bihar44e_eq.dat");
    }
    catch(FileNotFoundException exception)
    {
      System.out.println("pde_bihar44e_eq.dat  not opened");
    }
    catch(Exception exception) // catches all exceptions
    {
      System.out.println("pde_bihar44e_eq.java output-file error");
    } // end try
  } // end write_soln

  void write_eqn()
  {
    try
    {
      FileWriter fout = new FileWriter("pde_bihar44t_eq.eqn");
      BufferedWriter fileout = new BufferedWriter(fout);
      System.out.println("writing pde_bihar44t_eq.dat");
      fileout.write(nxyzt+" "+(nxyzt+1)+"\n");

      for(int i=0; i<nxyzt; i++)
      {
	for(int j=0; j<nxyzt; j++)
	{
	  fileout.write(kg[i][j]+" \n");
	} // end j
	fileout.write(fg[i]+" \n");
      } // end i
      fileout.close();
      System.out.println("finished writing pde_bihar44t_eq.dat");
    }
    catch(FileNotFoundException exception)
    {
      System.out.println("pde_bihar44t_eq.eqn  not opened");
    }
    catch(Exception exception) // catches all exceptions
    {
      System.out.println("pde_bihar44t_eq.java output-file error");
    } // end try
  } // end write_eqn

  void write_eqnij()
  {
    try
    {
      FileWriter fout = new FileWriter("pde_bihar44t_eqij.eqn");
      BufferedWriter fileout = new BufferedWriter(fout);
      System.out.println("writing pde_bihar44t_eqij.dat");
      fileout.write(nxyzt+" "+(nxyzt+1)+"\n");

      for(int i=0; i<nxyzt; i++)
      {
	for(int j=0; j<nxyzt; j++)
	{
	  if(kg[i][j]!=0.0)
	     fileout.write(i+" "+j+" "+kg[i][j]+" \n");
	} // end j
	fileout.write(i+" "+nxyzt+" "+fg[i]+" \n");
      } // end i
      fileout.close();
      System.out.println("finished writing pde_bihar44t_eqij.dat");
    }
    catch(FileNotFoundException exception)
    {
      System.out.println("pde_bihar44t_eqij.eqn  not opened");
    }
    catch(Exception exception) // catches all exceptions
    {
      System.out.println("pde_bihar44t_eq.java output-file error");
    } // end try
  } // end write_eqnij

  void check_soln()
                  // for all interior points
  {
    double x, y, z, t, err, maxerr, avgerr, rmserr;
    double sumerr, sumsqerr, r_deriv, a_deriv, eval_deriv;
    double t_deriv = 0.0;
    int cnt = 0;
    int i4;
    boolean check = false; // set for cases below

    System.out.println("check_soln values ");
    sumerr = 0.0;
    sumsqerr = 0.0;
    maxerr = 0.0;
    for(int i=1; i<nx-1; i++) // full equation 
    {
      x = xg[i];
      for(int ii=1; ii<ny-1; ii++)
      {
        y = yg[ii];
        for(int iii=1; iii<nz-1; iii++)
        {
          z = zg[iii];
          for(int iiii=1; iiii<nt-1; iiii++)
          {
            t = tg[iiii];
            i4 = s(i,ii,iii,iiii);

            if(debug && i==nx-2 && ii==ny-2 && iii==nz-2 && iiii==nt-2)
               check = true;
            if(check)
	    {
              System.out.println("checking on check against PDE i4="+i4);
              System.out.println("xg[i]="+xg[i]+", yg[ii]="+yg[ii]+
                   ", zg[iii]="+zg[iii]+", tg[iiii]="+tg[iiii]);
	    }

            // PDE 2 u(x,y,z,t) term 
            eval_deriv = 2.0*ug[i4];
            if(check) System.out.println("2.0*ug ="+(2.0*ug[i4])+", err="+
					 (eval_deriv-2.0*ub(x,y,z,t)));
            // t_deriv test derivative based on computed solution
            // a_deriv test derivative based on exact solution

            // PDE uxxxx term
            t_deriv = 0.0;           
            for(int j=1; j<nx-1; j++)
            {
	      t_deriv += cxxxx[i][j]*ug[s(j,ii,iii,iiii)];
              if(check) System.out.println("cxxxx["+i+"]["+j+"]="+
                                           cxxxx[i][j]+", ug[s("+j+
                                  ",ii,iii,iiii)]="+ug[s(j,ii,iii,iiii)]);
            } 
            t_deriv += cxxxx[i][0]*ub(xg[0],y,z,t);
            if(check) System.out.println("cxxxx["+i+"][0]="+
	              cxxxx[i][0]+", ub(xg[0],y,z,t)="+ub(xg[0],y,z,t));
            t_deriv += cxxxx[i][nx-1]*ub(xg[nx-1],y,z,t);
            if(check) System.out.println("cxxxx["+i+"][nx-1]="+
	      cxxxx[i][nx-1]+", ub(xg[nx-1],y,z,t)="+ub(xg[nx-1],y,z,t));
            eval_deriv += t_deriv;
            if(check)
	    {
              a_deriv = 0.0;           
              for(int j=1; j<nx-1; j++)
              {
		a_deriv += cxxxx[i][j]*ub(xg[j],y,z,t);
              } 
              a_deriv += cxxxx[i][0]*ub(xg[0],y,z,t);
              a_deriv += cxxxx[i][nx-1]*ub(xg[nx-1],y,z,t);
	      r_deriv = y*y*y*y*z*z*z*z*t*t*t*t*Math.exp(x*y*z*t);
	      err = r_deriv - t_deriv;
              System.out.println("uxxxx="+r_deriv+", eval="+
                                 t_deriv+", err="+err);
              System.out.println("uxxxx="+r_deriv+", a_deriv="+
                                 a_deriv+", err="+(r_deriv-a_deriv));
	    } // end uxxxx check

            // PDE uyyyy term 
            t_deriv = 0.0;
            for(int jj=1; jj<ny-1; jj++)
            {
	      t_deriv += cyyyy[ii][jj]*ug[s(i,jj,iii,iiii)];
              if(check) System.out.println("cyyyy["+ii+"]["+jj+"]="+
                                           cyyyy[ii][jj]+", ug[s(i,"+jj+
                                  ",iii,iiii)]="+ug[s(i,jj,iii,iiii)]);
            }
            t_deriv += cyyyy[ii][0]*ub(x,yg[0],z,t);
            if(check) System.out.println("cyyyy["+ii+"][0]="+
	     cyyyy[ii][0]+", ub(x,yg[0],z,t)="+ub(x,yg[0],z,t));
            t_deriv += cyyyy[ii][ny-1]*ub(x,yg[ny-1],z,t);
            if(check) System.out.println("cyyyy["+ii+"][ny-1]="+
	     cyyyy[ii][ny-1]+", ub(x,yg[ny-1],z,t)="+ub(x,yg[ny-1],z,t));
            eval_deriv += t_deriv;
            if(check)
	    {
              a_deriv = 0.0;           
              for(int jj=1; jj<ny-1; jj++)
              {
		a_deriv += cyyyy[ii][jj]*ub(x,yg[jj],z,t);
              } 
              a_deriv += cyyyy[ii][0]*ub(x,yg[0],z,t);
              a_deriv += cyyyy[ii][ny-1]*ub(x,yg[ny-1],z,t);
	      r_deriv = x*x*x*x*z*z*z*z*t*t*t*t*Math.exp(x*y*z*t);
	      err = r_deriv - t_deriv;
              System.out.println("uyyyy="+r_deriv+", eval="+
                                 t_deriv+", err="+err);
              System.out.println("uyyyy="+r_deriv+", a_deriv="+
                                 a_deriv+", err="+(r_deriv-a_deriv));
	    } // end uyyyy check

            // PDE uzzzz term 
            t_deriv = 0.0;           
            for(int jjj=1; jjj<nz-1; jjj++)
            {
	      t_deriv += czzzz[iii][jjj]*ug[s(i,ii,jjj,iiii)];
              if(check) System.out.println("czzzz["+iii+"]["+jjj+"]="+
                                           czzzz[iii][jjj]+", ug[s(i,ii,"+jjj+
                                  ",iiii)]="+ug[s(i,ii,jjj,iiii)]);
            } 
            t_deriv += czzzz[iii][0]*ub(x,y,zg[0],t);
            if(check) System.out.println("czzzz["+iii+"][0]="+
	     czzzz[iii][0]+", ub(x,y,zg[0],t)="+ub(x,y,zg[0],t));
            t_deriv += czzzz[iii][nz-1]*ub(x,y,zg[nz-1],t);
            if(check) System.out.println("czzzz["+iii+"][nz-1]="+
	     czzzz[iii][nz-1]+", ub(x,y,zg[nz-1],t)="+ub(x,y,zg[nz-1],t));
            eval_deriv += t_deriv;
            if(check)
	    {
              a_deriv = 0.0;           
              for(int jjj=1; jjj<nz-1; jjj++)
              {
		a_deriv += czzzz[iii][jjj]*ub(x,y,zg[jjj],t);
              } 
              a_deriv += czzzz[iii][0]*ub(x,y,zg[0],t);
              a_deriv += czzzz[iii][nz-1]*ub(x,y,zg[nz-1],t);
	      r_deriv = x*x*x*x*y*y*y*y*t*t*t*t*Math.exp(x*y*z*t);
	      err = r_deriv - t_deriv;
              System.out.println("uzzzz="+r_deriv+", eval="+
                                 t_deriv+", err="+err);
              System.out.println("uzzzz="+r_deriv+", a_deriv="+
                                 a_deriv+", err="+(r_deriv-a_deriv));
	    } // end uzzzz check

            // PDE utttt term 
            t_deriv = 0.0;           
            for(int jjjj=1; jjjj<nt-1; jjjj++)
            {
	      t_deriv += ctttt[iiii][jjjj]*ug[s(i,ii,iii,jjjj)];
              if(check) System.out.println("ctttt["+iiii+"]["+jjjj+"]="+
                                       ctttt[iiii][jjjj]+", ug[s(i,ii,iii,"+
                                  jjjj+")]="+ug[s(i,ii,iii,jjjj)]);
            }
            t_deriv += ctttt[iiii][0]*ub(x,y,z,tg[0]);
            if(check) System.out.println("ctttt["+iiii+"][0]="+
			ctttt[iiii][0]+", ub(x,y,z,tg[0])="+ub(x,y,z,tg[0]));
            t_deriv += ctttt[iiii][nt-1]*ub(x,y,z,tg[nt-1]);
            if(check) System.out.println("ctttt["+iiii+"][nt-1]="+
	      ctttt[iiii][nt-1]+", ub(x,y,z,tg[nz-1])="+ub(x,y,z,tg[nt-1]));
            eval_deriv += t_deriv;
            if(check)
	    {
              a_deriv = 0.0;           
              for(int jjjj=1; jjjj<nt-1; jjjj++)
              {
		a_deriv += ctttt[iiii][jjjj]*ub(x,y,z,tg[jjjj]);
              } 
              a_deriv += ctttt[iiii][0]*ub(x,y,z,tg[0]);
              a_deriv += ctttt[iiii][nt-1]*ub(x,y,z,tg[nt-1]);
	      r_deriv = x*x*x*x*y*y*y*y*z*z*z*z*Math.exp(x*y*z*t);
	      err = r_deriv - t_deriv;
              System.out.println("utttt="+r_deriv+", eval="+
                                 t_deriv+", err="+err);
              System.out.println("utttt="+r_deriv+", a_deriv="+
                                 a_deriv+", err="+(r_deriv-a_deriv));
	    } // end utttt check

            // PDE  2*uxx term 
            t_deriv = 0.0;
            for(int j=1; j<nx-1; j++)
            {
	      t_deriv += cxx[i][j]*ug[s(j,ii,iii,iiii)];
              if(check) System.out.println("cxx["+i+"]["+j+"]="+
                                           cxx[i][j]+", ug[s("+j+
                                  ",ii,iii,iiii)]="+ug[s(j,ii,iii,iiii)]);
            } // end j 
            t_deriv += cxx[i][0]*ub(xg[0],y,z,t);
            if(check) System.out.println("cxx["+i+"][0]="+
	              cxx[i][0]+", ub(xg[0],y,z,t)="+ub(xg[0],y,z,t));
            t_deriv += cxx[i][nx-1]*ub(xg[nx-1],y,z,t);
            if(check) System.out.println("cxx["+i+"][nx-1]="+
	      cxx[i][nx-1]+", ub(xg[nx-1],y,z,t)="+ub(xg[nx-1],y,z,t));
            eval_deriv += 2.0*t_deriv;
            if(check)
	    {
              a_deriv = 0.0;           
              for(int j=1; j<nx-1; j++)
              {
		a_deriv += cxx[i][j]*ub(xg[j],y,z,t);
              } 
              a_deriv += cxx[i][0]*ub(xg[0],y,z,t);
              a_deriv += cxx[i][nx-1]*ub(xg[nx-1],y,z,t);
	      r_deriv = y*y*z*z*t*t*Math.exp(x*y*z*t);
	      err = r_deriv - t_deriv;
              System.out.println("uxx="+r_deriv+", eval="+
                                 t_deriv+", err="+err);
              System.out.println("uxx="+r_deriv+", a_deriv="+
                                 a_deriv+", err="+(r_deriv-a_deriv));
	    } // end uxx check

            // PDE  2*uyy term 
            t_deriv = 0.0;           
            for(int jj=1; jj<ny-1; jj++)
            {
              t_deriv += cyy[ii][jj]*ug[s(i,jj,iii,iiii)];
            } // end jj 
            t_deriv += cyy[ii][0]*ub(x,yg[0],z,t);
            t_deriv += cyy[ii][ny-1]*ub(x,yg[ny-1],z,t);
            eval_deriv += 2.0*t_deriv;
            if(check)
	    {
              a_deriv = 0.0;           
              for(int jj=1; jj<ny-1; jj++)
              {
		  a_deriv += cyy[ii][jj]*ub(x,yg[jj],z,t);
              } 
              a_deriv += cyy[ii][0]*ub(x,yg[0],z,t);
              a_deriv += cyy[ii][ny-1]*ub(x,yg[ny-1],z,t);
	      r_deriv = x*x*z*z*t*t*Math.exp(x*y*z*t);
	      err = r_deriv - t_deriv;
              System.out.println("uyy="+r_deriv+", eval="+
                                 t_deriv+", err="+err);
              System.out.println("uyy="+r_deriv+", a_deriv="+
                                 a_deriv+", err="+(r_deriv-a_deriv));
	    } // end uyy check

            // PDE  2*uzz term 
            t_deriv = 0.0;           
            for(int jjj=1; jjj<nz-1; jjj++)
            {
              t_deriv += czz[iii][jjj]*ug[s(i,ii,jjj,iiii)];
            } // end jjj 
            t_deriv += czz[iii][0]*ub(x,y,zg[0],t);
            t_deriv += czz[iii][nz-1]*ub(x,y,zg[nz-1],t);
            eval_deriv += 2.0*t_deriv;
            if(check)
	    {
              a_deriv = 0.0;           
              for(int jjj=1; jjj<nz-1; jjj++)
              {
		  a_deriv += czz[iii][jjj]*ub(x,y,zg[jjj],t);
              } 
              a_deriv += czz[iii][0]*ub(x,y,zg[0],t);
              a_deriv += czz[iii][nz-1]*ub(x,y,zg[nz-1],t);
	      r_deriv = x*x*y*y*t*t*Math.exp(x*y*z*t);
	      err = r_deriv - t_deriv;
              System.out.println("uzz="+r_deriv+", eval="+
                                 t_deriv+", err="+err);
              System.out.println("uzz="+r_deriv+", a_deriv="+
                                 a_deriv+", err="+(r_deriv-a_deriv));
	    } // end uzz check

            // PDE  2*utt term 
            t_deriv = 0.0;           
            for(int jjjj=1; jjjj<nt-1; jjjj++)
            {
              t_deriv += ctt[iiii][jjjj]*ug[s(i,ii,iii,jjjj)];
            } // end jjjj 
            t_deriv += ctt[iiii][0]*ub(x,y,z,tg[0]);
            t_deriv += ctt[iiii][nt-1]*ub(x,y,z,tg[nt-1]);
            eval_deriv += 2.0*t_deriv;
            if(check)
	    {
              a_deriv = 0.0;           
              for(int jjjj=1; jjjj<nt-1; jjjj++)
              {
		a_deriv += ctt[iiii][jjjj]*ub(x,y,z,tg[jjjj]);
              } 
              a_deriv += ctt[iiii][0]*ub(x,y,z,tg[0]);
              a_deriv += ctt[iiii][nt-1]*ub(x,y,z,tg[nt-1]);
	      r_deriv = x*x*y*y*z*z*Math.exp(x*y*z*t);
	      err = r_deriv - t_deriv;
              System.out.println("utt="+r_deriv+", eval="+
                                 t_deriv+", err="+err);
              System.out.println("utt="+r_deriv+", a_deriv="+
                                 a_deriv+", err="+(r_deriv-a_deriv));
	    } // end utt check

            err = eval_deriv-F(x,y,z,t);
            if(check)
	    {
	      System.out.println("eval_deriv="+eval_deriv+", F(x,y,z,t)="+F(x,y,z,t));
	      System.out.println("err="+err);
	    }
            err = Math.abs(err);
            sumerr += err;
            sumsqerr += err*err;
            maxerr = Math.max(err,maxerr);
            cnt++;
	  } // end t 
        } // end z  
      } // end y 
    } // end x 
    rmserr = Math.sqrt(sumsqerr/(double)cnt);
    avgerr = sumerr/(double)cnt;
    System.out.println("maxerr="+maxerr+", rmserr="+rmserr+", avgerr="+avgerr);  
  } // end check_soln 

  void rderiv(int order, int npts, int point, double h, double c[])
  {
    nderiv D = new nderiv();
    int a[] = new int[npts]; // basic discretization from nderiv
    int b[] = new int[1];    // out parameter, denominator

    D.deriv(order,npts,point,a,b);
    for(int i=0; i<npts; i++)
    {
      c[i] = (double)a[i]/((double)b[0]);
      for(int j=0; j<order; j++) c[i] = c[i]/h;
    }
  } // end rderiv

  public static void main (String[] args)
  {
    pde_bihar44e_eq F2 = new pde_bihar44e_eq();
  
    F2.run(1.0, 1.0, 1.0, 1.0, 5, 5, 5, 5); // xmax, ymax, zmax, tmax,
                                            // nx, ny, nz, nt  
    F2.run(1.0, 1.0, 1.0, 1.0, 7, 5, 5, 5); // xmax, ymax, zmax, tmax,
                                            // nx, ny, nz, nt  
    F2.run(1.0, 1.0, 1.0, 1.0, 5, 7, 5, 5); // xmax, ymax, zmax, tmax,
                                            // nx, ny, nz, nt  
    F2.run(1.0, 1.0, 1.0, 1.0, 5, 5, 7, 5); // xmax, ymax, zmax, tmax,
                                            // nx, ny, nz, nt  
    F2.run(1.0, 1.0, 1.0, 1.0, 5, 5, 5, 7); // xmax, ymax, zmax, tmax,
                                            // nx, ny, nz, nt  
    F2.run(1.0, 1.0, 1.0, 1.0, 7, 7, 7, 7); // xmax, ymax, zmax, tmax,
                                            // nx, ny, nz, nt  
    try
    {
      Runtime runtime = Runtime.getRuntime();
      Process proc = runtime.exec("java -cp . plot4d < pde_bihar44e_eq.plot");
    }
    catch(Exception exception) // catches all exceptions
    {
      System.out.println("pde_bihar44e_eq.java runtime error");
    } // end try
  }
} // end class pde_bihar44e_eq