// pde3b_eq.java  3D non equal h boundary value PDE  non iterative    
//  u(x,y,z) = x^3  + y^3  +z^3  + 2x^2y + 3xy^2 + 4z^2x + 5y + 6  
//     du/dx = 3x^2 + 0    + 0   + 4xy   + 3y^2  + 4z^2  + 0  + 0  
// d^2u/dx^2 = 6x   + 0    + 0   + 4y    + 0     + 0     + 0  + 0  
//     du/dy = 0    + 3y^2 + 0   + 2x^2  + 6xy   + 0     + 5  + 0  
// d^2u/dy^2 = 0    + 6y   + 0   + 0     + 6x    + 0     + 0  + 0  
//     du/dz = 0    + 0    +3z^2 + 0     + 0     + 8zx   + 0  + 0  
// d^2u/dz^2 = 0    + 0    +6z   + 0     + 0     + 8x    + 0  + 0  
//                                                                 
// solve d^2u/dx^2 + d^2u/dy^2 + d^2u/dz^2 = c(x,y,z)              
// by second order method on cube -1,-1,-1 to 1.2,1,1              

// second order numerical difference                              
// d^2u/dx^2 + d^2u/dy^2 + d^2u/dz^2 =                            
// 1/hx^2(u(x-hx,y,z)+u(x+hx,y,z)-2u(x,y,z)) +                    
// 1/hy^2(u(x,y-hy,z)+u(x,y+hy,z)-2u(x,y,z)) +                    
// 1/hz^2(u(x,y,z-hz)+u(x,y,z+hz)-2u(x,y,z)) + O(h^2)             
//  = c(x,y,z) = 12x + 10y - 8z                                   

// solve system of linear equations for                           
// u(i-1,j,k)/hx^2 + u(i+1,j,k)/hx^2 +                            
// u(i,j-1,k)/hy^2 + u(i,j+1,k)/hy^2 +                            
// u(i,j,k-1)/hz^2 + u(i,j,k+1)/hz^2 -                            
// (2/hx^2 + 2/hy^2 + 2/hz^2)*u(i,j,k) = c(x,y,z)                 
// for i=1,nx-2 j=1,ny-2 k=1,nz-2 while substituting              
// u(i,j,k) boundary values i=0 or nx-1, j=0 or ny-1, k=0 or nz-1 
// x=xmin+i*hx, y=ymin+j*hz, z=zmin+k*hz                          

// create two dimensional array with                              
// (nx-2)*(ny-2)*(nz-2) columns, one for each equation            
// and (nx-2)*(ny-2)*(nz-2)+1 rows for variables plus constant    
// Initialize the matrix, see init_matrix in the source code      
// Then solve the system of linear equations to get the solution. 

public class pde3b_eq
{
  final int nx = 11;
  final int ny = 11;
  final int nz = 9;
  double us[] = new double[(nx-2)*(ny-2)*(nz-2)]; // solution being iterated        
  double ut[][] = new double[(nx-2)*(ny-2)*(nz-2)]
                            [(nx-2)*(ny-2)*(nz-2)+1]; // temporary for partial results 
                     // nx-2 by ny-2 by nz-2 internal, non boundary cells 
  final double xmin = -1.0;
  final double ymin = -1.0;
  final double zmin = -1.0;
  final double xmax =  1.2;
  final double ymax =  1.0;
  final double zmax =  1.0;
  double hx;   // x step = (xmax-xmin)/(nx-1) 
  double hy;   // y step = (ymax-ymin)/(ny-1) 
  double hz;   // z step = (zmax-zmin)/(nz-1)  
  double error;         // sum of absolute exact-computed         
  double avgerror;      // average error                          
  double maxerror;      // maximum cell error                     

