/* fempt_bihar2dps_la.c  Galerkin FEM biharmonic 4th order
 *                       using pthreads
 *                                                                   
 * solve  uxxxx(x,y) + uyyyy(x,y) + 2*uxxyy(x,y) = f(x,y)
 *
 *  f(x,y) := 4.0*sin(x)*sin(y)
 *
 *  in the rectangle  xmin< x <xmax and  ymin< y <ymax  ,
 *
 *     subject to boundary conditions of the first
 *     kind, Dirichlet
 *
 *  Dirichlet boundary conditions are computed using:
 *  ub(x,y) := sin(x)*sin(y)
 *
 *  uses laphi.c  simeq.c  deriv.c  gaulegf.c 
 */

/* Linux with glibc:
 *   _REENTRANT to grab thread-safe libraries
 *   _POSIX_SOURCE to get POSIX semantics
 */
#ifdef __linux__
#  define _REENTRANT
#  define _POSIX_SOURCE
#  define _P __P
#endif

#include <pthread.h>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <time.h>
#include "simeq.h"
#include "laphi.h"
#include "deriv.h"
#include "gaulegf.h"

#undef  abs
#define abs(a) ((a)<0.0?(-(a)):(a))
#undef  max
#define max(a,b) ((a)>(b)?(a):(b))

#define nx  9
#define ny  9
#define npx 9
#define npy 9

static double xmin = 0.0; /* problem parameters */
static double xmax = 3.14159; 
static double ymin = 0.0;
static double ymax = 3.14159;
static double start, now, prev; /* timing */

/* nx, ny, xg, yg needed by galk or galf, available at call*/
static double xg[nx]; /* X grig */
static double yg[ny]; /* Y grid */
static double hx, hy;
static double ug[nx*ny]; /* solution */
static double kg[nx*ny][nx*ny];
static double fg[nx*ny];
/* npx, npy  Gauss Legendre integration order */
static double xx[npx];
static double wx[npx];
static double yy[npy];
static double wy[npy];


/* function prototypes */
double F(double x, double y); /* RHS */
double ub(double x, double y);
/* Galerkin k function for this problem */
double galk(double x, double y, int i, int j, int ii, int jj); 
/* Galerkin f function for this problem */
double galf(double x, double y, int i, int ii); 
void write_soln();
void check_soln();
void *compute(void *arg);

