/* pde3_eq.c  just checking second order approximation      */
/*      u(x,y) = x^3  + y^3  + 2x^2y + 3xy^2 + 4x + 5y + 6  */
/*       du/dx = 3x^2 + 0    + 4xy   + 3y^2  + 4  + 0  + 0  */
/*   d^2u/dx^2 = 6x   + 0    + 4y    + 0     + 0  + 0  + 0  */
/*       du/dy = 0    + 3y^2 + 2x^2  + 6xy   + 0  + 5  + 0  */
/*   d^2u/dy^2 = 0    + 6y   + 0     + 6x    + 0  + 0  + 0  */
/*                                                          */
/* solve d^2u/dx^2 + d^2u/dy^2 = 12x + 10y = c(x,y)         */
/* by second order method on square -1,-1 to 1,1            */

/* second order numerical difference                        */
/* d^2u/dx^2 + d^2u/dy^2 =                                  */
/* 1/h^2(u(-1h,0)+u(1h,0)+u(0,-1h)+u(0,1h)-4u(0,0)) +O(h^2) */
/*  = c(x,y)                                                */
/* solve system of linear equations for                     */
/* u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1)-4u(i,j) = c(x,y)*h^2 */
/* see pde3_eq.txt for development                          */

#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 nn = 81; /* (nx-2)*(ny-2) */
static double ut[81][82]; /* simultaneous equations built here      */
   /* (nx-2)*(ny-2) by (nx-2)*(ny-2) internal, non boundary cells   */
   /* position is (i-1)+(nx-2)*(j-1) for internal solution u(i,j)   */
   /* [][(nx-2)*(ny-2)] is constant,y  as in  A x = y, solve for x  */
static double us[81];
static const double xmin = -1.0;
static const double ymin = -1.0;
static const double xmax =  1.0;
static const double ymax =  1.0;
static double h;   /* x  and  y step = (xmax-xmin)/(n-1) = (ymax-ymin) */
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 check */
{
  return x*x*x + y*y*y + 2.0*x*x*y + 3.0*x*y*y + 4.0*x + 5.0*y +6.0;
}

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


static void init_matrix()
{
  int i, j, ii, jj, k;
  
  k = (nx-2)*(ny-2); /* constant column subscript */
  
  /* zero internal cells */
  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, h=%g\n", h);          
  for(i=1; i<nx-1; i++) /* inside boundary */
  {
    for(j=1; j<ny-1; j++)
    {
      ii = (i-1)+(nx-2)*(j-1); /* row of matrix */
      printf("i=%d, j=%d, ii=%d, ", i, j, ii);
      ut[ii][ii] = -4.0;
      ut[ii][k] = h*h*c(xmin+i*h, ymin+j*h);
      if((j+1)<(ny-1))
      {
        jj = (i-1)+(nx-2)*((j-1)+1);
        ut[ii][jj] = 1.0;
        printf("jj+=%d, ", jj);
      }
      else
        ut[ii][k] = ut[ii][k] - u(xmin+i*h, ymin+(j+1)*h);
      if((j-1)>0)
      {
        jj = (i-1)+(nx-2)*((j-1)-1);
        ut[ii][jj] = 1.0;
        printf("jj-=%d, ", jj);
      }
      else
         ut[ii][k] = ut[ii][k] - u(xmin+i*h, ymin+(j-1)*h);
      if((i+1)<(nx-1))
      {
         jj = ((i-1)+1)+(nx-2)*(j-1);
         ut[ii][jj] = 1.0;
         printf("jji+=%d, ", jj);
      }
      else
         ut[ii][k] = ut[ii][k] - u(xmin+(i+1)*h, ymin+j*h);
      if((i-1)!=0)
      {
         jj = ((i-1)-1)+(nx-2)*(j-1);
         ut[ii][jj] = 1.0;
         printf("jji-=%d", jj);
      }
      else
         ut[ii][k] = ut[ii][k] - u(xmin+(i-1)*h, ymin+j*h);
      printf("\n");
    }
  }
}

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(xmin+i*h,ymin+j*h);
      err = abs(us[(i-1)+(nx-2)*(j-1)]-uexact);
      error = error + err;
      if(err>maxerror) maxerror=err;
    }
  }
}

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

static void printexact()
{
  int i, j;
  
  printf("\n exact solution\n");
  for(i=1; i<min(nx-1,11); i++)
  {
    for(j=1; j<min(ny-1,11); j++)
    {
      printf("%6.3f ", u(xmin+i*h,ymin+j*h));
    }
    printf("\n");
  }
  printf("\n");
}

static void printdiff()
{
  int i, j;
  
  printf("\n difference exact-computed solution\n");
  printf("\n"); 
  for(i=1; i<min(nx-1,10); i++)
  {
    for(j=1; j<min(ny-1,10); j++)
    {
      printf("%7.4f ", u(xmin+i*h,ymin+j*h)-us[(i-1)+(nx-2)*(j-1)]);
    }
    printf("\n");
  }
  printf("\n");
}

static void simeq(int n, int m, double B[81][82], double X[])
{
  
  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=0;
  
  printf("pde3_eq.c running\n");
  h = (xmax-xmin)/(nx-1); /* step size in x */
  /*= (ymax-ymin)/(ny-1);    step size in y */
  printf("xmin=%g, xmax=%g, nx=%d, h=%g \n", xmin, xmax, nx, h);
  printf("ymin=%g, ymax=%g, ny=%d, h=%g \n", ymin, ymax, ny, h);
  
  init_matrix();
  printf("matrix initialized\n");
  printf(" initial matrix, upper left \n");
  for(i=0; i<11; i++)
  {
    for(j=0; j<10; j++)
    {
      printf("%5.2f ",ut[i][j]);
    }
    printf("=%5.2f \n", ut[i][nn]);
  } 
  printf("\n initial matrix, lower right \n");
  for(i=nn-11; i<nn; i++)
  {
    for(j=nn-11; j<nn+1; j++)
    {
      printf("%5.2f ",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;
}