fft-benchmark/cuda_fft.c

#include <cuda.h>
#include <cuda_runtime.h>
#include <cuda_runtime_api.h>

#include "cuda_fft.h"

void loop_data_cuda (const char *vendor,
	                 CudaBenchmarkFunc func,
	                 int n_devices,
	                 OutputType outputType,
	                 TimeEntry *time_entries)
{
    Timer *timer;

    timer = timer_new ();

    for (int j = 0; j < n_devices; j++) {
        char vendor_name[50];
        int v_len = sprintf(vendor_name, "%s_%d", vendor, j);
        time_entries[j].lib_name = (char *)malloc(sizeof(char) * (v_len + 1));
        strcpy(time_entries[j].lib_name, vendor_name);

        time_entries[j].dim_entries = (DimEntry *)malloc(N_DIMS * sizeof(DimEntry));

        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];
            int num_entries = power_max - power_min + 1;

            time_entries[j].dim_entries[k].n_dims = dim;
            time_entries[j].dim_entries[k].sizes  = (unsigned int **)malloc(sizeof(unsigned int *) * num_entries);
            time_entries[j].dim_entries[k].times  = (double *)malloc(sizeof(double) * num_entries);
            time_entries[j].dim_entries[k].errors = (double *)malloc(sizeof(double) * num_entries);

            PRINT_DIM (dim);
            fflush (stdout);

            for (int m = power_min, i = 0; m <= power_max; m++, i++) {
                size_t size_bytes;
                cufftComplex *host_orig_mem;
                cufftComplex *host_result_mem;
                cufftComplex *dev_mem;
                cufftComplex *dev_out_mem;

                size_t side_size = pow(2,m);
                size_t size = pow(side_size,dim);

                size_bytes = size * sizeof (cufftComplex);
                host_orig_mem = (cufftComplex *)malloc(size_bytes);
                host_result_mem = (cufftComplex *)malloc(size_bytes);

                for (int l = 0; l < size; l++) {
                    host_orig_mem[l].x = rand() / ((float) RAND_MAX);
                    host_orig_mem[l].y = rand() / ((float) RAND_MAX);
                }

                CUDA_SAFE_CALL (cudaMalloc ((void **)&dev_mem, size_bytes));
                CUDA_SAFE_CALL (cudaMalloc ((void **)&dev_out_mem, size_bytes));

                PRINT_DIMS(dim, side_size);

                fflush (stdout);

                double time_sec;
                double sum;
                bool scale;

                printf (".");
                fflush (stdout);

                CUDA_SAFE_CALL (cudaMemcpy (dev_mem, host_orig_mem, size_bytes, cudaMemcpyHostToDevice));

                size_t fft_size[3] = { 1, 1, 1};
                for (int l = 0; l < dim; l++) {
                    fft_size[l] = side_size;
                }

                scale = func (dev_mem, dev_out_mem, dim, fft_size, N_RUNS, timer);

                /* Check precision */
                CUDA_SAFE_CALL (cudaMemcpy (host_result_mem, dev_out_mem, size_bytes, cudaMemcpyDeviceToHost));
                sum = sum_of_absolute_differences_complex (host_orig_mem, host_result_mem, size, scale);
        
                time_sec = timer_get_seconds (timer) / N_RUNS;

                time_entries[j].dim_entries[k].times[i] = get_measurements_with_format(outputType, size_bytes, time_sec);
                time_entries[j].dim_entries[k].errors[i] = sum / size;

                free (host_orig_mem);
                free (host_result_mem);

                CUDA_SAFE_CALL (cudaFree (dev_mem));
                CUDA_SAFE_CALL (cudaFree (dev_out_mem));
            }

            printf ("\n");
            fflush (stdout);
        }
    }

    printf ("\n");
    timer_destroy (timer);
}

bool compute_cuda_fft (cufftComplex *dev_mem,
	                   cufftComplex *out_mem,
	                   int n_dims,
	                   size_t *dims,
	                   int n_runs,
	                   Timer *timer)
{
    cufftHandle plan;

    int dim_sizes[3] = {1, 1, 1};

    switch (n_dims) {
        case 1:
            dim_sizes[0] = dims[0];
            break;
        case 2:
            dim_sizes[0] = dims[0];
            dim_sizes[1] = dims[1];
            break;
        case 3:
            dim_sizes[0] = dims[0];
            dim_sizes[1] = dims[1];
            dim_sizes[2] = dims[2];
            break;
        default:
            fprintf (stderr, "Unknown FFT dimensions\n");
            return true;
    }

    CUFFT_SAFE_CALL (cufftPlanMany(&plan, n_dims, dim_sizes, NULL, 1, 0, NULL, 1, 0, CUFFT_C2C, 1));

    timer_start (timer);

    for (int i = 0; i < n_runs; i++) {
        CUFFT_SAFE_CALL (cufftExecC2C (plan, (cufftComplex *)dev_mem, (cufftComplex *)out_mem, CUFFT_FORWARD));
        CUDA_SAFE_CALL (cudaDeviceSynchronize ());
    }

    timer_stop (timer);

    CUFFT_SAFE_CALL (cufftExecC2C (plan, (cufftComplex *)out_mem, (cufftComplex *)dev_mem, CUFFT_INVERSE));

    CUDA_SAFE_CALL (cudaDeviceSynchronize ());

    CUDA_SAFE_CALL (cudaMemcpy (out_mem, dev_mem, dim_sizes[0] * dim_sizes[1] * dim_sizes[2] * sizeof(cufftComplex), cudaMemcpyDeviceToDevice));   

    CUFFT_SAFE_CALL (cufftDestroy (plan));

    return true;
}

double sum_of_absolute_differences_complex (cufftComplex *a, cufftComplex *b, int n, bool scale)
{
    double sum = 0.0;

    for (int i = 0; i < n; i++) {
        sum += fabs (a[i].x - b[i].x / ((float)n));
        sum += fabs (a[i].y - b[i].y / ((float)n));
    }

    return sum;
}