/* pde_abc_eq.c  see WEB page for development, CMSC 455 supplemental  */
/*               The PDE to be solved is:                             */
/*        a1(x,y)*uxx(x,y) + b1(x,y)*uxy(x,y) + c1(x,y)*uyy(x,y) +    */
/*        d1(x,y)*ux(x,y) + e1(x,y)*uy(x,y) + f1(x,y)*u(x,y) = c(x,y) */
/*                                                                    */
/* Method is high order discretization                                */
/* The functions and other problem information is in the files:       */
/*    abc.h  and   abc.c                                              */ 
/*                             compile with abc.c and deriv.c         */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "abc.h"
#include "deriv.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))

/* nx and ny  are defined in abc.h */
static double us[(nx-2)*(ny-2)];  /* solution being computed */
static double uu[(nx-2)*(ny-2)];  /* solution if checking */
static double ut[(nx-2)*(ny-2)][(nx-2)*(ny-2)+1];
                                   /* equations to be solved  */
                /* nx-2 by ny-2  internal, non boundary cells */
/* all elements 0..nx-1 by 0..ny-1 includes boundary elements */
/* xmin, xmax, ymin, ymax are defined in  abc.h  */
static double hx; /* x step = (xmax-xmin)/(nx-1) */
static double hy; /* y step = (ymax-ymin)/(ny-1) */
static double hx2, hy2, hxy; /* hx^2, hy^2, hx*hy */
static int cs;    /* index of constant vector */
static double pdxx[nd][nd], pdxy[nd][nd], pdyy[nd][nd];
static double pdx[nd][nd], pdy[nd][nd], pd[nd][nd];
static int md;  /* mid point correction  nd/2 */
/* ifcheck defined in  abc.h  non zero for checking             */
static double error;         /* sum of absolute exact-computed  */
static double avgerror;      /* average error                   */
static double maxerror;      /* maximum cell error              */

static void check_soln();

static void clear()
{
  int i, j, ii;

  for(i=0; i<nx-2; i++) /* solution coordinates, compressed */
  {
    for (j=0; j<ny-2; j++)
    {
      ii = i+(nx-2)*j;
      us[ii] = 0.0;
      if(ifcheck>0) uu[ii] = u(xmin+(i+1)*hx,ymin+(j+1)*hy);
    }
  }
} /* end clear */

static void printus() /* print the computed solution */
{
  int i, j, k;
  
  printf("\n computed solution\n");
  printf("us[x,y] \n");
  for(i=1; i<nx-1; i++)
  {
    for(k=1; k<ny-1; k=k+8)
    {
      for(j=k; j<min(ny-1,k+8); j++)
      {
        printf("%8.2f ", us[(i-1)+(nx-2)*(j-1)]);
      }
      printf("\n");
    }
  }
  printf("\n");
} /* end printus */

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

static void printexact() /* typically not known */
{
  int i, j;
  
  printf("\n exact solution\n");
  printf("u(x,y) \n");
  for(i=1; i<nx-1; i++)
  {
    for(j=1; j<min(ny-1,9); j++)
    {
      printf("%8.2f ", u(xmin+i*hx, ymin+j*hy));
    }
    printf("\n");
  }
  if(min(ny,9)<ny) printf("Not all columns printed \n");
  printf("\n");
} /* end printexact */

static void printdiff() /* only when checking */
{
  int i, j;
  
  printf("\n difference exact-computed solution\n");
  printf("\n");
  printf("u(x,y)-us[x,y] \n");
  for(i=1; i<nx-1; i++)
  {
    for(j=1; j<min(ny-1,9); j++)
    {
      printf("%8.2f", u(xmin+i*hx, ymin+j*hy)-
             us[(i-1)+(nx-2)*(j-1)]);
    }
    printf("\n");
  }
  printf("\n");
} /* end printdiff */

