aboutsummaryrefslogtreecommitdiff
path: root/fftbench.c
diff options
context:
space:
mode:
Diffstat (limited to 'fftbench.c')
-rw-r--r--fftbench.c201
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);
-}