/* pde2sin_sparse_gmp.c    2d uniform grid boundary value PDE */
/*                     known solution for testing method  */
/*                     u(x,y) = sin(x-y)                  */
/*                                                        */
/*                     differential equation to solve     */
/*                     du/dx + du/dy=0                    */
/* uniform grid on rectangle 0 to 2Pi, 0 to 2Pi           */
/* set up and solve system of linear equations            */

/* maximum error measured against analytic solution
 * nx  ny  maxerror
 * 13  12  1.40e-4     same as nx=12, ny=13, they can not be equal
 * 15  14  4.38e-6
 * 17  16  1.02e-7
 * 19  18  1.84e-9
 * 21  20  2.70e-11
 * 23  22  3.23e-13
 * 25  24  3.23e-15
 * 31  30  1.26e-21
 */

#include <stdlib.h>
int digits = 100; /* external */
#include "mpf_math.h"
#include "mpf_deriv.h"
#include "mpf_sparse.h"

/* do NOT make nx=ny, this creates a singular matrix  */
/* including boundary at 0 and nx-1, 0 and ny-1    */
#define nx 13
#define ny 12
#define neqn ((nx-2)*(ny-2))
#define cs neqn
static mpf_t xmin, xmax, ymin, ymax, hx, hy;
static mpf_t twopi;

static mpf_t xg[nx];
static mpf_t yg[ny];
/* nx-2 by ny-2  internal, boundary cells */
static mpf_t cx[nx];
static mpf_t cy[ny];
static mpf_t coefdx, coefdy;
static mpf_t us[neqn]; /* solution vector */
/* static mpf_t ut[neqn][neqn+1];  includes RHS, now sparse B */
static struct sparse B;

static void u(mpf_t result, mpf_t x, mpf_t y)
{
  /* solution for boundary and optionalcheck */
  mpf_t xmy; /* x-y */

  mpf_init(xmy);
  mpf_sub(xmy, x, y);
  mpf_sin(result, xmy);
} /* end u */

static void c(mpf_t result, mpf_t x, mpf_t y)
{
  /* RHS of differential equation to solve */
  mpf_set_si(result, 0);
} /* end c */

static int S(int i, int j)
{
  if(i<1 || i>nx-1 || j<1 || j>ny-1) gmp_printf("BAD S i=%2d, j=%2d\n", i, j);
  return (i-1)*(ny-2)+(j-1); /* index into us and ut */
} /* end S */

static void print()
{
  int i, j;
  mpf_t error;
  mpf_t max_error;
  mpf_t avg_error;
  mpf_t val, tmp;


  mpf_init_set_si(error, 0);
  mpf_init_set_si(max_error, 0);
  mpf_init_set_si( avg_error, 0);
  mpf_init(val);
  mpf_init(tmp);
  gmp_printf("\n");
  gmp_printf("exact solution u, computed us, error\n");
  for(i=0; i<nx; i++)
  {
    for(j=0; j<ny; j++)
    {
      gmp_printf("xg[%2d]=%7.3Fg, ", i, xg[i]);
      gmp_printf("yg[%2d]=%7.3Fg, ", j, yg[j]);
      u(val, xg[i], yg[j]);
      gmp_printf("u=%7.3Fg, ", val);
      if(i==0 || i==nx-1 || j==0 || j==ny-1) /* boundary */
      {
        u(val, xg[i], yg[j]);
        mpf_set_ui(error, 0);
      }
      else
      {
        u(val, xg[i], yg[j]);
        mpf_sub(val, val, us[S(i,j)]);
        mpf_abs(error, val);
        mpf_add(avg_error, avg_error, error);
      }
      if(mpf_cmp(error, max_error)>0) mpf_set(max_error, error);
      gmp_printf("us=%Fg \n", val);
      gmp_printf("err=%Fg \n", error);
    }
  }
  mpf_div_ui(val, avg_error, neqn);
  gmp_printf("avg_error=%Fg \n", val);
  gmp_printf("max_error=%Fg \n", max_error);
} /* end print */

