#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdbool.h>
#include <string.h>
#include <math.h>
#include <limits.h>
#include <time.h>

#include "cleanbench.h"
#include "randnum.h"


/***********************
**  LU DECOMPOSITION  **
** (Linear Equations) **
************************
** These routines come from "Numerical Recipes in Pascal".
** Note that, as in the assignment algorithm, though we
** separately define LUARRAYROWS and LUARRAYCOLS, the two
** must be the same value (this routine depends on a square
** matrix).
*/

/*
** ARRAY_MAX
**
** This sets the upper limit on the number of arrays
** that the benchmark will attempt to build before
** flagging an error.  It is not a critical constant, and
** may be increased if your system has the horsepower.
*/
#define ARRAY_MAX 10000

#define LUARRAYROWS 101L
#define LUARRAYCOLS 101L

typedef struct
{       union
	{       double *p;
		double (*ap)[][LUARRAYCOLS];
	} ptrs;
} LUdblptr;

double *LUtempvv; /* FIXME: really? */

static clock_t DoLUIteration(double *a, double *b,
	double *abase, double *bbase,
	unsigned long num_arrays);
static void build_problem( double a[][LUARRAYCOLS],
	double b[LUARRAYROWS]);
static int ludcmp(double a[][LUARRAYCOLS],
	int n, int indx[], int *d);
static void lubksb(double a[][LUARRAYCOLS],
	int n, int indx[LUARRAYROWS],
	double b[LUARRAYROWS]);
static int lusolve(double a[][LUARRAYCOLS],
	int n, double b[LUARRAYROWS]);

/*********
** DoLU **
**********
** Perform the LU decomposition benchmark.
*/
double
DoLU(void)
{
        clock_t         total_time = 0;
        int             iterations = 0;
        double*         a = NULL;
        double*         b = NULL;
        double*         abase = NULL;
        double*         bbase = NULL;
        LUdblptr        ptra;
        static int      num_arrays = 0;
        static bool     is_adjusted = false;

        /*
        ** Our first step is to build a "solvable" problem.  This
        ** will become the "seed" set that all others will be
        ** derived from. (I.E., we'll simply copy these arrays
        ** into the others.
        */
        a = malloc(sizeof(double) * LUARRAYCOLS * LUARRAYROWS);
        b = malloc(sizeof(double) * LUARRAYROWS);

        /*
        ** We need to allocate a temp vector that is used by the LU
        ** algorithm.  This removes the allocation routine from the
        ** timing.
        */
        LUtempvv = malloc(sizeof(double)*LUARRAYROWS);

        /*
        ** Build a problem to be solved.
        */
        ptra.ptrs.p=a;                  /* Gotta coerce linear array to 2D array */
        build_problem(*ptra.ptrs.ap, b);

        /*
        ** Now that we have a problem built, see if we need to do
        ** auto-adjust.  If so, repeatedly call the DoLUIteration routine,
        ** increasing the number of solutions per iteration as you go.
        */
        if (is_adjusted == false) {
                is_adjusted = true;

                do {
                        ++num_arrays;

                        abase = realloc(abase, sizeof(double) * LUARRAYCOLS * LUARRAYROWS * (num_arrays + 1));

                        bbase = realloc(bbase, sizeof(double) * LUARRAYROWS * (num_arrays + 1));
                } while ((DoLUIteration(a, b, abase, bbase, num_arrays) <= MINIMUM_TICKS) && (num_arrays <= ARRAY_MAX));

                /*
                ** Were we able to do it?
                */
                if (num_arrays==0) {
                        fprintf(stderr, "FPU:LU -- Array limit reached\n");
                        free(a);
                        free(b);
                        free(abase);
                        free(bbase);
                        free(LUtempvv);
                        exit(1);
                }
        } else {
                /*
                ** Don't need to adjust -- just allocate the proper
                ** number of arrays and proceed.
                */
                abase = malloc(sizeof(double) * LUARRAYCOLS * LUARRAYROWS * num_arrays);

                bbase = malloc(sizeof(double) * LUARRAYROWS * num_arrays);
        }

        do {
                total_time += DoLUIteration(a, b, abase, bbase, num_arrays);
                iterations += num_arrays;
        } while(total_time < MINIMUM_SECONDS * CLOCKS_PER_SEC);

        free(a);
        free(b);
        free(abase);
        free(bbase);
        free(LUtempvv);

        return (double)(iterations * CLOCKS_PER_SEC) / (double)total_time;
}

