/* pde3b_eq.c  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                                   */

/* rewrite to collect  u(x,y,z)  term                             */
/* u(x-hx,y,z)/hx^2 + u(x+hx,y,z)/hx^2 +                          */
/* u(x,y-hy,z)/hy^2 + u(x,y+hy,z)/hy^2 +                          */
/* u(x,y,z-hz)/hz^2 + u(x,y,z+hz)/hz^2 -                          */
/* u(x,y,z)/(2/hx^2 + 2/hy^2 + 2/hz^2) = c(x,y,z)                 */

/* 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  -                           */
/* u(i,j,k)/(2/hx^2 + 2/hy^2 + 2/hz^2)= 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 */
/* using  x=xmin+i*hx, y=ymin+j*hz, z=zmin+k*hz                   */

/* create two dimensional array with                              */
/* (nx-2)*(ny-2)*(nz-2) rows, one for each equation               */
/* and (nx-2)*(ny-2)*(nz-2)+1 columns 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. */
/* using linear subscripts (i-1)+(nx-2)*(j-1)+(nx-2)*(ny-2)*(k-1) */

#include <stdio.h>
#include <stdlib.h>
#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))

static int nx = 11;
static int ny = 11;
static int nz = 9;
static double us[567]; /* solution being computed 9*9*7      */
static double ut[567][568]; /* temporary equations           */
        /* nx-2 by ny-2 by nz-2 internal, non boundary cells */
static const double xmin = -1.0;
static const double ymin = -1.0;
static const double zmin = -1.0;
static const double xmax =  1.2;
static const double ymax =  1.0;
static const double zmax =  1.0;
static double hx; /* x step = (xmax-xmin)/(nx-1) */
static double hy; /* y step = (ymax-ymin)/(ny-1) */
static double hz; /* z step = (zmax-zmin)/(nz-1) */
static double error;         /* sum of absolute exact-computed         */
static double avgerror;      /* average error                          */
static double maxerror;      /* maximum cell error                     */

static double u(double x, double y, double z) /* 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;
}

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

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

static void printus()
{
  int i, j, k;
  
  printf("\n computed solution\n");
  for(k=1; k<nz-1; k++)
  {
    printf("u[x,y,%d] plane \n", k);
    for(i=1; i<nx-1; i++)
    {
      for(j=1; j<min(ny-1,11); j++)
      {
        printf("%6.3f ", us[(i-1)+(nx-2)*(j-1)+(nx-2)*(ny-2)*(k-1)]);
      }
      printf("\n");
    }
  }
  printf("\n");
} /* end printus */

static void printexact() /* typically not known */
{
  int i, j, k;
  
  printf("\n exact solution\n");
  for(k=1; k<nz-1; k++)
  {
    printf("u[x,y,%d] plane \n", k);
    for(i=1; i<nx-1; i++)
    {
      for(j=1; j<min(ny-1,11); j++)
      {
        printf("%6.3f ", u(xmin+i*hx, ymin+j*hy, zmin+k*hz));
      }
      printf("\n");
    }
  }
  if(min(ny,11)<ny) printf("Not all columns printed \n");
  printf("\n");
} /* end printexact */

static void printdiff()
{
  int i, j, k;
  
  printf("\n difference exact-computed solution\n");
  printf("\n");
  for(k=1; k<nz-1; k++)
  {
    printf("u[x,y,%d] plane \n", k);
    for(i=1; i<nx-1; i++)
    {
      printf("       ");
      for(j=1; j<min(ny-1,10); j++)
      {
        printf("%7.4f", u(xmin+i*hx, ymin+j*hy, zmin+k*hz)-
               us[(i-1)+(nx-2)*(j-1)+(nx-2)*(ny-2)*(k-1)]);
      }
      printf("\n");
    }
    printf("\n");
  }
} /* end printdiff */

static 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;

  printf("internal cells zeroed \n");          
  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 */

static void simeq(int n, int m, double B[567][568], 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 */

int main(int argc, char * argv[])
{
  int i, j, k;
  
  printf("pde3b_eq.c running\n");
  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 */
  printf("xmin=%g, xmax=%g, nx=%d, hx=%g \n", xmin, xmax, nx, hx);
  printf("ymin=%g, ymax=%g, ny=%d, hy=%g \n", ymin, ymax, ny, hy);
  printf("zmin=%g, zmax=%g, nz=%d, hz=%g \n", zmin, zmax, nz, hz);

  init_matrix();
  printf("matrix initialized\n");
  printf(" initial matrix, left side, upper \n");
  for(i=0; i<11; i++)
  {
    for(j=0; j<11; j++)
    {
      printf("%5.2f ",ut[i][j]);
    }
    printf("\n");
  } 
  printf("\n initial matrix, right side, upper \n");
  for(i=0; i<11; i++)
  {
    for(j=557; j<568; j++)
    {
      printf("%5.2f ",ut[i][j]);
    }
    printf("\n");
  } 
  printf("\n");
  printf("\n initial matrix, right side, lower \n");
  for(i=557; i<567; i++)
  {
    for(j=557; j<568; j++)
    {
      printf("%5.2f ",ut[i][j]);
    }
    printf("\n");
  } 
  printf("\n");
  
  simeq((nx-2)*(ny-2)*(nz-2), (nx-2)*(ny-2)*(nz-2)+1, ut, us);
  printexact();
  printus();
  printdiff();

  check();
  avgerror = error/((nx-2)*(ny-2)*(nz-2));
  printf("total error=%g, average error=%g, max error=%g\n",
         error, avgerror, maxerror);
  return 0;
} /* end main of pde3b_eq.c */