/* time_mpf_tsimeqb.c  multiple precision threaded version */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include "udrnrt.h"
#include "mpf_tsimeq.h"

#define N 4097
static mpf_t A[N*N];
static mpf_t B[N*N];
static mpf_t X[N];
static int np=12;
int digits = 50; 

#undef abs
#define abs(x) ((x)<0.0?(-(x)):(x))

static void init_matrix(int n)
{
  int i, j, m;
  m = n+1;
  for(i=0; i<n; i++)
  {
    mpf_set_d(A[i*m+n], udrnrt());
    mpf_set(B[i*m+n], A[i*m+n]);
    for(j=0; j<n; j++)
    {
      mpf_set_d(A[i*m+j], udrnrt());
      mpf_set(B[i*m+j], A[i*m+j]);
    }
  }
  if(n<9)
  {
    for(i=0; i<n; i++)
    {
      for(j=0; j<n+1; j++)
	gmp_printf("%d,%d, %5.2Ff \n",i,j, A[i*m+j]);

    }
  }
}

int main(int argc, char * argv[])
{
  int i, j, k;
  int n, m;
  mpf_t t, total, p;
  double wall_start, wall_now, wall_next, cpu_start, cpu_now, cpu_next;
  
  gmp_printf("time_mpf_tsimeq.c running, np=%d, digits=%d \n", np, digits);
  mpf_set_default_prec(digits*3.32);
  mpf_init(t);
  mpf_init(total);
  mpf_init(p);
  for(i=0; i<N; i++)
  {
    mpf_init(X[i]);
    for(j=0; j<N; j++)
    {
      mpf_init(A[i*N+j]);
      mpf_init(B[i*N+j]);
    }
  }


  cpu_start = (double)clock()/(double)CLOCKS_PER_SEC;
  wall_start = (double)time(NULL);
  for(n=8; n<=N; n=n*2)
  {
    init_matrix(n); /* B and copy in A */
    cpu_now = (double)clock()/(double)CLOCKS_PER_SEC;
    wall_now = (double)time(NULL);
    gmp_printf("matrix initialized at wall %g, cpu %g\n",
           wall_now-wall_start, cpu_now-cpu_start);
    mpf_tsimeqb(np, n, B, X); /* B destroyed */
    cpu_next = (double)clock()/(double)CLOCKS_PER_SEC;
    wall_next = (double)time(NULL);
    mpf_set_ui(total, 0);
    m = n+1;
    for(i=0; i<n; i++)
    {
      mpf_set_ui(t, 0);
      for(j=0; j<n; j++)
      {
        mpf_mul(p, A[i*m+j], X[j]);
        mpf_add(t, t, p);
      }
      mpf_sub(p, t, A[i*m+n]);
      mpf_abs(p, p);
      mpf_add(total, total, p);
      if(n<9) gmp_printf("X[%d]=%Fg, t=%Fg, Y[%d]=%Fg \n",
                      i, X[i], t, i, A[i*m+n]);
    }
    gmp_printf("n=%d, total error=%Fg, mpf_tsimeqb took wall %g, cpu %g seconds\n",
           n, total, wall_next-wall_now, cpu_next-cpu_now);
  }
  gmp_printf("average number of processors used = %f \n",
         (cpu_next-cpu_now)/(wall_next-wall_now));
  gmp_printf("time_mpf_tsimeqb.c exiting, np=%d, digits=%d \n", np, digits);
  return 0;
} /* end main time_mpf_tsimeqb.c */