diff options
Diffstat (limited to 'fftbench.c')
| -rw-r--r-- | fftbench.c | 201 | 
1 files changed, 0 insertions, 201 deletions
| diff --git a/fftbench.c b/fftbench.c deleted file mode 100644 index 597c5693..00000000 --- a/fftbench.c +++ /dev/null @@ -1,201 +0,0 @@ -/* -    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); -} | 
