gpu_fftw.c.template

/* vi: set ft=c: */
#include "gpu_fftw.h"
#include "hello_fft/gpu_fft.h"
#include "hello_fft/mailbox.h"
#include <stdlib.h>
#include <fftw3.h>
#include <assert.h>
#include <string.h>

#ifndef S_SPLINT_S
#include <syslog.h>
#endif

/**********************
 *  Copy data to GPU  *
 **********************/

void SO_LOCAL copy_input(int N,PREFIX_complex* in,struct GPU_FFT* fft)
{
   GPU_FFT_COMPLEX* base;
   int i,j;
   int jobs=1;
   if (in==NULL) {
      say(LOG_ERR,"Received NULL pointer for input\n");
      return;
   }
   if (!strcmp("PREFIX","fftwf")) {
   for (j=0; j<jobs; j++) {
      base = fft->in + j*fft->step; // We can copy the memory, same type
      memcpy((void*) base, (void*) in, N*2*sizeof(float));
   }
   } else {
   for (j=0; j<jobs; j++) {
      base = fft->in + j*fft->step; // We need to convert float to bouble
      for (i=0; i<N; i++) {
         base[i][0] = (float) in[i][0];
         base[i][1] = (float) in[i][1];
      }
   }
   }
}

void SO_LOCAL copy_output(int N,PREFIX_complex* restrict out,struct GPU_FFT* restrict fft)
{
   GPU_FFT_COMPLEX* restrict base;
   int i,j;
   int jobs=1;
   if (out==NULL) {
      say(LOG_ERR,"Received NULL pointer for ouput\n");
      return;
   }
   if (!strcmp("PREFIX","fftwf")) {
      //Both arrays are float, can use memcpy
      // 2-3 times faster
      for (j=0; j<jobs; j++) {
         base = fft->out + j*fft->step; // output buffer
         memcpy((void*) out, (void*) base, N*2*sizeof(float));
      }
   } else {
      //Arrays of different types, have to copy one by one
      for (j=0; j<jobs; j++) {
         base = fft->out + j*fft->step; // output buffer
         for (i=0; i<N; i++) {
            out[i][0] = base[i][0];
            out[i][1] = base[i][1];
         }
   }

   }
}

/**********
 *  Plan  *
 **********/
typedef PREFIX_plan(*ptr_PREFIX_plan_dft_1d)
   (int n, PREFIX_complex *in, PREFIX_complex *out, int sign, unsigned flags);

/* PREFIX_plan here is really a pointer to GPU_FFT */
SO_LOCAL struct GPU_FFT* gpu_fftw_plan(int N,
      unsigned int log2_N,
      PREFIX_complex *in,
      PREFIX_complex *out,
      int sign,
      unsigned int flags)
{
   static int mb = -1;

   struct GPU_FFT* fft = NULL;
   if (mb==-1)
      mb=mbox_open();

   if (mb < 0)
      return NULL;

   int jobs = 1;
   int direction = (sign == FFTW_FORWARD) ? GPU_FFT_FWD:GPU_FFT_REV;
   int ret = gpu_fft_prepare(mb, log2_N, direction, jobs, &fft);

    switch(ret) {
        case -1:
           say(LOG_ERR,"Unable to enable V3D. Please check your firmware is up to date.\n");
           return NULL;
        case -2:
           say(LOG_ERR,"log2_N=%d not supported. This is a gpu_fftw bug. Please report it.\n", log2_N);
           return NULL;
        case -3:
           say(LOG_ERR,"GPU is out of memory.  Try a smaller batch or increase GPU memory.\n");
           return NULL;
        case -4:
           say(LOG_ERR,"Hardware error. Unable to map Videocore peripherals into ARM memory space.\n");
           return NULL;
    }
    return fft;
}