/******************
** DoLUIteration **
*******************
** Perform an iteration of the LU decomposition benchmark.
** An iteration refers to the repeated solution of several
** identical matrices.
*/
static clock_t
DoLUIteration(double *a,double *b, double *abase, double *bbase, unsigned long num_arrays)
{
        clock_t start, stop;
double *locabase;
double *locbbase;
LUdblptr ptra;  /* For converting ptr to 2D array */
unsigned long j,i;              /* Indexes */


/*
** Move the seed arrays (a & b) into the destination
** arrays;
*/
for(j=0;j<num_arrays;j++)
{       locabase=abase+j*LUARRAYROWS*LUARRAYCOLS;
	locbbase=bbase+j*LUARRAYROWS;
	for(i=0;i<LUARRAYROWS*LUARRAYCOLS;i++)
		*(locabase+i)=*(a+i);
	for(i=0;i<LUARRAYROWS;i++)
		*(locbbase+i)=*(b+i);
}

        start = clock();

for(i=0;i<num_arrays;i++)
{       locabase=abase+i*LUARRAYROWS*LUARRAYCOLS;
	locbbase=bbase+i*LUARRAYROWS;
	ptra.ptrs.p=locabase;
	lusolve(*ptra.ptrs.ap,LUARRAYROWS,locbbase);
}

        stop = clock();

        return stop - start;
}

/******************
** build_problem **
*******************
** Constructs a solvable set of linear equations.  It does this by
** creating an identity matrix, then loading the solution vector
** with random numbers.  After that, the identity matrix and
** solution vector are randomly "scrambled".  Scrambling is
** done by (a) randomly selecting a row and multiplying that
** row by a random number and (b) adding one randomly-selected
** row to another.
*/
static void build_problem(double a[][LUARRAYCOLS],
		double b[LUARRAYROWS])
{
long i,j,k,k1;  /* Indexes */
double rcon;     /* Random constant */

/*
** Reset random number generator
*/
/* randnum(13L); */
randnum((int32_t)13);

/*
** Build an identity matrix.
** We'll also use this as a chance to load the solution
** vector.
*/
for(i = 0; i < LUARRAYROWS; i++)
{       /* b[i]=(double)(abs_randwc(100L)+1L); */
	b[i]=(double)(abs_randwc((int32_t)100)+(int32_t)1);
	for(j = 0; j < LUARRAYROWS; j++)
		if(i==j)
		        /* a[i][j]=(double)(abs_randwc(1000L)+1L); */
			a[i][j]=(double)(abs_randwc((int32_t)1000)+(int32_t)1);
		else
			a[i][j]=(double)0.0;
}

/*
** Scramble.  Do this 8n times.  See comment above for
** a description of the scrambling process.
*/

for(i = 0; i < 8 * LUARRAYROWS; i++)
{
	/*
	** Pick a row and a random constant.  Multiply
	** all elements in the row by the constant.
	*/
 /*       k=abs_randwc((long)n);
	rcon=(double)(abs_randwc(20L)+1L);
	for(j=0;j<n;j++)
		a[k][j]=a[k][j]*rcon;
	b[k]=b[k]*rcon;
*/
	/*
	** Pick two random rows and add second to
	** first.  Note that we also occasionally multiply
	** by minus 1 so that we get a subtraction operation.
	*/
        /* k=abs_randwc((long)n); */
        /* k1=abs_randwc((long)n); */
	k=abs_randwc((int32_t)LUARRAYROWS);
	k1=abs_randwc((int32_t)LUARRAYROWS);
	if(k!=k1)
	{
		if(k<k1) rcon=(double)1.0;
			else rcon=(double)-1.0;
		for(j = 0; j < LUARRAYROWS; j++)
			a[k][j]+=a[k1][j]*rcon;;
		b[k]+=b[k1]*rcon;
	}
}
}


