/* pde4sin_nueqsp.c  4d uniform grid boundary value PDE */
/*                     known solution for testing method  */
/*                     u(x,y) = sin(x-y)                  */
/*                                                        */
/*                differential equation to solve          */
/*                du/dx + du/dy + du/dz + du/dt =0        */
/* uniform grid on rectangle 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
 */

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
int digits = 50; /* external */
#include "mpf_math.h"
#include "mpf_nuderiv.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 10
#define ny 9
#define nz 8
#define nt 7
#define neqn ((nx-2)*(ny-2)*(nz-2)*(nt-2))
#define cs neqn
static mpf_t xmin, xmax, ymin, ymax, zmin, zmax, tmin, tmax;
static mpf_t hx, hy, hz, ht;
static mpf_t twopi;

static mpf_t xg[nx];  /* index i */
static mpf_t yg[ny];  /* index ii */
static mpf_t zg[nz];  /* index iii */
static mpf_t tg[nt];  /* index iiii */
/* nx-2 by ny-2 by nz-2 by nt-2 internal, boundary cells */
static mpf_t cx[nx][nx];
static mpf_t cy[ny][ny];
static mpf_t cz[nz][nz];
static mpf_t ct[nt][nt];
static mpf_t coefdx, coefdy, coefdz, coefdt;
static mpf_t us[neqn]; /* solution vector */

/*static mpf_t ut[neqn][neqn+1]; includes RHS */
static struct sparse A;

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

  mpf_init(xmy);
  mpf_sub(xmy, x, y);
  mpf_add(xmy, xmy, z);
  mpf_sub(xmy, xmy, t);

  mpf_sin(result, xmy);
} /* end u */

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

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

