/
ffbench.c
372 lines (296 loc) · 10.2 KB
/
ffbench.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
// This small program computes a Fast Fourier Transform. It tests
// Valgrind's handling of FP operations. It is representative of all
// programs that do a lot of FP operations.
// Licensing: This program is closely based on the one of the same name from
// http://www.fourmilab.ch/. The front page of that site says:
//
// "Except for a few clearly-marked exceptions, all the material on this
// site is in the public domain and may be used in any manner without
// permission, restriction, attribution, or compensation."
/*
Two-dimensional FFT benchmark
Designed and implemented by John Walker in April of 1989.
This benchmark executes a specified number of passes (default
20) through a loop in which each iteration performs a fast
Fourier transform of a square matrix (default size 256x256) of
complex numbers (default precision double), followed by the
inverse transform. After all loop iterations are performed
the results are checked against known correct values.
This benchmark is intended for use on C implementations which
define "int" as 32 bits or longer and permit allocation and
direct addressing of arrays larger than one megabyte.
If CAPOUT is defined, the result after all iterations is
written as a CA Lab pattern file. This is intended for
debugging in case horribly wrong results are obtained on a
given machine.
Archival timings are run with the definitions below set as
follows: Float = double, Asize = 256, Passes = 20, CAPOUT not
defined.
Time (seconds) System
2393.93 Sun 3/260, SunOS 3.4, C, "-f68881 -O".
(John Walker).
1928 Macintosh IIx, MPW C 3.0, "-mc68020
-mc68881 -elems881 -m". (Hugh Hoover).
1636.1 Sun 4/110, "cc -O3 -lm". (Michael McClary).
The suspicion is that this is software
floating point.
1556.7 Macintosh II, A/UX, "cc -O -lm"
(Michael McClary).
1388.8 Sun 386i/250, SunOS 4.0.1 C
"-O /usr/lib/trig.il". (James Carrington).
1331.93 Sun 3/60, SunOS 4.0.1, C,
"-O4 -f68881 /usr/lib/libm.il"
(Bob Elman).
1204.0 Apollo Domain DN4000, C, "-cpu 3000 -opt 4".
(Sam Crupi).
1174.66 Compaq 386/25, SCO Xenix 386 C.
(Peter Shieh).
1068 Compaq 386/25, SCO Xenix 386,
Metaware High C. (Robert Wenig).
1064.0 Sun 3/80, SunOS 4.0.3 Beta C
"-O3 -f68881 /usr/lib/libm.il". (James Carrington).
1061.4 Compaq 386/25, SCO Xenix, High C 1.4.
(James Carrington).
1059.79 Compaq 386/25, 387/25, High C 1.4,
DOS|Extender 2.2, 387 inline code
generation. (Nathan Bender).
777.14 Compaq 386/25, IIT 3C87-25 (387 Compatible),
High C 1.5, DOS|Extender 2.2, 387 inline
code generation. (Nathan Bender).
751 Compaq DeskPro 386/33, High C 1.5 + DOS|Extender,
387 code generation. (James Carrington).
431.44 Compaq 386/25, Weitek 3167-25, DOS 3.31,
High C 1.4, DOS|Extender, Weitek code generation.
(Nathan Bender).
344.9 Compaq 486/25, Metaware High C 1.6, Phar Lap
DOS|Extender, in-line floating point. (Nathan
Bender).
324.2 Data General Motorola 88000, 16 Mhz, Gnu C.
323.1 Sun 4/280, C, "-O4". (Eric Hill).
254 Compaq SystemPro 486/33, High C 1.5 + DOS|Extender,
387 code generation. (James Carrington).
242.8 Silicon Graphics Personal IRIS, MIPS R2000A,
12.5 Mhz, "-O3" (highest level optimisation).
(Mike Zentner).
233.0 Sun SPARCStation 1, C, "-O4", SunOS 4.0.3.
(Nathan Bender).
187.30 DEC PMAX 3100, MIPS 2000 chip.
(Robert Wenig).
120.46 Sun SparcStation 2, C, "-O4", SunOS 4.1.1.
(John Walker).
120.21 DEC 3MAX, MIPS 3000, "-O4".
98.0 Intel i860 experimental environment,
OS/2, data caching disabled. (Kern
Sibbald).
34.9 Silicon Graphics Indigo², MIPS R4400,
175 Mhz, IRIX 5.2, "-O".
32.4 Pentium 133, Windows NT, Microsoft Visual
C++ 4.0.
17.25 Silicon Graphics Indigo², MIPS R4400,
175 Mhz, IRIX 6.5, "-O3".
14.10 Dell Dimension XPS R100, Pentium II 400 MHz,
Windows 98, Microsoft Visual C 5.0.
10.7 Hewlett-Packard Kayak XU 450Mhz Pentium II,
Microsoft Visual C++ 6.0, Windows NT 4.0sp3. (Nathan Bender).
5.09 Sun Ultra 2, UltraSPARC V9, 300 MHz, gcc -O3.
0.846 Dell Inspiron 9100, Pentium 4, 3.4 GHz, gcc -O3.
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
/* The program may be run with Float defined as either float or
double. With IEEE arithmetic, the same answers are generated for
either floating point mode. */
#define Float double /* Floating point type used in FFT */
#define Asize 256 /* Array edge size */
#define Passes 20 /* Number of FFT/Inverse passes */
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<=(b)?(a):(b))
/*
Multi-dimensional fast Fourier transform
Adapted from Press et al., "Numerical Recipes in C".
*/
#define SWAP(a,b) tempr=(a); (a)=(b); (b)=tempr
static void fourn(data, nn, ndim, isign)
Float data[];
int nn[], ndim, isign;
{
register int i1, i2, i3;
int i2rev, i3rev, ip1, ip2, ip3, ifp1, ifp2;
int ibit, idim, k1, k2, n, nprev, nrem, ntot;
Float tempi, tempr;
double theta, wi, wpi, wpr, wr, wtemp;
ntot = 1;
for (idim = 1; idim <= ndim; idim++)
ntot *= nn[idim];
nprev = 1;
for (idim = ndim; idim >= 1; idim--) {
n = nn[idim];
nrem = ntot / (n * nprev);
ip1 = nprev << 1;
ip2 = ip1 * n;
ip3 = ip2 * nrem;
i2rev = 1;
for (i2 = 1; i2 <= ip2; i2 += ip1) {
if (i2 < i2rev) {
for (i1 = i2; i1 <= i2 + ip1 - 2; i1 += 2) {
for (i3 = i1; i3 <= ip3; i3 += ip2) {
i3rev = i2rev + i3 - i2;
SWAP(data[i3], data[i3rev]);
SWAP(data[i3 + 1], data[i3rev + 1]);
}
}
}
ibit = ip2 >> 1;
while (ibit >= ip1 && i2rev > ibit) {
i2rev -= ibit;
ibit >>= 1;
}
i2rev += ibit;
}
ifp1 = ip1;
while (ifp1 < ip2) {
ifp2 = ifp1 << 1;
theta = isign * 6.28318530717959 / (ifp2 / ip1);
wtemp = sin(0.5 * theta);
wpr = -2.0 * wtemp * wtemp;
wpi = sin(theta);
wr = 1.0;
wi = 0.0;
for (i3 = 1; i3 <= ifp1; i3 += ip1) {
for (i1 = i3; i1 <= i3 + ip1 - 2; i1 += 2) {
for (i2 = i1; i2 <= ip3; i2 += ifp2) {
k1 = i2;
k2 = k1 + ifp1;
tempr = wr * data[k2] - wi * data[k2 + 1];
tempi = wr * data[k2 + 1] + wi * data[k2];
data[k2] = data[k1] - tempr;
data[k2 + 1] = data[k1 + 1] - tempi;
data[k1] += tempr;
data[k1 + 1] += tempi;
}
}
wr = (wtemp = wr) * wpr - wi * wpi + wr;
wi = wi * wpr + wtemp * wpi + wi;
}
ifp1 = ifp2;
}
nprev *= n;
}
}
#undef SWAP
int main()
{
int i, j, k, l, m, npasses = Passes, faedge;
Float *fdata /* , *fd */ ;
static int nsize[] = {0, 0, 0};
long fanum, fasize;
double mapbase, mapscale, /* x, */ rmin, rmax, imin, imax;
faedge = Asize; /* FFT array edge size */
fanum = faedge * faedge; /* Elements in FFT array */
fasize = ((fanum + 1) * 2 * sizeof(Float)); /* FFT array size */
nsize[1] = nsize[2] = faedge;
fdata = (Float *) malloc(fasize);
if (fdata == NULL) {
fprintf(stdout, "Can't allocate data array.\n");
exit(1);
}
/* Generate data array to process. */
#define Re(x,y) fdata[1 + (faedge * (x) + (y)) * 2]
#define Im(x,y) fdata[2 + (faedge * (x) + (y)) * 2]
memset(fdata, 0, fasize);
for (i = 0; i < faedge; i++) {
for (j = 0; j < faedge; j++) {
if (((i & 15) == 8) || ((j & 15) == 8))
Re(i, j) = 128.0;
}
}
for (i = 0; i < npasses; i++) {
/*printf("Pass %d\n", i);*/
/* Transform image to frequency domain. */
fourn(fdata, nsize, 2, 1);
/* Back-transform to image. */
fourn(fdata, nsize, 2, -1);
}
{
double r, ij, ar, ai;
rmin = 1e10; rmax = -1e10;
imin = 1e10; imax = -1e10;
ar = 0;
ai = 0;
for (i = 1; i <= fanum; i += 2) {
r = fdata[i];
ij = fdata[i + 1];
ar += r;
ai += ij;
rmin = min(r, rmin);
rmax = max(r, rmax);
imin = min(ij, imin);
imax = max(ij, imax);
}
#ifdef DEBUG
printf("Real min %.4g, max %.4g. Imaginary min %.4g, max %.4g.\n",
rmin, rmax, imin, imax);
printf("Average real %.4g, imaginary %.4g.\n",
ar / fanum, ai / fanum);
#endif
mapbase = rmin;
mapscale = 255 / (rmax - rmin);
}
/* See if we got the right answers. */
m = 0;
for (i = 0; i < faedge; i++) {
for (j = 0; j < faedge; j++) {
k = (Re(i, j) - mapbase) * mapscale;
l = (((i & 15) == 8) || ((j & 15) == 8)) ? 255 : 0;
if (k != l) {
m++;
fprintf(stdout,
"Wrong answer at (%d,%d)! Expected %d, got %d.\n",
i, j, l, k);
}
}
}
if (m == 0) {
fprintf(stdout, "%d passes. No errors in results.\n", npasses);
} else {
fprintf(stdout, "%d passes. %d errors in results.\n",
npasses, m);
}
#ifdef CAPOUT
/* Output the result of the transform as a CA Lab pattern
file for debugging. */
{
#define SCRX 322
#define SCRY 200
#define SCRN (SCRX * SCRY)
unsigned char patarr[SCRY][SCRX];
FILE *fp;
/* Map user external state numbers to internal state index */
#define UtoI(x) (((((x) >> 1) & 0x7F) | ((x) << 7)) & 0xFF)
/* Copy data from FFT buffer to map. */
memset(patarr, 0, sizeof patarr);
l = (SCRX - faedge) / 2;
m = (faedge > SCRY) ? 0 : ((SCRY - faedge) / 2);
for (i = 1; i < faedge; i++) {
for (j = 0; j < min(SCRY, faedge); j++) {
k = (Re(i, j) - mapbase) * mapscale;
patarr[j + m][i + l] = UtoI(k);
}
}
/* Dump pattern map to file. */
fp = fopen("fft.cap", "w");
if (fp == NULL) {
fprintf(stdout, "Cannot open output file.\n");
exit(0);
}
putc(':', fp);
putc(1, fp);
fwrite(patarr, SCRN, 1, fp);
putc(6, fp);
fclose(fp);
}
#endif
return 0;
}