/***********
** ludcmp **
************
** From the procedure of the same name in "Numerical Recipes in Pascal",
** by Press, Flannery, Tukolsky, and Vetterling.
** Given an nxn matrix a[], this routine replaces it by the LU
** decomposition of a rowwise permutation of itself.  a[] and n
** are input.  a[] is output, modified as follows:
**   --                       --
**  |  b(1,1) b(1,2) b(1,3)...  |
**  |  a(2,1) b(2,2) b(2,3)...  |
**  |  a(3,1) a(3,2) b(3,3)...  |
**  |  a(4,1) a(4,2) a(4,3)...  |
**  |  ...                      |
**   --                        --
**
** Where the b(i,j) elements form the upper triangular matrix of the
** LU decomposition, and the a(i,j) elements form the lower triangular
** elements.  The LU decomposition is calculated so that we don't
** need to store the a(i,i) elements (which would have laid along the
** diagonal and would have all been 1).
**
** indx[] is an output vector that records the row permutation
** effected by the partial pivoting; d is output as +/-1 depending
** on whether the number of row interchanges was even or odd,
** respectively.
** Returns 0 if matrix singular, else returns 1.
*/
static int ludcmp(double a[][LUARRAYCOLS],
		int n,
		int indx[],
		int *d)
{

double big;     /* Holds largest element value */
double sum;
double dum;     /* Holds dummy value */
int i,j,k;      /* Indexes */
int imax=0;     /* Holds max index value */
double tiny;    /* A really small number */

tiny=(double)1.0e-20;

*d=1;           /* No interchanges yet */

for(i=0;i<n;i++)
{       big=(double)0.0;
	for(j=0;j<n;j++)
		if((double)fabs(a[i][j]) > big)
			big=fabs(a[i][j]);
	/* Bail out on singular matrix */
	if(big == 0.0) return 0;
	LUtempvv[i]=1.0/big;
}

/*
** Crout's algorithm...loop over columns.
*/
for(j=0;j<n;j++)
{       if(j!=0)
		for(i=0;i<j;i++)
		{       sum=a[i][j];
			if(i!=0)
				for(k=0;k<i;k++)
					sum-=(a[i][k]*a[k][j]);
			a[i][j]=sum;
		}
	big = 0.0;
	for(i=j;i<n;i++)
	{       sum=a[i][j];
		if(j!=0)
			for(k=0;k<j;k++)
				sum-=a[i][k]*a[k][j];
		a[i][j]=sum;
		dum=LUtempvv[i]*fabs(sum);
		if(dum>=big)
		{       big=dum;
			imax=i;
		}
	}
	if(j!=imax)             /* Interchange rows if necessary */
	{       for(k=0;k<n;k++)
		{       dum=a[imax][k];
			a[imax][k]=a[j][k];
			a[j][k]=dum;
		}
		*d=-*d;         /* Change parity of d */
		dum=LUtempvv[imax];
		LUtempvv[imax]=LUtempvv[j]; /* Don't forget scale factor */
		LUtempvv[j]=dum;
	}
	indx[j]=imax;
	/*
	** If the pivot element is zero, the matrix is singular
	** (at least as far as the precision of the machine
	** is concerned.)  We'll take the original author's
	** recommendation and replace 0.0 with "tiny".
	*/
	if(a[j][j] == 0.0)
		a[j][j]=tiny;

	if(j!=(n-1))
	{       dum=1.0/a[j][j];
		for(i=j+1;i<n;i++)
			a[i][j]=a[i][j]*dum;
	}
}

return(1);
}

/***********
** lubksb **
************
** Also from "Numerical Recipes in Pascal".
** This routine solves the set of n linear equations A X = B.
** Here, a[][] is input, not as the matrix A, but as its
** LU decomposition, created by the routine ludcmp().
** Indx[] is input as the permutation vector returned by ludcmp().
**  b[] is input as the right-hand side an returns the
** solution vector X.
** a[], n, and indx are not modified by this routine and
** can be left in place for different values of b[].
** This routine takes into account the possibility that b will
** begin with many zero elements, so it is efficient for use in
** matrix inversion.
*/
static void lubksb( double a[][LUARRAYCOLS],
		int n,
		int indx[LUARRAYROWS],
		double b[LUARRAYROWS])
{

int i,j;        /* Indexes */
int ip;         /* "pointer" into indx */
int ii;
double sum;

/*
** When ii is set to a positive value, it will become
** the index of the first nonvanishing element of b[].
** We now do the forward substitution. The only wrinkle
** is to unscramble the permutation as we go.
*/
ii=-1;
for(i=0;i<n;i++)
{       ip=indx[i];
	sum=b[ip];
	b[ip]=b[i];
	if(ii!=-1)
		for(j=ii;j<i;j++)
			sum=sum-a[i][j]*b[j];
	else
		/*
		** If a nonzero element is encountered, we have
		** to do the sums in the loop above.
		*/
		if(sum != 0.0)
			ii=i;
	b[i]=sum;
}
/*
** Do backsubstitution
*/
for(i=(n-1);i>=0;i--)
{
	sum=b[i];
	if(i!=(n-1))
		for(j=(i+1);j<n;j++)
			sum=sum-a[i][j]*b[j];
	b[i]=sum/a[i][i];
}
return;
}

/************
** lusolve **
*************
** Solve a linear set of equations: A x = b
** Original matrix A will be destroyed by this operation.
** Returns 0 if matrix is singular, 1 otherwise.
*/
static int lusolve(double a[][LUARRAYCOLS],
		int n,
		double b[LUARRAYROWS])
{
int indx[LUARRAYROWS];
int d;
#ifdef DEBUG
int i,j;
#endif

if(ludcmp(a,n,indx,&d)==0) return(0);

/* Matrix not singular -- proceed */
lubksb(a,n,indx,b);

#ifdef DEBUG
printf("Solution:\n");
for(i=0;i<n;i++)
{
  for(j=0;j<n;j++){
  /*
    printf("%6.2f ",a[i][j]);
  */
  }
  printf("%6.2f\t",b[i]);
  /*
    printf("\n");
  */
}
printf("\n");
#endif

return(1);
}