/*
    fftbench.c

    Written by Scott Robert Ladd (scott@coyotegulch.com)
    No rights reserved. This is public domain software, for use by anyone.

    A number-crunching benchmark using LUP-decomposition to solve a large
    linear equation.

    The code herein is design for the purpose of testing computational
    performance; error handling is minimal.
    
    In fact, this is a weak implementation of the FFT; unfortunately, all
    of my really nifty FFTs are in commercial code, and I haven't had time
    to write a new FFT routine for this benchmark. I may add a Hartley
    transform to the seat, too.

    Actual benchmark results can be found at:
            http://www.coyotegulch.com

    Please do not use this information or algorithm in any way that might
    upset the balance of the universe or otherwise cause a disturbance in
    the space-time continuum.
*/

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

// embedded random number generator; ala Park and Miller
static long seed = 1325;
static const long IA = 16807;
static const long IM = 2147483647;
static const double AM = 4.65661287525E-10;
static const long IQ = 127773;
static const long IR = 2836;
static const long MASK = 123459876;

static double random_double()
{
    long k;
    double result;

    seed ^= MASK;
    k = seed / IQ;
    seed = IA * (seed - k * IQ) - IR * k;

    if (seed < 0)
	seed += IM;

    result = AM * seed;
    seed ^= MASK;

    return result;
}

static const int N = 800;
static const int NM1 = 799;	// N - 1
static const int NP1 = 801;	// N + 1

static int *lup_decompose(double **a)
{
    int i, j, k, k2, t;
    double p, temp;

    int *perm = (int *) malloc(sizeof(double) * N);

    for (i = 0; i < N; ++i)
	perm[i] = i;

    for (k = 0; k < NM1; ++k) {
	p = 0.0;

	for (i = k; i < N; ++i) {
	    temp = fabs(a[i][k]);

	    if (temp > p) {
		p = temp;
		k2 = i;
	    }
	}

	// check for invalid a
	if (p == 0.0)
	    return NULL;

	// exchange rows
	t = perm[k];
	perm[k] = perm[k2];
	perm[k2] = t;

	for (i = 0; i < N; ++i) {
	    temp = a[k][i];
	    a[k][i] = a[k2][i];
	    a[k2][i] = temp;
	}

	for (i = k + 1; i < N; ++i) {
	    a[i][k] /= a[k][k];

	    for (j = k + 1; j < N; ++j)
		a[i][j] -= a[i][k] * a[k][j];
	}
    }

    return perm;
}

static double *lup_solve(double **a, int *perm, double *b)
{
    int i, j, j2;
    double sum, u;

    double *y = (double *) malloc(sizeof(double) * N);
    double *x = (double *) malloc(sizeof(double) * N);

    for (i = 0; i < N; ++i) {
	y[i] = 0.0;
	x[i] = 0.0;
    }

    for (i = 0; i < N; ++i) {
	sum = 0.0;
	j2 = 0;

	for (j = 1; j <= i; ++j) {
	    sum += a[i][j2] * y[j2];
	    ++j2;
	}

	y[i] = b[perm[i]] - sum;
    }

    i = NM1;

    while (1) {
	sum = 0.0;
	u = a[i][i];

	for (j = i + 1; j < N; ++j)
	    sum += a[i][j] * x[j];

	x[i] = (y[i] - sum) / u;

	if (i == 0)
	    break;

	--i;
    }

    free(y);

    return x;
}

static double **a, *b, *r;
static int *p;

void fft_bench_init(void)
{
    int i, j;

    // generate test data            
    a = (double **) malloc(sizeof(double *) * N);

    for (i = 0; i < N; ++i) {
	a[i] = (double *) malloc(sizeof(double) * N);

	for (j = 0; j < N; ++j)
	    a[i][j] = random_double();
    }

    b = (double *) malloc(sizeof(double) * N);

    for (i = 0; i < N; ++i)
	b[i] = random_double();

}

void fft_bench_start(void)
{
    p = lup_decompose(a);
    r = lup_solve(a, p, b);
}

void fft_bench_finish(void)
{
    int i;
    
    // clean up
    for (i = 0; i < N; ++i)
	free(a[i]);

    free(a);
    free(b);
    free(p);
    free(r);
}