// pde47hn_eq.java  Biharmonic PDE in seven 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[i3], t[i4], u[i5],
//                        v[i6], w[i7]
//
// replace continuous derivatives with discrete derivatives
// solve for U(x,y,z,t,u,v,w) numerical approximation U
// A * U = F, solve simultaneous equations for U
//
 
import java.text.*;
import java.io.*;
// uses rderiv.java
// uses simeq.java

class pde47hn_eq
{
  int nx = 5;
  int ny = 5;
  int nz = 5;
  int nt = 5;
  int nu = 5;
  int nv = 5;
  int nw = 5;
  int nxyztuvw = (nx-2)*(ny-2)*(nz-2)*(nt-2)*(nu-2)*(nv-2)*(nw-2); // DOF
  double A[] = new double[nxyztuvw*(nxyztuvw+1)]; // last column is RHS 
  double U[] = new double[nxyztuvw]; // computed solution, no boundary
  double xg[] = new double[nx]; // coordinates
  double yg[] = new double[ny];
  double zg[] = new double[nz];
  double tg[] = new double[nt];
  double ug[] = new double[nu];
  double vg[] = new double[nv];
  double wg[] = new double[nw];
  double Ua[] = new double[nxyztuvw]; // analytic solution for checking, no bou
  double xmin = 0.0; // problem parameters 
  double xmax = 0.5;
  double ymin = 0.0;
  double ymax = 0.5;
  double zmin = 0.0;
  double zmax = 0.5;
  double tmin = 0.0;
  double tmax = 0.5;
  double umin = 0.0;
  double umax = 0.5;
  double vmin = 0.0;
  double vmax = 0.5;
  double wmin = 0.0;
  double wmax = 0.5;
  double hx, hy, hz, ht, hu, hv, hw;

