Changeset 19058 for trunk/psLib/src/fft/psImageFFT.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/psImageFFT.c (modified) (16 diffs)

trunk/psLib/src/fft/psImageFFT.c

@@ -6 +6 @@
 /// @author Robert DeSonia, MHPCC
 ///
-/// @version $Revision: 1.26 $ $Name: not supported by cvs2svn $
-/// @date $Date: 2008-04-08 01:25:42 $
+/// @version $Revision: 1.27 $ $Name: not supported by cvs2svn $
+/// @date $Date: 2008-08-14 03:23:13 $
 ///
 /// Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -20 +20 @@
 #include <complex.h>
 #include <fftw3.h>
+#include <pthread.h>
 
 #include "psAssert.h"
@@ -28 +29 @@
 #include "psImageStructManip.h"
 #include "psImageConvolve.h"
+#include "psThread.h"
+#include "psFFT.h"
 #include "psImageFFT.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
@@ -42 +55 @@
     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;
 
     int numCols = in->numCols;          // Number of columns
@@ -64 +81 @@
 
     // Do the FFT
+
+    FFTW_LOCK;
     fftwf_complex *out = fftwf_malloc((numCols/2 + 1) * numRows * sizeof(fftwf_complex)); // Output data
     fftwf_plan plan = fftwf_plan_dft_r2c_2d(numRows, numCols, input, out, FFTW_PLAN_RIGOR);
+    FFTW_UNLOCK;
+
     fftwf_execute(plan);
+
+    FFTW_LOCK;
     fftwf_destroy_plan(plan);
+    FFTW_UNLOCK;
+
     psFree(input);
 
@@ -87 +112 @@
         }
     }
+
+    FFTW_LOCK;
     fftwf_free(out);
+    FFTW_UNLOCK;
 
     return true;
@@ -101 +129 @@
     PS_ASSERT_PTR_NON_NULL(out, false);
 
+    bool threaded = psThreadPoolSize(); // Are we running threaded?
+
     int numCols = real->numCols;        // Number of columns
     int numRows = real->numRows;        // Number of rows
@@ -120 +150 @@
     }
 
-
     // Make sure the system-level wisdom information is imported.
+    FFTW_LOCK;
     if (!fftwWisdomImported) {
         fftwf_import_system_wisdom();
         fftwWisdomImported = true;
     }
+    FFTW_UNLOCK;
 
     // Stuff the real and imaginary parts in
     psF32 *target = psAlloc(2 * numCols * numRows * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); // Target for FFTW
+    FFTW_LOCK;
     fftwf_complex *in = fftwf_malloc(numCols * numRows * sizeof(fftwf_complex)); // Input data
+    FFTW_UNLOCK;
 
     for (int y = 0, index = 0; y < numRows; y++) {
@@ -145 +178 @@
 
     // Do the FFT
+    FFTW_LOCK;
     fftwf_plan plan = fftwf_plan_dft_c2r_2d(numRows, origCols, in, target, FFTW_PLAN_RIGOR);
+    FFTW_UNLOCK;
+
     fftwf_execute(plan);
+
+    FFTW_LOCK;
     fftwf_destroy_plan(plan);
     fftwf_free(in);
+    FFTW_UNLOCK;
 
     // Copy the target pixels into the output
@@ -288 +327 @@
     }
 
+    bool threaded = psThreadPoolSize(); // Are we running threaded?
+
     int numCols = in->numCols, numRows = in->numRows; // Size of image
     int xMin = kernel->xMin, xMax = kernel->xMax, yMin = kernel->yMin, yMax = kernel->yMax; // Kernel sizes
@@ -304 +345 @@
 
     // Create data array containing the padded image and padded kernel
+    FFTW_LOCK;
     psF32 *data = fftwf_malloc(2 * numPadded * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); // Data for FFTW
+    FFTW_UNLOCK;
     psF32 *dataPtr = data;              // Pointer into FFTW data
     psF32 **imageData = in->data.F32;   // Pointer into image data
@@ -317 +360 @@
     memset(dataPtr, 0, (paddedRows - numRows) * paddedCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
 
-#if 1
+#if 0
     {
         // Use this for inspecting the result of copying the image
@@ -379 +422 @@
     }
 
-#if 1
+#if 0
     {
         // Use this for inspecting the result of copying the kernel
@@ -409 +452 @@
     // Do the forward FFT
     // Note that the FFT images have different size from the input
+    FFTW_LOCK;
     fftwf_complex *fft = fftwf_malloc(2 * (paddedCols/2 + 1) * paddedRows * sizeof(fftwf_complex)); // FFT
+    FFTW_UNLOCK;
     int size[] = { paddedRows, paddedCols }; // Size of transforms
     int fftCols = paddedCols/2 + 1, fftRows = paddedRows; // Size of FFT images
     int fftPixels = fftCols * fftRows;  // Number of pixels in FFT image
+
+    FFTW_LOCK;
     fftwf_plan forward = fftwf_plan_many_dft_r2c(2, size, 2, data, NULL, 1, paddedCols * paddedRows,
                                                  fft, NULL, 1, fftPixels, FFTW_PLAN_RIGOR);
+    FFTW_UNLOCK;
+
     fftwf_execute(forward);
+
+    FFTW_LOCK;
     fftwf_destroy_plan(forward);
+    FFTW_UNLOCK;
 
     // Multiply the two transforms
@@ -434 +486 @@
 
     // Do the backward FFT
+    FFTW_LOCK;
     fftwf_plan backward = fftwf_plan_dft_c2r_2d(paddedRows, paddedCols, fft, data, FFTW_PLAN_RIGOR);
+    FFTW_UNLOCK;
+
     fftwf_execute(backward);
+
+    FFTW_LOCK;
     fftwf_destroy_plan(backward);
     fftwf_free(fft);
+    FFTW_UNLOCK;
 
     // Copy into the target, without the padding
@@ -446 +504 @@
         memcpy(*outData, dataPtr, goodBytes);
     }
+
+    FFTW_LOCK;
     fftwf_free(data);
+    FFTW_UNLOCK;
 
     return out;