test_main.c

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

/*
 * Copyright (c) 1997-1999, 2003 Massachusetts Institute of Technology
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

/*
 * test_main.c: driver for test programs (linked with fftw_test.c/rfftw_test.c)
 */

/* $Id: test_main.c,v 1.39 2003/03/16 23:43:46 stevenj Exp $ */
#include "fftw-int.h"
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <ctype.h>

#include "test_main.h"

#ifdef HAVE_GETOPT
#  ifdef HAVE_GETOPT_H
#    include <getopt.h>
#  elif defined(HAVE_UNISTD_H)
#    include <unistd.h>
#  endif
#endif

double mydrand(void)
{
     double d = rand();
     return (d / (double) RAND_MAX) - .5;
}

/* return random 0 or non-zero */
int coinflip(void)
{
     return (rand() & 8192);	/* higher-order bits are often more random 
				 * 
				 */
}

/* parse a string of the form N1xN2x... and return a size structure */
struct size parse_size(char *s)
{
     struct size sz;
     int n;

     sz.rank = 0;
     sz.is_nd = 0;

     if (*s == 'x' || *s == 'X' || *s == '*') {
	  ++s;
	  /* force ND transform of rank 1 */
	  sz.is_nd = 1;
     }
   accept_digit:
     n = 0;

     CHECK(isdigit(*s),
	   "invalid digit");

     while (isdigit(*s)) {
	  n = n * 10 + (*s - '0');
	  ++s;
     }

     sz.narray[sz.rank] = n;
     ++sz.rank;

     CHECK(sz.rank < MAX_CMDLINE_RANK,
	   "maximum rank exceeded");

     if (*s == 'x' || *s == 'X' || *s == '*') {
	  ++s;
	  goto accept_digit;
     }
     /* force ND transform if rank > 1 */
     if (sz.rank > 1)
	  sz.is_nd = 1;
     return sz;
}

/*******************
 * global variables
 *******************/
int verbose = 1;
int only_one_speed_test = 0;
int wisdom_flag = 0, measure_flag = FFTW_ESTIMATE;
int speed_flag = FFTW_MEASURE;
int no_vector_flag = 0;
int chk_mem_leak = 1;
int paranoid = 0;
int howmany_fields = 1;

/* maximum number of iterations to perform in "infinite"
   tests; default (0) means no limit: */
int max_iterations = 0;

/* When testing MPI stuff, only one process gets to do I/O: */
int io_okay = 1;
#define my_printf if (io_okay) printf
#define my_fprintf if (io_okay) fprintf
#define my_fflush if (io_okay) fflush

/*******************
 * procedures
 *******************/

/* smart time printer */
char *smart_sprint_time(double x)
{
     static char buf[128];

     if (x < 1.0E-6)
	  sprintf(buf, "%f ns", x * 1.0E9);
     else if (x < 1.0E-3)
	  sprintf(buf, "%f us", x * 1.0E6);
     else if (x < 1.0)
	  sprintf(buf, "%f ms", x * 1.0E3);
     else
	  sprintf(buf, "%f s", x);

     return buf;
}

/* greet the user */
/* jokes stolen from http://whereis.mit.edu/bin/map */
void please_wait(void)
{
     int i;
     const char *s[] =
     {
       "(while a large software vendor in Seattle takes over the world)",
       "(and remember, this is faster than Java)",
       "(and dream of faster computers)",
       "(checking the gravitational constant in your locale)",
       "(at least you are not on hold)",
       "(while X11 grows by another kilobyte)",
       "(while Windows NT reboots)",
       "(correcting for the phase of the moon)",
       "(your call is important to us)",
       "(while the Linux user-base doubles)",
       "(while you decide where you want to go tomorrow)",
       "(exorcising evil spirits)",
       "(while the C++ standard gains another page)",
     };
     int choices = sizeof(s) / sizeof(*s);

     i = rand() % choices;
     my_printf("Please wait %s.\n", s[i < 0 ? -i : i]);
}

void please_wait_forever(void)
{
     int i;
     const char *s[] =
     {
	  "(but it won't crash, either)",
	  "(at least in theory)",
	  "(please be patient)",
	  "(our next release will complete it more quickly)",
#ifdef HAVE_WIN32
	  "(by the way, Linux executes infinite loops faster)",
#endif
     };
     int choices = sizeof(s) / sizeof(*s);

     if (!max_iterations) {
	  i = rand() % choices;
	  my_printf("This test does not terminate %s.\n", s[i < 0 ? -i : i]);
     } else {
	  my_printf("This test will run for %d iterations.\n", max_iterations);
	  please_wait();
     }
}

/*************************************************
 * Speed tests
 *************************************************/

