Changeset 19058 for trunk/psLib/src/fft/psVectorFFT.c

Timestamp:

Aug 13, 2008, 5:23:13 PM (18 years ago)

Author:

Paul Price

Message:

Learned that of all the fftw_ calls, only fftw_execute is thread-safe: all other calls have to be wrapped in a mutex for thread safety.

File:

• trunk/psLib/src/fft/psVectorFFT.c (modified) (6 diffs)

trunk/psLib/src/fft/psVectorFFT.c

@@ -6 +6 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.39 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-02-09 00:45:51 $
+ *  @version $Revision: 1.40 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-08-14 03:23:13 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -27 +27 @@
 #include "psLogMsg.h"
 #include "psConstants.h"
+#include "psThread.h"
+#include "psFFT.h"
 #include "psVectorFFT.h"
+
+// Lock FFTW access
+#define FFTW_LOCK \
+if (threaded) { \
+    psFFTLock(); \
+}
+// Unlock FFTW access
+#define FFTW_UNLOCK \
+if (threaded) { \
+    psFFTUnlock(); \
+}
 
 #define FFTW_PLAN_RIGOR FFTW_ESTIMATE   // How rigorous the FFTW planning is
@@ -40 +53 @@
     PS_ASSERT_PTR_NON_NULL(imag, false);
 
+    bool threaded = psThreadPoolSize(); // Are we running threaded?
+
     // Make sure the system-level wisdom information is imported.
+    FFTW_LOCK;
     if (!fftwWisdomImported) {
         fftwf_import_system_wisdom();
         fftwWisdomImported = true;
     }
+    FFTW_UNLOCK;
 
     long num = in->n;                   // Number of elements
 
     // Do the FFT
+    FFTW_LOCK;
     fftwf_complex *out = fftwf_malloc((num/2 + 1) * sizeof(fftwf_complex)); // Output data
     fftwf_plan plan = fftwf_plan_dft_r2c_1d(num, in->data.F32, out, FFTW_PLAN_RIGOR);
+    FFTW_UNLOCK;
+
     fftwf_execute(plan);
+
+    FFTW_LOCK;
     fftwf_destroy_plan(plan);
+    FFTW_UNLOCK;
 
     // Pull the real and imaginary parts out
@@ -69 +92 @@
 #endif
     }
+
+    FFTW_LOCK;
     fftwf_free(out);
+    FFTW_UNLOCK;
 
     return true;
@@ -83 +109 @@
     PS_ASSERT_PTR_NON_NULL(out, false);
 
+    bool threaded = psThreadPoolSize(); // Are we running threaded?
+
     // Make sure the system-level wisdom information is imported.
+    FFTW_LOCK;
     if (!fftwWisdomImported) {
         fftwf_import_system_wisdom();
         fftwWisdomImported = true;
     }
+    FFTW_UNLOCK;
 
     long num = real->n;                 // Number of elements
 
     // Stuff the real and imaginary parts in
+    FFTW_LOCK;
     fftwf_complex *in = fftwf_malloc(num * sizeof(fftwf_complex)); // Input data
+    FFTW_UNLOCK;
     for (int i = 0; i < num; i++) {
 #if !defined(FFTW_NO_Complex) && defined(_Complex_I) && defined(complex) && defined(I)
@@ -107 +139 @@
     // Do the FFT
     *out = psVectorRecycle(*out, origNum, PS_TYPE_F32);
+    FFTW_LOCK;
     fftwf_plan plan = fftwf_plan_dft_c2r_1d(origNum, in, (*out)->data.F32, FFTW_PLAN_RIGOR);
+    FFTW_UNLOCK;
+
     fftwf_execute(plan);
+
+    FFTW_LOCK;
     fftwf_destroy_plan(plan);
-
     fftwf_free(in);
+    FFTW_UNLOCK;
 
     return true;