Browse Source

FFT: experimental implementation with test for complex-complex inplace forward FFT

benchmark-la
Christoph Ruegg 10 years ago
parent
commit
8639df292f
  1. 4
      src/NativeProviders/Linux/mkl_build.sh
  2. 33
      src/NativeProviders/MKL/fft.cpp
  3. 4
      src/NativeProviders/OSX/mkl_build.sh
  4. 1
      src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj
  5. 3
      src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj.filters
  6. 1
      src/Numerics/Numerics.csproj
  7. 2
      src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs
  8. 38
      src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs
  9. 38
      src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs
  10. 93
      src/Numerics/Providers/FourierTransform/Mkl/SafeNativeMethods.cs
  11. 81
      src/UnitTests/FourierTransformProviderTests/FourierTransformProviderTests.cs

4
src/NativeProviders/Linux/mkl_build.sh

@ -7,10 +7,10 @@ export OUT=../../../out/MKL/Linux
mkdir -p $OUT/x64
mkdir -p $OUT/x86
g++ -std=c++11 -D_M_X64 -DGCC -m64 --shared -fPIC -o $OUT/x64/MathNet.Numerics.MKL.dll -I$MKL/include -I../Common -I../MKL ../MKL/memory.c ../MKL/capabilities.cpp ../MKL/vector_functions.c ../Common/blas.c ../Common/lapack.cpp -Wl,--start-group $MKL/lib/intel64/libmkl_intel_lp64.a $MKL/lib/intel64/libmkl_intel_thread.a $MKL/lib/intel64/libmkl_core.a -Wl,--end-group -L$OPENMP/intel64 -liomp5 -lpthread -lm
g++ -std=c++11 -D_M_X64 -DGCC -m64 --shared -fPIC -o $OUT/x64/MathNet.Numerics.MKL.dll -I$MKL/include -I../Common -I../MKL ../MKL/memory.c ../MKL/capabilities.cpp ../MKL/vector_functions.c ../Common/blas.c ../Common/lapack.cpp ../MKL/fft.cpp -Wl,--start-group $MKL/lib/intel64/libmkl_intel_lp64.a $MKL/lib/intel64/libmkl_intel_thread.a $MKL/lib/intel64/libmkl_core.a -Wl,--end-group -L$OPENMP/intel64 -liomp5 -lpthread -lm
cp $OPENMP/intel64/libiomp5.so $OUT/x64/
g++ -std=c++11 -D_M_IX86 -DGCC -m32 --shared -fPIC -o $OUT/x86/MathNet.Numerics.MKL.dll -I$MKL/include -I../Common -I../MKL ../MKL/memory.c ../MKL/capabilities.cpp ../MKL/vector_functions.c ../Common/blas.c ../Common/lapack.cpp -Wl,--start-group $MKL/lib/ia32/libmkl_intel.a $MKL/lib/ia32/libmkl_intel_thread.a $MKL/lib/ia32/libmkl_core.a -Wl,--end-group -L$OPENMP/ia32 -liomp5 -lpthread -lm
g++ -std=c++11 -D_M_IX86 -DGCC -m32 --shared -fPIC -o $OUT/x86/MathNet.Numerics.MKL.dll -I$MKL/include -I../Common -I../MKL ../MKL/memory.c ../MKL/capabilities.cpp ../MKL/vector_functions.c ../Common/blas.c ../Common/lapack.cpp ../MKL/fft.cpp -Wl,--start-group $MKL/lib/ia32/libmkl_intel.a $MKL/lib/ia32/libmkl_intel_thread.a $MKL/lib/ia32/libmkl_core.a -Wl,--end-group -L$OPENMP/ia32 -liomp5 -lpthread -lm
cp $OPENMP/ia32/libiomp5.so $OUT/x86/

33
src/NativeProviders/MKL/fft.cpp

@ -0,0 +1,33 @@
#include "wrapper_common.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include "mkl_service.h"
#include "mkl_dfti.h"
extern "C" {
DLLEXPORT MKL_LONG z_fft_forward_inplace(MKL_LONG n, MKL_Complex16 x[])
{
MKL_LONG status = 0;
DFTI_DESCRIPTOR_HANDLE hand = 0;
status = DftiCreateDescriptor(&hand, DFTI_DOUBLE, DFTI_COMPLEX, 1, n);
if (0 != status) goto failed;
status = DftiCommitDescriptor(hand);
if (0 != status) goto failed;
status = DftiComputeForward(hand, x);
if (0 != status) goto failed;
cleanup:
DftiFreeDescriptor(&hand);
return status;
failed:
status = 1;
goto cleanup;
}
}

4
src/NativeProviders/OSX/mkl_build.sh