double mflops(double t, int N)
{
     return (5.0 * N * log((double) N) / (log(2.0) * t * 1.0e6));
}

void print_dims(struct size sz)
{
     int i;

     my_printf("%d", sz.narray[0]);
     for (i = 1; i < sz.rank; ++i)
	  my_printf("x%d", sz.narray[i]);
}

void test_speed(int n)
{
     int specific;

     please_wait();

     if (howmany_fields > 1)
	  WHEN_VERBOSE(1, my_printf("TIMING MULTIPLE-FIELD FFT: "
				 "howmany=%d, stride=%d, dist=%d\n\n",
				 howmany_fields, howmany_fields, 1));

     if (!only_one_speed_test) {
	  for (specific = 0; specific <= 1; ++specific) {
	       WHEN_VERBOSE(1,
 	      my_printf("SPEED TEST: n = %d, FFTW_FORWARD, out of place, %s\n",
			n, SPECIFICP(specific)));
	       test_speed_aux(n, FFTW_FORWARD, 0, specific);

	       WHEN_VERBOSE(1,
	       my_printf("SPEED TEST: n = %d, FFTW_FORWARD, in place, %s\n",
			n, SPECIFICP(specific)));
	       test_speed_aux(n, FFTW_FORWARD, FFTW_IN_PLACE, specific);

	       WHEN_VERBOSE(1,
	     my_printf("SPEED TEST: n = %d, FFTW_BACKWARD, out of place, %s\n",
		       n, SPECIFICP(specific)));
	       test_speed_aux(n, FFTW_BACKWARD, 0, specific);
	       
	       WHEN_VERBOSE(1,
	       my_printf("SPEED TEST: n = %d, FFTW_BACKWARD, in place, %s\n",
			 n, SPECIFICP(specific)));
	       test_speed_aux(n, FFTW_BACKWARD, FFTW_IN_PLACE, specific);
	  }
     }
     else {
	  WHEN_VERBOSE(1,
	     my_printf("SPEED TEST: n = %d, FFTW_FORWARD, out of place, %s\n",
			n, SPECIFICP(1)));
	       test_speed_aux(n, FFTW_FORWARD, 0, 1);
     }
}

void test_speed_nd(struct size sz)
{
     int specific;

     please_wait();

     if (howmany_fields > 1)
	  WHEN_VERBOSE(1, my_printf("TIMING MULTIPLE-FIELD FFT: "
				 "howmany=%d, stride=%d, dist=%d\n\n",
				 howmany_fields, howmany_fields, 1));

     if (!only_one_speed_test) {
	  for (specific = 0; specific <= 1; ++specific) {
	       my_printf("SPEED TEST: ");
	       WHEN_VERBOSE(1, print_dims(sz));
	       WHEN_VERBOSE(1, my_printf(", FFTW_FORWARD, in place, %s\n",
					 SPECIFICP(specific)));
	       test_speed_nd_aux(sz, FFTW_FORWARD,
				 FFTW_IN_PLACE, specific);
	       
	       WHEN_VERBOSE(1, my_printf("SPEED TEST: "));
	       print_dims(sz);
	       WHEN_VERBOSE(1, my_printf(", FFTW_BACKWARD, in place, %s\n",
					 SPECIFICP(specific)));
	       test_speed_nd_aux(sz, FFTW_BACKWARD, FFTW_IN_PLACE, specific);
	  }
     }
     else {
	  my_printf("SPEED TEST: ");
	  WHEN_VERBOSE(1, print_dims(sz));
	  WHEN_VERBOSE(1, my_printf(", FFTW_FORWARD, in place, %s\n",
				    SPECIFICP(1)));
	  test_speed_nd_aux(sz, FFTW_FORWARD,
			    FFTW_IN_PLACE, 1);
     }
}

/*************************************************
 * correctness tests
 *************************************************/

double compute_error_complex(fftw_complex * A, int astride,
			     fftw_complex * B, int bstride, int n)
{
     /* compute the relative error */
     double error = 0.0;
     int i;

     for (i = 0; i < n; ++i) {
	  double a;
	  double mag;
	  a = sqrt(SQR(c_re(A[i * astride]) - c_re(B[i * bstride])) +
		   SQR(c_im(A[i * astride]) - c_im(B[i * bstride])));
	  mag = 0.5 * (sqrt(SQR(c_re(A[i * astride]))
			    + SQR(c_im(A[i * astride]))) +
		       sqrt(SQR(c_re(B[i * bstride]))
			    + SQR(c_im(B[i * bstride])))) + TOLERANCE;

	  a /= mag;
	  if (a > error)
	       error = a;

#ifdef HAVE_ISNAN
	  CHECK(!isnan(a), "NaN in answer");
#endif
     }
     return error;
}