  public pde3b_eq()
  {
    int i, j;
    
    System.out.println("pde3b_eq.java running");
    hx = (xmax-xmin)/(nx-1); // step size in x 
    hy = (ymax-ymin)/(ny-1); // step size in y 
    hz = (zmax-zmin)/(nz-1); // step size in z 
    System.out.println("xmin="+xmin+ "xmax="+xmax+ "nx="+nx+ "hx="+hx);
    System.out.println("ymin="+ymin+ "ymax="+ymax+ "ny="+ny+ "hy="+hy);
    System.out.println("zmin="+zmin+ "zmax="+zmax+ "nz="+nz+ "hz="+hz);
  
    init_matrix();
    System.out.println("matrix initialized");
    System.out.println("initial matrix, left side, upper");
    for(i=0; i<4; i++)
    {
      for(j=0; j<4; j++)
      {
        System.out.print(" "+ut[i][j]);
      }
      System.out.println(" ");
    }
    System.out.println(" ");
    System.out.println("initial matrix, right side, upper");
    for(i=0; i<4; i++)
    {
      for(j=564; j<568; j++)
      {
        System.out.print(" "+ut[i][j]);
      }
      System.out.println(" ");
    } 
    System.out.println(" ");
    System.out.println("initial matrix, right side, lower");
    for(i=563; i<567; i++)
    {
      for(j=564; j<568; j++)
      {
        System.out.print(" "+ut[i][j]);
      }
      System.out.println(" ");
    } 
    System.out.println(" ");
  
  
    simeq((nx-2)*(ny-2)*(nz-2), ut, us);
    printexact();
    printus();
    printdiff();
    check();
    avgerror = error/((nx-2)*(ny-2)*(nz-2));
    System.out.println("error="+error+ "  avg="+avgerror+ "  max="+maxerror);
  } // end pde3b_eq

  double u(double x, double y, double z) // solution for boundary and check 
  {
    return x*x*x + y*y*y + z*z*z + 2.0*x*x*y + 3.0*x*y*y +
           4.0*z*z*x + 5.0*y +6.0;
  }

  double c(double x, double y, double z) // from solve 
  {
    return 20.0*x + 10.0*y + 6.0*z;
  }

  void check() // typically not known, instrumentation 
  {
    int i, j, k;
    double zexact, err;
  
    error = 0.0;
    maxerror = 0.0;
    for(i=1; i<nx-1; i++)
    {
      for(j=1; j<ny-1; j++)
      {
        for(k=1; k<nz-1; k++)
        {
          zexact = u(xmin+i*hx,ymin+j*hy,zmin+k*hz);
          err = Math.abs(us[(i-1)+(nx-2)*(j-1)+(nx-2)*(ny-2)*(k-1)]-zexact);
          error = error + err;
          if(err>maxerror) maxerror=err;
        }
      }
    }
  } // end check

  void printus()
  {
    int i, j, k;
  
    System.out.println(" computed solution");
    for(i=1; i<nx-1; i++)
    {
      for(j=1; j<ny-1; j++)
      {
        for(k=1; k<nz-1; k++)
        {
          System.out.println("us["+i+"]["+j+"]["+k+"]= "+
                              us[(i-1)+(nx-2)*(j-1)+(nx-2)*(ny-2)*(k-1)]);
        }
      }
    }
  } // end printus

  void printexact()
  {
    int i, j, k;
  
    System.out.println(" exact solution");
    for(i=0; i<nx; i++)
    {
      for(j=0; j<ny; j++)
      {
        for(k=0; k<nz; k++)
        {
          System.out.println("u["+i+"]["+j+"]["+k+"]="+
                             u(xmin+i*hx,ymin+j*hy,zmin+k*hz));
        }
      }
    }
  } // end printexact

  void printdiff()
  {
    int i, j, k;
  
    System.out.println(" difference exact-computed solution");
    for(i=1; i<nx-1; i++)
    {
      for(j=1; j<ny-1; j++)
      {
        for(k=1; k<nz-1; k++)
        {
          System.out.println("diff["+i+"]["+j+"]["+k+"]= "+
                     (u(xmin+i*hx,ymin+j*hy,zmin+k*hz)-
                     us[(i-1)+(nx-2)*(j-1)+(nx-2)*(ny-2)*(k-1)]));
        }
      }
    }
  } // end printdiff
  

