COMP: single-precision compilation with FFT (closes #678)

src/randomProcesses/fft/fft.C

@@ -29,20 +29,20 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace Foam
+Foam::tmp<Foam::complexField>
-{
+Foam::fft::realTransform1D(const scalarField& field)
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-Foam::tmp<Foam::complexField> fft::realTransform1D(const scalarField& field)
 {
     const label n = field.size();
     const label nBy2 = n/2;
 
     // Copy of input field for use by fftw
-    // - fftw requires non-const access to input and output...
+    // - require non-const access to input and output
-    scalar in[n], out[n];
+    // - use double to avoid additional libfftwf for single-precision
-    forAll(field, i)
+
+    List<double> in(n);
+    List<double> out(n);
+
+    for (label i=0; i < n; ++i)
     {
         in[i] = field[i];
     }
@@ -51,8 +51,8 @@ Foam::tmp<Foam::complexField> fft::realTransform1D(const scalarField& field)
     fftw_plan plan = fftw_plan_r2r_1d
     (
         n,
-        in,
+        in.data(),
-        out,
+        out.data(),
         FFTW_R2HC,
         FFTW_ESTIMATE
     );
@@ -77,7 +77,7 @@ Foam::tmp<Foam::complexField> fft::realTransform1D(const scalarField& field)
 }
 
 
-Foam::tmp<Foam::complexField> fft::realTransform1D
+Foam::tmp<Foam::complexField> Foam::fft::realTransform1D
 (
     const tmp<scalarField>& tfield
 )
@@ -88,7 +88,7 @@ Foam::tmp<Foam::complexField> fft::realTransform1D
 }
 
 
-void fft::transform
+void Foam::fft::transform
 (
     complexField& field,
     const UList<int>& nn,
@@ -110,7 +110,7 @@ void fft::transform
     }
 
     // Copy field into fftw containers
-    label N = field.size();
+    const label N = field.size();
     fftw_complex in[N], out[N];
 
     forAll(field, i)
@@ -128,7 +128,7 @@ void fft::transform
     //    fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
 
     // Generic 1..3-D plan
-    label rank = nn.size();
+    const label rank = nn.size();
     fftw_plan plan =
         fftw_plan_dft(rank, nn.begin(), in, out, dir, FFTW_ESTIMATE);
 
@@ -163,7 +163,7 @@ void fft::transform
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-tmp<complexField> fft::forwardTransform
+Foam::tmp<Foam::complexField> Foam::fft::forwardTransform
 (
     const tmp<complexField>& tfield,
     const UList<int>& nn
@@ -179,7 +179,7 @@ tmp<complexField> fft::forwardTransform
 }
 
 
-tmp<complexField> fft::reverseTransform
+Foam::tmp<Foam::complexField> Foam::fft::reverseTransform
 (
     const tmp<complexField>& tfield,
     const UList<int>& nn
@@ -195,7 +195,7 @@ tmp<complexField> fft::reverseTransform
 }
 
 
-tmp<complexVectorField> fft::forwardTransform
+Foam::tmp<Foam::complexVectorField> Foam::fft::forwardTransform
 (
     const tmp<complexVectorField>& tfield,
     const UList<int>& nn
@@ -224,7 +224,7 @@ tmp<complexVectorField> fft::forwardTransform
 }
 
 
-tmp<complexVectorField> fft::reverseTransform
+Foam::tmp<Foam::complexVectorField> Foam::fft::reverseTransform
 (
     const tmp<complexVectorField>& tfield,
     const UList<int>& nn
@@ -253,8 +253,4 @@ tmp<complexVectorField> fft::reverseTransform
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //

src/randomProcesses/noise/noiseFFT/noiseFFT.C

@@ -285,8 +285,8 @@ Foam::tmp<Foam::scalarField> Foam::noiseFFT::Pf
                 << abort(FatalError);
         }
 
-        List<scalar>& in = planInfo_.in;
+        List<double>& in = planInfo_.in;
-        const List<scalar>& out = planInfo_.out;
+        const List<double>& out = planInfo_.out;
         forAll(in, i)
         {
             in[i] = pn[i];
@@ -305,8 +305,8 @@ Foam::tmp<Foam::scalarField> Foam::noiseFFT::Pf
         result[0] = out[0];
         for (label i = 1; i <= nBy2; ++i)
         {
-            scalar re = out[i];
+            const auto re = out[i];
-            scalar im = out[n - i];
+            const auto im = out[n - i];
             result[i] = sqrt(re*re + im*im);
         }
 
@@ -443,7 +443,7 @@ Foam::graph Foam::noiseFFT::PSD(const graph& gPSDf) const
 Foam::graph Foam::noiseFFT::octaves
 (
     const graph& g,
-    const labelList& freqBandIDs,
+    const labelUList& freqBandIDs,
     bool integrate
 ) const
 {

src/randomProcesses/noise/noiseFFT/noiseFFT.H

@@ -68,13 +68,14 @@ class noiseFFT
 :
     public scalarField
 {
-    //- FFTW planner information
+    //- FFTW planner information.
+    //  Storage as double for use directly with FFTW.
     struct planInfo
     {
         bool active;
         label windowSize;
-        scalarList in;
+        List<double> in;
-        scalarList out;
+        List<double> out;
         fftw_plan plan;
     };
 
@@ -174,7 +175,7 @@ public:
         graph octaves
         (
             const graph& g,
-            const labelList& freqBandIDs,
+            const labelUList& freqBandIDs,
             bool integrate
         ) const;