/* test forever */
void test_all(void)
{
     int n;

     please_wait_forever();
     for (n = 1; !max_iterations || n <= max_iterations; ++n) {
	  test_correctness(n);
	  if (!(wisdom_flag & FFTW_USE_WISDOM) && chk_mem_leak)
	       fftw_check_memory_leaks();
     }
}

#define MAX_FACTOR 13

int rand_small_factors(int N)
{
     int f, n = 1;

     f = rand() % MAX_FACTOR + 1;

     while (n * f <= N) {
	  n *= f;
	  f = rand() % MAX_FACTOR + 1;
     }

     return n;
}

#define MAX_N 16384

struct size random_dims(int rank)
{
     int maxsize, dim;
     double maxsize_d;
     struct size sz;

     /* workaround to weird gcc warning */
     maxsize_d = pow((double) (rank == 1 ? MAX_N / 4 : MAX_N),
		     1.0 / (double) rank);
     maxsize = (int) maxsize_d;

     if (maxsize < 1)
	  maxsize = 1;

     sz.rank = rank;
     for (dim = 0; dim < rank; ++dim)
	  sz.narray[dim] = rand_small_factors(maxsize);

     return sz;
}

void test_random(void)
{
     static int counter = 0;
     struct size sz;

     if ((++counter) % 16 == 0) {
	  sz.rank = 1;
	  sz.narray[0] = rand() % (MAX_N / 16) + 1;
     } else {
	  sz = random_dims(1);
     }

     test_correctness(sz.narray[0]);
}

/*************************************************
 * multi-dimensional correctness tests
 *************************************************/

void testnd_correctness_both(struct size sz,
			     int alt_api, int specific, int force_buf)
{
     WHEN_VERBOSE(1,
		  my_printf("Testing nd correctness for size = ");
		  print_dims(sz);
		  my_printf("...");
		  my_fflush(stdout));

     if (alt_api)
	  WHEN_VERBOSE(1, my_printf("alt. api..."));
     if (specific)
	  WHEN_VERBOSE(1, my_printf("specific..."));
     if (force_buf)
	  WHEN_VERBOSE(1, my_printf("force buf..."));

     testnd_correctness(sz, FFTW_FORWARD, alt_api, specific, force_buf);
     testnd_correctness(sz, FFTW_BACKWARD, alt_api, specific, force_buf);

     WHEN_VERBOSE(1, my_printf("OK\n"));
}

void testnd_correctness_aux(struct size sz)
{
     int alt_api, specific, force_buf;

     for (alt_api = 0; alt_api <= 1; ++alt_api)
	  for (specific = 0; specific <= 1; ++specific)
	       for (force_buf = 0; force_buf <= 1; ++force_buf)
		    testnd_correctness_both(sz, alt_api, specific,
					    force_buf);
}

void testnd_correctness_square(int rank, int size)
{
     struct size sz;
     int alt_api, specific, force_buf;
     int i;

     sz.rank = rank;
     for (i = 0; i < rank; ++i)
	  sz.narray[i] = size;

     for (alt_api = 0; alt_api <= 1; ++alt_api)
	  for (specific = 0; specific <= 1; ++specific)
	       for (force_buf = 0; force_buf <= 1; ++force_buf)
		    testnd_correctness_both(sz, alt_api, 
					    specific, force_buf);

}

void testnd_random(int rank)
{
     struct size sz;

     sz = random_dims(rank);
     testnd_correctness_both(sz, coinflip(), coinflip(), coinflip());
}

/* loop forever */
void test_all_random(int rank)
{
     int counter;
     please_wait_forever();

     for (counter = 0; !max_iterations || counter < max_iterations; ++counter) {
	  if (rank > 0)
	       testnd_random(rank);
	  else if ((counter) % 2 == 0)
	       test_random();
	  else
	       testnd_random(rand() % MAX_RANK + 1);
     }

}

int pow2sqrt(int n)
/* return greatest power of two <= sqrt(n) */
{
     int s = 1;

     while (s * s * 4 <= n)
	  s *= 2;
     return s;
}

/* test forever */
void testnd_all(int rank)
{
     int n;

     please_wait_forever();

     for (n = 1; !max_iterations || n <= max_iterations; ++n)
	  testnd_correctness_square(rank, n);
}

fftw_direction random_dir(void)
{
     if (coinflip())
	  return FFTW_FORWARD;
     else
	  return FFTW_BACKWARD;
}

/*************************************************
 * timer tests
 *************************************************/

static int hack_sum_i;

void negative_time(void)
{
     my_fprintf(stderr,
	     "* PROBLEM: I measured a negative time interval.\n"
	     "* Please make sure you defined the timer correctly\n"
	     "* or contact [e-mail] for help.\n");
}