cannam@95: /* cannam@95: * Copyright (c) 2003, 2007-11 Matteo Frigo cannam@95: * Copyright (c) 2003, 2007-11 Massachusetts Institute of Technology cannam@95: * cannam@95: * This program is free software; you can redistribute it and/or modify cannam@95: * it under the terms of the GNU General Public License as published by cannam@95: * the Free Software Foundation; either version 2 of the License, or cannam@95: * (at your option) any later version. cannam@95: * cannam@95: * This program is distributed in the hope that it will be useful, cannam@95: * but WITHOUT ANY WARRANTY; without even the implied warranty of cannam@95: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the cannam@95: * GNU General Public License for more details. cannam@95: * cannam@95: * You should have received a copy of the GNU General Public License cannam@95: * along with this program; if not, write to the Free Software cannam@95: * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA cannam@95: * cannam@95: */ cannam@95: cannam@95: /* Lots of ugly duplication from verify-lib.c, plus lots of ugliness in cannam@95: general for all of the r2r variants...oh well, for now */ cannam@95: cannam@95: #include "verify.h" cannam@95: #include cannam@95: #include cannam@95: #include cannam@95: cannam@95: typedef struct { cannam@95: bench_problem *p; cannam@95: bench_tensor *probsz; cannam@95: bench_tensor *totalsz; cannam@95: bench_tensor *pckdsz; cannam@95: bench_tensor *pckdvecsz; cannam@95: } info; cannam@95: cannam@95: /* cannam@95: * Utility functions: cannam@95: */ cannam@95: cannam@95: static double dabs(double x) { return (x < 0.0) ? -x : x; } cannam@95: static double dmin(double x, double y) { return (x < y) ? x : y; } cannam@95: cannam@95: static double raerror(R *a, R *b, int n) cannam@95: { cannam@95: if (n > 0) { cannam@95: /* compute the relative Linf error */ cannam@95: double e = 0.0, mag = 0.0; cannam@95: int i; cannam@95: cannam@95: for (i = 0; i < n; ++i) { cannam@95: e = dmax(e, dabs(a[i] - b[i])); cannam@95: mag = dmax(mag, dmin(dabs(a[i]), dabs(b[i]))); cannam@95: } cannam@95: if (dabs(mag) < 1e-14 && dabs(e) < 1e-14) cannam@95: e = 0.0; cannam@95: else cannam@95: e /= mag; cannam@95: cannam@95: #ifdef HAVE_ISNAN cannam@95: BENCH_ASSERT(!isnan(e)); cannam@95: #endif cannam@95: return e; cannam@95: } else cannam@95: return 0.0; cannam@95: } cannam@95: cannam@95: #define by2pi(m, n) ((K2PI * (m)) / (n)) cannam@95: cannam@95: /* cannam@95: * Improve accuracy by reducing x to range [0..1/8] cannam@95: * before multiplication by 2 * PI. cannam@95: */ cannam@95: cannam@95: static trigreal bench_sincos(trigreal m, trigreal n, int sinp) cannam@95: { cannam@95: /* waiting for C to get tail recursion... */ cannam@95: trigreal half_n = n * 0.5; cannam@95: trigreal quarter_n = half_n * 0.5; cannam@95: trigreal eighth_n = quarter_n * 0.5; cannam@95: trigreal sgn = 1.0; cannam@95: cannam@95: if (sinp) goto sin; cannam@95: cos: cannam@95: if (m < 0) { m = -m; /* goto cos; */ } cannam@95: if (m > half_n) { m = n - m; goto cos; } cannam@95: if (m > eighth_n) { m = quarter_n - m; goto sin; } cannam@95: return sgn * COS(by2pi(m, n)); cannam@95: cannam@95: msin: cannam@95: sgn = -sgn; cannam@95: sin: cannam@95: if (m < 0) { m = -m; goto msin; } cannam@95: if (m > half_n) { m = n - m; goto msin; } cannam@95: if (m > eighth_n) { m = quarter_n - m; goto cos; } cannam@95: return sgn * SIN(by2pi(m, n)); cannam@95: } cannam@95: cannam@95: static trigreal cos2pi(int m, int n) cannam@95: { cannam@95: return bench_sincos((trigreal)m, (trigreal)n, 0); cannam@95: } cannam@95: cannam@95: static trigreal sin2pi(int m, int n) cannam@95: { cannam@95: return bench_sincos((trigreal)m, (trigreal)n, 1); cannam@95: } cannam@95: cannam@95: static trigreal cos00(int i, int j, int n) cannam@95: { cannam@95: return cos2pi(i * j, n); cannam@95: } cannam@95: cannam@95: static trigreal cos01(int i, int j, int n) cannam@95: { cannam@95: return cos00(i, 2*j + 1, 2*n); cannam@95: } cannam@95: cannam@95: static trigreal cos10(int i, int j, int n) cannam@95: { cannam@95: return cos00(2*i + 1, j, 2*n); cannam@95: } cannam@95: cannam@95: static trigreal cos11(int i, int j, int n) cannam@95: { cannam@95: return cos00(2*i + 1, 2*j + 1, 4*n); cannam@95: } cannam@95: cannam@95: static trigreal sin00(int i, int j, int n) cannam@95: { cannam@95: return sin2pi(i * j, n); cannam@95: } cannam@95: cannam@95: static trigreal sin01(int i, int j, int n) cannam@95: { cannam@95: return sin00(i, 2*j + 1, 2*n); cannam@95: } cannam@95: cannam@95: static trigreal sin10(int i, int j, int n) cannam@95: { cannam@95: return sin00(2*i + 1, j, 2*n); cannam@95: } cannam@95: cannam@95: static trigreal sin11(int i, int j, int n) cannam@95: { cannam@95: return sin00(2*i + 1, 2*j + 1, 4*n); cannam@95: } cannam@95: cannam@95: static trigreal realhalf(int i, int j, int n) cannam@95: { cannam@95: UNUSED(i); cannam@95: if (j <= n - j) cannam@95: return 1.0; cannam@95: else cannam@95: return 0.0; cannam@95: } cannam@95: cannam@95: static trigreal coshalf(int i, int j, int n) cannam@95: { cannam@95: if (j <= n - j) cannam@95: return cos00(i, j, n); cannam@95: else cannam@95: return cos00(i, n - j, n); cannam@95: } cannam@95: cannam@95: static trigreal unity(int i, int j, int n) cannam@95: { cannam@95: UNUSED(i); cannam@95: UNUSED(j); cannam@95: UNUSED(n); cannam@95: return 1.0; cannam@95: } cannam@95: cannam@95: typedef trigreal (*trigfun)(int, int, int); cannam@95: cannam@95: static void rarand(R *a, int n) cannam@95: { cannam@95: int i; cannam@95: cannam@95: /* generate random inputs */ cannam@95: for (i = 0; i < n; ++i) { cannam@95: a[i] = mydrand(); cannam@95: } cannam@95: } cannam@95: cannam@95: /* C = A + B */ cannam@95: static void raadd(R *c, R *a, R *b, int n) cannam@95: { cannam@95: int i; cannam@95: cannam@95: for (i = 0; i < n; ++i) { cannam@95: c[i] = a[i] + b[i]; cannam@95: } cannam@95: } cannam@95: cannam@95: /* C = A - B */ cannam@95: static void rasub(R *c, R *a, R *b, int n) cannam@95: { cannam@95: int i; cannam@95: cannam@95: for (i = 0; i < n; ++i) { cannam@95: c[i] = a[i] - b[i]; cannam@95: } cannam@95: } cannam@95: cannam@95: /* B = rotate left A + rotate right A */ cannam@95: static void rarolr(R *b, R *a, int n, int nb, int na, cannam@95: r2r_kind_t k) cannam@95: { cannam@95: int isL0 = 0, isL1 = 0, isR0 = 0, isR1 = 0; cannam@95: int i, ib, ia; cannam@95: cannam@95: for (ib = 0; ib < nb; ++ib) { cannam@95: for (i = 0; i < n - 1; ++i) cannam@95: for (ia = 0; ia < na; ++ia) cannam@95: b[(ib * n + i) * na + ia] = cannam@95: a[(ib * n + i + 1) * na + ia]; cannam@95: cannam@95: /* ugly switch to do boundary conditions for various r2r types */ cannam@95: switch (k) { cannam@95: /* periodic boundaries */ cannam@95: case R2R_DHT: cannam@95: case R2R_R2HC: cannam@95: for (ia = 0; ia < na; ++ia) { cannam@95: b[(ib * n + n - 1) * na + ia] = cannam@95: a[(ib * n + 0) * na + ia]; cannam@95: b[(ib * n + 0) * na + ia] += cannam@95: a[(ib * n + n - 1) * na + ia]; cannam@95: } cannam@95: break; cannam@95: cannam@95: case R2R_HC2R: /* ugh (hermitian halfcomplex boundaries) */ cannam@95: if (n > 2) { cannam@95: if (n % 2 == 0) cannam@95: for (ia = 0; ia < na; ++ia) { cannam@95: b[(ib * n + n - 1) * na + ia] = 0.0; cannam@95: b[(ib * n + 0) * na + ia] += cannam@95: a[(ib * n + 1) * na + ia]; cannam@95: b[(ib * n + n/2) * na + ia] += cannam@95: + a[(ib * n + n/2 - 1) * na + ia] cannam@95: - a[(ib * n + n/2 + 1) * na + ia]; cannam@95: b[(ib * n + n/2 + 1) * na + ia] += cannam@95: - a[(ib * n + n/2) * na + ia]; cannam@95: } cannam@95: else cannam@95: for (ia = 0; ia < na; ++ia) { cannam@95: b[(ib * n + n - 1) * na + ia] = 0.0; cannam@95: b[(ib * n + 0) * na + ia] += cannam@95: a[(ib * n + 1) * na + ia]; cannam@95: b[(ib * n + n/2) * na + ia] += cannam@95: + a[(ib * n + n/2) * na + ia] cannam@95: - a[(ib * n + n/2 + 1) * na + ia]; cannam@95: b[(ib * n + n/2 + 1) * na + ia] += cannam@95: - a[(ib * n + n/2 + 1) * na + ia] cannam@95: - a[(ib * n + n/2) * na + ia]; cannam@95: } cannam@95: } else /* n <= 2 */ { cannam@95: for (ia = 0; ia < na; ++ia) { cannam@95: b[(ib * n + n - 1) * na + ia] = cannam@95: a[(ib * n + 0) * na + ia]; cannam@95: b[(ib * n + 0) * na + ia] += cannam@95: a[(ib * n + n - 1) * na + ia]; cannam@95: } cannam@95: } cannam@95: break; cannam@95: cannam@95: /* various even/odd boundary conditions */ cannam@95: case R2R_REDFT00: cannam@95: isL1 = isR1 = 1; cannam@95: goto mirrors; cannam@95: case R2R_REDFT01: cannam@95: isL1 = 1; cannam@95: goto mirrors; cannam@95: case R2R_REDFT10: cannam@95: isL0 = isR0 = 1; cannam@95: goto mirrors; cannam@95: case R2R_REDFT11: cannam@95: isL0 = 1; cannam@95: isR0 = -1; cannam@95: goto mirrors; cannam@95: case R2R_RODFT00: cannam@95: goto mirrors; cannam@95: case R2R_RODFT01: cannam@95: isR1 = 1; cannam@95: goto mirrors; cannam@95: case R2R_RODFT10: cannam@95: isL0 = isR0 = -1; cannam@95: goto mirrors; cannam@95: case R2R_RODFT11: cannam@95: isL0 = -1; cannam@95: isR0 = 1; cannam@95: goto mirrors; cannam@95: cannam@95: mirrors: cannam@95: cannam@95: for (ia = 0; ia < na; ++ia) cannam@95: b[(ib * n + n - 1) * na + ia] = cannam@95: isR0 * a[(ib * n + n - 1) * na + ia] cannam@95: + (n > 1 ? isR1 * a[(ib * n + n - 2) * na + ia] cannam@95: : 0); cannam@95: cannam@95: for (ia = 0; ia < na; ++ia) cannam@95: b[(ib * n) * na + ia] += cannam@95: isL0 * a[(ib * n) * na + ia] cannam@95: + (n > 1 ? isL1 * a[(ib * n + 1) * na + ia] : 0); cannam@95: cannam@95: } cannam@95: cannam@95: for (i = 1; i < n; ++i) cannam@95: for (ia = 0; ia < na; ++ia) cannam@95: b[(ib * n + i) * na + ia] += cannam@95: a[(ib * n + i - 1) * na + ia]; cannam@95: } cannam@95: } cannam@95: cannam@95: static void raphase_shift(R *b, R *a, int n, int nb, int na, cannam@95: int n0, int k0, trigfun t) cannam@95: { cannam@95: int j, jb, ja; cannam@95: cannam@95: for (jb = 0; jb < nb; ++jb) cannam@95: for (j = 0; j < n; ++j) { cannam@95: trigreal c = 2.0 * t(1, j + k0, n0); cannam@95: cannam@95: for (ja = 0; ja < na; ++ja) { cannam@95: int k = (jb * n + j) * na + ja; cannam@95: b[k] = a[k] * c; cannam@95: } cannam@95: } cannam@95: } cannam@95: cannam@95: /* A = alpha * A (real, in place) */ cannam@95: static void rascale(R *a, R alpha, int n) cannam@95: { cannam@95: int i; cannam@95: cannam@95: for (i = 0; i < n; ++i) { cannam@95: a[i] *= alpha; cannam@95: } cannam@95: } cannam@95: cannam@95: /* cannam@95: * compute rdft: cannam@95: */ cannam@95: cannam@95: /* copy real A into real B, using output stride of A and input stride of B */ cannam@95: typedef struct { cannam@95: dotens2_closure k; cannam@95: R *ra; cannam@95: R *rb; cannam@95: } cpyr_closure; cannam@95: cannam@95: static void cpyr0(dotens2_closure *k_, cannam@95: int indxa, int ondxa, int indxb, int ondxb) cannam@95: { cannam@95: cpyr_closure *k = (cpyr_closure *)k_; cannam@95: k->rb[indxb] = k->ra[ondxa]; cannam@95: UNUSED(indxa); UNUSED(ondxb); cannam@95: } cannam@95: cannam@95: static void cpyr(R *ra, bench_tensor *sza, R *rb, bench_tensor *szb) cannam@95: { cannam@95: cpyr_closure k; cannam@95: k.k.apply = cpyr0; cannam@95: k.ra = ra; k.rb = rb; cannam@95: bench_dotens2(sza, szb, &k.k); cannam@95: } cannam@95: cannam@95: static void dofft(info *nfo, R *in, R *out) cannam@95: { cannam@95: cpyr(in, nfo->pckdsz, (R *) nfo->p->in, nfo->totalsz); cannam@95: after_problem_rcopy_from(nfo->p, (bench_real *)nfo->p->in); cannam@95: doit(1, nfo->p); cannam@95: after_problem_rcopy_to(nfo->p, (bench_real *)nfo->p->out); cannam@95: cpyr((R *) nfo->p->out, nfo->totalsz, out, nfo->pckdsz); cannam@95: } cannam@95: cannam@95: static double racmp(R *a, R *b, int n, const char *test, double tol) cannam@95: { cannam@95: double d = raerror(a, b, n); cannam@95: if (d > tol) { cannam@95: ovtpvt_err("Found relative error %e (%s)\n", d, test); cannam@95: { cannam@95: int i, N; cannam@95: N = n > 300 && verbose <= 2 ? 300 : n; cannam@95: for (i = 0; i < N; ++i) cannam@95: ovtpvt_err("%8d %16.12f %16.12f\n", i, cannam@95: (double) a[i], cannam@95: (double) b[i]); cannam@95: } cannam@95: bench_exit(EXIT_FAILURE); cannam@95: } cannam@95: return d; cannam@95: } cannam@95: cannam@95: /***********************************************************************/ cannam@95: cannam@95: typedef struct { cannam@95: int n; /* physical size */ cannam@95: int n0; /* "logical" transform size */ cannam@95: int i0, k0; /* shifts of input/output */ cannam@95: trigfun ti, ts; /* impulse/shift trig functions */ cannam@95: } dim_stuff; cannam@95: cannam@95: static void impulse_response(int rnk, dim_stuff *d, R impulse_amp, cannam@95: R *A, int N) cannam@95: { cannam@95: if (rnk == 0) cannam@95: A[0] = impulse_amp; cannam@95: else { cannam@95: int i; cannam@95: N /= d->n; cannam@95: for (i = 0; i < d->n; ++i) { cannam@95: impulse_response(rnk - 1, d + 1, cannam@95: impulse_amp * d->ti(d->i0, d->k0 + i, d->n0), cannam@95: A + i * N, N); cannam@95: } cannam@95: } cannam@95: } cannam@95: cannam@95: /***************************************************************************/ cannam@95: cannam@95: /* cannam@95: * Implementation of the FFT tester described in cannam@95: * cannam@95: * Funda Ergün. Testing multivariate linear functions: Overcoming the cannam@95: * generator bottleneck. In Proceedings of the Twenty-Seventh Annual cannam@95: * ACM Symposium on the Theory of Computing, pages 407-416, Las Vegas, cannam@95: * Nevada, 29 May--1 June 1995. cannam@95: * cannam@95: * Also: F. Ergun, S. R. Kumar, and D. Sivakumar, "Self-testing without cannam@95: * the generator bottleneck," SIAM J. on Computing 29 (5), 1630-51 (2000). cannam@95: */ cannam@95: cannam@95: static double rlinear(int n, info *nfo, R *inA, R *inB, R *inC, R *outA, cannam@95: R *outB, R *outC, R *tmp, int rounds, double tol) cannam@95: { cannam@95: double e = 0.0; cannam@95: int j; cannam@95: cannam@95: for (j = 0; j < rounds; ++j) { cannam@95: R alpha, beta; cannam@95: alpha = mydrand(); cannam@95: beta = mydrand(); cannam@95: rarand(inA, n); cannam@95: rarand(inB, n); cannam@95: dofft(nfo, inA, outA); cannam@95: dofft(nfo, inB, outB); cannam@95: cannam@95: rascale(outA, alpha, n); cannam@95: rascale(outB, beta, n); cannam@95: raadd(tmp, outA, outB, n); cannam@95: rascale(inA, alpha, n); cannam@95: rascale(inB, beta, n); cannam@95: raadd(inC, inA, inB, n); cannam@95: dofft(nfo, inC, outC); cannam@95: cannam@95: e = dmax(e, racmp(outC, tmp, n, "linear", tol)); cannam@95: } cannam@95: return e; cannam@95: } cannam@95: cannam@95: static double rimpulse(dim_stuff *d, R impulse_amp, cannam@95: int n, int vecn, info *nfo, cannam@95: R *inA, R *inB, R *inC, cannam@95: R *outA, R *outB, R *outC, cannam@95: R *tmp, int rounds, double tol) cannam@95: { cannam@95: double e = 0.0; cannam@95: int N = n * vecn; cannam@95: int i; cannam@95: int j; cannam@95: cannam@95: /* test 2: check that the unit impulse is transformed properly */ cannam@95: cannam@95: for (i = 0; i < N; ++i) { cannam@95: /* pls */ cannam@95: inA[i] = 0.0; cannam@95: } cannam@95: for (i = 0; i < vecn; ++i) { cannam@95: inA[i * n] = (i+1) / (double)(vecn+1); cannam@95: cannam@95: /* transform of the pls */ cannam@95: impulse_response(nfo->probsz->rnk, d, impulse_amp * inA[i * n], cannam@95: outA + i * n, n); cannam@95: } cannam@95: cannam@95: dofft(nfo, inA, tmp); cannam@95: e = dmax(e, racmp(tmp, outA, N, "impulse 1", tol)); cannam@95: cannam@95: for (j = 0; j < rounds; ++j) { cannam@95: rarand(inB, N); cannam@95: rasub(inC, inA, inB, N); cannam@95: dofft(nfo, inB, outB); cannam@95: dofft(nfo, inC, outC); cannam@95: raadd(tmp, outB, outC, N); cannam@95: e = dmax(e, racmp(tmp, outA, N, "impulse", tol)); cannam@95: } cannam@95: return e; cannam@95: } cannam@95: cannam@95: static double t_shift(int n, int vecn, info *nfo, cannam@95: R *inA, R *inB, R *outA, R *outB, R *tmp, cannam@95: int rounds, double tol, cannam@95: dim_stuff *d) cannam@95: { cannam@95: double e = 0.0; cannam@95: int nb, na, dim, N = n * vecn; cannam@95: int i, j; cannam@95: bench_tensor *sz = nfo->probsz; cannam@95: cannam@95: /* test 3: check the time-shift property */ cannam@95: /* the paper performs more tests, but this code should be fine too */ cannam@95: cannam@95: nb = 1; cannam@95: na = n; cannam@95: cannam@95: /* check shifts across all SZ dimensions */ cannam@95: for (dim = 0; dim < sz->rnk; ++dim) { cannam@95: int ncur = sz->dims[dim].n; cannam@95: cannam@95: na /= ncur; cannam@95: cannam@95: for (j = 0; j < rounds; ++j) { cannam@95: rarand(inA, N); cannam@95: cannam@95: for (i = 0; i < vecn; ++i) { cannam@95: rarolr(inB + i * n, inA + i*n, ncur, nb,na, cannam@95: nfo->p->k[dim]); cannam@95: } cannam@95: dofft(nfo, inA, outA); cannam@95: dofft(nfo, inB, outB); cannam@95: for (i = 0; i < vecn; ++i) cannam@95: raphase_shift(tmp + i * n, outA + i * n, ncur, cannam@95: nb, na, d[dim].n0, d[dim].k0, d[dim].ts); cannam@95: e = dmax(e, racmp(tmp, outB, N, "time shift", tol)); cannam@95: } cannam@95: cannam@95: nb *= ncur; cannam@95: } cannam@95: return e; cannam@95: } cannam@95: cannam@95: /***********************************************************************/ cannam@95: cannam@95: void verify_r2r(bench_problem *p, int rounds, double tol, errors *e) cannam@95: { cannam@95: R *inA, *inB, *inC, *outA, *outB, *outC, *tmp; cannam@95: info nfo; cannam@95: int n, vecn, N; cannam@95: double impulse_amp = 1.0; cannam@95: dim_stuff *d; cannam@95: int i; cannam@95: cannam@95: if (rounds == 0) cannam@95: rounds = 20; /* default value */ cannam@95: cannam@95: n = tensor_sz(p->sz); cannam@95: vecn = tensor_sz(p->vecsz); cannam@95: N = n * vecn; cannam@95: cannam@95: d = (dim_stuff *) bench_malloc(sizeof(dim_stuff) * p->sz->rnk); cannam@95: for (i = 0; i < p->sz->rnk; ++i) { cannam@95: int n0, i0, k0; cannam@95: trigfun ti, ts; cannam@95: cannam@95: d[i].n = n0 = p->sz->dims[i].n; cannam@95: if (p->k[i] > R2R_DHT) cannam@95: n0 = 2 * (n0 + (p->k[i] == R2R_REDFT00 ? -1 : cannam@95: (p->k[i] == R2R_RODFT00 ? 1 : 0))); cannam@95: cannam@95: switch (p->k[i]) { cannam@95: case R2R_R2HC: cannam@95: i0 = k0 = 0; cannam@95: ti = realhalf; cannam@95: ts = coshalf; cannam@95: break; cannam@95: case R2R_DHT: cannam@95: i0 = k0 = 0; cannam@95: ti = unity; cannam@95: ts = cos00; cannam@95: break; cannam@95: case R2R_HC2R: cannam@95: i0 = k0 = 0; cannam@95: ti = unity; cannam@95: ts = cos00; cannam@95: break; cannam@95: case R2R_REDFT00: cannam@95: i0 = k0 = 0; cannam@95: ti = ts = cos00; cannam@95: break; cannam@95: case R2R_REDFT01: cannam@95: i0 = k0 = 0; cannam@95: ti = ts = cos01; cannam@95: break; cannam@95: case R2R_REDFT10: cannam@95: i0 = k0 = 0; cannam@95: ti = cos10; impulse_amp *= 2.0; cannam@95: ts = cos00; cannam@95: break; cannam@95: case R2R_REDFT11: cannam@95: i0 = k0 = 0; cannam@95: ti = cos11; impulse_amp *= 2.0; cannam@95: ts = cos01; cannam@95: break; cannam@95: case R2R_RODFT00: cannam@95: i0 = k0 = 1; cannam@95: ti = sin00; impulse_amp *= 2.0; cannam@95: ts = cos00; cannam@95: break; cannam@95: case R2R_RODFT01: cannam@95: i0 = 1; k0 = 0; cannam@95: ti = sin01; impulse_amp *= n == 1 ? 1.0 : 2.0; cannam@95: ts = cos01; cannam@95: break; cannam@95: case R2R_RODFT10: cannam@95: i0 = 0; k0 = 1; cannam@95: ti = sin10; impulse_amp *= 2.0; cannam@95: ts = cos00; cannam@95: break; cannam@95: case R2R_RODFT11: cannam@95: i0 = k0 = 0; cannam@95: ti = sin11; impulse_amp *= 2.0; cannam@95: ts = cos01; cannam@95: break; cannam@95: default: cannam@95: BENCH_ASSERT(0); cannam@95: return; cannam@95: } cannam@95: cannam@95: d[i].n0 = n0; cannam@95: d[i].i0 = i0; cannam@95: d[i].k0 = k0; cannam@95: d[i].ti = ti; cannam@95: d[i].ts = ts; cannam@95: } cannam@95: cannam@95: cannam@95: inA = (R *) bench_malloc(N * sizeof(R)); cannam@95: inB = (R *) bench_malloc(N * sizeof(R)); cannam@95: inC = (R *) bench_malloc(N * sizeof(R)); cannam@95: outA = (R *) bench_malloc(N * sizeof(R)); cannam@95: outB = (R *) bench_malloc(N * sizeof(R)); cannam@95: outC = (R *) bench_malloc(N * sizeof(R)); cannam@95: tmp = (R *) bench_malloc(N * sizeof(R)); cannam@95: cannam@95: nfo.p = p; cannam@95: nfo.probsz = p->sz; cannam@95: nfo.totalsz = tensor_append(p->vecsz, nfo.probsz); cannam@95: nfo.pckdsz = verify_pack(nfo.totalsz, 1); cannam@95: nfo.pckdvecsz = verify_pack(p->vecsz, tensor_sz(nfo.probsz)); cannam@95: cannam@95: e->i = rimpulse(d, impulse_amp, n, vecn, &nfo, cannam@95: inA, inB, inC, outA, outB, outC, tmp, rounds, tol); cannam@95: e->l = rlinear(N, &nfo, inA, inB, inC, outA, outB, outC, tmp, rounds,tol); cannam@95: e->s = t_shift(n, vecn, &nfo, inA, inB, outA, outB, tmp, cannam@95: rounds, tol, d); cannam@95: cannam@95: /* grr, verify-lib.c:preserves_input() only works for complex */ cannam@95: if (!p->in_place && !p->destroy_input) { cannam@95: bench_tensor *totalsz_swap, *pckdsz_swap; cannam@95: totalsz_swap = tensor_copy_swapio(nfo.totalsz); cannam@95: pckdsz_swap = tensor_copy_swapio(nfo.pckdsz); cannam@95: cannam@95: for (i = 0; i < rounds; ++i) { cannam@95: rarand(inA, N); cannam@95: dofft(&nfo, inA, outB); cannam@95: cpyr((R *) nfo.p->in, totalsz_swap, inB, pckdsz_swap); cannam@95: racmp(inB, inA, N, "preserves_input", 0.0); cannam@95: } cannam@95: cannam@95: tensor_destroy(totalsz_swap); cannam@95: tensor_destroy(pckdsz_swap); cannam@95: } cannam@95: cannam@95: tensor_destroy(nfo.totalsz); cannam@95: tensor_destroy(nfo.pckdsz); cannam@95: tensor_destroy(nfo.pckdvecsz); cannam@95: bench_free(tmp); cannam@95: bench_free(outC); cannam@95: bench_free(outB); cannam@95: bench_free(outA); cannam@95: bench_free(inC); cannam@95: bench_free(inB); cannam@95: bench_free(inA); cannam@95: bench_free(d); cannam@95: } cannam@95: cannam@95: cannam@95: typedef struct { cannam@95: dofft_closure k; cannam@95: bench_problem *p; cannam@95: int n0; cannam@95: } dofft_r2r_closure; cannam@95: cannam@95: static void cpyr1(int n, R *in, int is, R *out, int os, R scale) cannam@95: { cannam@95: int i; cannam@95: for (i = 0; i < n; ++i) cannam@95: out[i * os] = in[i * is] * scale; cannam@95: } cannam@95: cannam@95: static void mke00(C *a, int n, int c) cannam@95: { cannam@95: int i; cannam@95: for (i = 1; i + i < n; ++i) cannam@95: a[n - i][c] = a[i][c]; cannam@95: } cannam@95: cannam@95: static void mkre00(C *a, int n) cannam@95: { cannam@95: mkreal(a, n); cannam@95: mke00(a, n, 0); cannam@95: } cannam@95: cannam@95: static void mkimag(C *a, int n) cannam@95: { cannam@95: int i; cannam@95: for (i = 0; i < n; ++i) cannam@95: c_re(a[i]) = 0.0; cannam@95: } cannam@95: cannam@95: static void mko00(C *a, int n, int c) cannam@95: { cannam@95: int i; cannam@95: a[0][c] = 0.0; cannam@95: for (i = 1; i + i < n; ++i) cannam@95: a[n - i][c] = -a[i][c]; cannam@95: if (i + i == n) cannam@95: a[i][c] = 0.0; cannam@95: } cannam@95: cannam@95: static void mkro00(C *a, int n) cannam@95: { cannam@95: mkreal(a, n); cannam@95: mko00(a, n, 0); cannam@95: } cannam@95: cannam@95: static void mkio00(C *a, int n) cannam@95: { cannam@95: mkimag(a, n); cannam@95: mko00(a, n, 1); cannam@95: } cannam@95: cannam@95: static void mkre01(C *a, int n) /* n should be be multiple of 4 */ cannam@95: { cannam@95: R a0; cannam@95: a0 = c_re(a[0]); cannam@95: mko00(a, n/2, 0); cannam@95: c_re(a[n/2]) = -(c_re(a[0]) = a0); cannam@95: mkre00(a, n); cannam@95: } cannam@95: cannam@95: static void mkro01(C *a, int n) /* n should be be multiple of 4 */ cannam@95: { cannam@95: c_re(a[0]) = c_im(a[0]) = 0.0; cannam@95: mkre00(a, n/2); cannam@95: mkro00(a, n); cannam@95: } cannam@95: cannam@95: static void mkoddonly(C *a, int n) cannam@95: { cannam@95: int i; cannam@95: for (i = 0; i < n; i += 2) cannam@95: c_re(a[i]) = c_im(a[i]) = 0.0; cannam@95: } cannam@95: cannam@95: static void mkre10(C *a, int n) cannam@95: { cannam@95: mkoddonly(a, n); cannam@95: mkre00(a, n); cannam@95: } cannam@95: cannam@95: static void mkio10(C *a, int n) cannam@95: { cannam@95: mkoddonly(a, n); cannam@95: mkio00(a, n); cannam@95: } cannam@95: cannam@95: static void mkre11(C *a, int n) cannam@95: { cannam@95: mkoddonly(a, n); cannam@95: mko00(a, n/2, 0); cannam@95: mkre00(a, n); cannam@95: } cannam@95: cannam@95: static void mkro11(C *a, int n) cannam@95: { cannam@95: mkoddonly(a, n); cannam@95: mkre00(a, n/2); cannam@95: mkro00(a, n); cannam@95: } cannam@95: cannam@95: static void mkio11(C *a, int n) cannam@95: { cannam@95: mkoddonly(a, n); cannam@95: mke00(a, n/2, 1); cannam@95: mkio00(a, n); cannam@95: } cannam@95: cannam@95: static void r2r_apply(dofft_closure *k_, bench_complex *in, bench_complex *out) cannam@95: { cannam@95: dofft_r2r_closure *k = (dofft_r2r_closure *)k_; cannam@95: bench_problem *p = k->p; cannam@95: bench_real *ri, *ro; cannam@95: int n, is, os; cannam@95: cannam@95: n = p->sz->dims[0].n; cannam@95: is = p->sz->dims[0].is; cannam@95: os = p->sz->dims[0].os; cannam@95: cannam@95: ri = (bench_real *) p->in; cannam@95: ro = (bench_real *) p->out; cannam@95: cannam@95: switch (p->k[0]) { cannam@95: case R2R_R2HC: cannam@95: cpyr1(n, &c_re(in[0]), 2, ri, is, 1.0); cannam@95: break; cannam@95: case R2R_HC2R: cannam@95: cpyr1(n/2 + 1, &c_re(in[0]), 2, ri, is, 1.0); cannam@95: cpyr1((n+1)/2 - 1, &c_im(in[n-1]), -2, ri + is*(n-1), -is, 1.0); cannam@95: break; cannam@95: case R2R_REDFT00: cannam@95: cpyr1(n, &c_re(in[0]), 2, ri, is, 1.0); cannam@95: break; cannam@95: case R2R_RODFT00: cannam@95: cpyr1(n, &c_re(in[1]), 2, ri, is, 1.0); cannam@95: break; cannam@95: case R2R_REDFT01: cannam@95: cpyr1(n, &c_re(in[0]), 2, ri, is, 1.0); cannam@95: break; cannam@95: case R2R_REDFT10: cannam@95: cpyr1(n, &c_re(in[1]), 4, ri, is, 1.0); cannam@95: break; cannam@95: case R2R_RODFT01: cannam@95: cpyr1(n, &c_re(in[1]), 2, ri, is, 1.0); cannam@95: break; cannam@95: case R2R_RODFT10: cannam@95: cpyr1(n, &c_im(in[1]), 4, ri, is, 1.0); cannam@95: break; cannam@95: case R2R_REDFT11: cannam@95: cpyr1(n, &c_re(in[1]), 4, ri, is, 1.0); cannam@95: break; cannam@95: case R2R_RODFT11: cannam@95: cpyr1(n, &c_re(in[1]), 4, ri, is, 1.0); cannam@95: break; cannam@95: default: cannam@95: BENCH_ASSERT(0); /* not yet implemented */ cannam@95: } cannam@95: cannam@95: after_problem_rcopy_from(p, ri); cannam@95: doit(1, p); cannam@95: after_problem_rcopy_to(p, ro); cannam@95: cannam@95: switch (p->k[0]) { cannam@95: case R2R_R2HC: cannam@95: if (k->k.recopy_input) cannam@95: cpyr1(n, ri, is, &c_re(in[0]), 2, 1.0); cannam@95: cpyr1(n/2 + 1, ro, os, &c_re(out[0]), 2, 1.0); cannam@95: cpyr1((n+1)/2 - 1, ro + os*(n-1), -os, &c_im(out[1]), 2, 1.0); cannam@95: c_im(out[0]) = 0.0; cannam@95: if (n % 2 == 0) cannam@95: c_im(out[n/2]) = 0.0; cannam@95: mkhermitian1(out, n); cannam@95: break; cannam@95: case R2R_HC2R: cannam@95: if (k->k.recopy_input) { cannam@95: cpyr1(n/2 + 1, ri, is, &c_re(in[0]), 2, 1.0); cannam@95: cpyr1((n+1)/2 - 1, ri + is*(n-1), -is, &c_im(in[1]), 2,1.0); cannam@95: } cannam@95: cpyr1(n, ro, os, &c_re(out[0]), 2, 1.0); cannam@95: mkreal(out, n); cannam@95: break; cannam@95: case R2R_REDFT00: cannam@95: if (k->k.recopy_input) cannam@95: cpyr1(n, ri, is, &c_re(in[0]), 2, 1.0); cannam@95: cpyr1(n, ro, os, &c_re(out[0]), 2, 1.0); cannam@95: mkre00(out, k->n0); cannam@95: break; cannam@95: case R2R_RODFT00: cannam@95: if (k->k.recopy_input) cannam@95: cpyr1(n, ri, is, &c_im(in[1]), 2, -1.0); cannam@95: cpyr1(n, ro, os, &c_im(out[1]), 2, -1.0); cannam@95: mkio00(out, k->n0); cannam@95: break; cannam@95: case R2R_REDFT01: cannam@95: if (k->k.recopy_input) cannam@95: cpyr1(n, ri, is, &c_re(in[0]), 2, 1.0); cannam@95: cpyr1(n, ro, os, &c_re(out[1]), 4, 2.0); cannam@95: mkre10(out, k->n0); cannam@95: break; cannam@95: case R2R_REDFT10: cannam@95: if (k->k.recopy_input) cannam@95: cpyr1(n, ri, is, &c_re(in[1]), 4, 2.0); cannam@95: cpyr1(n, ro, os, &c_re(out[0]), 2, 1.0); cannam@95: mkre01(out, k->n0); cannam@95: break; cannam@95: case R2R_RODFT01: cannam@95: if (k->k.recopy_input) cannam@95: cpyr1(n, ri, is, &c_re(in[1]), 2, 1.0); cannam@95: cpyr1(n, ro, os, &c_im(out[1]), 4, -2.0); cannam@95: mkio10(out, k->n0); cannam@95: break; cannam@95: case R2R_RODFT10: cannam@95: if (k->k.recopy_input) cannam@95: cpyr1(n, ri, is, &c_im(in[1]), 4, -2.0); cannam@95: cpyr1(n, ro, os, &c_re(out[1]), 2, 1.0); cannam@95: mkro01(out, k->n0); cannam@95: break; cannam@95: case R2R_REDFT11: cannam@95: if (k->k.recopy_input) cannam@95: cpyr1(n, ri, is, &c_re(in[1]), 4, 2.0); cannam@95: cpyr1(n, ro, os, &c_re(out[1]), 4, 2.0); cannam@95: mkre11(out, k->n0); cannam@95: break; cannam@95: case R2R_RODFT11: cannam@95: if (k->k.recopy_input) cannam@95: cpyr1(n, ri, is, &c_im(in[1]), 4, -2.0); cannam@95: cpyr1(n, ro, os, &c_im(out[1]), 4, -2.0); cannam@95: mkio11(out, k->n0); cannam@95: break; cannam@95: default: cannam@95: BENCH_ASSERT(0); /* not yet implemented */ cannam@95: } cannam@95: } cannam@95: cannam@95: void accuracy_r2r(bench_problem *p, int rounds, int impulse_rounds, cannam@95: double t[6]) cannam@95: { cannam@95: dofft_r2r_closure k; cannam@95: int n, n0 = 1; cannam@95: C *a, *b; cannam@95: aconstrain constrain = 0; cannam@95: cannam@95: BENCH_ASSERT(p->kind == PROBLEM_R2R); cannam@95: BENCH_ASSERT(p->sz->rnk == 1); cannam@95: BENCH_ASSERT(p->vecsz->rnk == 0); cannam@95: cannam@95: k.k.apply = r2r_apply; cannam@95: k.k.recopy_input = 0; cannam@95: k.p = p; cannam@95: n = tensor_sz(p->sz); cannam@95: cannam@95: switch (p->k[0]) { cannam@95: case R2R_R2HC: constrain = mkreal; n0 = n; break; cannam@95: case R2R_HC2R: constrain = mkhermitian1; n0 = n; break; cannam@95: case R2R_REDFT00: constrain = mkre00; n0 = 2*(n-1); break; cannam@95: case R2R_RODFT00: constrain = mkro00; n0 = 2*(n+1); break; cannam@95: case R2R_REDFT01: constrain = mkre01; n0 = 4*n; break; cannam@95: case R2R_REDFT10: constrain = mkre10; n0 = 4*n; break; cannam@95: case R2R_RODFT01: constrain = mkro01; n0 = 4*n; break; cannam@95: case R2R_RODFT10: constrain = mkio10; n0 = 4*n; break; cannam@95: case R2R_REDFT11: constrain = mkre11; n0 = 8*n; break; cannam@95: case R2R_RODFT11: constrain = mkro11; n0 = 8*n; break; cannam@95: default: BENCH_ASSERT(0); /* not yet implemented */ cannam@95: } cannam@95: k.n0 = n0; cannam@95: cannam@95: a = (C *) bench_malloc(n0 * sizeof(C)); cannam@95: b = (C *) bench_malloc(n0 * sizeof(C)); cannam@95: accuracy_test(&k.k, constrain, -1, n0, a, b, rounds, impulse_rounds, t); cannam@95: bench_free(b); cannam@95: bench_free(a); cannam@95: }