rfftw_test.c

FFTW, a collection of fast C routines to compute the Discrete Fourier Transform in one or more dime

     rfftwnd_destroy_plan(ip);

     fftw_free(out2);
     fftw_free(in2);
     fftw_free(in1);
}

void testnd_correctness(struct size sz, fftw_direction dir,
			int alt_api, int specific, int force_buf)
{
     fftwnd_plan validated_plan;

     if (dir != FFTW_FORWARD)
	  return;
     if (force_buf)
	  return;

     validated_plan = fftwnd_create_plan(sz.rank, sz.narray, 
					 dir, measure_flag | wisdom_flag);
     CHECK(validated_plan != NULL, "can't create plan");

     testnd_out_of_place(sz.rank, sz.narray, validated_plan);
     testnd_in_place(sz.rank, sz.narray,
		     validated_plan, alt_api, specific);

     fftwnd_destroy_plan(validated_plan);
}

/*************************************************
 * planner tests
 *************************************************/

void test_planner(int rank)
{
     /* 
      * create and destroy many plans, at random.  Check the
      * garbage-collecting allocator of twiddle factors 
      */
     int i, dim;
     int r, s;
     fftw_plan p[PLANNER_TEST_SIZE];
     fftwnd_plan pnd[PLANNER_TEST_SIZE];
     int *narr, maxdim;

     chk_mem_leak = 0;
     verbose--;

     please_wait();
     if (rank < 1)
	  rank = 1;

     narr = (int *) fftw_malloc(rank * sizeof(int));

     maxdim = (int) pow(8192.0, 1.0/rank);

     for (i = 0; i < PLANNER_TEST_SIZE; ++i) {
	  p[i] = (fftw_plan) 0;
	  pnd[i] = (fftwnd_plan) 0;
     }

     for (i = 0; i < PLANNER_TEST_SIZE * PLANNER_TEST_SIZE; ++i) {
	  r = rand();
	  if (r < 0)
	       r = -r;
	  r = r % PLANNER_TEST_SIZE;

	  for (dim = 0; dim < rank; ++dim) {
	       do {
		    s = rand();
		    if (s < 0)
			 s = -s;
		    s = s % maxdim + 1;
	       } while (s == 0);
	       narr[dim] = s;
	  }

	  if (rank == 1) {
	       if (p[r])
		    rfftw_destroy_plan(p[r]);

	       p[r] = rfftw_create_plan(narr[0], random_dir(), measure_flag |
					wisdom_flag);
	       if (paranoid && narr[0] < 200)
		    test_correctness(narr[0]);
	  }
	  if (pnd[r])
	       rfftwnd_destroy_plan(pnd[r]);

	  pnd[r] = rfftwnd_create_plan(rank, narr,
				       random_dir(), measure_flag |
				       wisdom_flag);

	  if (i % (PLANNER_TEST_SIZE * PLANNER_TEST_SIZE / 20) == 0) {
	       WHEN_VERBOSE(0, printf("test planner: so far so good\n"));
	       WHEN_VERBOSE(0, printf("test planner: iteration %d out of %d\n",
			      i, PLANNER_TEST_SIZE * PLANNER_TEST_SIZE));
	  }
     }

     for (i = 0; i < PLANNER_TEST_SIZE; ++i) {
	  if (p[i])
	       rfftw_destroy_plan(p[i]);
	  if (pnd[i])
	       rfftwnd_destroy_plan(pnd[i]);
     }

     fftw_free(narr);
     verbose++;
     chk_mem_leak = 1;
}


/*************************************************
 * Ergun's test for real->complex transforms
 *************************************************/
static void rfill_random(fftw_real *a, int n)
{
     int i;

     for (i = 0; i < n; ++i) {
	  a[i] = DRAND();
     }
}

static void rarray_copy(fftw_real *out, fftw_real *in, int n)
{
     int i;

     for (i = 0; i < n; ++i)
	  out[i] = in[i];
}

/* C = A + B */
void rarray_add(fftw_real *C, fftw_real *A, fftw_real *B, int n)
{
     int i;

     for (i = 0; i < n; ++i) {
	  C[i] = A[i] + B[i];
     }
}

/* C = A - B */
void rarray_sub(fftw_real *C, fftw_real *A, fftw_real *B, int n)
{
     int i;

     for (i = 0; i < n; ++i) {
	  C[i] = A[i] - B[i];
     }
}

