From 8639df292fe96ee06dc8912503d7c9210d4380a0 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Fri, 7 Oct 2016 14:52:21 +0200 Subject: [PATCH] FFT: experimental implementation with test for complex-complex inplace forward FFT --- src/NativeProviders/Linux/mkl_build.sh | 4 +- src/NativeProviders/MKL/fft.cpp | 33 +++++++ src/NativeProviders/OSX/mkl_build.sh | 4 +- .../Windows/MKL/MKLWrapper.vcxproj | 1 + .../Windows/MKL/MKLWrapper.vcxproj.filters | 3 + src/Numerics/Numerics.csproj | 1 + .../IFourierTransformProvider.cs | 2 +- .../ManagedFourierTransformProvider.cs | 38 +++++++- .../Mkl/MklFourierTransformProvider.cs | 38 +++++++- .../FourierTransform/Mkl/SafeNativeMethods.cs | 93 +++++++++++++++++++ .../FourierTransformProviderTests.cs | 81 ++++++++++++++++ 11 files changed, 284 insertions(+), 14 deletions(-) create mode 100644 src/NativeProviders/MKL/fft.cpp create mode 100644 src/Numerics/Providers/FourierTransform/Mkl/SafeNativeMethods.cs create mode 100644 src/UnitTests/FourierTransformProviderTests/FourierTransformProviderTests.cs diff --git a/src/NativeProviders/Linux/mkl_build.sh b/src/NativeProviders/Linux/mkl_build.sh index b485a4db..3f497dac 100644 --- a/src/NativeProviders/Linux/mkl_build.sh +++ b/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/ diff --git a/src/NativeProviders/MKL/fft.cpp b/src/NativeProviders/MKL/fft.cpp new file mode 100644 index 00000000..8fd0c091 --- /dev/null +++ b/src/NativeProviders/MKL/fft.cpp @@ -0,0 +1,33 @@ +#include "wrapper_common.h" + +#include +#include +#include +#include +#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; + } +} diff --git a/src/NativeProviders/OSX/mkl_build.sh b/src/NativeProviders/OSX/mkl_build.sh index 5f5d5559..15483cbd 100755 --- a/src/NativeProviders/OSX/mkl_build.sh +++ b/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/ diff --git a/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj b/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj index f14c9fb1..58f87706 100644 --- a/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj +++ b/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj @@ -294,6 +294,7 @@ + diff --git a/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj.filters b/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj.filters index ce59f97e..336ae134 100644 --- a/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj.filters +++ b/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj.filters @@ -33,6 +33,9 @@ Source Files + + Source Files + diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 333f4d05..601d96ac 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -175,6 +175,7 @@ + diff --git a/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs index 28256d78..ba52fc2f 100644 --- a/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs @@ -1,4 +1,4 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics diff --git a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs index 58c3320e..14a37198 100644 --- a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs @@ -1,4 +1,32 @@ -using System.Numerics; +// +// 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. +// + +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); diff --git a/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs index 9f99097c..1c472ffe 100644 --- a/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs @@ -1,7 +1,32 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; +// +// 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. +// + +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); + } } } diff --git a/src/Numerics/Providers/FourierTransform/Mkl/SafeNativeMethods.cs b/src/Numerics/Providers/FourierTransform/Mkl/SafeNativeMethods.cs new file mode 100644 index 00000000..60078c92 --- /dev/null +++ b/src/Numerics/Providers/FourierTransform/Mkl/SafeNativeMethods.cs @@ -0,0 +1,93 @@ +// +// 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. +// + +#if NATIVE + +using System.Numerics; +using System.Runtime.InteropServices; +using System.Security; + +namespace MathNet.Numerics.Providers.FourierTransform.Mkl +{ + /// + /// P/Invoke methods to the native math libraries. + /// + [SuppressUnmanagedCodeSecurity] + [SecurityCritical] + internal static class SafeNativeMethods + { + // ReSharper disable InconsistentNaming + + /// + /// Name of the native DLL. + /// + 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 diff --git a/src/UnitTests/FourierTransformProviderTests/FourierTransformProviderTests.cs b/src/UnitTests/FourierTransformProviderTests/FourierTransformProviderTests.cs new file mode 100644 index 00000000..068b7e02 --- /dev/null +++ b/src/UnitTests/FourierTransformProviderTests/FourierTransformProviderTests.cs @@ -0,0 +1,81 @@ +// +// 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. +// + +using System; +using NUnit.Framework; + +namespace MathNet.Numerics.UnitTests.FourierTransformProviderTests +{ +#if NOSYSNUMERICS + using Complex = Numerics.Complex; +#else + using Complex = System.Numerics.Complex; +#endif + + /// + /// Base class for linear algebra provider tests. + /// + [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"); + } + } + } + } +}