   void init_matrix()
  {
    int i, j, k, ii, jj, kk, cs;
    double hx2, hy2, hz2;

    hx2 = hx*hx;
    hy2 = hy*hy;
    hz2 = hz*hz;
    cs = (nx-2)*(ny-2)*(nz-2); // constant column subscript 
  
    // zero internal cells 
    for(i=1; i<nx-1; i++)
      for(j=1; j<ny-1; j++)
        for(k=1; k<nz-1; k++)
         for(ii=1; ii<nx-1; ii++)
           for(jj=1; jj<ny-1; jj++)
             for(kk=1; kk<nz-1; kk++) 
               ut[(i-1)+(nx-2)*(j-1)+(nx-2)*(ny-2)*(k-1)]
                 [(ii-1)+(nx-2)*(jj-1)+(nx-2)*(ny-2)*(kk-1)] = 0.0;

    System.out.println("internal cells zeroed");          
    for(i=1; i<nx-1; i++) // inside boundary 
    {
      for(j=1; j<ny-1; j++)
      {
        for(k=1; k<nz-1; k++)
        {
          ii = (i-1)+(nx-2)*(j-1)+(nx-2)*(ny-2)*(k-1); // column of matrix 
          ut[ii][ii] = -(2.0/hx2 + 2.0/hy2 + 2.0/hz2);
          ut[ii][cs] = c(xmin+i*hx, ymin+j*hy, zmin+k*hz);
          if((i+1)<(nx-1))
            ut[ii][((i-1)+1)+(nx-2)*(j-1)+(nx-2)*(ny-2)*(k-1)] = 1.0/hx2;
          else
            ut[ii][cs] = ut[ii][cs] - u(xmin+(i+1)*hx, ymin+j*hy, zmin+k*hz)/hx2;
          if((i-1)!=0)
            ut[ii][((i-1)-1)+(nx-2)*(j-1)+(nx-2)*(ny-2)*(k-1)] = 1.0/hx2;
          else
            ut[ii][cs] = ut[ii][cs] - u(xmin+(i-1)*hx, ymin+j*hy, zmin+k*hz)/hx2;
          if((j+1)<(ny-1))
            ut[ii][(i-1)+(nx-2)*((j-1)+1)+(nx-2)*(ny-2)*(k-1)] = 1.0/hy2;
          else
            ut[ii][cs] = ut[ii][cs] - u(xmin+i*hx, ymin+(j+1)*hy, zmin+k*hz)/hy2;
          if((j-1)!=0)
            ut[ii][(i-1)+(nx-2)*((j-1)-1)+(nx-2)*(ny-2)*(k-1)] = 1.0/hy2;
          else
            ut[ii][cs] = ut[ii][cs] - u(xmin+i*hx, ymin+(j-1)*hy, zmin+k*hz)/hy2;
          if((k+1)<(nz-1))
            ut[ii][(i-1)+(nx-2)*(j-1)+(nx-2)*(ny-2)*((k-1)+1)] = 1.0/hz2;
          else
            ut[ii][cs] = ut[ii][cs] - u(xmin+i*hx, ymin+j*hy, zmin+(k+1)*hz)/hz2;
          if((k-1)!=0)
            ut[ii][(i-1)+(nx-2)*(j-1)+(nx-2)*(ny-2)*((k-1)-1)] = 1.0/hz2;
          else
            ut[ii][cs] = ut[ii][cs] - u(xmin+i*hx, ymin+j*hy, zmin+(k-1)*hz)/hz2;
        }
      }
    }
  } // end init_matrix 

  void simeq(int n, double B[][], double X[])
  { // tailored version for this problem 
    int row[] = new int[n];  // row interchange indices 
    int hold, i_pivot;       // pivot indices 
    double pivot;            // pivot element value 
    double abs_pivot;
    int i, j, k;

    // set up row  interchange vectors 
    for(k=0; k<n; k++) row[k] = k;
  
    // begin main reduction loop 
    for(k=0; k<n; k++)
    {
      // find largest element for pivot 
      pivot = B[row[k]][k];
      abs_pivot = Math.abs(pivot);
      i_pivot = k;
      for(i=k; i<n; i++)
      {
        if(Math.abs(B[row[i]][k]) > abs_pivot)
        {
          i_pivot = i;
          pivot = B[row[i]][k];
          abs_pivot = Math.abs(pivot);
        }
      }

      // have pivot, interchange row pointers 
      hold = row[k];
      row[k] = row[i_pivot];
      row[i_pivot] = hold;
 
      // check for near singular 
      if( abs_pivot < 1.0E-10 )
      {
        for(j=k+1; j<n+1; j++)
        {
          B[row[k]][j] = 0.0;
        }
        System.out.println("redundant row (singular) "+row[k]);
      } // end singular, delete row 
      else
      {
        // reduce about pivot 
        for(j=k+1; j<n+1; j++)
        {
          B[row[k]][j] = B[row[k]][j] / B[row[k]][k];
        }
        // inner reduction loop 
        for(i=0; i<n; i++)
        {
          if(i != k)
          {
            for(j=k+1; j<n+1; j++)
            {
              B[row[i]][j] = B[row[i]][j] - B[row[i]][k] * B[row[k]][j];
            }
          }
        }
      } // finished inner reduction 
    }

    // end of main reduction loop 
    // build  x  for return, unscrambling rows 
    for(i=0; i<n; i++)
    {
      X[i] = B[row[i]][n];
    }
  } // end simeq 

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