static void print_coef()
{
  int i,j;

  printf("pdxx     j=0     j=1     j=2     j=3     j=4  y\n");
  for(i=0; i<nd; i++)
  {
    printf("i=%d, ", i);
    for(j=0; j<nd; j++)
    {
      printf("%7.2f ", pdxx[i][j]);
    }
    printf("\n");
  }
  printf("x\n\n");

  printf("pdxy     j=0     j=1     j=2     j=3     j=4  y\n");
  for(i=0; i<nd; i++)
  {
    printf("i=%d, ", i);
    for(j=0; j<nd; j++)
    {
      printf("%7.2f ", pdxy[i][j]);
    }
    printf("\n");
  }
  printf("\n");

  printf("pdyy     j=0     j=1     j=2     j=3     j=4  y\n");
  for(i=0; i<nd; i++)
  {
    printf("i=%d, ", i);
    for(j=0; j<nd; j++)
    {
      printf("%7.2f ", pdyy[i][j]);
    }
    printf("\n");
  }
  printf("\n");

  printf("pdx      j=0     j=1     j=2     j=3     j=4  y\n");
  for(i=0; i<nd; i++)
  {
    printf("i=%d, ", i);
    for(j=0; j<nd; j++)
    {
      printf("%7.2f ", pdx[i][j]);
    }
    printf("\n");
  }
  printf("\n");

  printf("pdy      j=0     j=1     j=2     j=3     j=4  y\n");
  for(i=0; i<nd; i++)
  {
    printf("i=%d, ", i);
    for(j=0; j<nd; j++)
    {
      printf("%7.2f ", pdy[i][j]);
    }
    printf("\n");
  }
  printf("\n");

} /* end print_coef */

static void build_coef(int ix, int iy, int ic)
{
  int i,j;
  int a[25], bb;

  if((nd&1)==0 || nd<3)
    { printf("nd must be odd and >1 , bye.\n"); exit(1);}
  if(nd!=5) { printf("init_matrix needs enhancement, bye.\n"); exit(1);}

  /* build discretization matrices */
  if(ifcheck) printf("build_coef for point ix=%d, iy=%d \n", ix, iy);
  for(i=0; i<nd; i++)
  {
    for(j=0; j<nd; j++)
    {
      pdxx[i][j]=0.0; /* coefficient of a1(x,y) */
      pdxy[i][j]=0.0; /* coefficient of b1(x,y) */
      pdyy[i][j]=0.0; /* coefficient of c1(x,y) */
      pdx [i][j]=0.0; /* coefficient of d1(x,y) */
      pdy [i][j]=0.0; /* coefficient of e1(x,y) */
      pd  [i][j]=0.0; /* coefficient of f1(x,y) */
    }
  }
  pd[ix][iy] = 1.0;
  deriv(2, nd, ix, a, &bb); /* for pdxx */
  if(ifcheck>3) printf("deriv pdxx returns bb=%d, aa=%d %d %d %d %d\n",
		     bb, a[0], a[1], a[2], a[3], a[4]);
  for(i=0; i<nd; i++)
    pdxx[i][iy] = (double)a[i]/((double)bb*hx2);

  deriv(2, nd, iy, a, &bb); /* for pdyy */
  if(ifcheck>3) printf("deriv pdyy returns bb=%d, aa=%d %d %d %d %d \n",
		     bb, a[0], a[1], a[2], a[3], a[4]);
  for(j=0; j<nd; j++)
    pdyy[ix][j] = (double)a[j]/((double)bb*hy2);

  deriv(1, nd, ix, a, &bb); /* for pdx */
  if(ifcheck>3) printf("deriv pdx returns bb=%d, aa=%d %d %d %d %d \n",
		      bb, a[0], a[1], a[2], a[3], a[4]);
  for(i=0; i<nd; i++)
    pdx[i][iy] = (double)a[i]/((double)bb*hx);

  deriv(1, nd, iy, a, &bb); /* for pdy */
  if(ifcheck>3) printf("deriv pdy returns bb=%d, aa=%d %d %d %d %d \n",
		      bb, a[0], a[1], a[2], a[3], a[4]);
  for(j=0; j<nd; j++)
    pdy[ix][j] = (double)a[j]/((double)bb*hy);

  for(i=0; i<nd; i++)
    for(j=0; j<nd; j++)
      pdxy[i][j] = pdx[i][iy]*pdy[ix][j];

  if(ifcheck>2) print_coef();
} /* end build_coef */

static void check_row(int row)
{
  int i, j, ii;
  double sum = 0.0;

  for(i=0; i<nx-2; i++)
  {
    for(j=0; j<ny-2; j++)
    {
      ii = i+(nx-2)*j;
      sum += ut[row][ii]*uu[ii];
    }
  }
  sum -= ut[row][cs];
  if(abs(sum)>0.001) printf("row=%d is bad err=%g \n", row, sum);
} /* end check_row */

