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