int main(int argc, char * argv[])
{
    int nxy = nx*ny;
    double x, y;
    double Ua[nx*ny];
    int stat;
    int worker;
    int ids[nx];                          /* worker id number */
    pthread_t threads[nx];                /* holds thread info */
    int errcode;                          /* holds pthread error code */
    int *status;                          /* holds return code */
    double pnow, wnow;                    /* time CPU and wall */
    int i, j, ii, jj;
    double val, maxerr, err, avgerr;


    printf("fempt_bihar2dps_la.c running.\n");
    now = (double)time(NULL);
    start = now;
    prev = now;
    printf("solve  uxxxx(x,y) + 2*uxxyy(x,y) + uyyyy(x,y) = f(x,y)\n");
    printf(" \n");
    printf("f(x,y) := 4.0*sin(x)*sin(y) \n");
    printf(" \n");
    printf("in the rectangle   xmin< x <xmax  and  ymin< y <ymax,\n");
    printf(" \n");
    printf("   subject to boundary conditions of the first kind, Dirichlet\n");
    printf("   ub(x,y) := sin(x)*sin(y) \n");
    printf("xmin=%7.4f, xmax=%7.4f, ymin=%7.4f, ymax=%7.4f\n",
           xmin, xmax, ymin, ymax);
    printf("nx=%d, ny=%d\n",nx, ny);
    printf("x grid and analytic solution at ymin \n");
    hx = (xmax-xmin)/(double)(nx-1);
    for(i=0; i<nx; i++)
    {
      xg[i] = xmin+i*hx;
      printf("i=%d, xg[i]=%f, ub(xg[i],ymin)=%f\n",
             i, xg[i], ub(xg[i],ymin));
    }  
    printf(" \n");
    printf("y grid and analytic solution at xmin \n");
    hy = (ymax-ymin)/(double)(ny-1);
    for(ii=0; ii<ny; ii++)
    {
      yg[ii] = ymin+ii*hy;
      printf("ii=%d, yg[ii]=%f, ub(xmin,yg[ii])=%f \n",
             ii, yg[ii], ub(xmin, yg[ii]));
    }  
    printf(" \n");
    /* put solution in Ua vector */ 
    for(i=0; i<nx; i++)
    {
      x = xmin+(double)i*hx;
      for(ii=0; ii<ny; ii++)
      {
	y = ymin+(double)ii*hy;
        Ua[i*ny+ii] = ub(x,y);
      }
    }
    printf(" \n");

    /* initialize k and f */ 
    for(i=0; i<nx; i++)
    {
      for(j=0; j<nx; j++)
      {
        for(ii=0; ii<ny; ii++)
        {
          for(jj=0; jj<ny; jj++)
          {
            kg[(i*ny+ii)][(j*ny+jj)] = 0.0;
	  }
        }
      }
    }
    for(i=0; i<nx; i++)
    {
      for(ii=0; ii<ny; ii++)
      {
        fg[i*ny+ii] = 0.0;
      }
    }
    /* set boundary conditions, four sides */ 
    for(i=0; i<nx; i++)
    {
      x = xg[i];
      ii=0;
      y = yg[ii];
      kg[(i*ny+ii)][(i*ny+ii)] = 1.0;
      fg[i*ny+ii] = ub(x,y);
      ii=ny-1;
      y = yg[ii];
      kg[(i*ny+ii)][(i*ny+ii)] = 1.0;
      fg[i*ny+ii] = ub(x,y);
    }
    for(ii=0; ii<ny; ii++)
    {
      y = yg[ii];
      i=0;
      x = xg[i];
      kg[(i*ny+ii)][(i*ny+ii)] = 1.0;
      fg[i*ny+ii] = ub(x,y);
      i=nx-1;
      x = xg[i];
      kg[(i*ny+ii)][(i*ny+ii)] = 1.0;
      fg[i*ny+ii] = ub(x,y);
    }

    printf("calling gauleg xmin=%f, xmax=%f, npx=%d \n",
           xmin, xmax, npx);
    gaulegf(xmin, xmax, xx, wx, npx); /* xx and wx set up for integrations */ 
    printf("xx[1]=%f, xx[2]=%f, wx[1]=%f, wx[2]=%f \n",
           xx[1],xx[2],wx[1],wx[2]);
    printf("calling gauleg ymin=%f, ymax=%f, npy=%d \n",
           ymin,ymax,npy);
    gaulegf(ymin, ymax, yy, wy, npy); /* yy and wy set up for integrations */
    printf("yy[1]=%f, yy[2]=%f, wy[1]=%f, wy[2]=%f \n",
           yy[1],yy[2],wy[1],wy[2]);
    i=1; ii=1;  j=1; jj=1;
    val = galk(xx[2],yy[2],i,j,ii,jj);
    printf("galk(xx[2],yy[2],i,j,ii,jj)=%g\n",val);
    val = galf(xx[2],yy[2],i,ii);
    printf("galf(xx[2],yy[2],i,ii)=%g\n",val);

    now = (double)time(NULL);
    fprintf(stderr, "time to initialize = %f seconds \n", now-prev);
    printf("time to initialize = %f seconds \n", now-prev);
    prev = now;

    /* this is the main thread's control code, no significant computation */
    printf("starting threads \n");

    pnow = (double)clock()/(double)CLOCKS_PER_SEC;
    wnow = (double)time(NULL);

    /* create the threads */
    for (worker=1; worker<nx-1; worker++) /* outter loop */
    {
      ids[worker]=worker;  /* we are passing an address, 'worker' may change */
      if (errcode=pthread_create(&threads[worker], /* thread struct          */
		       NULL,                    /* default thread attributes */
		       compute,                 /* start thread              */
		       &ids[worker]))           /* arg to thread             */
      {
        printf("worker &d failed create: %s\n", worker, strerror(errcode));
        fflush(stdout);
        exit(1);
      }
    }
    /* wait for all threads to exit */
    printf("wait for workers to build their part of stiffness matrix \n");
    fflush(stdout);
    for (worker=1; worker<nx-1; worker++)
    {
      if (errcode=pthread_join(threads[worker],(void *) &status))
      { 
        printf("worker &d failed join: %s\n", worker, strerror(errcode));
        fflush(stdout);
        /* exit(1); */
      }
      /* check thread's exit status and release its resources */
      if (*status != worker)
      {
        printf("worker %d terminated abnormally \n", worker);
        fflush(stdout);
        /* exit(1); */
      }
    }
    pnow = (double)clock()/(double)CLOCKS_PER_SEC - pnow;
    printf("total CPU time is %f seconds \n", pnow);
    fflush(stdout);
    wnow = (double)time(NULL) - wnow;
    printf("wall time is %f seconds \n", wnow);
    printf("average number of processors used = %f \n", now/wnow);

    printf("k computed stiffness matrix, see above \n ");
    printf("f computed forcing function, see above \n ");

    now = (double)time(NULL);
    fprintf(stderr, "time for stiffness and forcing = %f seconds \n", now-prev);
    printf("time for stiffness and forcing = %f seconds \n", now-prev);
    prev = now;
    simeq(nx*ny, (double *)kg, fg, ug);
    now = (double)time(NULL);
    fprintf(stderr, "time for solve = %f seconds \n", now-prev);
    printf("time for solve = %f seconds \n", now-prev);
    prev = now;

    write_soln(); /* to file x y soln */

    check_soln(); /* independent of nowing analytic solution */

    avgerr = 0.0;
    maxerr = 0.0;
    printf("ug computed Galerkin, Ua analytic, error \n");
    for(i=0; i<nx; i++)
    {
      for(ii=0; ii<ny; ii++)
      {
        err = abs(ug[i*ny+ii]-Ua[i*ny+ii]);
        if(err>maxerr) maxerr = err;
        avgerr = avgerr + err;
        /* printf("ug[%d,%d]=%f, Ua=%f, err=%g \n", */
	/*       i,ii,ug[i*ny+ii],Ua[i*ny+ii],ug[i*ny+ii]-Ua[i*ny+ii]); */
      }
    }
    printf("xmax=%f, ymax=%f, nx=%d, ny=%d, npx=%d, npy=%d \n",
           xmax,ymax,nx,ny,npx,npy);
    printf("                  maxerr=%g, avgerr=%g \n",
           maxerr,avgerr/(double)(nx*ny));
    printf("\n");

  now = (double)time(NULL);
  fprintf(stderr, "total wall time = %f seconds \n", now-start);
  printf("total wall time = %f seconds \n", now-start);
  printf("calling gnuplot \n");
  fflush(stdout);

  stat = execl("/usr/bin/gnuplot","-persist","fempt_bihar2dps_la_c.plot", NULL);
  /* only returns upon error, no return on sucess */
  printf("returned from gnuplot with stat= %d \n", stat);
  return 0;
} /* end fempt_bihar2dps_la main */