static void build_row(int i, int j, int bi, int bj)
{
  /* i,j are in element coordinates, ir,jr are all element coord */
  int ki, kj, ii, jj, ij, ir, jr, ia, ja;
  double  p, a1c, b1c, c1c, d1c, e1c, f1c, cc;

  ir = i;
  jr = j;
  ii = (i-1)+(nx-2)*(j-1); /* row of matrix */
      if(ifcheck>2)
      {
         printf("working row %d, i=%d, j=%d, bi=%d, bj=%d  \n",
                             ii, i, j, bi, bj);
         /* printf("x=%g, y=%g \n", xmin+ir*hx, ymin+jr*hy); */
      }
      a1c = a1(xmin+ir*hx,ymin+jr*hy);
      b1c = b1(xmin+ir*hx,ymin+jr*hy);
      c1c = c1(xmin+ir*hx,ymin+jr*hy);
      d1c = d1(xmin+ir*hx,ymin+jr*hy);
      e1c = e1(xmin+ir*hx,ymin+jr*hy);
      f1c = f1(xmin+ir*hx,ymin+jr*hy);
      cc  =  c(xmin+ir*hx,ymin+jr*hy);
      if(ifcheck>5)
      {
        printf("a1c=%7.3f, b1c=%7.3f, c1c=%7.3f \n", a1c, b1c, c1c);
        printf("d1c=%7.3f, d1c=%7.3f, cc =%7.3f \n", d1c, e1c, cc);
      }
      for(ki=0; ki<nd; ki++)
      {
        for(kj=0; kj<nd; kj++)
	{
          ia = ir-bi+ki-md; /* real i,j in all element coordinates */
          ja = jr-bj+kj-md; /* may be boundary */
          ij = (ia-1)+(nx-2)*(ja-1); /* in ut coordinates */
          p =  a1c*pdxx[ki][kj];
          p += b1c*pdxy[ki][kj];
          p += c1c*pdyy[ki][kj];
          p += d1c* pdx[ki][kj];
          p += e1c* pdy[ki][kj];
          p += f1c*  pd[ki][kj];
          if(ia<0)    printf("***-error in ia=%d \n", ia);
          if(ia>nx-1) printf("***+error in ia=%d \n", ia);
          if(ja<0)    printf("***-error in ja=%d \n", ja);
          if(ja>ny-1) printf("***+error in ja=%d \n", ja);
          if(ia==0)  /* boundary value tests */
	  {
            ut[ii][cs] -= p * u(xmin+(ia)*hx, ymin+(ja)*hy);
          }
          else if(ia==nx-1)
            ut[ii][cs] -= p * u(xmin+(ia)*hx, ymin+(ja)*hy);
          else if(ja==0)
            ut[ii][cs] -= p * u(xmin+(ia)*hx, ymin+(ja)*hy);
          else if(ja==ny-1)
            ut[ii][cs] -= p * u(xmin+(ia)*hx, ymin+(ja)*hy);
          else /* non boundary case */
	  {
            ut[ii][ij] = p;
	    if(ifcheck>5) printf("ij=%d, p=u[ii][ij]=%g \n", ij, ut[ii][ij]);
	  } 
	} 
      }
      ut[ii][cs] += cc; 
      printf("ii=%d, ut[ii][ii]=%g, ut[ii][cs]=%g \n\n",
             ii, ut[ii][ii], ut[ii][cs]);
      if(ifcheck>0) check_row(ii);
} /* end build_row */

static void init_matrix()
{
  int i, j, ii, jj;

  /* zero internal cells, the matrix does not include boundary 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 \n");          
  build_coef(md, md, md);
  for(i=2; i<nx-2; i++) /* compute standard central cells */
  {
    for(j=2; j<ny-2; j++)
    {
      build_row(i, j, 0, 0);
    } /* end j loop */
  } /* end i loop */
  printf("standard central cells generated \n");

  /* do four sides, one element in. special cases */
  build_coef(md-1, md, md);
  for(j=2; j<ny-2; j++) /* do i=1 */
  {
    i = 1;
    if(ifcheck>2) printf("working row edge i=%d, j=%d \n", i, j);
    build_row(i, j, -1, 0);
  } /* end j loop on edge */

  build_coef(md+1, md, md);
  for(j=2; j<ny-2; j++) /* do  i=nx-2 */
  {
    i = nx-2;
    if(ifcheck>2) printf("working row edge i=%d, j=%d \n", i, j);
    build_row(i, j, 1, 0);
  } /* end j loop on edge */

  build_coef(md, md-1, md);
  for(i=2; i<nx-2; i++) /* do j=1 */
  {
    j = 1;
    if(ifcheck>2) printf("working row edge i=%d, j=%d \n", i, j);
    build_row(i, j, 0 , -1);
  } /* end i loop on edge */

  build_coef(md, md+1, md);
  for(i=2; i<nx-2; i++) /* do  j=ny-2 */
  {
    j = ny-2;
    if(ifcheck>2) printf("working row edge i=%d, j=%d \n", i, j);
    build_row(i, j, 0 , 1);
  } /* end i loop on edge */
  printf("four sides, one element in, special case, generated \n");
  /* now do four corners */
  i=1;
  j=1;
  build_coef(md-1, md-1, md);
  build_row(i, j, -1, -1);
  i=1;
  j=ny-2;
  build_coef(md-1, md+1, md);
  build_row(i, j, -1, 1);
  i=nx-2;
  j=1;
  build_coef(md+1, md-1, md);
  build_row(i, j, 1, -1);
  i=nx-2;
  j=ny-2;
  build_coef(md+1, md+1, md);
  build_row(i, j, 1, 1);
} /* end init_matrix */

