/* pdenu24_eq.c  2D non uniform grid boundary value PDE              */
/*               known solution to test method                       */
/*               u(x,y) = sin(x)+cos(y)                              */
/*                                                                   */
/*               differential equation to solve                     */
/* du/dy + d^4u/dx^4 + d^2u/dx^2 - du/dx  =  -sin(y)-cos(x)          */
/*                                                                   */
/* non uniform grid on rectangle 0.0,0.1 to 1.0,1.0                */
/*                                                                   */
/* 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))

#define nx 17          /* includes boundary */
#define ny 17          /* includes boundary */

/*
static double xg[nx] = {0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0};
static double yg[ny] = {0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0};
*/
static double xg[nx]; /* includes boundary */
static double yg[ny];
static double cx[nx]; /* nuderiv coefficients */
static double cxx[nx];
static double cxxxx[nx];
static double cy[ny];
static double coefdy, coefdx, coefdxx, coefdxxxx;

                            /* nx-2 by ny-2  internal, non boundary cells */
static double us[(nx-2)*(ny-2)];                  /* solution being computed */
static double ut[(nx-2)*(ny-2)][(nx-2)*(ny-2)+1]; /* 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 sin(x)+cos(y);
}

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

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, err=%g \n",
             i, xg[i], j, yg[j], us[(i-1)+(nx-2)*(j-1)],
             us[(i-1)+(nx-2)*(j-1)]-u(xg[i],yg[j]));
    }
  }
  printf("\n");
} /* end printus */


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 */
        /* coefdy * d u(x,y)/dy */
        coefdy = 1.0;
        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];
	}

        /* coefdx * d u(x,y)/dx */
        coefdx = -1.0;
        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];
	}

        /* 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];
	}

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

        /* 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[(nx-2)*(ny-2)][(nx-2)*(ny-2)+1],
                  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;
  double dx, dy;
  
  printf("pdenu24_eq.c running\n");
  xg[0] = 0.0;
  xg[nx-1] = 1.0;
  yg[0] = 0.0;
  yg[ny-1] = 1.0;
  dx = (xg[nx-1]-xg[0])/(double)(nx-1);
  dy = (yg[ny-1]-yg[0])/(double)(ny-1);
  for(i=1; i<nx; i++) {xg[i]=xg[0]+i*dx; printf("xg[%d]=%f\n",i,xg[i]);}
  for(j=1; j<ny; j++) {yg[j]=yg[0]+j*dy; printf("yg[%d]=%f\n",j,yg[j]);}

  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[1],yg[1])=%g \n", u(xg[1],yg[1]));
  printf("c(xg[1],yg[1])=%g \n", c(xg[1],yg[1]));

  init_matrix();
  printf("matrix initialized\n");
  
  simeq((nx-2)*(ny-2), (nx-2)*(ny-2)+1, ut, us);
  printus();
  printdiff();

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