PREFIX_plan SO_EXPORT PREFIX_plan_dft_1d(int n, PREFIX_complex *in, PREFIX_complex *out, int sign, unsigned flags)
{
   static unsigned char fprint[FINGERPRINTSZ] = FINGERPRINT;
   static ptr_PREFIX_plan_dft_1d fftw_plan_fun=NULL;
   static struct GPU_FFTW_PLAN mplan_arr[MPLAN_ARRSIZE];
   struct GPU_FFTW_PLAN* metaplan;
   static int mplan_n = 0;
   static char* env=NULL;

   /* This is triggered by the destructor function */
   if (n==-1) { 
      /* Clean unreleased plans */
      for (int i = mplan_n - 1; i >= 0; --i)
         PREFIX_destroy_plan((PREFIX_plan) &mplan_arr[i]);
      return NULL;
   }

   if (fftw_plan_fun==NULL)
      fftw_plan_fun=orig_func("PREFIX_plan_dft_1d",PREFIX_plan_dft_1d);

   metaplan = NULL;
   if (mplan_n<MPLAN_ARRSIZE)
      metaplan = &mplan_arr[mplan_n++];
   else {
      say(LOG_ERR,"Sorry, ran out of space for plans: max %d plans",MPLAN_ARRSIZE);
      goto use_fftw3;
   }

   memcpy((void*)metaplan->fprint,(void*)fprint,FINGERPRINTSZ);
   metaplan->count_ptr = &mplan_n;
   metaplan->n=n;
   metaplan->out=(void*) out;
   metaplan->gpu=false;
   env = getenv("GPU_FFTW_DISABLE");
   if (env!=NULL)
      goto use_fftw3;

   /* Use FFTW if:
    *     - n is not a power of 2 */
   unsigned int u = (unsigned int) n;
   if ( u >=256 && IS_2_PWR(u) ) {
      u = log2u(u);
      metaplan->plan=(void*)gpu_fftw_plan(n,u,in,out,sign,flags);
      metaplan->gpu=true;
      if (metaplan->plan==NULL)
         goto use_fftw3;
      else {
         copy_input(n,in,(struct GPU_FFT*) metaplan->plan);
         gpu_active(true);
         say(LOG_INFO,"plan #%d - gpu fft plan running (N=%d,%s)\n",mplan_n,n,
               (sign==FFTW_FORWARD) ? "forward":"backward");
         return (PREFIX_plan) metaplan;
      }
   } else {

use_fftw3:
      say(LOG_INFO,"%s","Falling back to fftw3\n");
      PREFIX_plan plan = (*fftw_plan_fun)(n,in,out,sign,flags);
      metaplan->gpu=false;
      metaplan->plan=(void*)plan;
      gpu_active(false);
      return (PREFIX_plan) metaplan;
   }
}

/*************
 *  Execute  *
 *************/
typedef void(*ptr_PREFIX_execute)(const PREFIX_plan);

void SO_EXPORT PREFIX_execute(const PREFIX_plan plan)
{
   static ptr_PREFIX_execute fftw_execute_fun=NULL;
   struct GPU_FFTW_PLAN* metaplan = (struct GPU_FFTW_PLAN*) plan;

   if (fftw_execute_fun==NULL)
      fftw_execute_fun=orig_func("PREFIX_execute",PREFIX_execute);

   if (!fingerprint_ok((void*)metaplan)) {
      say(LOG_ERR,"Falling back to fftw3 - fingerprint mismatch in execute\n");
      (*fftw_execute_fun)(plan);
      return;
   }

   if ( metaplan->gpu ) {
      gpu_fft_execute( (struct GPU_FFT*) metaplan->plan);
      copy_output(metaplan->n,
            (PREFIX_complex*) metaplan->out,
            (struct GPU_FFT*) metaplan->plan);
   }
   else
      (*fftw_execute_fun)((PREFIX_plan) metaplan->plan);
}

/*************
 *  Destroy  *
 *************/
typedef void(*ptr_PREFIX_destroy_plan)(PREFIX_plan);

void SO_EXPORT PREFIX_destroy_plan(PREFIX_plan plan)
{
   static ptr_PREFIX_destroy_plan fftw_destroy_plan_fun=NULL;
   struct GPU_FFTW_PLAN* metaplan = (struct GPU_FFTW_PLAN*) plan;

   if (fftw_destroy_plan_fun==NULL)
      fftw_destroy_plan_fun=orig_func("PREFIX_destroy_plan",PREFIX_destroy_plan);

   if (!fingerprint_ok((void*)metaplan)) {
      say(LOG_ERR,"Falling back to fftw3 - fingerprint mismatch in destroy_plan\n");
      (*fftw_destroy_plan_fun)(plan);
      return;
   }

   (*metaplan->count_ptr)--;

   if (metaplan->gpu) {
      say(LOG_INFO,"plan #%d - released (N=%d)\n",*metaplan->count_ptr + 1,metaplan->n);
      gpu_fft_release( (struct GPU_FFT*) metaplan->plan );
   }
   else
      (*fftw_destroy_plan_fun)((PREFIX_plan) metaplan->plan);
}

void SO_LOCAL __attribute__ ((destructor)) clean_undestroyed_plans()
{
   PREFIX_plan_dft_1d(-1, NULL, NULL, 0, 0);
}

/*==========================================================================*
 *                      FORTRAN API
 *==========================================================================*
 * SUFFIX: fortran functions have a _ or a __ at the end
 *         depending on the compiler, so we provide both versions
 * PARAMS: function parameters in fortran are always a pointer, and
 *         from what I see have no return value */

/*****************
 *  Execute      *
 *****************/

/* This is done by looking at the source code in:
 * https://github.com/FFTW/fftw3/blob/cde4559ba9b822166cb88a84a0994fdb83a2061c/api/f77funcs.h
 *
 * We see for example, that the execute function does more than simply call the C
 * equivalent, so we need to override that function;

 * This is a partial implementation, functions are missing, but the basic 1d
 * plan, execute and destroy should work, b/c plan and destroy simply
 * call the equivalent C function, they don't need to be overriden  */
void SO_EXPORT FPREFIXfftw_execute_(const PREFIX_plan* plan)
{
   PREFIX_execute(*plan);
}