  pde47hn_eq()
  {
    double x, y, z, t, u, v, w;
    int    i,ii,i3,i4,i5,i6,i7,
           j,jj,j3,j4,j5,j6,j7;
    int k, cs;
    // coefficients of terms in first equation 
    // derivative point arrays, dynamic in this version 
    double cxx[] = new double[nx];
    double cyy[] = new double[ny];
    double czz[] = new double[nz];
    double ctt[] = new double[nt];
    double cuu[] = new double[nu];
    double cvv[] = new double[nv]; 
    double cww[] = new double[nw]; 
    double cxxxx[] = new double[nx];
    double cyyyy[] = new double[ny];
    double czzzz[] = new double[nz];
    double ctttt[] = new double[nt];
    double cuuuu[] = new double[nu];
    double cvvvv[] = new double[nv];
    double cwwww[] = new double[nw];
    double val, err, avgerr, maxerr;
    double time_start;
    double time_now;
    double time_last;

    System.out.println("pde47hn_eq.java running ");
    System.out.println("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) +");
    System.out.println("      Utttt(x,y,z,t,u,v,w) + Uuuuu(x,y,z,t,u,v,w) + Uvvvv(x,y,z,t,u,v,w) +");
    System.out.println("      Uwwww(x,y,z,t,u,v,w) +");
    System.out.println("      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) +");
    System.out.println("      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) +");
    System.out.println("      2 Uww(x,y,z,t,u,v,w) +");
    System.out.println("      8   U(x,y,z,t,u,v,w) = f(x,y,z,t,u,v,w) ");
    System.out.println(" "); 
    System.out.println("xmin<=x<=xmax ymin<=y<=ymax zmin<=z<=zmax ");
    System.out.println("tmin<=t<=tmax umin<=u<=umax vmin<=v<=vmax ");
    System.out.println("wmin<=w<=wmax    Boundaries ");
    System.out.println("Analytic solution U(x,y,z,t,u,v,w) =  sin(x+y+z+t+u+v+w) ");
    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("umin="+umin+", umax="+umax+
                       ", vmin="+vmin+", vmax="+vmax);
    System.out.println("wmin="+wmin+", wmax="+wmax);
    System.out.println("nx="+nx+", ny="+ny+", nz="+nz+", nt="+nt+
                     ", nu="+nu+", nv="+nv+", nw="+nw);
    time_start = (double)System.currentTimeMillis()/1000.0;
    System.out.println("time start="+time_start+" seconds");
    System.out.println("DOF, degrees of freedom, nxyztuvw="+nxyztuvw);

    time_last = time_start;

    System.out.println(" ");
    System.out.println("x grid and analytic solution at minimums ");
    hx = (xmax-xmin)/(double)(nx-1);
    i = 0;
    xg[i] = xmin;
    System.out.println("i="+i+", Ua("+xg[i]+","+ymin+","+zmin+
                       ","+tmin+","+umin+","+vmin+","+wmin+")="+
                       Ub(xg[i],ymin,zmin,tmin,umin,vmin,wmin)); 
    for(i=1; i<nx; i++)
    {
      xg[i] = xmin+(double)(i)*hx;
      System.out.println("i="+i+", Ua("+xg[i]+","+ymin+","+zmin+
                         ","+tmin+","+umin+","+umin+","+wmin+")="+
                         Ub(xg[i],ymin,zmin,tmin,umin,vmin,wmin)); 
    } // i  x   
    System.out.println(" ");
    System.out.println("y grid and analytic solution at minimums ");
    hy = (ymax-ymin)/(double)(ny-1);
    ii = 0;
    yg[ii] = ymin;
    System.out.println("ii="+ii+", Ua("+xmin+","+yg[ii]+","+zmin+
                     ","+tmin+","+umin+","+vmin+","+wmin+")="+
		       Ub(xmin,yg[ii],zmin,tmin,umin,vmin,wmin)); 
    for(ii=1; ii<ny; ii++)
    {
      yg[ii] = ymin+(double)(ii)*hy;
      System.out.println("ii="+ii+", Ua("+xmin+","+yg[ii]+","+zmin+
                         ","+tmin+","+umin+","+vmin+","+wmin+")="+
                         Ub(xmin,yg[ii],zmin,tmin,umin,vmin,wmin)); 
    } // ii  y   
    System.out.println(" ");
    System.out.println("z grid and analytic solution at minimums ");
    hz = (zmax-zmin)/(double)(nz-1);
    i3 = 0;
    zg[i3] = zmin;
    System.out.println("i3="+i3+", Ua("+xmin+","+ymin+","+zg[i3]+
                     ","+tmin+","+umin+","+vmin+","+wmin+")="+
		       Ub(xmin,ymin,zg[i3],tmin,umin,vmin,wmin)); 
    for(i3=1; i3<nz; i3++)
    {
      zg[i3] = zmin+(double)(i3)*hz;
      System.out.println("i3="+i3+", Ua("+xmin+","+ymin+","+zg[i3]+
                         ","+tmin+","+umin+","+vmin+","+wmin+")="+
                         Ub(xmin,ymin,zg[i3],tmin,umin,vmin,wmin)); 
    } // i3  z   
    System.out.println(" ");
    System.out.println("t grid and analytic solution at minimums ");
    ht = (tmax-tmin)/(double)(nt-1);
    i4 = 0;
    tg[i4] = tmin;
    System.out.println("i4="+i4+", Ua("+xmin+","+ymin+","+zmin+
                       ","+tg[i4]+","+umin+","+vmin+","+wmin+")="+
                       Ub(xmin,ymin,zmin,tg[i4],umin,vmin,wmin)); 
    for(i4=1; i4<nt; i4++)
    {
      tg[i4] = tmin+(double)(i4)*ht;
      System.out.println("i4="+i4+", Ua("+xmin+","+ymin+","+zmin+
                         ","+tg[i4]+","+umin+","+vmin+","+wmin+")="+
                         Ub(xmin,ymin,zmin,tg[i4],umin,vmin,wmin)); 
    } // i4  t 
    System.out.println(" ");
    System.out.println("u grid and analytic solution at minimums ");
    hu = (umax-umin)/(double)(nu-1);
    i5 = 0;
    ug[i5] = umin;
    System.out.println("i5="+i5+", Ua("+xmin+","+ymin+","+zmin+
                       ","+tmin+","+ug[i5]+","+vmin+","+wmin+")="+
                       Ub(xmin,ymin,zmin,tmin,ug[i5],vmin,wmin)); 
    for(i5=1; i5<nu; i5++)
    {
      ug[i5] = umin+(double)(i5)*hu;
      System.out.println("i5="+i5+", Ua("+xmin+","+ymin+","+zmin+
                         ","+tmin+","+ug[i5]+","+vmin+","+wmin+")="+
                         Ub(xmin,ymin,zmin,tmin,ug[i5],vmin,wmin)); 
    } // i5  u 
    System.out.println(" ");
    System.out.println("v grid and analytic solution at minimums ");
    hv = (vmax-vmin)/(double)(nv-1);
    i6 = 0;
    vg[i6] = vmin;
    System.out.println("i6="+i6+", Ua("+xmin+","+ymin+","+zmin+
                     ","+tmin+","+umin+","+vg[i6]+","+wmin+")="+
		     Ub(xmin,ymin,zmin,tmin,umin,vg[i6],wmin)); 
    for(i6=1; i6<nv; i6++)
    {
      vg[i6] = vmin+(double)(i6)*hv;
      System.out.println("i6="+i6+", Ua("+xmin+","+ymin+","+zmin+
                         ","+tmin+","+umin+","+vg[i6]+","+wmin+")="+
                         Ub(xmin,ymin,zmin,tmin,umin,vg[i6],wmin)); 
    } // i6  v 
    System.out.println(" ");
    System.out.println("w grid and analytic solution at minimums ");
hw = (wmax-wmin)/(double)(nw-1);
    i7 = 0;
    wg[i7] = wmin;
    System.out.println("i7="+i7+", Ua("+xmin+","+ymin+","+zmin+
                     ","+tmin+","+umin+","+vmin+","+wg[i7]+")="+
		       Ub(xmin,ymin,zmin,tmin,umin,vmin,wg[i7])); 
    for(i7=1; i7<nw; i7++)
    {
      wg[i7] = wmin+(double)(i7)*hw;
      System.out.println("i7="+i7+", Ua("+xmin+","+ymin+","+zmin+
                       ","+tmin+","+umin+","+vmin+","+wg[i7]+")="+
		       Ub(xmin,ymin,zmin,tmin,umin,vmin,wg[i7])); 
    } // i7  w 
    System.out.println(" ");

    System.out.println("check rderiv(4,nv,0,hv,cvvvv ");
    new rderiv(4,nv,0,hv,cvvvv);
    for(i=0; i<nv; i++) System.out.println("cvvvv["+i+"]="+cvvvv[i]);
    System.out.println("check rderiv(4,nv,nv-1,hv,cvvvv ");
    new rderiv(4,nv,nv-1,hv,cvvvv);
    for(i=0; i<nv; i++) System.out.println("cvvvv["+i+"]="+cvvvv[i]);
    System.out.println("");

  System.out.println("put solution in Ua vector ");
  // 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(i3=1; i3<nz-1; i3++)
      {
        z = zg[i3];
        for(i4=1; i4<nt-1; i4++)
        {
          t = tg[i4];
          for(i5=1; i5<nu-1; i5++)
          {
            u = ug[i5];
            for(i6=1; i6<nv-1; i6++)
            {
              v = vg[i6];
              for(i7=1; i7<nw-1; i7++)
              {
                w = wg[i7];
                Ua[s(i,ii,i3,i4,i5,i6,i7)] = 
                   Ub(x,y,z,t,u,v,w);
              } // i7 
            } // i6 
          } // i5 
        } // i4 
      } // i3 
    } // ii 
  } // i 
  time_now = (double)System.currentTimeMillis()/1000.0;
  System.out.println("time now ="+(time_now-time_start)+
                     " seconds, time for previous section ="+
                     (time_now-time_last)+" seconds");
  time_last = time_now;
  System.out.println("");

  System.out.println("initialize A ");
  // initialize A 
  for(i=1; i<nx-1; i++)
  {
    for(ii=1; ii<ny-1; ii++)
    {
      for(i3=1; i3<nz-1; i3++)
      {
        for(i4=1; i4<nt-1; i4++)
        {
          for(i5=1; i5<nu-1; i5++)
          {
            for(i6=1; i6<nv-1; i6++)
            {
              for(i7=1; i7<nw-1; i7++)
              {
                for(j=1; j<nx-1; j++)
                {
                  for(jj=1; jj<ny-1; jj++)
                  {
                    for(j3=1; j3<nz-1; j3++)
                    {
                      for(j4=1; j4<nt-1; j4++)
                      {
                        for(j5=1; j5<nu-1; j5++)
                        {
                          for(j6=1; j6<nv-1; j6++)
                          {
                            for(j7=1; j7<nw-1; j7++)
                            {
			      A[ss(i,ii,i3,i4,i5,i6,i7,
				   j,jj,j3,j4,j5,j6,j7)] = 0.0;
	                    } // j7 
	                  } // j6 
	                } // j5 
	              } // j4 
	            } // j3 
	          } // jj 
                } // j 
              } // i7 
            } // i6 
          } // i5 
        } // i4 
      } // i3 
    } // ii 
  } // i 

  time_now = (double)System.currentTimeMillis()/1000.0;
  System.out.println("time now ="+(time_now-time_start)+
                     " seconds, time for previous section ="+
                     (time_now-time_last)+" seconds");
  time_last = time_now;
  System.out.println("");

  // compute entries in A matrix, inside boundary 
  System.out.println("compute A matrix ");
  cs = nxyztuvw; // constant RHS column 
  val = 0.0;
  for(i=1; i<nx-1; i++)
    {
    for(ii=1; ii<ny-1; ii++)
    {
      for(i3=1; i3<nz-1; i3++)
      {
        for(i4=1; i4<nt-1; i4++)
        {
          for(i5=1; i5<nu-1; i5++)
          {
            for(i6=1; i6<nv-1; i6++)
            { 
              for(i7=1; i7<nw-1; i7++)
              { 
		  k = s(i,ii,i3,i4,i5,i6,i7); // row index 

                // for each term in first equation 

	        // for d^4U/dx^4 
                new rderiv(4,nx,i,hx,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[i3],tg[i4],ug[i5],
                         vg[i6],wg[i7]);
		  else A[sk(k,j,ii,i3,i4,i5,i6,i7)] += val;
	        }

	        // for d^4U/dy^4 
                new rderiv(4,ny,ii,hy,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[i3],tg[i4],ug[i5],
                         vg[i6],wg[i7]);
                  else A[sk(k,i,jj,i3,i4,i5,i6,i7)] += val;
	        }

	        // for d^4U/dz^4 
                new rderiv(4,nz,i3,hz,czzzz);
                for(j3=0; j3<nz; j3++)
                {
                  val = czzzz[j3];
                  if(j3==0 || j3==nz-1)
		    A[k*(nxyztuvw+1)+cs] -= val*
		      Ub(xg[i],yg[ii],zg[j3],tg[i4],ug[i5],
                         tg[i6],wg[i7]);
                  else A[sk(k,i,ii,j3,i4,i5,i6,i7)] += val;
	        }

	        // for d^4U/dt^4 
                new rderiv(4,nt,i4,ht,ctttt);
                for(j4=0; j4<nt; j4++)
                {
                  val = ctttt[j4];
                  if(j4==0 || j4==nt-1)
	            A[k*(nxyztuvw+1)+cs] -=  val*
		      Ub(xg[i],yg[ii],zg[i3],tg[j4],ug[i5],
                         vg[i6],wg[i7]);
                  else A[sk(k,i,ii,i3,j4,i5,i6,i7)] += val;
	        }

	        // for d^4U/du^4 
                new rderiv(4,nu,i5,hu,cuuuu);
                for(j5=0; j5<nu; j5++)
                {
                  val = cuuuu[j5];
                  if(j5==0 || j5==nu-1)
	            A[k*(nxyztuvw+1)+cs] -= val*
		      Ub(xg[i],yg[ii],zg[i3],tg[i4],ug[j5],
                         vg[i6],wg[i7]);
                  else A[sk(k,i,ii,i3,i4,j5,i6,i7)] += val;
	        }

	        // for d^4U/dv^4 
                new rderiv(4,nv,i6,hv,cvvvv);
                for(j6=0; j6<nv; j6++)
                {
                  val = cvvvv[j6];
                  if(j6==0 || j6==nv-1)
	            A[k*(nxyztuvw+1)+cs] -= val*
		      Ub(xg[i],yg[ii],zg[i3],tg[i4],ug[i5],
                         vg[j6],wg[i7]);
                  else A[sk(k,i,ii,i3,i4,i5,j6,i7)] += val;
	        }

	        // for d^4U/dw^4 
                new rderiv(4,nw,i7,hw,cwwww);
                for(j7=0; j7<nw; j7++)
                {
                  val = cwwww[j7];
                  if(j7==0 || j7==nw-1)
	            A[k*(nxyztuvw+1)+cs] -= val*
		      Ub(xg[i],yg[ii],zg[i3],tg[i4],ug[i5],
                         vg[i6],wg[j7]);
                  else A[sk(k,i,ii,i3,i4,i5,i6,j7)] += val;
	        }


	        // for 2 d^2U/dx^2 
                new rderiv(2,nx,ii,hx,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[j],yg[ii],zg[i3],tg[i4],ug[i5],
                         vg[i6],wg[i7]);
                  else A[sk(k,j,ii,i3,i4,i5,i6,i7)] += val;
	        }

	        // for 2 d^2U/dy^2 
                new rderiv(2,ny,ii,hy,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[i3],tg[i4],ug[i5],
                         vg[i6],wg[i7]);
                  else A[sk(k,i,jj,i3,i4,i5,i6,i7)] += val;
	        }

	        // for 2 d^2U/dz^2 
                new rderiv(2,nz,i3,hz,czz);
                for(j3=0; j3<nz; j3++)
                {
                  val = 2.0*czz[j3];
                  if(j3==0 || j3==nz-1)
	            A[k*(nxyztuvw+1)+cs] -= val*
		      Ub(xg[i],yg[ii],zg[j3],tg[i4],ug[i5],
                      vg[i6],wg[i7]);
                  else A[sk(k,i,ii,j3,i4,i5,i6,i7)] += val;
	        }

	        // for 2 d^2U/dt^2 
                new rderiv(2,nt,i4,ht,ctt);
                for(j4=0; j4<nt; j4++)
                {
                  val = 2.0*ctt[j4];
                  if(j4==0 || j4==nt-1)
	            A[k*(nxyztuvw+1)+cs] -=  val*
		      Ub(xg[i],yg[ii],zg[i3],tg[j4],ug[i5],
                         vg[i6],wg[i7]);
                  else A[sk(k,i,ii,i3,j4,i5,i6,i7)] += val;
	        }

	        // for 2 d^2U/du^2 
                new rderiv(2,nu,i5,hu,cuu);
                for(j5=0; j5<nu; j5++)
                { 
                  val = 2.0*cuu[j5];
                  if(j5==0 || j5==nu-1)
	            A[k*(nxyztuvw+1)+cs] -= val*
		      Ub(xg[i],yg[ii],zg[i3],tg[i4],ug[j5],
                         vg[i6],wg[i7]);
                  else A[sk(k,i,ii,i3,i4,j5,i6,i7)] += val;
	        }

	        // for 2 d^2U/dv^2 
                new rderiv(2,nv,i6,hv,cvv);
                for(j6=0; j6<nv; j6++)
                {
                  val = 2.0*cvv[j6];
                  if(j6==0 || j6==nv-1)
	            A[k*(nxyztuvw+1)+cs] -= val*
		      Ub(xg[i],yg[ii],zg[i3],tg[i4],ug[i5],
                         vg[j6],wg[i7]);
                  else A[sk(k,i,ii,i3,i4,i5,j6,i7)] += val;
	        }

	        // for 2 d^2U/dw^2 
                new rderiv(2,nw,i7,hw,cww);
                for(j7=0; j7<nw; j7++)
                {
                  val = 2.0*cww[j7];
                  if(j7==0 || j7==nw-1)
	            A[k*(nxyztuvw+1)+cs] -= val*
		      Ub(xg[i],yg[ii],zg[i3],tg[i4],ug[i5],
                         vg[i6],wg[j7]);
                  else A[sk(k,i,ii,i3,i4,i5,i6,j7)] += val;
	        }

	        // for 8 u(x,y,z,t,u,v,w) 
                A[sk(k,i,ii,i3,i4,i5,i6,i7)] += 8.0;

                // for f(x,y,z,t,u,v,w) RHS 
	        A[k*(nxyztuvw+1)+cs] += f(xg[i],yg[ii],zg[i3],
                                    tg[i4],ug[i5],vg[i6],wg[i7]);

                // end terms 

	      } // end i7 loop 
	    } // end i6 loop 
	  } // end i5 loop 
	} // end i4 loop 
      } // end i3 loop 
    } // end ii loop 
  } // end i loop 

    System.out.println("A computed stiffness matrix ");
    time_now = (double)System.currentTimeMillis()/1000.0;
    System.out.println("time now = "+(time_now-time_start)+
                       " seconds, time for previous section ="+
                       (time_now-time_last)+" seconds");
    time_last = time_now;
    System.out.println("");

    System.out.println("solve "+nxyztuvw+" equations in "+nxyztuvw+" unknowns ");

    new simeq(nxyztuvw, A, U); // A destroyed 

    System.out.println("equations solved ");
    time_now = (double)System.currentTimeMillis()/1000.0;
    System.out.println("time now = "+(time_now-time_start)+
                       " seconds, time for previous section ="+
                       (time_now-time_last)+" seconds");
    time_last = time_now;
    System.out.println("");
  
    System.out.println("check PDE against computed solution");
    check_soln(U);
    write_soln(U);
    System.out.println("");

    avgerr = 0.0;
    maxerr = 0.0;
    System.out.println("          U computed, Ua analytic, error ");
    for(i=1; i<nx-1; i++)
    {
      for(ii=1; ii<ny-1; ii++)
      {
        for(i3=1; i3<nz-1; i3++)
        {
          for(i4=1; i4<nt-1; i4++)
          {
            for(i5=1; i5<nu-1; i5++)
            {
              for(i6=1; i6<nv-1; i6++)
              {
                for(i7=1; i7<nw-1; i7++)
                {
		    err = Math.abs(U[s(i,ii,i3,i4,i5,i6,i7)]-
				   Ua[s(i,ii,i3,i4,i5,i6,i7)]);
                if(err>maxerr) maxerr = err;
                avgerr = avgerr + err;
                System.out.println("ug["+i+","+ii+","+i3+","+i4+","+
                                   i5+","+i6+","+i7+"]="+
                                   U[s(i,ii,i3,i4,i5,i6,i7)]+
		                   ", Ua="+Ua[s(i,ii,i3,i4,i5,i6,i7)]+
                                   ", err="+err);
                } // i7 
              } // i6 
            } // i5 
          } // i4 
        } // i3 
      } // ii 
    } // i 
    time_now = (double)System.currentTimeMillis()/1000.0;
    System.out.println("total CPU time="+(time_now-time_start)+" seconds");
    System.out.println("");

    System.out.println("                                  maxerr="+maxerr+
                       ", avgerr="+(avgerr/(double)nxyztuvw));
    System.out.println(" ");
    System.out.println("end pde47hn_eq.java ");
  } // end pde47hn_eq

  // utility functions
  int s(int i, int ii, int i3, int i4, int i5,
        int i6, int i7)
       // for x,y,z,t,u,v,w  in loops  for(int i=0; i<nx; i++) i!=0, i!=ny-1 
       // etc., no boundary
  {
    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) +
        (i3-1)*(nt-2)*(nu-2)*(nv-2)*(nw-2) + 
        (i4-1)*(nu-2)*(nv-2)*(nw-2) +
        (i5-1)*(nv-2)*(nw-2) + 
        (i6-1)*(nw-2) + 
        (i7-1);
    if(k<0 || k>=nxyztuvw) System.out.println("bad s("+i+","+ii+","+i3+
		","+i4+","+i5+","+i6+","+i7+") ="+k);  
    return k;
  } // end s 

  int sk(int k, int i, int ii, int i3, int i4, int i5,
         int i6, int i7)   // 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) +
           (i3-1)*(nt-2)*(nu-2)*(nv-2)*(nw-2) + 
	   (i4-1)*(nu-2)*(nv-2)*(nw-2) +
	   (i5-1)*(nv-2)*(nw-2) + 
           (i6-1)*(nw-2) + 
           (i7-1);
  } // end sk 

  int ss(int i, int ii, int i3, int i4, int i5, int i6, int i7,
         int j, int jj, int j3, int j4, int j5, int j6, int j7)
  {
    return s(i,ii,i3,i4,i5,i6,i7)*(nxyztuvw+1) + 
	   s(j,jj,j3,j4,j5,j6,j7);
  } // end ss 

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

  double f(double x, double y, double z, double t, double u,
            double v, double w)     //analytic solution and boundaries
  {
    return Math.sin(x+y+z+t+u+v+w);
  } // end f

  double eval_derivative(int xord, int yord, int zord, int tord, int uord,
                         int vord, int word, int i, int ii, int i3, int i4,
                         int i5, int i6, int i7, double U[])
  {
    int j, jj, j3, j4, j5, j6, j7;
    double cx[] = new double[nx];
    double cy[] = new double[ny];
    double cz[] = new double[nz];
    double ct[] = new double[nt];
    double cu[] = new double[nu];
    double cv[] = new double[nv];
    double cw[] = new double[nw];
    double discrete;
    double x, y, z, t, u, v, w;

    new rderiv(xord, nx, i, hx, cx);
    x = xg[i];
    new rderiv(yord, ny, ii, hy, cy);
    y = yg[ii];
    new rderiv(zord, nz, i3, hz, cz);
    z = zg[i3];
    new rderiv(tord, nt, i4, ht, ct);
    t = tg[i4];
    new rderiv(uord, nu, i5, hu, cu);
    u = ug[i5];
    new rderiv(vord, nv, i6, hv, cv);
    v = vg[i6];
    new rderiv(word, nw, i7, hw, cw);
    w = wg[i6];
    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,i3,i4,i5,i6,i7)];
    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,i3,i4,i5,i6,i7)];
    else if(xord==0 && yord==0 && zord!=0 && tord==0 && uord==0 &&
            vord==0 && word==0)
      for(j3=0; j3<nz; j3++)
	  if(j3==0 || j3==nz-1) discrete += cz[j3]*Ub(x,y,zg[j3],t,u,v,w);
	  else discrete += cz[j3]*U[s(i,ii,j3,i4,i5,i6,i6)];
    else if(xord==0 && yord==0 && zord==0 && tord!=0 && uord==0 &&
            vord==0 && word==0)
      for(j4=0; j4<nt; j4++)
	  if(j4==0 || j4==nt-1) discrete += ct[j4]*
                                                Ub(x,y,z,tg[j4],u,v,w);
	  else discrete += ct[j4]*U[s(i,ii,i3,j4,i5,i6,i7)];
    else if(xord==0 && yord==0 && zord==0 && tord==0 && uord!=0 &&
            vord==0 && word==0)
      for(j5=0; j5<nu; j5++)
        if(j5==0 || j5==nu-1) discrete += cu[j5]*
					        Ub(x,y,z,t,ug[j5],v,w);
        else discrete += cu[j5]*U[s(i,ii,i3,i4,j5,i6,i7)];
    else if(xord==0 && yord==0 && zord==0 && tord==0 && uord==0 &&
            vord!=0 && word==0)
      for(j6=0; j6<nv; j6++)
        if(j6==0 || j6==nv-1) discrete += cv[j6]*
					  Ub(x,y,z,t,u,vg[j6],w);
        else discrete += ct[j6]*U[s(i,ii,i3,i4,i5,j6,i7)];
    else if(xord==0 && yord==0 && zord==0 && tord==0 && uord==0 &&
            vord==0 && word!=0)
      for(j7=0; j7<nw; j7++)
        if(j7==0 || j7==nw-1) discrete += cw[j7]*
					  Ub(x,y,z,t,u,v,wg[j7]);
        else discrete += cw[j7]*U[s(i,ii,i3,i4,i5,i6,j7)];
    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, i3, i4, i5, i6, i7;
    int ib=0, iib=0, i3b=0, i4b=0, i5b=0, i6b=0, i7b=0;

    smaxerr = 0.0;
    for(i=1; i<nx-1; i++)
    {
      for(ii=1; ii<ny-1; ii++)
      {
        for(i3=1; i3<nz-1; i3++)
        {
          for(i4=1; i4<nt-1; i4++)
          {
            for(i5=1; i5<nu-1; i5++)
            {
              for(i6=1; i6<nv-1; i6++)
              {
                for(i7=1; i7<nw-1; i7++)
                {
      val = 0.0;
      // Uxxxx(x,y,z,t,u,v,w) 
      val += eval_derivative(4,0,0,0,0,0,0,i,ii,i3,i4,i5,i6,i7,U);
      // + Uyyyy(x,y,z,t,u,v,w) 
      val += eval_derivative(0,4,0,0,0,0,0,i,ii,i3,i4,i5,i6,i7,U);
      // + Uzzzz(x,y,z,t,u,v,w) 
      val += eval_derivative(0,0,4,0,0,0,0,i,ii,i3,i4,i5,i6,i7,U);
      // + Utttt(x,y,z,t,u,v,w) 
      val += eval_derivative(0,0,0,4,0,0,0,i,ii,i3,i4,i5,i6,i7,U);
      // + Uuuuu(x,y,z,t,u,v,w) 
      val += eval_derivative(0,0,0,0,4,0,0,i,ii,i3,i4,i5,i6,i7,U);
      // + Uvvvv(x,y,z,t,u,v,w) 
      val += eval_derivative(0,0,0,0,4,0,0,i,ii,i3,i4,i5,i6,i7,U);
      // + Uwwww(x,y,z,t,u,v,w) 
      val += eval_derivative(0,0,0,0,0,0,4,i,ii,i3,i4,i5,i6,i7,U);
      // + 2 Uxx(x,y,z,t,u,v,w) 
      val += 2.0*eval_derivative(2,0,0,0,0,0,0,i,ii,i3,i4,i5,i6,i7,U);
      // + 2 Uyy(x,y,z,t,u,v,w) 
      val += 2.0*eval_derivative(0,2,0,0,0,0,0,i,ii,i3,i4,i5,i6,i7,U);
      // + 2 Uzz(x,y,z,t,u,v,w) 
      val += 2.0*eval_derivative(0,0,2,0,0,0,0,i,ii,i3,i4,i5,i6,i7,U);
      // + 2 Utt(x,y,z,t,u,v,w) 
      val += 2.0*eval_derivative(0,0,0,2,0,0,0,i,ii,i3,i4,i5,i6,i7,U);
      // + 2 Uuu(x,y,z,t,u,v,w) 
      val += 2.0*eval_derivative(0,0,0,0,2,0,0,i,ii,i3,i4,i5,i6,i7,U);
      // + 2 Uvv(x,y,z,t,u,v,w) 
      val += 2.0*eval_derivative(0,0,0,0,0,2,0,i,ii,i3,i4,i5,i6,i7,U);
      // + 2 Uww(x,y,z,t,u,v,w) 
      val += 2.0*eval_derivative(0,0,0,0,0,0,2,i,ii,i3,i4,i5,i6,i7,U);
      // + 8 u(x,y,z,t,u,v,w) 
      val += 8.0*Ub(xg[i],yg[ii],zg[i3],tg[i4],ug[i5],vg[i6],wg[i7]);
      // -       F(x,y,z,t,u,v,w)
      val -=  f(xg[i],yg[ii],zg[i3],tg[i4],ug[i5],vg[i6],wg[i7]);
      err = Math.abs(val);
      if(err>smaxerr) 
      {
        smaxerr = err;
        ib = i;
        iib = ii;
        i3b = i3;
        i4b = i4;
        i5b = i5;
        i6b = i6;
        i7b = i7;
      }
               } // end i7 
             } // end i6 
            } // end i5 
          } // end i4 
        } // end i3 
      } // end ii 
    } // end i 
    System.out.println("check_soln maxerr="+smaxerr);
    System.out.println("at("+ib+","+iib+","+i3b+","+i4b+
                       ","+i5b+","+i6b+","+i7b+")");
  } // end check_soln 

  void write_soln(double U[]) // for plot7d_gl 
  {
    double x, y, z, t, u, v, w;
    try
    {
      FileWriter fout = new FileWriter("pde47hn_eq_java.dat");
      BufferedWriter fileout = new BufferedWriter(fout);
      System.out.println("writing pde47hn_eq_java.dat");
      fileout.write("7 "+nx+" "+ny+" "+nz+" "+nt+" "+nu+" "+nv+" "+nw+"\n");
      for(int i=0; i<nx; i++)
      {
	x = xg[i];
        for(int ii=0; ii<ny; ii++)
        {
          y = yg[ii];
          for(int i3=0; i3<nz; i3++)
          {
            z = zg[i3];
            for(int i4=0; i4<nt; i4++)
            {
              t = tg[i4];
              for(int i5=0; i5<nu; i5++)
              {
                u = ug[i5];
                for(int i6=0; i6<nv; i6++)
                {
                  v = vg[i6];
                  for(int i7=0; i7<nw; i7++)
                  {
                    w = wg[i7];
		    if(i ==0 || i ==nx-1 || ii==0 || ii==ny-1 ||
                       i3==0 || i3==nz-1 || i4==0 || i4==nt-1 ||
                       i5==0 || i5==nu-1 || i6==0 || i6==nv-1 ||
                       i7==0 || i7==nw-1)
                      fileout.write(" "+x+" "+y+" "+z+" "+t+" "+u+
		   		    " "+v+" "+w+" "+Ub(x,y,z,t,u,v,w)+" \n");
                    else
                      fileout.write(" "+x+" "+y+" "+z+" "+t+" "+u+" "+v+
		        " "+w+" "+
                        U[s(i,ii,i3,i4,i5,i6,i7)]+
                                  " \n");
	          } // end i7
	        } // end i6
	      } // end i5
	    } // end i4
	  } // end i3
        } // end ii
      } // end i

      fileout.close();
      System.out.println("finished writing pde47hn_eq_java.dat");
    }
    catch(FileNotFoundException exception)
    {
      System.out.println("pde47hn_eq_java.dat  not opened");
    }
    catch(Exception exception) // catches all exceptions
    {
      System.out.println("pde47hn_eq_java.java output-file error");
    } // end try
  } // end write_soln
 

  public static void main (String[] args)
  {
    new pde47hn_eq();
  } // end main
} // end pde47hn_eq.java