Try to partly integrate multi-dim support into FFTW

benchmark.c

@@ -55,13 +55,19 @@ exit(EXIT_FAILURE); \
 
 #include "timer.h"
 
+#define N_DIMS 3
+#define N_RUNS 4
+#define INITIAL_SIZE 8
+
 const int N_ARRAYS = 6;
-const int N_DIMS = 3;
-const int N_RUNS = 4;
 
-const int INITIAL_SIZE = 8;
+const int DIMS[N_DIMS] = {1, 2, 3};
+const int N_DIM_ARRAYS[N_DIMS] = {5, 4, 4};
+const int N_POWERS_INTERVALS[N_DIMS][2] = {{5, 11}, {8, 12}, {7, 8}};
 
 #define UPDATE_SIZE(size) size *= 8;
+#define PRINT_DIM_SIZE(side_size,dim) { \
+    printf(" %zu", side_size); while (dim != 1) { printf("x%zu", side_size);dim--; } printf("."); }
 
 #ifdef HAVE_OPENCL
 #define OCL_CHECK_ERROR(error) { \
@@ -261,11 +267,107 @@ static void
 loop_data_fftw (double *times, double *errors, FILE *fp)
 {
     Timer *timer;
-    
-    size_t size = INITIAL_SIZE;
+
     timer = timer_new ();
     fprintf (fp, "FFTW_bw FFTW_err ");
 
+    for (int k = 0; k < N_DIMS; k++) {
+        int dim = DIMS[k];
+        int power_min = N_POWERS_INTERVALS[k][0];
+        int power_max = N_POWERS_INTERVALS[k][1];
+
+        printf ("%dD:", dim);
+        fflush (stdout);
+
+        for (int m = power_min; m <= power_max; m++) {
+            fftw_complex *host_orig_mem;
+            fftw_complex *host_result_mem;
+            fftw_complex *host_immediate_mem;
+            fftw_plan plan;
+            fftw_plan inverse_plan;
+            //double time;
+            //double mflops;
+            double sum = 0.0;
+
+            size_t side_size = pow(2,m);
+            size_t size = pow(side_size,dim);
+
+            host_orig_mem = fftw_malloc (sizeof (fftw_complex) * size);
+            host_immediate_mem = fftw_malloc (sizeof (fftw_complex) * size);
+            host_result_mem = fftw_malloc (sizeof (fftw_complex) * size);
+
+            switch (dim) {
+                case 1:
+                    plan = fftw_plan_dft_1d (side_size, 
+                                             host_orig_mem, host_immediate_mem,
+                                             FFTW_FORWARD, FFTW_ESTIMATE);
+                    inverse_plan = fftw_plan_dft_1d (side_size, 
+                                                     host_immediate_mem,
+                                                     host_result_mem,
+                                                     FFTW_BACKWARD, FFTW_ESTIMATE);
+                    break;
+                case 2:
+                    plan = fftw_plan_dft_2d (side_size, side_size,
+                                             host_orig_mem, host_immediate_mem,
+                                             FFTW_FORWARD, FFTW_ESTIMATE);
+                    inverse_plan = fftw_plan_dft_2d (side_size, side_size,
+                                                     host_immediate_mem, host_result_mem,
+                                                     FFTW_BACKWARD, FFTW_ESTIMATE);
+                    break;
+                case 3:
+                    plan = fftw_plan_dft_3d (side_size, side_size, side_size,
+                                             host_orig_mem, host_immediate_mem,
+                                             FFTW_FORWARD, FFTW_ESTIMATE);
+                    inverse_plan = fftw_plan_dft_3d (side_size, side_size, side_size,
+                                                     host_immediate_mem, host_result_mem,
+                                                     FFTW_BACKWARD, FFTW_ESTIMATE);
+                    break;
+                default:
+                    fprintf (stderr, "Unknown FFT dimensions\n");
+                    return;
+            }
+
+            for (int j = 0; j < size; j++) {
+                host_orig_mem[j][0] = rand() / ((double) RAND_MAX);
+                host_orig_mem[j][1] = rand() / ((double) RAND_MAX);
+            }
+
+            if (dim == 1) printf (" %zu.", side_size);
+            else if (dim == 2) printf (" %zux%zu.", side_size, side_size);
+            else printf (" %zux%zux%zu.", side_size, side_size, side_size);
+            fflush (stdout);
+
+            timer_start (timer);
+
+            for (int j = 0; j < N_RUNS; j++) {
+                fftw_execute (plan);
+            }
+
+            timer_stop (timer);
+
+            /* Check precision */
+            fftw_execute (inverse_plan);
+
+            for (int j = 0; j < size; j++) {
+                sum += fabs (host_result_mem[j][0] / size - host_orig_mem[j][0]);
+                sum += fabs (host_result_mem[j][1] / size - host_orig_mem[j][1]);
+            }
+
+            //time = timer_get_seconds (timer) / N_RUNS / 1000.0;
+            //mflops = 5 * size * log (size) / log (2) / time;
+            //times[i] = mflops;
+            //errors[i] = sum / size;
+
+            fftw_destroy_plan (inverse_plan);
+            fftw_destroy_plan (plan);
+            fftw_free (host_orig_mem);
+            fftw_free (host_immediate_mem);
+            fftw_free (host_result_mem);
+        }
+        printf ("\n");
+        fflush (stdout);
+    }
+    /*
     for (int i = 0; i < N_ARRAYS; i++) {
             fftw_complex *host_orig_mem;
             fftw_complex *host_result_mem;
@@ -300,7 +402,6 @@ loop_data_fftw (double *times, double *errors, FILE *fp)
 
             timer_stop (timer);
 
-            /* Check precision */
             inverse_plan = fftw_plan_dft_1d (size, host_immediate_mem, host_result_mem, FFTW_BACKWARD, FFTW_ESTIMATE);
             fftw_execute (inverse_plan);
 
@@ -320,6 +421,7 @@ loop_data_fftw (double *times, double *errors, FILE *fp)
             fftw_free (host_immediate_mem);
             fftw_free (host_result_mem);
     }
+    */
     printf ("\n");
     timer_destroy (timer);
 }