static void check_row(int row) /* only for checking */
{
  /* row is in ut solution coordinates */
  int i, j, ii;
  mpf_t sum, val, tmp, bval;

  mpf_init_set_si(sum, 0);
  mpf_init(val);
  mpf_init(tmp);
  mpf_init(bval);

  for(i=1; i<nx-1; i++)
  {
    for(j=1; j<ny-1; j++)
    {
      ii = S(i,j);
      u(val, xg[i], yg[j]); /* u is check solution */
      sparseGet(bval, B, row, ii);
      mpf_mul(val, val, bval);
      mpf_add(sum, sum, val); 
    }
  }
  sparseGet(bval, B, row, cs);
  mpf_sub(sum, sum, bval);
  mpf_abs(val, sum);
  mpf_set_d(tmp, 0.001);
  if(mpf_cmp(sum,tmp)>0)
    gmp_printf("BAD row=%4d, err=%Fg\n", row, sum);
} /* end check_row */

static void init_matrix()
{
  int i, j, ii, jj;
  int k, kj, ik;
  mpf_t val;

  mpf_init(val);
  /* cs is RHS, constant column subscript */
  /* zero internal cells, not with sparse */
  for(i=0; i<neqn; i++) mpf_init(us[i]);

  for(i=1; i<nx-1; i++) /* inside boundary */
  {
    for(j=1; j<ny-1; j++)
    {
      ii = S(i,j); /* equation */
      /* for each term */

      /* coefdx * d u(x,y)/dx */
      mpf_init_set_si(coefdx, 1);
      mpf_rderiv(1,nx,i,hx,cx);
      for(k=0; k<nx; k++)
      {
        mpf_mul(cx[k], cx[k], coefdx);
        if(k==0 || k==nx-1)
	{
          u(val, xg[k], yg[j]);
          mpf_mul(val, val, cx[k]);
          /* mpf_sub(ut[ii][cs], ut[ii][cs], val); */
          mpf_neg(val, val);
          sparseAdd(B, ii, cs, val);
	}
        else
	{
          kj = S(k,j);
          /* mpf_add(ut[ii][kj], ut[ii][kj], cx[k]); */
          sparseAdd(B, ii, kj, cx[k]);
        }
      } /* k */
      /* coefdy * d u(x,y)/dy */
      mpf_init_set_si(coefdy, 1);
      mpf_rderiv(1,ny,j,hy,cy);
      for(k=0; k<ny; k++)
      {
        mpf_mul(cy[k], cy[k], coefdy);
        if(k==0 || k==ny-1)
	{
          u(val, xg[i],yg[k]);
          mpf_mul(val, val, cy[k]);
	  /* mpf_sub(ut[ii][cs], ut[ii][cs], val); */
          mpf_neg(val, val);
          sparseAdd(B, ii, cs, val); 
	}
        else
	{
          ik = S(i,k);
          /* mpf_add(ut[ii][ik], ut[ii][ik], cy[k]); */
          sparseAdd(B, ii, ik, cy[k]);
        }
      } /* k */
      /* RHS constant */
      /* ut[ii][cs] = ut[ii][cs] + c(xg[i],yg[j]); */
      /* end terms */
      check_row(ii);
    } /* j */
  } /* i */
} /* end init_matrix */

static void check_soln()  /* check when solution unknown */
{
  int i, ii, j, jj;
  mpf_t val, smaxerr, err, x, y, ux, uy;
  /*  ux(x,y) + uy(x,y) = c(x,y) */
  gmp_printf("check_soln against PDE\n");
  mpf_init_set_si(smaxerr, 0);
  mpf_init(val);
  mpf_init(err);
  mpf_init(x);
  mpf_init(y);
  mpf_init(ux);
  mpf_init(uy);

  for(i=1; i<nx-1; i++) /* through dof equations */
  {
    mpf_rderiv(1,nx,i,hx,cx);
    mpf_set(x, xg[i]);
    for(ii=1; ii<ny-1; ii++)
    {
      mpf_rderiv(1,ny,ii,hy,cy);
      mpf_set(y, yg[ii]);
      c(val, x, y);
      mpf_neg(err, val); /* add other side of equation terms */
      mpf_set_si(ux, 0);
      for(j=0; j<nx; j++) /* compute numerical derivative at (i,ii) */
      {
        if(j==0 || j==nx-1)
	{
          u(val, xg[j], yg[ii]); /* boundary */
	  mpf_mul(val, val, cx[j]);
	  mpf_add(ux, ux, val);
	}
        else
	{
          mpf_mul(val, cx[j], us[S(j,ii)]); 
          mpf_add(ux, ux, val); 
	}
      }
      mpf_set_si(uy, 0);
      for(jj=0; jj<ny; jj++)
      {
        if(jj==0 || jj==ny-1)
	{
          u(val, xg[i], yg[jj]); /* boundary */
          mpf_mul(val, val, cy[jj]); 
	  mpf_add(uy, uy, val);
	}
        else
	{
          mpf_mul(val, cy[jj], us[S(i,jj)]); 
          mpf_add(uy, uy, val); 
	}
      }
      mpf_add(err, err, ux);
      mpf_add(err, err, uy);
      mpf_abs(val, err);
      if(mpf_cmp(val, smaxerr)>0) mpf_set(smaxerr, val);
    } /* ii Y */
  } /* i X */
  gmp_printf("check_soln max error=%Fg \n", smaxerr);
} /* end Check_Soln */