/* B = rotate left A */
void rarray_rol(fftw_real *B, fftw_real *A,
	       int n, int n_before, int n_after)
{
     int i, ib, ia;

     for (ib = 0; ib < n_before; ++ib) {
	  for (i = 0; i < n - 1; ++i)
	       for (ia = 0; ia < n_after; ++ia)
		    B[(ib * n + i) * n_after + ia] =
			A[(ib * n + i + 1) * n_after + ia];

	  for (ia = 0; ia < n_after; ++ia)
	       B[(ib * n + n - 1) * n_after + ia] = A[ib * n * n_after + ia];
     }
}

/* A = alpha * B  (out of place) */
void rarray_scale(fftw_real *A, fftw_real *B, fftw_real alpha, int n)
{
     int i;

     for (i = 0; i < n; ++i) {
	  A[i] = B[i] * alpha;
     }
}

void rarray_compare(fftw_real *A, fftw_real *B, int n)
{
     double d = compute_error(A, 1, B, 1, n);
     if (d > TOLERANCE) {
	  fflush(stdout);
	  fprintf(stderr, "Found relative error %e\n", d);
	  fftw_die("failure in Ergun's verification procedure\n");
     }
}

/*
 * guaranteed out-of-place transform.  Does the necessary
 * copying if the plan is in-place.
 */
static void rfftw_out_of_place(fftw_plan plan, int n,
			      fftw_real *in, fftw_real *out)
{
     if (plan->flags & FFTW_IN_PLACE) {
	  rarray_copy(out, in, n);
	  rfftw(plan, 1, out, 1, n, (fftw_real *)0, 1, n);
     } else {
	  rfftw(plan, 1, in, 1, n, out, 1, n);
     }
}

/*
 * This is a real (as opposed to complex) variation of the FFT tester
 * described in
 *
 * Funda Erg黱. Testing multivariate linear functions: Overcoming the
 * generator bottleneck. In Proceedings of the Twenty-Seventh Annual
 * ACM Symposium on the Theory of Computing, pages 407-416, Las Vegas,
 * Nevada, 29 May--1 June 1995.
 */
