// pde23_eqn2.c  3D uniform grid boundary value PDE              
//              user defined via web using pdedef232.h pdedef232.c
// set up and solve system of linear equations                       
// create two dimensional array with                                 
// (nx-2)*(ny-2)*(nz-2) rows, one for each equation                  

#include "nuderiv.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "pdedef232.h"  // defines  nx, ny, nz

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

              // not including boundary
#define nxyz  (nx-2)*(ny-2)*(nz-2)
static double xg[nx];
static double yg[ny];
static double zg[nz];
                        // nx-2 by ny-2 by nz-2 internal, non boundary cells 
static double cxx[nx];       // nuderiv coefficients 
static double cyy[ny];
static double czz[nz];
static double cxy[nx][ny];
static double cxz[nx][nz];
static double cyz[ny][nz];
static double cx[nx];
static double cy[ny];
static double cz[nz];
static double coefdxx, coefdyy, coefdzz, 
              coefdxdy, coefdxdz, coefdydz,
              coefdx, coefdy, coefdz, coefu;

static double us[nxyz];         // solution being computed 9*8*7 
static double ut[nxyz][nxyz+1]; // temporary equations 
static double error;            // sum of absolute exact-computed         
static double avgerror;         // average error                          
static double maxerror;         // maximum cell error                     


int S(int i, int j, int k) // index into us and row into ut
{
  return (i-1)+(nx-2)*(j-1)+(nx-2)*(ny-2)*(k-1);
} // end S

static void check() // typically not known, instrumentation 
{
  int i, j, k;
  double uexact, 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++)
      {
        uexact = ub(xg[i], yg[j],zg[k]);
        err = abs(us[S(i,j,k)]-uexact);
        error = error + err;
        if(err>maxerror) maxerror=err;
      }
    }
  }
} // end check 

static void printdiff()
{
  int i, j, k;
  
  printf("difference exact-computed solution, partial\n");
  k = 1;
  for(i=1; i<nx-1; i++)
  {
    for(j=1; j<ny-1; j++)
    {
      for(k=1; k<nz-1; k++)
      {
        printf("xg[%d]=%f,yg[%d]=%f, zg[%d]=%f, diff(x,y,z)=%g \n",
               i, xg[i], j, yg[j], k, zg[k], ub(xg[i],yg[j],zg[k])
               -us[S(i,j,k)]);
      }
    }
  }
  printf("\n");
} // end printdiff 