/* PDE problem definition functions  f and ub */

double F(double x, double y)
{
  return 4.0*sin(x)*sin(y);
}

double ub(double x, double y)
{
  return sin(x)*sin(y);
}

/* Galerkin k stiffness function for this problem      */ 
/* solve  uxxxx(x,y) + 2*uxxyy(x,y) + uyyyy(x,y)  = f(x,y) */
double galk(double x, double y, int i, int j, int ii, int jj)
{           
  return ( phipppp(x,j,nx-1,xg)*     phi(y,jj,ny-1,yg) +
        2.0* phipp(x,j,nx-1,xg)*   phipp(y,jj,ny-1,yg) +
               phi(x,j,nx-1,xg)* phipppp(y,jj,ny-1,yg) )*
               phi(x,i,nx-1,xg)*     phi(y,ii,ny-1,yg);
} /* end galk used by integration */ 

/* Galerkin f function for this problem */
double galf(double x, double y, int i, int ii) 
{
  return F(x,y)*phi(x,i,nx-1,xg)*phi(y,ii,ny-1,yg);
} /* end galf used by integration */ 

void write_soln()
{
  FILE * fout;
  int i, ii;

  fout = fopen("fempt_bihar2dps_la_c.dat", "w");
  if(fout==NULL)
  {
    printf("file fempt_bihar2dps_la_c.dat can not be opened for writing\n");
    return;
  }
  printf("writing fempt_bihar2dps_la_c.dat \n");
  for(i=0; i<nx; i++)
  {
    for(ii=0; ii<ny; ii++)
    {
      fprintf(fout," %f %f %f \n",xg[i],yg[ii],ug[i*ny+ii]);
    } /* end y */
    fprintf(fout," \n");
  } /* end x */
  fflush(fout);
  fclose(fout);
  printf("finished writing fempt_bihar2dps_la_c.dat \n");
} /* end write_soln */