@ -6,10 +6,10 @@ export OUT=../../../out/MKL/OSX
mkdir -p $OUT/x64
mkdir -p $OUT/x86
clang++ -std=c++11 -D_M_X64 -DGCC -m64 --shared -fPIC -o $OUT/x64/MathNet.Numerics.MKL.dll -I$MKL/include -I../Common -I../MKL ../MKL/memory.c ../MKL/capabilities.cpp ../MKL/vector_functions.c ../Common/blas.c ../Common/lapack.cpp $MKL/lib/libmkl_intel_lp64.a $MKL/lib/libmkl_core.a $MKL/lib/libmkl_intel_thread.a -L$OPENMP -liomp5 -lpthread -lm
clang++ -std=c++11 -D_M_X64 -DGCC -m64 --shared -fPIC -o $OUT/x64/MathNet.Numerics.MKL.dll -I$MKL/include -I../Common -I../MKL ../MKL/memory.c ../MKL/capabilities.cpp ../MKL/vector_functions.c ../Common/blas.c ../Common/lapack.cpp ../MKL/fft.cpp $MKL/lib/libmkl_intel_lp64.a $MKL/lib/libmkl_core.a $MKL/lib/libmkl_intel_thread.a -L$OPENMP -liomp5 -lpthread -lm
cp $OPENMP/libiomp5.dylib $OUT/x64/
clang++ -std=c++11 -D_M_IX86 -DGCC -m32 --shared -fPIC -o $OUT/x86/MathNet.Numerics.MKL.dll -I$MKL/include -I../Common -I../MKL ../MKL/memory.c ../MKL/capabilities.cpp ../MKL/vector_functions.c ../Common/blas.c ../Common/lapack.cpp $MKL/lib/libmkl_intel_lp64.a $MKL/lib/libmkl_core.a $MKL/lib/libmkl_intel_thread.a -L$OPENMP -liomp5 -lpthread -lm
clang++ -std=c++11 -D_M_IX86 -DGCC -m32 --shared -fPIC -o $OUT/x86/MathNet.Numerics.MKL.dll -I$MKL/include -I../Common -I../MKL ../MKL/memory.c ../MKL/capabilities.cpp ../MKL/vector_functions.c ../Common/blas.c ../Common/lapack.cpp ../MKL/fft.cpp $MKL/lib/libmkl_intel_lp64.a $MKL/lib/libmkl_core.a $MKL/lib/libmkl_intel_thread.a -L$OPENMP -liomp5 -lpthread -lm
cp $OPENMP/libiomp5.dylib $OUT/x86/

1
src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj

@ -294,6 +294,7 @@
<ClCompile Include="..\..\MKL\capabilities.cpp" />
<ClCompile Include="..\..\MKL\memory.c" />
<ClCompile Include="..\..\MKL\vector_functions.c" />
<ClCompile Include="..\..\MKL\fft.cpp" />
</ItemGroup>
<ItemGroup>
<ResourceCompile Include="..\..\MKL\resource.rc" />

3
src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj.filters

@ -33,6 +33,9 @@
<ClCompile Include="..\..\Common\blas.c">
<Filter>Source Files</Filter>
</ClCompile>
<ClCompile Include="fft.cpp">
<Filter>Source Files</Filter>
</ClCompile>
</ItemGroup>
<ItemGroup>
<ResourceCompile Include="..\..\MKL\resource.rc">

1
src/Numerics/Numerics.csproj

@ -175,6 +175,7 @@
<Compile Include="Providers\FourierTransform\IFourierTransformProvider.cs" />
<Compile Include="Providers\FourierTransform\ManagedFourierTransformProvider.cs" />
<Compile Include="Providers\FourierTransform\Mkl\MklFourierTransformProvider.cs" />
<Compile Include="Providers\FourierTransform\Mkl\SafeNativeMethods.cs" />
<Compile Include="Providers\LinearAlgebra\Acml\AcmlLinearAlgebraProvider.Complex.cs" />
<Compile Include="Providers\LinearAlgebra\Acml\AcmlLinearAlgebraProvider.Complex32.cs" />
<Compile Include="Providers\LinearAlgebra\Acml\AcmlLinearAlgebraProvider.Double.cs" />

2
src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs

@ -1,4 +1,4 @@
// <copyright file="ILinearAlgebraProvider.cs" company="Math.NET">
// <copyright file="IFourierTransformProvider.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics

38
src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs

