// pde47h_nu.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[iii], t[iiii], u[iiiii],
//                        v[iiiiii], w[iiiiiii]
//
// 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 nuderiv.java
// uses simeq.java

class pde47h_nu
{
  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;

  pde47h_nu()
  {
    double x, y, z, t, u, v, w;
    int i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,
        j,jj,jjj,jjjj,jjjjj,jjjjjj,jjjjjjj;
    int k, cs;
    // coefficients of terms in first equation 
    // derivative point arrays, dynamic in this version 
    double cxx[] = 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("pde47h_nu.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 = (double)(nx*(nx-1))/2.0;
    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] = xg[i-1]+(xmax-xmin)*(double)(nx-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 = (double)(ny*(ny-1))/2.0;
    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] = yg[ii-1]+(ymax-ymin)*(double)(ny-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 = (double)(nz*(nz-1))/2.0;
    iii = 0;
    zg[iii] = zmin;
    System.out.println("iii="+iii+", Ua("+xmin+","+ymin+","+zg[iii]+
                     ","+tmin+","+umin+","+vmin+","+wmin+")="+
		       Ub(xmin,ymin,zg[iii],tmin,umin,vmin,wmin)); 
    for(iii=1; iii<nz; iii++)
    {
      zg[iii] = zg[iii-1]+(zmax-zmin)*(double)(nz-iii)/hz;
      System.out.println("iii="+iii+", Ua("+xmin+","+ymin+","+zg[iii]+
                         ","+tmin+","+umin+","+vmin+","+wmin+")="+
                         Ub(xmin,ymin,zg[iii],tmin,umin,vmin,wmin)); 
    } // iii  z   
    System.out.println(" ");
    System.out.println("t grid and analytic solution at minimums ");
    ht = (double)(nt*(nt-1))/2.0;
    iiii = 0;
    tg[iiii] = tmin;
    System.out.println("iiii="+iiii+", Ua("+xmin+","+ymin+","+zmin+
                       ","+tg[iiii]+","+umin+","+vmin+","+wmin+")="+
                       Ub(xmin,ymin,zmin,tg[iiii],umin,vmin,wmin)); 
    for(iiii=1; iiii<nt; iiii++)
    {
      tg[iiii] = tg[iiii-1]+(tmax-tmin)*(double)(nt-iiii)/ht;
      System.out.println("iiii="+iiii+", Ua("+xmin+","+ymin+","+zmin+
                         ","+tg[iiii]+","+umin+","+vmin+","+wmin+")="+
                         Ub(xmin,ymin,zmin,tg[iiii],umin,vmin,wmin)); 
    } // iiii  t 
    System.out.println(" ");
    System.out.println("u grid and analytic solution at minimums ");
    hu = (double)(nu*(nu-1))/2.0;
    iiiii = 0;
    ug[iiiii] = umin;
    System.out.println("iiiii="+iiiii+", Ua("+xmin+","+ymin+","+zmin+
                       ","+tmin+","+ug[iiiii]+","+vmin+","+wmin+")="+
                       Ub(xmin,ymin,zmin,tmin,ug[iiiii],vmin,wmin)); 
    for(iiiii=1; iiiii<nu; iiiii++)
    {
      ug[iiiii] = ug[iiiii-1]+(umax-umin)*(double)(nu-iiiii)/hu;
      System.out.println("iiiii="+iiiii+", Ua("+xmin+","+ymin+","+zmin+
                         ","+tmin+","+ug[iiiii]+","+vmin+","+wmin+")="+
                         Ub(xmin,ymin,zmin,tmin,ug[iiiii],vmin,wmin)); 
    } // iiiii  u 
    System.out.println(" ");
    System.out.println("v grid and analytic solution at minimums ");
    hv = (double)(nv*(nv-1))/2.0;
    iiiiii = 0;
    vg[iiiiii] = vmin;
    System.out.println("iiiiii="+iiiiii+", Ua("+xmin+","+ymin+","+zmin+
                     ","+tmin+","+umin+","+vg[iiiiii]+","+wmin+")="+
		     Ub(xmin,ymin,zmin,tmin,umin,vg[iiiiii],wmin)); 
    for(iiiiii=1; iiiiii<nv; iiiiii++)
    {
      vg[iiiiii] = vg[iiiiii-1]+(vmax-vmin)*(double)(nv-iiiiii)/hv;
      System.out.println("iiiiii="+iiiiii+", Ua("+xmin+","+ymin+","+zmin+
                         ","+tmin+","+umin+","+vg[iiiiii]+","+wmin+")="+
                         Ub(xmin,ymin,zmin,tmin,umin,vg[iiiiii],wmin)); 
    } // iiiiii  v 
    System.out.println(" ");
    System.out.println("w grid and analytic solution at minimums ");
    hw = (double)(nw*(nw-1))/2.0;
    iiiiiii = 0;
    wg[iiiiiii] = wmin;
    System.out.println("iiiiiii="+iiiiiii+", Ua("+xmin+","+ymin+","+zmin+
                     ","+tmin+","+umin+","+vmin+","+wg[iiiiiii]+")="+
		       Ub(xmin,ymin,zmin,tmin,umin,vmin,wg[iiiiiii])); 
    for(iiiiiii=1; iiiiiii<nw; iiiiiii++)
    {
      wg[iiiiiii] = wg[iiiiiii-1]+(wmax-wmin)*(double)(nw-iiiiiii)/hw;
      System.out.println("iiiiiii="+iiiiiii+", Ua("+xmin+","+ymin+","+zmin+
                       ","+tmin+","+umin+","+vmin+","+wg[iiiiiii]+")="+
		       Ub(xmin,ymin,zmin,tmin,umin,vmin,wg[iiiiiii])); 
    } // iiiiiii  w 
    System.out.println(" ");

    System.out.println("check nuderiv(4,nv,0,vg,cvvvv ");
    new nuderiv(4,nv,0,vg,cvvvv);
    for(i=0; i<nv; i++) System.out.println("cvvvv["+i+"]="+cvvvv[i]);
    System.out.println("check nuderiv(4,nv,nv-1,vg,cvvvv ");
    new nuderiv(4,nv,nv-1,vg,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(iii=1; iii<nz-1; iii++)
      {
        z = zg[iii];
        for(iiii=1; iiii<nt-1; iiii++)
        {
          t = tg[iiii];
          for(iiiii=1; iiiii<nu-1; iiiii++)
          {
            u = ug[iiiii];
            for(iiiiii=1; iiiiii<nv-1; iiiiii++)
            {
              v = vg[iiiiii];
              for(iiiiiii=1; iiiiiii<nw-1; iiiiiii++)
              {
                w = wg[iiiiiii];
                Ua[s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)] = 
                   Ub(x,y,z,t,u,v,w);
              } // iiiiiii 
            } // iiiiii 
          } // iiiii 
        } // iiii 
      } // iii 
    } // ii 
  } // i 
  time_now = (double)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(iii=1; iii<nz-1; iii++)
      {
        for(iiii=1; iiii<nt-1; iiii++)
        {
          for(iiiii=1; iiiii<nu-1; iiiii++)
          {
            for(iiiiii=1; iiiiii<nv-1; iiiiii++)
            {
              for(iiiiiii=1; iiiiiii<nw-1; iiiiiii++)
              {
                for(j=1; j<nx-1; j++)
                {
                  for(jj=1; jj<ny-1; jj++)
                  {
                    for(jjj=1; jjj<nz-1; jjj++)
                    {
                      for(jjjj=1; jjjj<nt-1; jjjj++)
                      {
                        for(jjjjj=1; jjjjj<nu-1; jjjjj++)
                        {
                          for(jjjjjj=1; jjjjjj<nv-1; jjjjjj++)
                          {
                            for(jjjjjjj=1; jjjjjjj<nw-1; jjjjjjj++)
                            {
			      A[ss(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii,
				   j,jj,jjj,jjjj,jjjjj,jjjjjj,jjjjjjj)] = 0.0;
	                    } // jjjjjjj 
	                  } // jjjjjj 
	                } // jjjjj 
	              } // jjjj 
	            } // jjj 
	          } // jj 
                } // j 
              } // iiiiiii 
            } // iiiiii 
          } // iiiii 
        } // iiii 
      } // iii 
    } // ii 
  } // i 

  time_now = (double)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(iii=1; iii<nz-1; iii++)
      {
        for(iiii=1; iiii<nt-1; iiii++)
        {
          for(iiiii=1; iiiii<nu-1; iiiii++)
          {
            for(iiiiii=1; iiiiii<nv-1; iiiiii++)
            { 
              for(iiiiiii=1; iiiiiii<nw-1; iiiiiii++)
              { 
		  k = s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii); // row index 

                // for each term in first equation 

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

                // end terms 

	      } // end iiiiiii loop 
	    } // end iiiiii loop 
	  } // end iiiii loop 
	} // end iiii loop 
      } // end iii 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(iii=1; iii<nz-1; iii++)
        {
          for(iiii=1; iiii<nt-1; iiii++)
          {
            for(iiiii=1; iiiii<nu-1; iiiii++)
            {
              for(iiiiii=1; iiiiii<nv-1; iiiiii++)
              {
                for(iiiiiii=1; iiiiiii<nw-1; iiiiiii++)
                {
		    err = Math.abs(U[s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)]-
				   Ua[s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)]);
                if(err>maxerr) maxerr = err;
                avgerr = avgerr + err;
                System.out.println("ug["+i+","+ii+","+iii+","+iiii+","+
                                   iiiii+","+iiiiii+","+iiiiiii+"]="+
                                   U[s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)]+
		                   ", Ua="+Ua[s(i,ii,iii,iiii,iiiii,iiiiii,iiiiiii)]+
                                   ", err="+err);
                } // iiiiiii 
              } // iiiiii 
            } // iiiii 
          } // iiii 
        } // iii 
      } // 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 pde47h_nu.java ");
  } // end pde47h_nu

  // utility functions
  int s(int i, int ii, int iii, int iiii, int iiiii,
        int iiiiii, int iiiiiii)
       // 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) +
        (iii-1)*(nt-2)*(nu-2)*(nv-2)*(nw-2) + 
        (iiii-1)*(nu-2)*(nv-2)*(nw-2) +
        (iiiii-1)*(nv-2)*(nw-2) + 
        (iiiiii-1)*(nw-2) + 
        (iiiiiii-1);
    if(k<0 || k>=nxyztuvw) System.out.println("bad s("+i+","+ii+","+iii+
		","+iiii+","+iiiii+","+iiiiii+","+iiiiiii+") ="+k);  
    return k;
  } // end s 

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

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

  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 iii, int iiii,
                         int iiiii, int iiiiii, int iiiiiii, double U[])
  {
    int j, jj, jjj, jjjj, jjjjj, jjjjjj, jjjjjjj;
    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 nuderiv(vord, nx, i, xg, cx);
    x = xg[i];
    new nuderiv(word, ny, ii, yg, cy);
    y = yg[ii];
    new nuderiv(xord, nz, iii, zg, cz);
    z = zg[iii];
    new nuderiv(yord, nt, iiii, tg, ct);
    t = tg[iiii];
    new nuderiv(zord, nu, iiiii, ug, cu);
    u = ug[iiiii];
    new nuderiv(tord, nv, iiiiii, vg, cv);
    v = vg[iiiiii];
    new nuderiv(tord, nw, iiiiiii, wg, cw);
    w = wg[iiiiii];
    discrete = 0.0;


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

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

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

  void write_soln(double U[]) // for plot7d_gl 
  {
    double x, y, z, t, u, v, w;
    try
    {
      FileWriter fout = new FileWriter("pde47h_nu_java.dat");
      BufferedWriter fileout = new BufferedWriter(fout);
      System.out.println("writing pde47h_nu_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 iii=0; iii<nz; iii++)
          {
            z = zg[iii];
            for(int iiii=0; iiii<nt; iiii++)
            {
              t = tg[iiii];
              for(int iiiii=0; iiiii<nu; iiiii++)
              {
                u = ug[iiiii];
                for(int iiiiii=0; iiiiii<nv; iiiiii++)
                {
                  v = vg[iiiiii];
                  for(int iiiiiii=0; iiiiiii<nw; iiiiiii++)
                  {
                    w = wg[iiiiiii];
		    if(i==0 || i==nv-1 || ii==0 || ii==nw-1 ||
                       iii==0 || iii==nx-1 || iiii==0 || iiii==ny-1 ||
                       iiiii==0 || iiiii==nz-1 || iiiiii==0 || iiiiii==nt-1 ||
                       iiiiiii==0 || iiiiiii==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,iii,iiii,iiiii,iiiiii,iiiiiii)]+
                                  " \n");
	          } // end iiiiiii
	        } // end iiiiii
	      } // end iiiii
	    } // end iiii
	  } // end iii
        } // end ii
      } // end i

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

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