static void print()
{
  int i, ii, iii, iiii;
  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=1; i<nx-1; i++)
  {
    for(ii=1; ii<ny-1; ii++)
    {
      for(iii=1; iii<nz-1; iii++)
      {
        for(iiii=1; iiii<nt-1; iiii++)
        {
          gmp_printf("xg[%2d]=%7.3Fg, ", i, xg[i]);
          gmp_printf("yg[%2d]=%7.3Fg, ", ii, yg[ii]);
          gmp_printf("zg[%2d]=%7.3Fg, ", iii, zg[iii]);
          gmp_printf("tg[%2d]=%7.3Fg, ", iiii, tg[iiii]);
          u(val, xg[i], yg[ii], zg[iii], tg[iiii]);
          gmp_printf("u=%7.3Fg, ", val);
          mpf_sub(val, val, us[S(i,ii,iii,iiii)]);
          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 dsetup()
{
  int i, ii, iii, iiii;

  for(i=0; i<nx; i++)
  {
    mpf_nuderiv(1,nx,i,xg,cx[i]); /* makes cx[about_point][j] */
  }
  for(ii=0; ii<ny; ii++)
  {
    mpf_nuderiv(1,ny,ii,yg,cy[ii]);
  }
  for(iii=0; iii<nz; iii++)
  {
    mpf_nuderiv(1,nz,iii,zg,cz[iii]);
  }
  for(iiii=0; iiii<nt; iiii++)
  {
    mpf_nuderiv(1,nt,iiii,tg,ct[iiii]);
  }
} /* end dsetup */


static void check_row(int row) /* only for checking, when solution known */
{
  /* row is in ut solution coordinates */
  int i, ii, iii, iiii, kk;
  mpf_t sum, val, tmp, utval;

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

  for(i=1; i<nx-1; i++)
  {
    for(ii=1; ii<ny-1; ii++)
    {
      for(iii=1; iii<nz-1; iii++)
      {
        for(iiii=1; iiii<nt-1; iiii++)
        {
          kk = S(i,ii,iii,iiii);
          u(val, xg[i], yg[ii], zg[iii], tg[iiii]); /* u is check solution */
          sparseGet(utval, A, row, kk);
          mpf_mul(val, val, utval);
          mpf_add(sum, sum, val);
	}
      } 
    }
  }
  sparseGet(utval, A, row, cs);
  mpf_sub(sum, sum, utval);
  mpf_abs(val, sum);
  mpf_set_d(tmp, 0.01);
  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, iii, jjj, iiii, jjjj;
  int k, kk, ki, kii, kiii, kiiii;
  mpf_t val, utval;

  mpf_init(val);
  mpf_init(utval);
  /* cs is RHS, constant column subscript */

  /* zero internal cells */
  sparseInit(&A, neqn);

  gmp_printf("internal cells zeroed \n");
  for(i=1; i<nx-1; i++) /* inside boundary */
  {
    for(ii=1; ii<ny-1; ii++)
    {
      for(iii=1; iii<nz-1; iii++)
      {
        for(iiii=1; iiii<nt-1; iiii++)
        {
          kk = S(i,ii,iii,iiii); /* equation */
          /* for each term */

          /* coefdx * d u(x,y,z,t)/dx */
          mpf_set_si(coefdx, 1);
          for(k=0; k<nx; k++)
          {
            mpf_mul(cx[i][k], cx[i][k], coefdx);
            if(k==0 || k==nx-1) /* boundary */
	    {
              u(val, xg[k], yg[ii], zg[iii], tg[iiii]);
              mpf_mul(val, val, cx[i][k]);
              /* mpf_sub(ut[kk][cs], ut[kk][cs], val); */
              mpf_neg(val,val);
              sparseAdd(A, kk, cs, val);
	    }
            else
	    {
              ki = S(k,ii,iii,iiii);
              /* mpf_add(ut[kk][ki], ut[kk][ki], cx[i][k]); */
              sparseAdd(A, kk, ki, cx[i][k]);
            }
          } /* k */

          /* coefdy * d u(x,y,z,t)/dy */
          mpf_set_si(coefdy, 1);
          for(k=0; k<ny; k++)
          {
            mpf_mul(cy[ii][k], cy[ii][k], coefdy);
            if(k==0 || k==ny-1) /* boundary */
	    {
              u(val, xg[i],yg[k],zg[iii],tg[iiii]);
              mpf_mul(val, val, cy[ii][k]);
	      /* mpf_sub(ut[kk][cs], ut[kk][cs], val); */
              mpf_neg(val,val);
              sparseAdd(A, kk, cs, val); 
	    }
            else
	    {
              kii = S(i,k,iii,iiii);
              /* mpf_add(ut[kk][kii], ut[kk][kii], cy[ii][k]); */
              sparseAdd(A, kk, kii, cy[ii][k]);
            }
          } /* k */

          /* coefdz * d u(x,y,z,t)/dz */
          mpf_set_si(coefdz, 1);
          for(k=0; k<nz; k++)
          {
            mpf_mul(cz[iii][k], cz[iii][k], coefdz);
            if(k==0 || k==nz-1) /* boundary */
	    {
              u(val, xg[i],yg[ii],zg[k],tg[iiii]);
              mpf_mul(val, val, cz[iii][k]);
	      /* mpf_sub(ut[kk][cs], ut[kk][cs], val); */
              mpf_neg(val,val);
              sparseAdd(A, kk, cs, val);
	    }
            else
	    {
              kiii = S(i,ii,k,iiii);
              /* mpf_add(ut[kk][kiii], ut[kk][kiii], cz[iii][k]); */
              sparseAdd(A, kk, kiii, cz[iii][k]);
            }
          } /* k */

          /* coefdt * d u(x,y,z,t)/dt */
          mpf_set_si(coefdt, 1);
          for(k=0; k<nt; k++)
          {
            mpf_mul(ct[iiii][k], ct[iiii][k], coefdt);
            if(k==0 || k==nt-1) /* boundary */
	    {
              u(val, xg[i],yg[ii],zg[iii],tg[k]);
              mpf_mul(val, val, ct[iiii][k]);
	      /* mpf_sub(ut[kk][cs], ut[kk][cs], val); */
              mpf_neg(val,val);
              sparseAdd(A, kk, cs, val);
	    }
            else
	    {
              kiiii = S(i,ii,iii,k);
              /* mpf_add(ut[kk][kiiii], ut[kk][kiiii], ct[iiii][k]); */
              sparseAdd(A, kk, kiiii, ct[iiii][k]);
            }
          } /* k */

          /* RHS constant */
          /* ut[kk][cs] = ut[kk][cs] + c(xg[i],yg[ii],zg[iii],tg[iiii]); */
          /* end terms */
          check_row(kk);
        } /* iiii */
      } /* iii */
    } /* ii */
  } /* i */
} /* end init_matrix */

static void check_soln()  /* check when solution unknown */
{
  int i, ii, iii, iiii, j, jj, jjj, jjjj;
  mpf_t val, smaxerr, err, x, y, z, t, ux, uy, uz, ut;
  /*  ux(x,y,z,t) + uy(x,y,z,t) + uz(x,y,z,t) + ut(x,y,z,t) = c(x,y) */
  gmp_printf("check_soln against PDE\n");
  fflush(stdin);
  mpf_init_set_si(smaxerr, 0);
  mpf_init(val);
  mpf_init(err);
  mpf_init(x);
  mpf_init(y);
  mpf_init(z);
  mpf_init(t);
  mpf_init(ux);
  mpf_init(uy);
  mpf_init(uz);
  mpf_init(ut);

  for(i=1; i<nx-1; i++) /* through dof equations */
  {
    mpf_set(x, xg[i]);
    for(ii=1; ii<ny-1; ii++)
    {
      mpf_set(y, yg[ii]);
      for(iii=1; iii<nz-1; iii++) 
      {
        mpf_set(x, zg[iii]);
        for(iiii=1; iiii<nt-1; iiii++)
        {
          mpf_set(t, tg[iiii]);

	  c(val, x, y, z, t);
          mpf_neg(err, val); /* add other side of equation terms */
          mpf_set_si(ux, 0);
          /* compute numerical derivative at (i,ii,iii,iiii) */
          for(j=0; j<nx; j++)
          {
            if(j==0 || j==nx-1)
	    {
              u(val, xg[j], yg[ii], zg[iii], tg[iiii]); /* boundary */
	      mpf_mul(val, val, cx[i][j]);
	      mpf_add(ux, ux, val);
	    }
            else
	    {
              mpf_mul(val, cx[i][j], us[S(j,ii,iii,iiii)]); 
              mpf_add(ux, ux, val); 
	    }
          } /* j for ux */

          mpf_set_si(uy, 0);
          for(jj=0; jj<ny; jj++)
          {
            if(jj==0 || jj==ny-1)
	    {
              u(val, xg[i], yg[jj], zg[iii], tg[iiii]); /* boundary */
              mpf_mul(val, val, cy[ii][jj]); 
	      mpf_add(uy, uy, val);
	    }
            else
	    {
              mpf_mul(val, cy[ii][jj], us[S(i,jj,iii,iiii)]); 
              mpf_add(uy, uy, val); 
	    }
          }  /* jj for uy */

          mpf_set_si(uz, 0);
          for(jjj=0; jjj<nz; jjj++)
          {
            if(jjj==0 || jjj==nz-1)
	    {
              u(val, xg[i], yg[ii], zg[jjj], tg[iiii]); /* boundary */
              mpf_mul(val, val, cz[iii][jjj]); 
	      mpf_add(uz, uz, val);
	    }
            else
	    {
              mpf_mul(val, cz[iii][jjj], us[S(i,ii,jjj,iiii)]); 
              mpf_add(uz, uz, val); 
	    }
          }  /* jjj for uz */

          mpf_set_si(ut, 0);
          for(jjjj=0; jjjj<nt; jjjj++)
          {
            if(jjjj==0 || jjjj==nt-1)
	    {
              u(val, xg[i], yg[ii], zg[iii], tg[jjjj]); /* boundary */
              mpf_mul(val, val, ct[iiii][jjjj]); 
	      mpf_add(ut, ut, val);
	    }
            else
	    {
              mpf_mul(val, ct[iiii][jjjj], us[S(i,ii,iii,jjjj)]); 
              mpf_add(ut, ut, val); 
	    }
          }  /* jjjj for ut */

          mpf_add(err, err, ux);
          mpf_add(err, err, uy);
          mpf_add(err, err, uz);
          mpf_add(err, err, ut);
          mpf_abs(val, err);
          if(mpf_cmp(val, smaxerr)>0) mpf_set(smaxerr, val);
        } /* iiii T */
      } /* iii Z */
    } /* ii Y */
  } /* i X */
  gmp_printf("check_soln max error=%Fg \n", smaxerr);
} /* end Check_Soln */

int main(int argc, char *argv[])
{
  int i, ii, iii, iiii, j, jj, jjj, jjjj;
  mpf_t tmp, val;
  double time_start, time_now;

  /* pde4sin_nueq  */
  mpf_set_default_prec(digits*3.32);
  time_start = (double)clock()/(double)CLOCKS_PER_SEC;
  gmp_printf("pde4sin_nueqsp.c  running, digits=%d\n", digits);
  gmp_printf("differential equation to solve\n");
  gmp_printf("du/dx + du/dy + du/dz + du/dt= c(x,y,z,t)\n");
  gmp_printf("c(x,y,z,t)=0\n");
  gmp_printf("uniform grid on rectangle 0 to 2Pi in all dimensions\n");
  gmp_printf("known Solution, for testing method\n");
  gmp_printf("u(x,y,z,t) = sin(x-y+z-t)\n");
  gmp_printf("do not make nx=ny or any pair, this creates a singular matrix\n");

  mpf_init(tmp);
  mpf_init(val);
  mpf_init(coefdx);
  mpf_init(coefdy);
  mpf_init(coefdz);
  mpf_init(coefdt);

  mpf_init(twopi);
  mpf_halfpi(twopi); /* two pi parallel */
  mpf_init_set_si(xmin, 0);
  mpf_init_set(xmax, twopi);
  mpf_init_set_si(ymin, 0);
  mpf_init_set(ymax, twopi);
  mpf_init_set_si(zmin, 0);
  mpf_init_set(zmax, twopi);
  mpf_init_set_si(tmin, 0);
  mpf_init_set(tmax, 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);
  mpf_init(hz);
  mpf_sub(hz, zmax, zmin);
  mpf_div_ui(hz, hz, nz-1);
  mpf_init(ht);
  mpf_sub(ht, tmax, tmin);
  mpf_div_ui(ht, ht, nt-1);

  for(i=0; i<nx; i++)
  {
    mpf_init(xg[i]);
    for(j=0; j<nx; j++) mpf_init(cx[i][j]);
    mpf_mul_ui(tmp, hx, i);
    mpf_add(xg[i], xmin, tmp);
    gmp_printf("xg[%2d]=%Fg\n", i, xg[i]);
  }
  for(ii=0; ii<ny; ii++)
  {
    mpf_init(yg[ii]);
    for(jj=0; jj<ny; jj++) mpf_init(cy[ii][jj]);
    mpf_mul_ui(tmp, hy, ii);
    mpf_add(yg[ii], ymin, tmp);
    gmp_printf("yg[%2d]=%Fg\n", ii, yg[ii]);
  }
  for(iii=0; iii<nz; iii++)
  {
    mpf_init(zg[iii]);
    for(jjj=0; jjj<nz; jjj++) mpf_init(cz[iii][jjj]);
    mpf_mul_ui(tmp, hz, iii);
    mpf_add(zg[iii], zmin, tmp);
    gmp_printf("zg[%2d]=%Fg\n", iii, zg[iii]);
  }
  for(iiii=0; iiii<nt; iiii++)
  {
    mpf_init(tg[iiii]);
    for(jjjj=0; jjjj<nt; jjjj++) mpf_init(ct[iiii][jjjj]);
    mpf_mul_ui(tmp, ht, iiii);
    mpf_add(tg[iiii], tmin, tmp);
    gmp_printf("tg[%2d]=%Fg\n", iiii, tg[iiii]);
  }
  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);
  gmp_printf("zmin=%Fg \n", zmin);
  gmp_printf("zmax=%Fg \n", zmax);
  gmp_printf("hz=%Fg \n", hz);
  gmp_printf("nz=%3d \n", nz);
  gmp_printf("tmin=%Fg \n", tmin);
  gmp_printf("tmax=%Fg \n", tmax);
  gmp_printf("ht=%Fg \n",ht);
  gmp_printf("nt=%3d\n", nt);
  time_now = (double)clock()/(double)CLOCKS_PER_SEC - time_start;
  gmp_printf("initial printing, at %f seconds\n", time_now);

  dsetup();
  time_now = (double)clock()/(double)CLOCKS_PER_SEC - time_start;
  gmp_printf("derivatives set up, at %f seconds\n", time_now);

  init_matrix();
  time_now = (double)clock()/(double)CLOCKS_PER_SEC - time_start;
  gmp_printf("matrix initialized, at %f seconds\n", time_now);

  /* debug stuff */
  /* gmp_printf("matrix B includes RHS\n");
     for(i=0; i<neqn; i++)
       for(j=0; j<cs; j++)
         gmp_printf("i=%2d, j=%2d, B[i][j]=%Fg\n", i, j, ut[i][j]);
  */
  /* mpf_simeqb(neqn, (mpf_t *)ut, us); */
  for(i=0; i<neqn; i++) mpf_init(us[i]);
  sparseSimeq(A, us);
  time_now = (double)clock()/(double)CLOCKS_PER_SEC - time_start;
  gmp_printf("equations solved, at %f seconds\n", time_now);
  fflush(stdin);
  for(i=0; i<neqn; i++) gmp_printf("us[%d]=%Fg\n",i,us[i]);

  check_soln();
  time_now = (double)clock()/(double)CLOCKS_PER_SEC - time_start;
  gmp_printf("solution checked, at %f seconds\n", time_now);

  print();
  gmp_printf("sparse used %d of %d\n", sparseGetCellCount, neqn*(neqn+1));
  time_now = (double)clock()/(double)CLOCKS_PER_SEC - time_start;
  gmp_printf("pde4sin_nueqsp.c finished, at %f seconds\n", time_now);
  return 0;
} /* end pde4sin_nueqsp.c */