@ -1,4 +1,32 @@
using System.Numerics;
// <copyright file="ManagedFourierTransformProvider.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://mathnet.opensourcedotnet.info
//
// Copyright (c) 2009-2016 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
using System.Numerics;
using MathNet.Numerics.IntegralTransforms;
namespace MathNet.Numerics.Providers.FourierTransform
@ -9,17 +37,17 @@ namespace MathNet.Numerics.Providers.FourierTransform
{
}
public void ForwardInplace(Complex[] complex)
public virtual void ForwardInplace(Complex[] complex)
{
Fourier.BluesteinForward(complex, FourierOptions.Default);
}
public void BackwardInplace(Complex[] complex)
public virtual void BackwardInplace(Complex[] complex)
{
Fourier.BluesteinInverse(complex, FourierOptions.Default);
}
public Complex[] Forward(Complex[] complexTimeSpace)
public virtual Complex[] Forward(Complex[] complexTimeSpace)
{
Complex[] work = new Complex[complexTimeSpace.Length];
complexTimeSpace.Copy(work);
@ -27,7 +55,7 @@ namespace MathNet.Numerics.Providers.FourierTransform
return work;
}
public Complex[] Backward(Complex[] complexFrequenceSpace)
public virtual Complex[] Backward(Complex[] complexFrequenceSpace)
{
Complex[] work = new Complex[complexFrequenceSpace.Length];
complexFrequenceSpace.Copy(work);

38
src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs

@ -1,7 +1,32 @@
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
// <copyright file="MklFourierTransformProvider.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://mathnet.opensourcedotnet.info
//
// Copyright (c) 2009-2016 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
using System.Numerics;
namespace MathNet.Numerics.Providers.FourierTransform.Mkl
{
@ -10,5 +35,10 @@ namespace MathNet.Numerics.Providers.FourierTransform.Mkl
public override void InitializeVerify()
{
}
public override void ForwardInplace(Complex[] complex)
{
SafeNativeMethods.z_fft_forward_inplace(complex.Length, complex);
}
}
}

93
src/Numerics/Providers/FourierTransform/Mkl/SafeNativeMethods.cs

@ -0,0 +1,93 @@
// <copyright file="SafeNativeMethods.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://mathnet.opensourcedotnet.info
//
// Copyright (c) 2009-2016 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
#if NATIVE
using System.Numerics;
using System.Runtime.InteropServices;
using System.Security;
namespace MathNet.Numerics.Providers.FourierTransform.Mkl
{
/// <summary>
/// P/Invoke methods to the native math libraries.
/// </summary>
[SuppressUnmanagedCodeSecurity]
[SecurityCritical]
internal static class SafeNativeMethods
{
// ReSharper disable InconsistentNaming
/// <summary>
/// Name of the native DLL.
/// </summary>
const string _DllName = "MathNet.Numerics.MKL.dll";
internal static string DllName { get { return _DllName; } }
[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern int query_capability(int capability);
[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern void set_consistency_mode(int mode);
[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern void set_vml_mode(uint mode);
[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern void set_max_threads(int num_threads);
#region Memory
[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern void free_buffers();
[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern void thread_free_buffers();
[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern int disable_fast_mm();
[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern long mem_stat([Out]out int allocatedBuffers);
[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern long peak_mem_usage(int mode);
#endregion Memory
#region FFT
[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern long z_fft_forward_inplace(long n, [In, Out] Complex[] x);
#endregion FFT
// ReSharper restore InconsistentNaming
}
}
#endif

81
src/UnitTests/FourierTransformProviderTests/FourierTransformProviderTests.cs

@ -0,0 +1,81 @@
// <copyright file="FourierTransformProviderTests.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
//
// Copyright (c) 2009-2016 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
using System;
using NUnit.Framework;
namespace MathNet.Numerics.UnitTests.FourierTransformProviderTests
{
#if NOSYSNUMERICS
using Complex = Numerics.Complex;
#else
using Complex = System.Numerics.Complex;
#endif
/// <summary>
/// Base class for linear algebra provider tests.
/// </summary>
[TestFixture, Category("LAProvider")]
public class LinearAlgebraProviderTests
{
[Test]
public void ForwardInplace()
{
var samples = Generate.PeriodicMap(16, w => new Complex(Math.Sin(w), 0), 16, 1.0, Constants.Pi2);
var spectrum = new Complex[samples.Length];
// real-odd transforms to imaginary odd
samples.Copy(spectrum);
Control.FourierTransformProvider.ForwardInplace(spectrum);
// all real components must be zero
foreach (var c in spectrum)
{
Assert.AreEqual(0, c.Real, 1e-12, "real");
}
// all imaginary components except second and last musth be zero
for (var i = 0; i < spectrum.Length; i++)
{
if (i == 1)
{
Assert.AreEqual(-8, spectrum[i].Imaginary, 1e-12, "imag second");
}
else if (i == spectrum.Length - 1)
{
Assert.AreEqual(8, spectrum[i].Imaginary, 1e-12, "imag last");
}
else
{
Assert.AreEqual(0, spectrum[i].Imaginary, 1e-12, "imag");
}
}
}
}
}
Loading…
Cancel
Save