/* pdenu23_eq.c  3D non uniform grid boundary value PDE              */
/*               known solution to test method                       */
/*               u(x,y,z) = x^3+2y^3+3z^3+4x^2y+5zy^2+6zx^2+7xy+8z+9 */
/*                                                                   */
/*               differential equation to solve                     */
/* d^2u/dx^2 + 2*d^2u/dy^2 + 3*d^2u/dz^2 + 4*d^2u/dxdy + 5*d^u/dxdz +*/
/* 6*d^2u/dydz + 7*du/dx + 8*du/dy + 9*du/dz + 10*u =                */
/* 10x^3+20y^3+30z^3+40x^2y+50y^2z+60z^2x+53x^2+93y^2+123z^2+126xy+  */
/* 80yz+108zx+130x+141y+214z+190                                     */
/* non uniform grid on cude (-1.0,-1.1,-1.2) to (1.2,1.3,1.4)        */
/*                                                                   */
/* 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>
#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;          /* includes boundary */
static int ny = 10;          /* includes boundary */
static int nz = 9;           /* includes boundary */
static double xg[11] = {-1.0, -0.9, -0.8, -0.6, -0.3, 0.0,
			 0.35, 0.65, 0.85, 1.0, 1.2};
static double yg[10] = {-1.1, -1.0, -0.9, -0.6, -0.2,
                         0.3,  0.8,  1.0,  1.2,  1.3};
static double zg[10] = {-1.2, -1.0, -0.9, -0.5, 0.0,
                         0.5,  0.9,  1.2,  1.4};
                        /* nx-2 by ny-2 by nz-2 internal, non boundary cells */
static double cxx[11];       /* nuderiv coefficients */
static double cyy[10];
static double czz[9];
static double cxy[11][10];
static double cxz[11][9];
static double cyz[10][9];
static double cx[11];
static double cy[10];
static double cz[9];
static double coefdxx, coefdyy, coefdzz, 
              coefdxdy, coefdxdz, coefdydz,
              coefdx, coefdy, coefdz, coefu;

static double us[504];       /* solution being computed 9*8*7 */
static double ut[504][505];  /* temporary equations */
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 optional check */
{
  return x*x*x + 2.0*y*y*y + 3.0*z*z*z + 4.0*x*x*y + 5.0*y*y*z + 6.0*z*z*x + 
         7.0*x*y + 8.0*z + 9.0;
}

static double c(double x, double y, double z)
                              /* RHS of differential equation to solve */
{
  return 10.0*x*x*x + 20.0*y*y*y + 30.0*z*z*z + 40.0*x*x*y + 50.0*y*y*z + 
    60.0*z*z*x + 53.0*x*x + 93.0*y*y + 123.0*z*z + 126.0*x*y + 108.0*x*z +
    80.0*y*z + 130.0*x + 141.0*y + 214.0*z + 190.0; 
         
}

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 = u(xg[i], yg[j],zg[k]);
        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("computed solution\n");
  for(i=1; i<nx-1; i+=nx-3)
  {
    for(j=1; j<ny-1; j+=ny-3)
    {
      for(k=1; k<nz-1; k+=nz-3)
      {
        printf("xg[%d]=%f,yg[%d]=%f, zg[%d]=%f, us=%g \n",
               i, xg[i], j, yg[j], k, zg[k], 
               us[(i-1)+(nx-2)*(j-1)+(nx-2)*(ny-2)*(k-1)]);
      }
    }
  }
  printf("\n");
} /* end printus */

static void printexact() /* typically not known */
{
  int i, j, k;
  
  printf("exact solution\n");
  for(i=1; i<nx-1; i+=nx-3)
  {
    for(j=1; j<ny-1; j+=ny-3)
    {
      for(k=1; k<nz-1; k+=nz-3)
      {
        printf("xg[%d]=%f,yg[%d]=%f, zg[%d]=%f, u(x,y,z)=%g \n",
               i, xg[i], j, yg[j], k, zg[k], u(xg[i],yg[j],zg[k]));
      }
    }
  }
  printf("\n");
} /* end printexact */

