|
|
// 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<EFBFBD>, MIPS R4400, |
|
|
175 Mhz, IRIX 5.2, "-O". |
|
|
|
|
|
32.4 Pentium 133, Windows NT, Microsoft Visual |
|
|
C++ 4.0. |
|
|
|
|
|
17.25 Silicon Graphics Indigo<EFBFBD>, 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; |
|
|
}
|
|
|
|
