/* pdenu_eq.c  2D non uniform grid boundary value PDE  non iterative  */
/*             known solution to test method                          */
/*             u(x,y) = exp(x*y)*sin(x+y)                             */
/*                                                                    */
/*             differential equation to solve                         */
/* d^2u/dx^2 + 2*d^2u/dy^2 - 3*d^2u/dxdy + y*du/dx - x*du/dy = c(x,y) */
/* exp(x*y)*sin(x+y)*(x*x + 2*y*y -3*x*y -3)                          */
/* non uniform grid on rectangle -1.0,-1.1 to 1.2,1.3                 */
/*                                                                    */
/* set up and solve system of linear equations                        */
/*                                                                    */
/* create two dimensional array with                                  */
/* (nx-2)*(ny-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 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};
                             /* nx-2 by ny-2  internal, non boundary cells */
static double cxx[11];       /* nuderiv coefficients */
static double cyy[10];
static double cxy[11][10];
static double cx[11];
static double cy[10];
static double coefdxx, coefdyy, coefdxdy, coefdx, coefdy, coefu;

static double us[72];        /* solution being computed 9*8 */
static double ut[72][73];    /* 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) /* for boundary and optional check */
{
  return exp(x*y)*sin(x+y);
}

static double c(double x, double y) /* RHS of differential equation to solve */
{
  return exp(x*y)*sin(x+y)*(x*x + 2.0*y*y -3.0*x*y -3.0);
}

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

static void printus()
{
  int i, j;
  
  printf("computed solution\n");
  for(i=1; i<nx-1; i++)
  {
    for(j=1; j<ny-1; j++)
    {
      printf("xg[%d]=%f,yg[%d]=%f, us=%g \n",
             i, xg[i], j, xg[j], us[(i-1)+(nx-2)*(j-1)]);
    }
  }
  printf("\n");
} /* end printus */

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

static void printdiff()
{
  int i, j, k;
  
  printf("difference exact-computed solution\n");
  for(i=1; i<nx-1; i++)
  {
    printf(" ");
    for(j=1; j<min(ny-1,10); j++)
    {
      printf("%8.3f", u(xg[i],yg[j])-us[(i-1)+(nx-2)*(j-1)]);
    }
    printf("\n");
  }
  printf("\n");
} /* end printdiff */

static void init_matrix()
{
  int i, j, k, ii, jj, ik, kj, kx, ky, kxky, cs;

  cs = (nx-2)*(ny-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(ii=1; ii<nx-1; ii++)
       for(jj=1; jj<ny-1; jj++)
         ut[(i-1)+(nx-2)*(j-1)][(ii-1)+(nx-2)*(jj-1)] = 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 */
    {
      ii = (i-1)+(nx-2)*(j-1); /* row of matrix */
      /* for each term */
        /* coefdxx * d^2 u(x,y)/dx^2 */
        coefdxx = 1.0;
        nuderiv(2,nx,i,xg,cxx);
        for(k=0; k<nx; k++)
	{
          cxx[k] = cxx[k] * coefdxx;
          kj = (k-1)+(nx-2)*(j-1);
	  if(k==0 || k==nx-1)
               ut[ii][cs] = ut[ii][cs] - cxx[k]*u(xg[k],yg[j]);
	  else ut[ii][kj] = ut[ii][kj] + cxx[k];
	}

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

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

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

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

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

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

static void simeq(int n, int m, double B[72][73], 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("pdenu_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("u(xg[0],yg[0])=%g \n", u(xg[0],yg[0]));
  printf("c(xg[0],yg[0])=%g \n", c(xg[0],yg[0]));

  init_matrix();
  printf("matrix initialized\n");
  printf("initial matrix, left side, upper \n");
  for(i=0; i<9; i++)
  {
    for(j=0; j<9; j++)
    {
      printf("%8.3f",ut[i][j]);
    }
    printf("\n");
  } 
  printf("initial matrix, right side, upper \n");
  for(i=0; i<9; i++)
  {
    for(j=64; j<73; j++)
    {
      printf("%8.3f",ut[i][j]);
    }
    printf("\n");
  } 
  printf("initial matrix, left side, lower \n");
  for(i=63; i<72; i++)
  {
    for(j=0; j<9; j++)
    {
      printf("%8.3f",ut[i][j]);
    }
    printf("\n");
  } 
  printf("initial matrix, right side, lower \n");
  for(i=63; i<72; i++)
  {
    for(j=64; j<73; j++)
    {
      printf("%8.3f",ut[i][j]);
    }
    printf("\n");
  } 
  printf("\n");
  
  simeq((nx-2)*(ny-2), (nx-2)*(ny-2)+1, ut, us);
  printexact();
  printus();
  printdiff();

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