← Back to branch summary

~noskcaj/ubuntu/trusty/hardinfo/ftbfs

« back to all changes in this revision

Viewing changes to fftbench.c

  • Committer: Bazaar Package Importer
  • Author(s): Agney Lopes Roth Ferraz
  • Date: 2009-03-28 22:55:02 UTC
  • mfrom: (3.1.2 squeeze)
  • Revision ID: james.westby@ubuntu.com-20090328225502-p4bnvi8q6hr95cij
Tags: 0.5c-1
http://bugs.debian.org/517591
http://bugs.debian.org/511237
http://bugs.debian.org/457703
http://bugs.debian.org/519256
http://bugs.debian.org/449250
http://bugs.debian.org/457820
http://bugs.debian.org/497758
New upstream version. 
(Closes: #517591, #511237, #457703, #519256, #449250, #457820, #497758) 

Show diffs side-by-side

added added

removed removed

Lines of Context:
fftbench.c
 
1
/*
 
2
    fftbench.c
 
3
 
 
4
    Written by Scott Robert Ladd (scott@coyotegulch.com)
 
5
    No rights reserved. This is public domain software, for use by anyone.
 
6
 
 
7
    A number-crunching benchmark using LUP-decomposition to solve a large
 
8
    linear equation.
 
9
 
 
10
    The code herein is design for the purpose of testing computational
 
11
    performance; error handling is minimal.
 
12
    
 
13
    In fact, this is a weak implementation of the FFT; unfortunately, all
 
14
    of my really nifty FFTs are in commercial code, and I haven't had time
 
15
    to write a new FFT routine for this benchmark. I may add a Hartley
 
16
    transform to the seat, too.
 
17
 
 
18
    Actual benchmark results can be found at:
 
19
            http://www.coyotegulch.com
 
20
 
 
21
    Please do not use this information or algorithm in any way that might
 
22
    upset the balance of the universe or otherwise cause a disturbance in
 
23
    the space-time continuum.
 
24
*/
 
25
 
 
26
#include <time.h>
 
27
#include <string.h>
 
28
#include <stdlib.h>
 
29
#include <math.h>
 
30
#include <stdbool.h>
 
31
#include <stdio.h>
 
32
 
 
33
// embedded random number generator; ala Park and Miller
 
34
static long seed = 1325;
 
35
static const long IA = 16807;
 
36
static const long IM = 2147483647;
 
37
static const double AM = 4.65661287525E-10;
 
38
static const long IQ = 127773;
 
39
static const long IR = 2836;
 
40
static const long MASK = 123459876;
 
41
 
 
42
static double random_double()
 
43
{
 
44
    long k;
 
45
    double result;
 
46
 
 
47
    seed ^= MASK;
 
48
    k = seed / IQ;
 
49
    seed = IA * (seed - k * IQ) - IR * k;
 
50
 
 
51
    if (seed < 0)
 
52
        seed += IM;
 
53
 
 
54
    result = AM * seed;
 
55
    seed ^= MASK;
 
56
 
 
57
    return result;
 
58
}
 
59
 
 
60
static const int N = 800;
 
61
static const int NM1 = 799;     // N - 1
 
62
static const int NP1 = 801;     // N + 1
 
63
 
 
64
static int *lup_decompose(double **a)
 
65
{
 
66
    int i, j, k, k2, t;
 
67
    double p, temp;
 
68
 
 
69
    int *perm = (int *) malloc(sizeof(double) * N);
 
70
 
 
71
    for (i = 0; i < N; ++i)
 
72
        perm[i] = i;
 
73
 
 
74
    for (k = 0; k < NM1; ++k) {
 
75
        p = 0.0;
 
76
 
 
77
        for (i = k; i < N; ++i) {
 
78
            temp = fabs(a[i][k]);
 
79
 
 
80
            if (temp > p) {
 
81
                p = temp;
 
82
                k2 = i;
 
83
            }
 
84
        }
 
85
 
 
86
        // check for invalid a
 
87
        if (p == 0.0)
 
88
            return NULL;
 
89
 
 
90
        // exchange rows
 
91
        t = perm[k];
 
92
        perm[k] = perm[k2];
 
93
        perm[k2] = t;
 
94
 
 
95
        for (i = 0; i < N; ++i) {
 
96
            temp = a[k][i];
 
97
            a[k][i] = a[k2][i];
 
98
            a[k2][i] = temp;
 
99
        }
 
100
 
 
101
        for (i = k + 1; i < N; ++i) {
 
102
            a[i][k] /= a[k][k];
 
103
 
 
104
            for (j = k + 1; j < N; ++j)
 
105
                a[i][j] -= a[i][k] * a[k][j];
 
106
        }
 
107
    }
 
108
 
 
109
    return perm;
 
110
}
 
111
 
 
112
static double *lup_solve(double **a, int *perm, double *b)
 
113
{
 
114
    int i, j, j2;
 
115
    double sum, u;
 
116
 
 
117
    double *y = (double *) malloc(sizeof(double) * N);
 
118
    double *x = (double *) malloc(sizeof(double) * N);
 
119
 
 
120
    for (i = 0; i < N; ++i) {
 
121
        y[i] = 0.0;
 
122
        x[i] = 0.0;
 
123
    }
 
124
 
 
125
    for (i = 0; i < N; ++i) {
 
126
        sum = 0.0;
 
127
        j2 = 0;
 
128
 
 
129
        for (j = 1; j <= i; ++j) {
 
130
            sum += a[i][j2] * y[j2];
 
131
            ++j2;
 
132
        }
 
133
 
 
134
        y[i] = b[perm[i]] - sum;
 
135
    }
 
136
 
 
137
    i = NM1;
 
138
 
 
139
    while (1) {
 
140
        sum = 0.0;
 
141
        u = a[i][i];
 
142
 
 
143
        for (j = i + 1; j < N; ++j)
 
144
            sum += a[i][j] * x[j];
 
145
 
 
146
        x[i] = (y[i] - sum) / u;
 
147
 
 
148
        if (i == 0)
 
149
            break;
 
150
 
 
151
        --i;
 
152
    }
 
153
 
 
154
    free(y);
 
155
 
 
156
    return x;
 
157
}
 
158
 
 
159
static double **a, *b, *r;
 
160
static int *p;
 
161
 
 
162
void fft_bench_init(void)
 
163
{
 
164
    int i, j;
 
165
 
 
166
    // generate test data            
 
167
    a = (double **) malloc(sizeof(double *) * N);
 
168
 
 
169
    for (i = 0; i < N; ++i) {
 
170
        a[i] = (double *) malloc(sizeof(double) * N);
 
171
 
 
172
        for (j = 0; j < N; ++j)
 
173
            a[i][j] = random_double();
 
174
    }
 
175
 
 
176
    b = (double *) malloc(sizeof(double) * N);
 
177
 
 
178
    for (i = 0; i < N; ++i)
 
179
        b[i] = random_double();
 
180
 
 
181
}
 
182
 
 
183
void fft_bench_start(void)
 
184
{
 
185
    p = lup_decompose(a);
 
186
    r = lup_solve(a, p, b);
 
187
}
 
188
 
 
189
void fft_bench_finish(void)
 
190
{
 
191
    int i;
 
192
    
 
193
    // clean up
 
194
    for (i = 0; i < N; ++i)
 
195
        free(a[i]);
 
196
 
 
197
    free(a);
 
198
    free(b);
 
199
    free(p);
 
200
    free(r);
 
201
}