int main(int argc, char *argv[])
{
  int i, j;
  mpf_t tmp, val;

  /* pde2sin_sparse  */
  gmp_printf("pde2sin_sparse_gmp.c   running\n");
  gmp_printf("differential equation to solve\n");
  gmp_printf("du/dx + du/dy = c(x,y)\n");
  gmp_printf("c(x,y)=0\n");
  gmp_printf("uniform grid on rectangle 0,2Pi to 0,2Pi\n");
  gmp_printf("known Solution, for testing method\n");
  gmp_printf("u(x,y) = sin(x-y)\n");
  gmp_printf("do not make nx=ny, this creates a singular matrix\n");

  mpf_init(tmp);
  mpf_init(val);
  mpf_init(twopi);
  mpf_twopi(twopi);
  mpf_init_set_si(xmin, 0);
  mpf_init_set(xmax, twopi);
  mpf_init_set_si(ymin, 0);
  mpf_init_set(ymax, twopi);
  mpf_init(hx);
  mpf_sub(hx, xmax, xmin);
  mpf_div_ui(hx, hx, nx-1);
  mpf_init(hy);
  mpf_sub(hy, ymax, ymin);
  mpf_div_ui(hy, hy, ny-1);

  for(i=0; i<nx; i++)
  {
    mpf_init(xg[i]);
    mpf_init(cx[i]);
    mpf_mul_ui(tmp, hx, i);
    mpf_add(xg[i], xmin, tmp);
    gmp_printf("xg[%2d]=%Fg\n", i, xg[i]);
  }
  for(j=0; j<ny; j++)
  {
    mpf_init(yg[j]);
    mpf_init(cy[j]);
    mpf_mul_ui(tmp, hy, j);
    mpf_add(yg[j], ymin, tmp);
    gmp_printf("yg[%2d]=%Fg\n", j, yg[j]);
  }
  gmp_printf("xmin=%Fg \n", xmin);
  gmp_printf("xmax=%Fg \n", xmax);
  gmp_printf("hx=%Fg \n", hx);
  gmp_printf("nx=%3d \n", nx);
  gmp_printf("ymin=%Fg \n", ymin);
  gmp_printf("ymax=%Fg \n", ymax);
  gmp_printf("hy=%Fg \n",hy);
  gmp_printf("ny=%3d\n", ny);
  u(val, xg[0], yg[0]);
  gmp_printf("u(xg[0],yg[0])=%Fg \n", val);
  c(val, xg[0], yg[0]);
  gmp_printf("c(xg[0],yg[0])=%Fg \n", val);

  /* debug stuff */
  gmp_printf("debug print\n");
  mpf_rderiv(1,nx,2,hx,cx);
  for(i=0; i<nx; i++) gmp_printf("rderiv 1,2,%2d cx[%2d]=%Fg\n", nx, i, cx[i]);
  mpf_rderiv(1,ny,2,hy,cy);
  for(j=0; j<ny; j++) gmp_printf("rderiv 1,2,%2d cy[%2d]=%Fg\n", ny, j, cy[j]);

  sparseInit(&B, neqn);
  init_matrix();
  gmp_printf("matrix initialized\n");

  /* debug stuff */
  gmp_printf("matrix B includes RHS\n");
  sparseWrite_All(B);

  sparseSimeq(B, us); /* returned us */
  check_soln();
  print();
  gmp_printf("pde2sin_sparse_gmp.c finished\n");
  return 0;
} /* end pde2sin_sparse_gmp.c */