void check_soln()
{
               /* for all interior points */
  double x, y, err, maxerr, avgerr, rmserr;
  double sumerr, sumsqerr, eval_deriv;
  double cxx[nx][nx];
  double cxxxx[nx][nx];
  double cyy[ny][ny];
  double cyyyy[ny][ny];
  double tcxxyy[ny];
  int cnt = 0;
  int i, ii, j, jj, k, m;

  for(i=1; i<nx-1; i++)
  {
    rderiv(2,nx,i,hx,cxx[i]);
    rderiv(4,nx,i,hx,cxxxx[i]);
  }
  for(ii=1; ii<ny-1; ii++)
  {
    rderiv(2,ny,ii,hy,cyy[ii]);
    rderiv(4,ny,ii,hy,cyyyy[ii]);
  }

  sumerr = 0.0;
  sumsqerr = 0.0;
  maxerr = 0.0;
  for(i=1; i<nx-1; i++) /* full equation */
  {
    x = xmin + (double)i*hx;
    for(j=1; j<ny-1; j++)
    {
      y = ymin + (double)j*hy;
      eval_deriv = 0.0; /* no PDE u term */
        
      for(k=0; k<nx; k++)
      {
	eval_deriv += cxxxx[i][k]*ug[k*ny+j];   /* PDE  uxxxx term */
      } 
      for(k=0; k<ny; k++)
      {
	eval_deriv += cyyyy[j][k]*ug[i*ny+k];    /* PDE uyyyy term */
      }
      /*  + 2*uxxyy */
      for(jj=0; jj<ny; jj++) /* walk jj in y, computing uxx */
      {
	tcxxyy[jj] = 0.0; /* numeric  uxx part of uxxyy */
        for(k=0; k<nx; k++)
        {
	  tcxxyy[jj] += cxx[i][k]*ug[k*ny+jj];
	}
      }
      for(k=0; k<ny; k++)  /* uyy part of uxxyy */
      {
        eval_deriv += 2.0*cyy[j][k]*tcxxyy[k];  /* PDE uxxyy term */
      }

      err = eval_deriv-F(x,y);
      err = abs(err);
      sumerr += err;
      sumsqerr += err*err;
      maxerr = max(err,maxerr);
      cnt++;  
    } /* end y */
  } /* end x */
  rmserr = sqrt(sumsqerr/(double)cnt);
  avgerr = sumerr/(double)cnt;
  printf("check_soln values \n");
  printf("maxerr=%g, rmserr=%g, avgerr=%g \n",maxerr,rmserr,avgerr);
} /* end check_soln */ 


/* this is the main computation thread */
void *compute(void *arg)
{
  int i=*(int *)arg;
  double val;
  int j, ii, jj, px, py;

  printf("compute %d: running  \n", i);
  fflush(stdout); /* a good idea inside a thread */
  /* compute entries in stiffness matrix */ 
  /*  for(i=1; i<nx-1; i++) thread */
  /*  {                     thread */
    for(j=0; j<nx; j++)
    {
      for(ii=1; ii<ny-1; ii++)
      {
        for(jj=0; jj<ny; jj++)
        {
          /* compute integral for k(i,j) */  
          val = 0.0;
          for(px=1; px<=npx; px++)
          {
            for(py=1; py<=npy; py++)
            {
              val = val + wx[px]*wy[py]*galk(xx[px],yy[py],i,j,ii,jj);
            }
          }
          kg[(i*ny+ii)][(j*ny+jj)] = val;
        } /* end jj loop */
      } /* end ii loop */
    } /* end j loop */
  /* }  end i loop  thread */

  printf("compute %d: mid  \n", i);
  fflush(stdout); /* a good idea inside a thread */

  /* compute integral for f(i) */  
  /*for(i=1; i<nx-1; i++) thread */
  /*{                     thread */
    for(ii=1; ii<ny-1; ii++)
    {
      val = 0.0;
      for(px=1; px<=npx; px++)
      {
        for(py=1; py<=npy; py++)
        {
          val = val + wx[px]*wy[py]*galf(xx[px],yy[py],i,ii);
        }
      }
      fg[i*ny+ii] = val;
    }
  /* } thread */
  printf("compute %d: ending  \n", i);
  fflush(stdout); /* a good idea inside a thread */
  return arg;
} /* end compute */

/* end fempt_bihar2dps_la.c */