static void printdiff()
{
  int i, j, k;
  
  printf("difference exact-computed solution\n");
  k = 1;
  for(i=1; i<nx-1; i+=nx-3)
  {
    for(j=1; j<ny-1; j+=ny-3)
    {
      for(k=1; k<nz-1; k+=nz-3)
      {
        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], u(xg[i],yg[j],zg[k])
               -us[(i-1)+(nx-2)*(j-1)+(nx-2)*(ny-2)*(k-1)]);
      }
    }
  }
  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;

  cs = (nx-2)*(ny-2)*(nz-2); /* 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[(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;
        ut[(i-1)+(nx-2)*(j-1)+(nx-2)*(ny-2)*(k-1)][cs] = 0.0;
      }
    }
  }
  printf("internal cells zeroed \n");          
  for(i=1; i<nx-1; i++) /* inside boundary */
  {
    for(j=1; j<ny-1; j++) /* inside boundary */
    {
      for(k=1; k<nz-1; k++) /* inside boundary */
      {
        ii = (i-1)+(nx-2)*(j-1)+(nx-2)*(ny-2)*(k-1); /* row of matrix */
        /* for each term */
          /* coefdxx * d^2 u(x,y,z)/dx^ */
          coefdxx = 1.0;
          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]*u(xg[m],yg[j],zg[k]);
	    else ut[ii][mjk] = ut[ii][mjk] + cxx[m];
	  }

          /* coefdyy * d^2 u(x,y,z)/dy^2 */
          coefdyy = 2.0;
          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]*u(xg[i],yg[m],zg[k]);
	    else ut[ii][imk] = ut[ii][imk] + cyy[m];
 	  }

          /* coefdzz * d^2 u(x,y,z)/dz^2 */
          coefdzz = 3.0;
          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]*u(xg[i],yg[j],zg[m]);
	    else ut[ii][ijm] = ut[ii][ijm] + czz[m];
 	  }

          /* coefdxdy * d^2 u(x,y,z)/dxdy */
          coefdxdy = 4.0;
          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]*u(xg[kx],yg[ky],zg[k]);
	      else ut[ii][kxky] = ut[ii][kxky] + cxy[kx][ky];
            }
	  }

          /* coefdxdz * d^2 u(x,y,z)/dxdz */
          coefdxdz = 5.0;
          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]*u(xg[kx],yg[j],zg[kz]);
	      else ut[ii][kxkz] = ut[ii][kxkz] + cxz[kx][kz];
            }
	  }

          /* coefdydz * d^2 u(x,y,z)/dydz */
          coefdydz = 6.0;
          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]*u(xg[i],yg[ky],zg[kz]);
	      else ut[ii][kykz] = ut[ii][kykz] + cyz[ky][kz];
            }
	  }

          /* coefdx * d u(x,y,z)/dx */
          coefdx = 7.0;
          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]*u(xg[m],yg[j],zg[k]);
	    else ut[ii][mjk] = ut[ii][mjk] + cx[m];
	  }

          /* coefdy * d u(x,y,z)/dy */
          coefdy = 8.0;
          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]*u(xg[i],yg[m],zg[k]);
	    else ut[ii][imk] = ut[ii][imk] + cy[m];
	  }

          /* coefdz * d u(x,y,z)/dz */
          coefdz = 9.0;
          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]*u(xg[i],yg[j],zg[m]);
	    else ut[ii][ijm] = ut[ii][ijm] + cz[m];
	  }

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

          /* RHS constant */
          ut[ii][cs] = ut[ii][cs] + c(xg[i],yg[j],zg[k]);
        /* end terms */
      }
    }
  }
} /* end init_matrix */

static void simeq(int n, int m, double B[504][505], 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("pdenu23_eq.c running\n");
  printf("xmin=%g, xmax=%g, nx=%d \n", xg[0], xg[nx-1], nx);
  printf("ymin=%g, ymax=%g, ny=%d \n", yg[0], yg[ny-1], ny);
  printf("zmin=%g, zmax=%g, nz=%d \n", zg[0], zg[nz-1], nz);
  printf("u(xg[1],yg[1],zg[1])=%g \n",u(xg[1],yg[1],zg[1]));
  printf("c(xg[1],yg[1],zg[1])=%g \n",c(xg[1],yg[1],zg[1]));

  init_matrix();
  printf("matrix initialized\n");
  printf("initial matrix, left side, upper \n");
  for(i=0; i<7; i++)
  {
    for(j=0; j<7; j++)
    {
      printf("%10.3f ",ut[i][j]);
    }
    printf("\n");
  } 
  printf("initial matrix, right side, upper \n");
  for(i=0; i<7; i++)
  {
    for(j=498; j<505; j++)
    {
      printf("%10.3f ",ut[i][j]);
    }
    printf("\n");
  } 
  printf("initial matrix, left side, lower \n");
  for(i=497; i<504; i++)
  {
    for(j=0; j<7; j++)
    {
      printf("%10.3f ",ut[i][j]);
    }
    printf("\n");
  } 
  printf("initial matrix, right side, lower \n");
  for(i=497; i<504; i++)
  {
    for(j=498; j<505; j++)
    {
      printf("%10.3f ",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 pdenu23_eq.c */