void test_ergun(int n, fftw_direction dir, fftw_plan plan)
{
     fftw_real *inA, *inB, *inC, *outA, *outB, *outC;
     fftw_real *inA1, *inB1;
     fftw_real *tmp;
     int i;
     int rounds = 20;
     FFTW_TRIG_REAL twopin = FFTW_K2PI / (FFTW_TRIG_REAL) n;

     inA = (fftw_real *) fftw_malloc(n * sizeof(fftw_real));
     inB = (fftw_real *) fftw_malloc(n * sizeof(fftw_real));
     inA1 = (fftw_real *) fftw_malloc(n * sizeof(fftw_real));
     inB1 = (fftw_real *) fftw_malloc(n * sizeof(fftw_real));
     inC = (fftw_real *) fftw_malloc(n * sizeof(fftw_real));
     outA = (fftw_real *) fftw_malloc(n * sizeof(fftw_real));
     outB = (fftw_real *) fftw_malloc(n * sizeof(fftw_real));
     outC = (fftw_real *) fftw_malloc(n * sizeof(fftw_real));
     tmp = (fftw_real *) fftw_malloc(n * sizeof(fftw_real));

     WHEN_VERBOSE(2,
		  printf("Validating plan, n = %d, dir = %s\n", n,
			 dir == FFTW_REAL_TO_COMPLEX ? 
			 "REAL_TO_COMPLEX" : "COMPLEX_TO_REAL"));

     /* test 1: check linearity */
     for (i = 0; i < rounds; ++i) {
	  fftw_real alpha, beta;
	  alpha = DRAND();
	  beta = DRAND();
	  rfill_random(inA, n);
	  rfill_random(inB, n);
	  rarray_scale(inA1, inA, alpha, n);
	  rarray_scale(inB1, inB, beta, n);
	  rarray_add(inC, inA1, inB1, n);
	  rfftw_out_of_place(plan, n, inA, outA);
	  rfftw_out_of_place(plan, n, inB, outB);
	  rarray_scale(outA, outA, alpha, n);
	  rarray_scale(outB, outB, beta, n);
	  rarray_add(tmp, outA, outB, n);
	  rfftw_out_of_place(plan, n, inC, outC);
	  rarray_compare(outC, tmp, n);
     }

     /* test 2: check that the unit impulse is transformed properly */
     for (i = 0; i < n; ++i) {
	  /* impulse */
	  inA[i] = 0.0;
	  
	  /* transform of the impulse */
	  if (2 * i <= n)
	       outA[i] = 1.0;
	  else
	       outA[i] = 0.0;
     }
     inA[0] = 1.0;

     if (dir == FFTW_REAL_TO_COMPLEX) {
	  for (i = 0; i < rounds; ++i) {
	       rfill_random(inB, n);
	       rarray_sub(inC, inA, inB, n);
	       rfftw_out_of_place(plan, n, inB, outB);
	       rfftw_out_of_place(plan, n, inC, outC);
	       rarray_add(tmp, outB, outC, n);
	       rarray_compare(tmp, outA, n);
	  }
     } else {
	  for (i = 0; i < rounds; ++i) {
	       rfill_random(outB, n);
	       rarray_sub(outC, outA, outB, n);
	       rfftw_out_of_place(plan, n, outB, inB);
	       rfftw_out_of_place(plan, n, outC, inC);
	       rarray_add(tmp, inB, inC, n);
	       rarray_scale(tmp, tmp, 1.0 / ((double) n), n);
	       rarray_compare(tmp, inA, n);
	  }
     }

     /* test 3: check the time-shift property */
     /* the paper performs more tests, but this code should be fine too */
     if (dir == FFTW_REAL_TO_COMPLEX) {
	  for (i = 0; i < rounds; ++i) {
	       int j;

	       rfill_random(inA, n);
	       rarray_rol(inB, inA, n, 1, 1);
	       rfftw_out_of_place(plan, n, inA, outA);
	       rfftw_out_of_place(plan, n, inB, outB);
	       tmp[0] = outB[0];
	       for (j = 1; 2 * j < n; ++j) {
		    FFTW_TRIG_REAL s = dir * FFTW_TRIG_SIN(j * twopin);
		    FFTW_TRIG_REAL c = FFTW_TRIG_COS(j * twopin);
		    tmp[j] = outB[j] * c - outB[n - j] * s;
		    tmp[n - j] = outB[j] * s + outB[n - j] * c;
	       }
	       if (2 * j == n)
		    tmp[j] = -outB[j];

	       rarray_compare(tmp, outA, n);
	  }
     } else {
	  for (i = 0; i < rounds; ++i) {
	       int j;

	       rfill_random(inA, n);
	       inB[0] = inA[0];
	       for (j = 1; 2 * j < n; ++j) {
		    FFTW_TRIG_REAL s = dir * FFTW_TRIG_SIN(j * twopin);
		    FFTW_TRIG_REAL c = FFTW_TRIG_COS(j * twopin);
		    inB[j] = inA[j] * c - inA[n - j] * s;
		    inB[n - j] = inA[j] * s + inA[n - j] * c;
	       }
	       if (2 * j == n)
		    inB[j] = -inA[j];

	       rfftw_out_of_place(plan, n, inA, outA);
	       rfftw_out_of_place(plan, n, inB, outB);	       
	       rarray_rol(tmp, outA, n, 1, 1);
	       rarray_compare(tmp, outB, n);
	  }
     }

     WHEN_VERBOSE(2, printf("Validation done\n"));

     fftw_free(tmp);
     fftw_free(outC);
     fftw_free(outB);
     fftw_free(outA);
     fftw_free(inC);
     fftw_free(inB1);
     fftw_free(inA1);
     fftw_free(inB);
     fftw_free(inA);
}

static void rfftw_plan_hook_function(fftw_plan p)
{
     WHEN_VERBOSE(3, printf("Validating tentative plan\n"));
     WHEN_VERBOSE(3, fftw_print_plan(p));
     test_ergun(p->n, p->dir, p);
}

/*************************************************
 * test initialization
 *************************************************/

void test_init(int *argc, char ***argv)
{
}

void test_finish(void)
{
}

void enter_paranoid_mode(void)
{
     rfftw_plan_hook = rfftw_plan_hook_function;
}

int get_option(int argc, char **argv, char *argval, int argval_maxlen)
{
     return default_get_option(argc,argv,argval,argval_maxlen);
}