static void init_matrix()
{
  int i, j, k, m, ii, jj, kk, imk, mjk, ijm, kx, ky, kz ,kxky, kxkz, kykz, cs;
  double x, y, z;

  cs = nxyz; // constant column subscript 
  
  // zero internal cells (zero based row and column)
  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[S(i,j,k)][S(ii,jj,kk)] = 0.0;
        ut[S(i,j,k)][cs] = 0.0;
      }
    }
  }
  printf("internal cells zeroed \n");          
  for(i=1; i<nx-1; i++) // inside boundary 
  {
    x = xg[i];
    for(j=1; j<ny-1; j++) // inside boundary 
    {
      y = yg[j];
      for(k=1; k<nz-1; k++) // inside boundary 
      {
        z = zg[k];
        ii = S(i,j,k); // row of matrix 
        // for each term 
          // coefdxx * d^2 ub(x,y,z)/dx^ 
	  coefdxx = ccxx(x,y,z);
          nuderiv(2,nx,i,xg,cxx);
          for(m=0; m<nx; m++)
  	  {
            cxx[m] = cxx[m] * coefdxx;
            mjk = (m-1)+(nx-2)*(j-1)+(nx-2)*(ny-2)*(k-1);
	    if(m==0 || m==nx-1)
                 ut[ii][cs]  = ut[ii][cs]  - cxx[m]*ub(xg[m],yg[j],zg[k]);
	    else ut[ii][mjk] = ut[ii][mjk] + cxx[m];
	  }

          // coefdyy * d^2 ub(x,y,z)/dy^2 
          coefdyy = ccyy(x,y,z);
          nuderiv(2,ny,j,yg,cyy);
          for(m=0; m<ny; m++)
 	  {
            cyy[m] = cyy[m] * coefdyy;
            imk = (i-1)+(nx-2)*(m-1)+(nx-2)*(ny-2)*(k-1);
	    if(m==0 || m==ny-1)
                 ut[ii][cs]  = ut[ii][cs]  - cyy[m]*ub(xg[i],yg[m],zg[k]);
	    else ut[ii][imk] = ut[ii][imk] + cyy[m];
 	  }

          // coefdzz * d^2 ub(x,y,z)/dz^2 
          coefdzz = cczz(x,y,z);
          nuderiv(2,nz,k,zg,czz);
          for(m=0; m<nz; m++)
 	  {
            czz[m] = czz[m] * coefdzz;
            ijm = (i-1)+(nx-2)*(j-1)+(nx-2)*(ny-2)*(m-1);
	    if(m==0 || m==nz-1)
                 ut[ii][cs]  = ut[ii][cs]  - czz[m]*ub(xg[i],yg[j],zg[m]);
	    else ut[ii][ijm] = ut[ii][ijm] + czz[m];
 	  }

          // coefdxdy * d^2 ub(x,y,z)/dxdy 
          coefdxdy = ccxy(x,y,z);
          nuderiv(1,nx,i,xg,cx);
          nuderiv(1,ny,j,yg,cy);
          for(kx=0; kx<nx; kx++)
	  {
            for(ky=0; ky<ny; ky++)
	    {
              cxy[kx][ky] = cx[kx] * cy[ky] * coefdxdy;
              kxky = (kx-1)+(nx-2)*(ky-1)+(nx-2)*(ny-2)*(k-1);
	      if(kx==0 || kx==nx-1 || ky==0 || ky==ny-1)
                   ut[ii][cs]   = ut[ii][cs]   - cxy[kx][ky]*ub(xg[kx],yg[ky],zg[k]);
	      else ut[ii][kxky] = ut[ii][kxky] + cxy[kx][ky];
            }
	  }

          // coefdxdz * d^2 ub(x,y,z)/dxdz 
          coefdxdz = ccxz(x,y,z);
          nuderiv(1,nx,i,xg,cx);
          nuderiv(1,nz,k,zg,cz);
          for(kx=0; kx<nx; kx++)
	  {
            for(kz=0; kz<nz; kz++)
	    {
              cxz[kx][kz] = cx[kx] * cz[kz] * coefdxdz;
              kxkz = (kx-1)+(nx-2)*(j-1)+(nx-2)*(ny-2)*(kz-1);
	      if(kx==0 || kx==nx-1 || kz==0 || kz==nz-1)
                   ut[ii][cs]   = ut[ii][cs]   - cxz[kx][kz]*ub(xg[kx],yg[j],zg[kz]);
	      else ut[ii][kxkz] = ut[ii][kxkz] + cxz[kx][kz];
            }
	  }

          // coefdydz * d^2 ub(x,y,z)/dydz 
          coefdydz = ccyz(x,y,z);
          nuderiv(1,ny,j,yg,cy);
          nuderiv(1,nz,k,zg,cz);
          for(ky=0; ky<ny; ky++)
	  {
            for(kz=0; kz<nz; kz++)
	    {
              cyz[ky][kz] = cy[ky] * cz[kz] * coefdydz;
              kykz = (i-1)+(nx-2)*(ky-1)+(nx-2)*(ny-2)*(kz-1);
	      if(ky==0 || ky==ny-1 || kz==0 || kz==nz-1)
                   ut[ii][cs]   = ut[ii][cs]   - cyz[ky][kz]*ub(xg[i],yg[ky],zg[kz]);
	      else ut[ii][kykz] = ut[ii][kykz] + cyz[ky][kz];
            }
	  }

          // coefdx * d ub(x,y,z)/dx 
          coefdx = ccx(x,y,z);
          nuderiv(1,nx,i,xg,cx);
          for(m=0; m<nx; m++)
	  {
            cx[m] = cx[m] * coefdx;
            mjk = (m-1)+(nx-2)*(j-1)+(nx-2)*(ny-2)*(k-1);
	    if(m==0 || m==nx-1)
                 ut[ii][cs]  = ut[ii][cs]  - cx[m]*ub(xg[m],yg[j],zg[k]);
	    else ut[ii][mjk] = ut[ii][mjk] + cx[m];
	  }

          // coefdy * d ub(x,y,z)/dy 
          coefdy = ccy(x,y,z);
          nuderiv(1,ny,j,yg,cy);
          for(m=0; m<ny; m++)
	  {
            cy[m] = cy[m] * coefdy;
            imk = (i-1)+(nx-2)*(m-1)+(nx-2)*(ny-2)*(k-1);
	    if(m==0 || m==ny-1)
                 ut[ii][cs]  = ut[ii][cs]  - cy[m]*ub(xg[i],yg[m],zg[k]);
	    else ut[ii][imk] = ut[ii][imk] + cy[m];
	  }

          // coefdz * d ub(x,y,z)/dz 
          coefdz = ccz(x,y,z);
          nuderiv(1,nz,k,zg,cz);
          for(m=0; m<nz; m++)
	  {
            cz[m] = cz[m] * coefdz;
            ijm = (i-1)+(nx-2)*(j-1)+(nx-2)*(ny-2)*(m-1);
	    if(m==0 || m==nz-1)
                 ut[ii][cs]  = ut[ii][cs]  - cz[m]*ub(xg[i],yg[j],zg[m]);
	    else ut[ii][ijm] = ut[ii][ijm] + cz[m];
	  }

          // coefu * ub(x,y,z) 
          coefu = ccc(x,y,z);
          ut[ii][ii] = ut[ii][ii] + 1.0*coefu;

          // RHS constant 
          ut[ii][cs] = ut[ii][cs] + f(x,y,z);
        // end terms 
      }
    }
  }
} // end init_matrix 

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

  row = (int *)calloc(n, sizeof(int));
  // 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 = abs(pivot);
    i_pivot = k;
    for(i=k; i<n; i++)
    {
      if(abs(B[row[i]][k]) > abs_pivot)
      {
        i_pivot = i;
        pivot = B[row[i]][k];
        abs_pivot = 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;
      }
      printf("redundant row (singular) %d \n", 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 

static void  write_soln3(double U[], char file_name[])
{
  FILE * Fout;
  int i, j, k;

  Fout = fopen(file_name, "w");
  if(Fout==NULL)
  {
    printf("ERROR unable to open %s for writing \n", file_name);
    exit(1);
  }
  printf("writing %s\n", file_name);
  fprintf(Fout, "3 %d %d %d \n", nx, ny, nz);
  for(i=0; i<nx; i++)
  {
    for(j=0; j<ny; j++)
    {
      for(k=0; k<nz; k++)
      {
	if(i==0 || i==nx-1 || j==0 || j==ny-1 || k==0 || k==nz-1)
	{
	  fprintf(Fout," %f %f %f %f \n",
		  xg[i], yg[j], zg[k], ub(xg[i],yg[j],zg[k])); // boundary 
	}
	else
	{
	  fprintf(Fout,"%f %f %f %f \n",
	          xg[i],yg[j],zg[k],U[S(i,j,k)]); // computed solution 
	}
      } // k 
    } // j 
  } // i 
  fclose(Fout);
  printf("finished writing %s\n", file_name);
}  // end write_soln3 

int main(int argc, char * argv[])
{
  int i, j, k, ii, jj, kk, cs;
  double xmin, xmax, ymin, ymax, zmin, zmax;
  double hx, hy, hz;
  double coef[9];

  printf("pde23_eqn2.c running\n");
  printf("calling coefs in pdedef232.c \n");
  coefs(coef);
  xmin = coef[0];
  xmax = coef[1];
  ymin = coef[3];
  ymax = coef[4];
  zmin = coef[6];
  zmax = coef[7];

  printf("xmin=%g, xmax=%g, nx=%d \n", xmin, xmax, nx);
  printf("ymin=%g, ymax=%g, ny=%d \n", ymin, ymax, ny);
  printf("zmin=%g, zmax=%g, nz=%d \n", zmin, zmax, nz);
  printf("number of equations nxyz=%d \n", nxyz);
  cs = nxyz;

  hx = (xmax-xmin)/(nx-1);
  printf("hx=%f \n", hx);
  for(i=0; i<nx; i++)
  {
    xg[i] = xmin+i*hx;
    printf("xg[%2d]=%f \n", i, xg[i]);
  }

  hy = (ymax-ymin)/(ny-1);
  printf("hy=%f \n", hy);
  for(j=0; j<ny; j++)
  {
    yg[j] = ymin+j*hy;
    printf("yg[%2d]=%f \n", j, yg[j]);
  }

  hz = (zmax-zmin)/(nz-1);
  printf("hz=%f \n", hz);
  for(k=0; k<nz; k++)
  {
    zg[k] = zmin+k*hz;
    printf("zg[%2d]=%f \n", k, zg[k]);
  }

  printf("ub(xg[1],yg[1],zg[1])=%g \n",ub(xg[1],yg[1],zg[1]));
  printf("f(xg[1],yg[1],zg[1])=%g \n",f(xg[1],yg[1],zg[1]));

  init_matrix();
  printf("matrix initialized\n");
  
  simeq(nxyz, nxyz+1, ut, us);
  printdiff();
  write_soln3(us, "pde23_eqn2_c.dat");

  check();
  avgerror = error/((double)nxyz);
  printf("total error=%g, average error=%g, max error=%g\n",
         error, avgerror, maxerror);
  printf("pde23_eqn2.c finished \n");
  return 0;
} // end main of pde23_eqn2.c