static void print_mat()
{
  int i, j, k, ii;

  printf(" initial matrix \n");
  for(i=0; i<(nx-2)*(ny-2); i++)
  {
    printf("i=%d, j=0 .. %d \n", i, (nx-2)*(ny-2));
    for(k=0; k<=(nx-2)*(ny-2); k=k+nx-2)
    {
      for(j=0; j<nx-2 && j+k<=(nx-2)*(ny-2); j++)
      {
        printf("%10.3f ",ut[i][j+k]);
      }
      printf("\n");
    }
    printf("\n");
  } 
  printf("\n");
} /* end print_mat */

static void check_soln() /* against PDE, added later, used rderiv and xg,yg */
{
               /* for all interior points */
  double x, y, err, maxerr, avgerr, rmserr;
  double sumerr, sumsqerr;
  double cx[nx][nx];
  double cxx[nx][nx];
  double cy[ny][ny];
  double cyy[ny][ny];
  double tcxy[ny];
  double dxx, dxy, dyy, dx, dy; /* a1, b1, c1, d1, e1 coefficients */ 
  double xg[nx]; /* may be available */
  double yg[ny]; 
  double ug[nx][ny]; /* full, including boundaries, copied from us[] */
  int cnt = 0;
  int i, ii, j, jj, k;

  if(ifcheck>5) printf("check_soln setup \n");
  for(i=0; i<nx; i++) xg[i] = xmin + (double)i*hx;
  for(j=0; j<ny; j++) yg[j] = ymin + (double)j*hy;
  /* transform us[] into ug, add boundaries, make code simpler */
  for(i=1; i<nx-1; i++) /* get from us[] */
  {
    for(j=1; j<ny-1; j++)
    {
      ii = (i-1)+(nx-2)*(j-1); /* ii for  us[ii] computed solution at i,j */
      ug[i][j] = us[ii];
      if(ifcheck>5) printf("ug[%d][%d]=%g \n", i, j, ug[i][j]);
    }
  }
  for(i=0; i<nx; i++) /* fill boundaries, loops below easier */
  {
    ug[i][0]    = u(xg[i],yg[0]);
    ug[i][ny-1] = u(xg[i],yg[ny-1]);
  }
  for(j=0; j<ny; j++) /* fill boundaries */
  {
    ug[0][j]    = u(xg[0],   yg[j]);
    ug[nx-1][j] = u(xg[nx-1],yg[j]);
  }

  for(i=0; i<nx; i++)
  {
    rderiv(1,nx,i,hx,cx[i]); /* fills row of derivative matrix */
    rderiv(2,nx,i,hx,cxx[i]);
  }
  for(j=0; j<ny; j++)
  {
    rderiv(1,ny,j,hy,cy[j]);
    rderiv(2,ny,j,hy,cyy[j]);
  }
  if(ifcheck>5)
  {
    printf("cx=\n");
    for(i=0; i<nx; i++)
    {
      for(k=0; k<nx; k++)
      {
        printf("%8.3f ",cx[i][k]);
      }
      printf("\n");
    }
    printf("\n");
    printf("cxx=\n");
    for(i=0; i<nx; i++)
    {
      for(k=0; k<nx; k++)
      {
        printf("%8.3f ",cxx[i][k]);
      }
      printf("\n");
    }
    printf("\n");
    printf("cy=\n");
    for(j=0; j<ny; j++)
    {
      for(k=0; k<ny; k++)
      {
        printf("%8.3f ",cy[j][k]);
      }
      printf("\n");
    }
    printf("\n");
    printf("cyy=\n");
    for(j=0; j<ny; j++)
    {
      for(k=0; k<ny; k++)
      {
        printf("%8.3f ",cyy[j][k]);
      }
      printf("\n");
    }
    printf("\n");
  } /* end ifcheck */

  sumerr = 0.0;
  sumsqerr = 0.0;
  maxerr = 0.0;
  for(i=1; i<nx-1; i++) /* check full equation, at every point */
  {
    x = xg[i];
    for(j=1; j<ny-1; j++)
    {
      y = yg[j];
      /* compute derivatives, plug into PDE, subtract RHS, should be zero */
      dxx = 0.0;
      dx  = 0.0;
      for(k=0; k<nx; k++)
      {
        /* PDE  a1(x,y)*uxx(x,y) term */
	dxx += cxx[i][k]*ug[k][j];
        /* PDE  d1(x,y)*ux(x,y) term */
	dx += cx[i][k]*ug[k][j];
      } 
      dyy = 0.0;
      dy  = 0.0;
      for(k=0; k<ny; k++)
      {
        /* PDE c1(x,y)*uyy(x,y) term */
	dyy += cyy[j][k]*ug[i][k];
        /* PDE e1(x,y)*uy(x,y) term */
	dy += cy[j][k]*ug[i][k];
      }
      /* PDE b1(x,y)*uxy(x,y) term */
      dxy = 0.0;
      for(jj=0; jj<ny; jj++) /* walk jj in y, computing uxx */
      {
	tcxy[jj] = 0.0; /* numeric  ux part of uxy */
        for(k=0; k<nx; k++)
        {
	  tcxy[jj] += cx[i][k]*ug[k][jj];
	}
      }
      for(k=0; k<ny; k++)  /* uy part of uxy */
      {
        dxy += cy[j][k]*tcxy[k];  /* PDE uxy term */
      }
      /* check */
      if(ifcheck>5)
      {
        printf("numerical dxx=%g, dxy=%g,dyy=%g\n", dxx, dxy, dyy);
        printf("analytic  dxx=%g, dxy=%g,dyy=%g\n",
	        6.0*x+6.0*y, 6.0*x+8.0*y+5.0, 12.0*y+8.0*x);
        printf("numerical dx=%g , dy=%g ,u=%g\n", dx, dy, ug[i][j]);
        printf("analytic  dx=%g , dy=%g ,u=%g\n",
               3.0*x*x+6.0*x*y+4.0*y*y+5.0*y+6.0,
               6.0*y*y+3.0*x*x+8.0*x*y+5.0*x+7.0,
               u(x,y));
      }
      err = a1(xg[i],yg[j])*dxx +
	    b1(xg[i],yg[j])*dxy +
            c1(xg[i],yg[j])*dyy +
            d1(xg[i],yg[j])*dx  +
            e1(xg[i],yg[j])*dy  +
            f1(xg[i],yg[j])*ug[i][j] - /* PDE u term at i,j*/
	     c(xg[i],yg[j]);           /* RHS */
      err = abs(err);
      sumerr += err;
      sumsqerr += err*err;
      maxerr = max(err,maxerr);
      cnt++;  
    } /* end j */
  } /* end i */
  rmserr = sqrt(sumsqerr/(double)cnt);
  avgerr = sumerr/(double)cnt;
  printf("check_soln values against PDE \n");
  printf("maxerr=%g, rmserr=%g, avgerr=%g \n",maxerr,rmserr,avgerr);
} /* end check_soln */ 

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, ii, jj;
  
  printf("pde_abc_eq.c running with user files  abc.h abc.c \n");
  md = nd/2; /* mid point correction */
  hx = (xmax-xmin)/(nx-1); /* step size in x */
  hy = (ymax-ymin)/(ny-1); /* step size in y */
  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);
  hx2 = hx*hx;
  hy2 = hy*hy;
  hxy = hx*hy;
  cs = (nx-2)*(ny-2); /* constant column subscript */

  clear();
  init_matrix();
  printf("matrix initialized\n");
  if(ifcheck) print_mat();

  simeq((nx-2)*(ny-2), (nx-2)*(ny-2)+1, ut, us);

  printus();
  check_soln(); /* do for testing code, actual below in else part */

  if(ifcheck)
  {
    printexact(); /* us vector is the solution */
    printdiff();
    check();
    avgerror = error/((nx-2)*(ny-2));
    printf("check against known analytic solution \n");
    printf("total error=%g, average error=%g, max error=%g\n",
           error, avgerror, maxerror);
  }
  else
  {
    /* use derivative computation on solution and check against PDE */
    check_soln();
  }
  return 0;
} /* end main of pde_abc_eq.c */