From 7f5b414bbee7ce2fb1d635d255a88ca58689330d Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Fri, 7 Oct 2016 13:39:27 +0200 Subject: [PATCH 01/17] FFT: fft provider mvp interface --- src/Numerics/Numerics.csproj | 2 + .../IFourierTransformProvider.cs | 45 +++++++++++++++++++ .../ManagedFourierTransformProvider.cs | 34 ++++++++++++++ 3 files changed, 81 insertions(+) create mode 100644 src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs create mode 100644 src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 0063697a..d224559c 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -172,6 +172,8 @@ True Resources.resx + + diff --git a/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs new file mode 100644 index 00000000..7c6e5c77 --- /dev/null +++ b/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs @@ -0,0 +1,45 @@ +// +// 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. +// + +namespace MathNet.Numerics.Providers.FourierTransform +{ + +#if !NOSYSNUMERICS + using Complex = System.Numerics.Complex; +#endif + + public interface IFourierTransformProvider + { + void ForwardInplace(Complex[] complex); + void BackwardInplace(Complex[] complex); + + Complex[] Forward(Complex[] complexTimeSpace); + Complex[] Backward(Complex[] complexFrequenceSpace); + } +} diff --git a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs new file mode 100644 index 00000000..f4afed72 --- /dev/null +++ b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs @@ -0,0 +1,34 @@ +using System.Numerics; +using MathNet.Numerics.IntegralTransforms; + +namespace MathNet.Numerics.Providers.FourierTransform +{ + public class ManagedFourierTransformProvider : IFourierTransformProvider + { + public void ForwardInplace(Complex[] complex) + { + Fourier.BluesteinForward(complex, FourierOptions.Default); + } + + public void BackwardInplace(Complex[] complex) + { + Fourier.BluesteinInverse(complex, FourierOptions.Default); + } + + public Complex[] Forward(Complex[] complexTimeSpace) + { + Complex[] work = new Complex[complexTimeSpace.Length]; + complexTimeSpace.Copy(work); + ForwardInplace(work); + return work; + } + + public Complex[] Backward(Complex[] complexFrequenceSpace) + { + Complex[] work = new Complex[complexFrequenceSpace.Length]; + complexFrequenceSpace.Copy(work); + BackwardInplace(work); + return work; + } + } +} From 0486971c4aa01d57392d74f682a694571fec1d0a Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Fri, 7 Oct 2016 14:05:17 +0200 Subject: [PATCH 02/17] FFT: wire provider to Control class --- src/Numerics/Control.cs | 44 +++++++++++++++++-- src/Numerics/Numerics.csproj | 1 + .../IFourierTransformProvider.cs | 5 +++ .../ManagedFourierTransformProvider.cs | 4 ++ .../Mkl/MklFourierTransformProvider.cs | 14 ++++++ 5 files changed, 65 insertions(+), 3 deletions(-) create mode 100644 src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs diff --git a/src/Numerics/Control.cs b/src/Numerics/Control.cs index e7f9fd23..dfd8fe32 100644 --- a/src/Numerics/Control.cs +++ b/src/Numerics/Control.cs @@ -30,6 +30,7 @@ using MathNet.Numerics.Providers.LinearAlgebra; using System; using System.Threading.Tasks; +using MathNet.Numerics.Providers.FourierTransform; namespace MathNet.Numerics { @@ -45,6 +46,7 @@ namespace MathNet.Numerics static int _parallelizeOrder; static int _parallelizeElements; static ILinearAlgebraProvider _linearAlgebraProvider; + static IFourierTransformProvider _fourierTransformProvider; static readonly object _staticLock = new object(); static Control() @@ -66,11 +68,11 @@ namespace MathNet.Numerics TaskScheduler = TaskScheduler.Default; } - private static void InitializeDefaultLinearAlgebraProvider() + private static void InitializeDefaultProviders() { lock (_staticLock) { - if (_linearAlgebraProvider == null) + if (_linearAlgebraProvider == null || _fourierTransformProvider == null) { #if NATIVE try @@ -113,6 +115,7 @@ namespace MathNet.Numerics public static void UseManaged() { LinearAlgebraProvider = new ManagedLinearAlgebraProvider(); + FourierTransformProvider = new ManagedFourierTransformProvider(); } #if NATIVE @@ -124,6 +127,7 @@ namespace MathNet.Numerics public static void UseNativeMKL() { LinearAlgebraProvider = new Providers.LinearAlgebra.Mkl.MklLinearAlgebraProvider(); + FourierTransformProvider = new Providers.FourierTransform.Mkl.MklFourierTransformProvider(); } /// @@ -137,6 +141,7 @@ namespace MathNet.Numerics Providers.LinearAlgebra.Mkl.MklAccuracy accuracy = Providers.LinearAlgebra.Mkl.MklAccuracy.High) { LinearAlgebraProvider = new Providers.LinearAlgebra.Mkl.MklLinearAlgebraProvider(consistency, precision, accuracy); + FourierTransformProvider = new Providers.FourierTransform.Mkl.MklFourierTransformProvider(); } /// @@ -158,6 +163,10 @@ namespace MathNet.Numerics public static void UseNativeCUDA() { LinearAlgebraProvider = new Providers.LinearAlgebra.Cuda.CudaLinearAlgebraProvider(); + if (_fourierTransformProvider == null) + { + FourierTransformProvider = new ManagedFourierTransformProvider(); + } } /// @@ -179,6 +188,10 @@ namespace MathNet.Numerics public static void UseNativeOpenBLAS() { LinearAlgebraProvider = new Providers.LinearAlgebra.OpenBlas.OpenBlasLinearAlgebraProvider(); + if (_fourierTransformProvider == null) + { + FourierTransformProvider = new ManagedFourierTransformProvider(); + } } /// @@ -226,6 +239,7 @@ namespace MathNet.Numerics ThreadSafeRandomNumberGenerators = false; LinearAlgebraProvider.InitializeVerify(); + FourierTransformProvider.InitializeVerify(); } public static void UseMultiThreading() @@ -234,6 +248,7 @@ namespace MathNet.Numerics ThreadSafeRandomNumberGenerators = true; LinearAlgebraProvider.InitializeVerify(); + FourierTransformProvider.InitializeVerify(); } /// @@ -266,7 +281,7 @@ namespace MathNet.Numerics get { if (_linearAlgebraProvider == null) - InitializeDefaultLinearAlgebraProvider(); + InitializeDefaultProviders(); return _linearAlgebraProvider; } @@ -279,6 +294,28 @@ namespace MathNet.Numerics } } + /// + /// Gets or sets the fourier transform provider. Consider to use UseNativeMKL or UseManaged instead. + /// + /// The linear algebra provider. + public static IFourierTransformProvider FourierTransformProvider + { + get + { + if (_fourierTransformProvider == null) + InitializeDefaultProviders(); + + return _fourierTransformProvider; + } + set + { + value.InitializeVerify(); + + // only actually set if verification did not throw + _fourierTransformProvider = value; + } + } + /// /// Gets or sets a value indicating how many parallel worker threads shall be used /// when parallelization is applicable. @@ -293,6 +330,7 @@ namespace MathNet.Numerics // Reinitialize providers: LinearAlgebraProvider.InitializeVerify(); + FourierTransformProvider.InitializeVerify(); } } diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index d224559c..333f4d05 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -174,6 +174,7 @@ + diff --git a/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs index 7c6e5c77..28256d78 100644 --- a/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs @@ -36,6 +36,11 @@ namespace MathNet.Numerics.Providers.FourierTransform public interface IFourierTransformProvider { + /// + /// Initialize and verify that the provided is indeed available. If not, fall back to alternatives like the managed provider + /// + void InitializeVerify(); + void ForwardInplace(Complex[] complex); void BackwardInplace(Complex[] complex); diff --git a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs index f4afed72..58c3320e 100644 --- a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs @@ -5,6 +5,10 @@ namespace MathNet.Numerics.Providers.FourierTransform { public class ManagedFourierTransformProvider : IFourierTransformProvider { + public virtual void InitializeVerify() + { + } + public void ForwardInplace(Complex[] complex) { Fourier.BluesteinForward(complex, FourierOptions.Default); diff --git a/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs new file mode 100644 index 00000000..9f99097c --- /dev/null +++ b/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs @@ -0,0 +1,14 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +namespace MathNet.Numerics.Providers.FourierTransform.Mkl +{ + public class MklFourierTransformProvider : ManagedFourierTransformProvider + { + public override void InitializeVerify() + { + } + } +} From 8639df292fe96ee06dc8912503d7c9210d4380a0 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Fri, 7 Oct 2016 14:52:21 +0200 Subject: [PATCH 03/17] 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"); + } + } + } + } +} From 2ab73ac781250a10388e9f57506d7c88b0cb0a7c Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Fri, 7 Oct 2016 15:13:10 +0200 Subject: [PATCH 04/17] FFT: verify parseval theorem --- .../FourierTransformProviderTests.cs | 30 ++++++++++++++++++- 1 file changed, 29 insertions(+), 1 deletion(-) diff --git a/src/UnitTests/FourierTransformProviderTests/FourierTransformProviderTests.cs b/src/UnitTests/FourierTransformProviderTests/FourierTransformProviderTests.cs index 068b7e02..588768b5 100644 --- a/src/UnitTests/FourierTransformProviderTests/FourierTransformProviderTests.cs +++ b/src/UnitTests/FourierTransformProviderTests/FourierTransformProviderTests.cs @@ -28,6 +28,8 @@ // using System; +using MathNet.Numerics.Distributions; +using MathNet.Numerics.Statistics; using NUnit.Framework; namespace MathNet.Numerics.UnitTests.FourierTransformProviderTests @@ -45,7 +47,7 @@ namespace MathNet.Numerics.UnitTests.FourierTransformProviderTests public class LinearAlgebraProviderTests { [Test] - public void ForwardInplace() + public void ForwardInplaceRealSine() { var samples = Generate.PeriodicMap(16, w => new Complex(Math.Sin(w), 0), 16, 1.0, Constants.Pi2); var spectrum = new Complex[samples.Length]; @@ -77,5 +79,31 @@ namespace MathNet.Numerics.UnitTests.FourierTransformProviderTests } } } + + [TestCase(0x1000)] + [TestCase(0x7FF)] + public void ForwardInplaceParsevalTheorem(int count) + { + var samples = Generate.RandomComplex(count, GetUniform(1)); + var timeSpaceEnergy = Generate.Map(samples, s => s.MagnitudeSquared()).Mean(); + + var work = new Complex[samples.Length]; + samples.Copy(work); + + Control.FourierTransformProvider.ForwardInplace(work); + + var frequencySpaceEnergy = Generate.Map(work, s => s.MagnitudeSquared()).Mean(); + + // TODO: normalize scaling - this should instead be controllable, not needed by default + frequencySpaceEnergy /= count; + + Assert.AreEqual(timeSpaceEnergy, frequencySpaceEnergy, 1e-12); + } + + IContinuousDistribution GetUniform(int seed) + { + return new ContinuousUniform(-1, 1, new System.Random(seed)); + } + } } From 31c140cf8008f946b92e66e02f39675f559d9bab Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Fri, 7 Oct 2016 15:58:34 +0200 Subject: [PATCH 05/17] FFT: MKL backward transformation --- src/NativeProviders/MKL/fft.cpp | 51 +++++++++++++------ .../Mkl/MklFourierTransformProvider.cs | 27 ++++++++-- .../FourierTransform/Mkl/SafeNativeMethods.cs | 9 ++++ 3 files changed, 69 insertions(+), 18 deletions(-) diff --git a/src/NativeProviders/MKL/fft.cpp b/src/NativeProviders/MKL/fft.cpp index 8fd0c091..c35985b6 100644 --- a/src/NativeProviders/MKL/fft.cpp +++ b/src/NativeProviders/MKL/fft.cpp @@ -7,27 +7,48 @@ #include "mkl_service.h" #include "mkl_dfti.h" +template +inline MKL_LONG fft_inplace(MKL_LONG n, Data x[], DFTI_CONFIG_VALUE precision, DFTI_CONFIG_VALUE domain, FFT fft) +{ + MKL_LONG status = 0; + DFTI_DESCRIPTOR_HANDLE descriptor = 0; + status = DftiCreateDescriptor(&descriptor, precision, domain, 1, n); + if (0 != status) goto failed; + + status = DftiCommitDescriptor(descriptor); + if (0 != status) goto failed; + + status = fft(descriptor, x); + if (0 != status) goto failed; + +cleanup: + DftiFreeDescriptor(&descriptor); + return status; + +failed: + status = 1; + goto cleanup; +} + 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; + return fft_inplace(n, x, DFTI_DOUBLE, DFTI_COMPLEX, DftiComputeForward); + } - status = DftiComputeForward(hand, x); - if (0 != status) goto failed; + DLLEXPORT MKL_LONG c_fft_forward_inplace(MKL_LONG n, MKL_Complex8 x[]) + { + return fft_inplace(n, x, DFTI_SINGLE, DFTI_COMPLEX, DftiComputeForward); + } - cleanup: - DftiFreeDescriptor(&hand); - return status; + DLLEXPORT MKL_LONG z_fft_backward_inplace(MKL_LONG n, MKL_Complex16 x[]) + { + return fft_inplace(n, x, DFTI_DOUBLE, DFTI_COMPLEX, DftiComputeBackward); + } - failed: - status = 1; - goto cleanup; + DLLEXPORT MKL_LONG c_fft_backward_inplace(MKL_LONG n, MKL_Complex8 x[]) + { + return fft_inplace(n, x, DFTI_SINGLE, DFTI_COMPLEX, DftiComputeBackward); } } diff --git a/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs index 1c472ffe..49f57311 100644 --- a/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs @@ -30,15 +30,36 @@ using System.Numerics; namespace MathNet.Numerics.Providers.FourierTransform.Mkl { - public class MklFourierTransformProvider : ManagedFourierTransformProvider + public class MklFourierTransformProvider : IFourierTransformProvider { - public override void InitializeVerify() + public void InitializeVerify() { } - public override void ForwardInplace(Complex[] complex) + public void ForwardInplace(Complex[] complex) { SafeNativeMethods.z_fft_forward_inplace(complex.Length, complex); } + + public void BackwardInplace(Complex[] complex) + { + SafeNativeMethods.z_fft_backward_inplace(complex.Length, complex); + } + + public Complex[] Forward(Complex[] complexTimeSpace) + { + Complex[] work = new Complex[complexTimeSpace.Length]; + complexTimeSpace.Copy(work); + ForwardInplace(work); + return work; + } + + public Complex[] Backward(Complex[] complexFrequenceSpace) + { + Complex[] work = new Complex[complexFrequenceSpace.Length]; + complexFrequenceSpace.Copy(work); + BackwardInplace(work); + return work; + } } } diff --git a/src/Numerics/Providers/FourierTransform/Mkl/SafeNativeMethods.cs b/src/Numerics/Providers/FourierTransform/Mkl/SafeNativeMethods.cs index 60078c92..27cfcb5a 100644 --- a/src/Numerics/Providers/FourierTransform/Mkl/SafeNativeMethods.cs +++ b/src/Numerics/Providers/FourierTransform/Mkl/SafeNativeMethods.cs @@ -84,6 +84,15 @@ namespace MathNet.Numerics.Providers.FourierTransform.Mkl [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] internal static extern long z_fft_forward_inplace(long n, [In, Out] Complex[] x); + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern long c_fft_forward_inplace(long n, [In, Out] Complex32[] x); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern long z_fft_backward_inplace(long n, [In, Out] Complex[] x); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern long c_fft_backward_inplace(long n, [In, Out] Complex32[] x); + #endregion FFT // ReSharper restore InconsistentNaming From 3b92b38afdd34654510bd81cd4df2cecde3f4655 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Fri, 7 Oct 2016 16:44:28 +0200 Subject: [PATCH 06/17] FFT: forward and backward scaling support, drop manual correction in parseval test --- src/NativeProviders/MKL/fft.cpp | 26 +++++++----- .../IFourierTransformProvider.cs | 15 +++++-- .../ManagedFourierTransformProvider.cs | 30 ++++++++++---- .../Mkl/MklFourierTransformProvider.cs | 41 +++++++++++++++---- .../FourierTransform/Mkl/SafeNativeMethods.cs | 8 ++-- .../FourierTransformProviderTests.cs | 14 ++----- 6 files changed, 90 insertions(+), 44 deletions(-) diff --git a/src/NativeProviders/MKL/fft.cpp b/src/NativeProviders/MKL/fft.cpp index c35985b6..bc92dfb4 100644 --- a/src/NativeProviders/MKL/fft.cpp +++ b/src/NativeProviders/MKL/fft.cpp @@ -7,14 +7,20 @@ #include "mkl_service.h" #include "mkl_dfti.h" -template -inline MKL_LONG fft_inplace(MKL_LONG n, Data x[], DFTI_CONFIG_VALUE precision, DFTI_CONFIG_VALUE domain, FFT fft) +template +inline MKL_LONG fft_1d_inplace(const MKL_LONG n, Data x[], const Precision forward_scale, const Precision backward_scale, const DFTI_CONFIG_VALUE precision, const DFTI_CONFIG_VALUE domain, FFT fft) { MKL_LONG status = 0; DFTI_DESCRIPTOR_HANDLE descriptor = 0; status = DftiCreateDescriptor(&descriptor, precision, domain, 1, n); if (0 != status) goto failed; + status = DftiSetValue(descriptor, DFTI_FORWARD_SCALE, forward_scale); + if (0 != status) goto failed; + + status = DftiSetValue(descriptor, DFTI_BACKWARD_SCALE, backward_scale); + if (0 != status) goto failed; + status = DftiCommitDescriptor(descriptor); if (0 != status) goto failed; @@ -32,23 +38,23 @@ failed: extern "C" { - DLLEXPORT MKL_LONG z_fft_forward_inplace(MKL_LONG n, MKL_Complex16 x[]) + DLLEXPORT MKL_LONG z_fft_forward_inplace(const MKL_LONG n, const double scaling, MKL_Complex16 x[]) { - return fft_inplace(n, x, DFTI_DOUBLE, DFTI_COMPLEX, DftiComputeForward); + return fft_1d_inplace(n, x, scaling, 1.0, DFTI_DOUBLE, DFTI_COMPLEX, DftiComputeForward); } - DLLEXPORT MKL_LONG c_fft_forward_inplace(MKL_LONG n, MKL_Complex8 x[]) + DLLEXPORT MKL_LONG c_fft_forward_inplace(const MKL_LONG n, const float scaling, MKL_Complex8 x[]) { - return fft_inplace(n, x, DFTI_SINGLE, DFTI_COMPLEX, DftiComputeForward); + return fft_1d_inplace(n, x, scaling, 1.0f, DFTI_SINGLE, DFTI_COMPLEX, DftiComputeForward); } - DLLEXPORT MKL_LONG z_fft_backward_inplace(MKL_LONG n, MKL_Complex16 x[]) + DLLEXPORT MKL_LONG z_fft_backward_inplace(const MKL_LONG n, const double scaling, MKL_Complex16 x[]) { - return fft_inplace(n, x, DFTI_DOUBLE, DFTI_COMPLEX, DftiComputeBackward); + return fft_1d_inplace(n, x, 1.0, scaling, DFTI_DOUBLE, DFTI_COMPLEX, DftiComputeBackward); } - DLLEXPORT MKL_LONG c_fft_backward_inplace(MKL_LONG n, MKL_Complex8 x[]) + DLLEXPORT MKL_LONG c_fft_backward_inplace(const MKL_LONG n, const float scaling, MKL_Complex8 x[]) { - return fft_inplace(n, x, DFTI_SINGLE, DFTI_COMPLEX, DftiComputeBackward); + return fft_1d_inplace(n, x, 1.0f, scaling, DFTI_SINGLE, DFTI_COMPLEX, DftiComputeBackward); } } diff --git a/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs index ba52fc2f..2e4dea88 100644 --- a/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs @@ -34,6 +34,13 @@ namespace MathNet.Numerics.Providers.FourierTransform using Complex = System.Numerics.Complex; #endif + public enum FourierTransformScaling : int + { + NoScaling = 0, + SymmetricScaling = 1, + AsymmetricScaling = 2 + } + public interface IFourierTransformProvider { /// @@ -41,10 +48,10 @@ namespace MathNet.Numerics.Providers.FourierTransform /// void InitializeVerify(); - void ForwardInplace(Complex[] complex); - void BackwardInplace(Complex[] complex); + void ForwardInplace(Complex[] complex, FourierTransformScaling scaling); + void BackwardInplace(Complex[] complex, FourierTransformScaling scaling); - Complex[] Forward(Complex[] complexTimeSpace); - Complex[] Backward(Complex[] complexFrequenceSpace); + Complex[] Forward(Complex[] complexTimeSpace, FourierTransformScaling scaling); + Complex[] Backward(Complex[] complexFrequenceSpace, FourierTransformScaling scaling); } } diff --git a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs index 14a37198..696f8f95 100644 --- a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs @@ -37,30 +37,44 @@ namespace MathNet.Numerics.Providers.FourierTransform { } - public virtual void ForwardInplace(Complex[] complex) + public virtual void ForwardInplace(Complex[] complex, FourierTransformScaling scaling) { - Fourier.BluesteinForward(complex, FourierOptions.Default); + Fourier.BluesteinForward(complex, Options(scaling)); } - public virtual void BackwardInplace(Complex[] complex) + public virtual void BackwardInplace(Complex[] complex, FourierTransformScaling scaling) { - Fourier.BluesteinInverse(complex, FourierOptions.Default); + Fourier.BluesteinInverse(complex, Options(scaling)); } - public virtual Complex[] Forward(Complex[] complexTimeSpace) + public virtual Complex[] Forward(Complex[] complexTimeSpace, FourierTransformScaling scaling) { Complex[] work = new Complex[complexTimeSpace.Length]; complexTimeSpace.Copy(work); - ForwardInplace(work); + ForwardInplace(work, scaling); return work; } - public virtual Complex[] Backward(Complex[] complexFrequenceSpace) + public virtual Complex[] Backward(Complex[] complexFrequenceSpace, FourierTransformScaling scaling) { Complex[] work = new Complex[complexFrequenceSpace.Length]; complexFrequenceSpace.Copy(work); - BackwardInplace(work); + BackwardInplace(work, scaling); return work; } + + private FourierOptions Options(FourierTransformScaling scaling) + { + switch (scaling) + { + case FourierTransformScaling.NoScaling: + return FourierOptions.NoScaling; + case FourierTransformScaling.AsymmetricScaling: + return FourierOptions.AsymmetricScaling; + case FourierTransformScaling.SymmetricScaling: + default: + return FourierOptions.Default; + } + } } } diff --git a/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs index 49f57311..c3c66484 100644 --- a/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs @@ -26,6 +26,7 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System; using System.Numerics; namespace MathNet.Numerics.Providers.FourierTransform.Mkl @@ -36,30 +37,54 @@ namespace MathNet.Numerics.Providers.FourierTransform.Mkl { } - public void ForwardInplace(Complex[] complex) + public void ForwardInplace(Complex[] complex, FourierTransformScaling scaling) { - SafeNativeMethods.z_fft_forward_inplace(complex.Length, complex); + SafeNativeMethods.z_fft_forward_inplace(complex.Length, ForwardScaling(scaling, complex.Length), complex); } - public void BackwardInplace(Complex[] complex) + public void BackwardInplace(Complex[] complex, FourierTransformScaling scaling) { - SafeNativeMethods.z_fft_backward_inplace(complex.Length, complex); + SafeNativeMethods.z_fft_backward_inplace(complex.Length, BackwardScaling(scaling, complex.Length), complex); } - public Complex[] Forward(Complex[] complexTimeSpace) + public Complex[] Forward(Complex[] complexTimeSpace, FourierTransformScaling scaling) { Complex[] work = new Complex[complexTimeSpace.Length]; complexTimeSpace.Copy(work); - ForwardInplace(work); + ForwardInplace(work, scaling); return work; } - public Complex[] Backward(Complex[] complexFrequenceSpace) + public Complex[] Backward(Complex[] complexFrequenceSpace, FourierTransformScaling scaling) { Complex[] work = new Complex[complexFrequenceSpace.Length]; complexFrequenceSpace.Copy(work); - BackwardInplace(work); + BackwardInplace(work, scaling); return work; } + + private double ForwardScaling(FourierTransformScaling scaling, int length) + { + switch (scaling) + { + case FourierTransformScaling.SymmetricScaling: + return Math.Sqrt(1.0/length); + default: + return 1.0; + } + } + + private double BackwardScaling(FourierTransformScaling scaling, int length) + { + switch (scaling) + { + case FourierTransformScaling.SymmetricScaling: + return Math.Sqrt(1.0/length); + case FourierTransformScaling.AsymmetricScaling: + return 1.0/length; + default: + return 1.0; + } + } } } diff --git a/src/Numerics/Providers/FourierTransform/Mkl/SafeNativeMethods.cs b/src/Numerics/Providers/FourierTransform/Mkl/SafeNativeMethods.cs index 27cfcb5a..02be1430 100644 --- a/src/Numerics/Providers/FourierTransform/Mkl/SafeNativeMethods.cs +++ b/src/Numerics/Providers/FourierTransform/Mkl/SafeNativeMethods.cs @@ -82,16 +82,16 @@ namespace MathNet.Numerics.Providers.FourierTransform.Mkl #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); + internal static extern long z_fft_forward_inplace(long n, double scaling, [In, Out] Complex[] x); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern long c_fft_forward_inplace(long n, [In, Out] Complex32[] x); + internal static extern long c_fft_forward_inplace(long n, float scaling, [In, Out] Complex32[] x); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern long z_fft_backward_inplace(long n, [In, Out] Complex[] x); + internal static extern long z_fft_backward_inplace(long n, double scaling, [In, Out] Complex[] x); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern long c_fft_backward_inplace(long n, [In, Out] Complex32[] x); + internal static extern long c_fft_backward_inplace(long n, float scaling, [In, Out] Complex32[] x); #endregion FFT diff --git a/src/UnitTests/FourierTransformProviderTests/FourierTransformProviderTests.cs b/src/UnitTests/FourierTransformProviderTests/FourierTransformProviderTests.cs index 588768b5..91e67507 100644 --- a/src/UnitTests/FourierTransformProviderTests/FourierTransformProviderTests.cs +++ b/src/UnitTests/FourierTransformProviderTests/FourierTransformProviderTests.cs @@ -29,6 +29,7 @@ using System; using MathNet.Numerics.Distributions; +using MathNet.Numerics.Providers.FourierTransform; using MathNet.Numerics.Statistics; using NUnit.Framework; @@ -54,7 +55,7 @@ namespace MathNet.Numerics.UnitTests.FourierTransformProviderTests // real-odd transforms to imaginary odd samples.Copy(spectrum); - Control.FourierTransformProvider.ForwardInplace(spectrum); + Control.FourierTransformProvider.ForwardInplace(spectrum, FourierTransformScaling.AsymmetricScaling); // all real components must be zero foreach (var c in spectrum) @@ -87,15 +88,8 @@ namespace MathNet.Numerics.UnitTests.FourierTransformProviderTests var samples = Generate.RandomComplex(count, GetUniform(1)); var timeSpaceEnergy = Generate.Map(samples, s => s.MagnitudeSquared()).Mean(); - var work = new Complex[samples.Length]; - samples.Copy(work); - - Control.FourierTransformProvider.ForwardInplace(work); - - var frequencySpaceEnergy = Generate.Map(work, s => s.MagnitudeSquared()).Mean(); - - // TODO: normalize scaling - this should instead be controllable, not needed by default - frequencySpaceEnergy /= count; + Control.FourierTransformProvider.ForwardInplace(samples, FourierTransformScaling.SymmetricScaling); + var frequencySpaceEnergy = Generate.Map(samples, s => s.MagnitudeSquared()).Mean(); Assert.AreEqual(timeSpaceEnergy, frequencySpaceEnergy, 1e-12); } From df3264c1739c2aaf5f30df36371eb3b24770f114 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Fri, 7 Oct 2016 16:55:27 +0200 Subject: [PATCH 07/17] FFT: do not mask MKL error status --- src/NativeProviders/MKL/fft.cpp | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/src/NativeProviders/MKL/fft.cpp b/src/NativeProviders/MKL/fft.cpp index bc92dfb4..44dc474e 100644 --- a/src/NativeProviders/MKL/fft.cpp +++ b/src/NativeProviders/MKL/fft.cpp @@ -13,27 +13,23 @@ inline MKL_LONG fft_1d_inplace(const MKL_LONG n, Data x[], const Precision forwa MKL_LONG status = 0; DFTI_DESCRIPTOR_HANDLE descriptor = 0; status = DftiCreateDescriptor(&descriptor, precision, domain, 1, n); - if (0 != status) goto failed; + if (0 != status) goto cleanup; status = DftiSetValue(descriptor, DFTI_FORWARD_SCALE, forward_scale); - if (0 != status) goto failed; + if (0 != status) goto cleanup; status = DftiSetValue(descriptor, DFTI_BACKWARD_SCALE, backward_scale); - if (0 != status) goto failed; + if (0 != status) goto cleanup; status = DftiCommitDescriptor(descriptor); - if (0 != status) goto failed; + if (0 != status) goto cleanup; status = fft(descriptor, x); - if (0 != status) goto failed; + if (0 != status) goto cleanup; cleanup: DftiFreeDescriptor(&descriptor); return status; - -failed: - status = 1; - goto cleanup; } extern "C" { From 6a42fb7c83648e689041c89c33b553dc47063a07 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 8 Oct 2016 08:44:43 +0200 Subject: [PATCH 08/17] FFT: fix PCL builds, cleanup --- MathNet.Numerics.All.sln | 4 ++-- src/NativeProviders/MKL/fft.cpp | 6 ++---- src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj.filters | 2 +- .../FourierTransform/ManagedFourierTransformProvider.cs | 6 +++++- .../FourierTransform/Mkl/MklFourierTransformProvider.cs | 4 ++++ 5 files changed, 14 insertions(+), 8 deletions(-) diff --git a/MathNet.Numerics.All.sln b/MathNet.Numerics.All.sln index ec7daa56..e8d2ac77 100644 --- a/MathNet.Numerics.All.sln +++ b/MathNet.Numerics.All.sln @@ -1,7 +1,7 @@  Microsoft Visual Studio Solution File, Format Version 12.00 -# Visual Studio 2013 -VisualStudioVersion = 12.0.31101.0 +# Visual Studio 14 +VisualStudioVersion = 14.0.25420.1 MinimumVisualStudioVersion = 10.0.40219.1 Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Readme", "Readme", "{C2F37492-38AE-4186-8A7F-17B0B080942C}" ProjectSection(SolutionItems) = preProject diff --git a/src/NativeProviders/MKL/fft.cpp b/src/NativeProviders/MKL/fft.cpp index 44dc474e..b4603d71 100644 --- a/src/NativeProviders/MKL/fft.cpp +++ b/src/NativeProviders/MKL/fft.cpp @@ -4,14 +4,13 @@ #include #include #include -#include "mkl_service.h" #include "mkl_dfti.h" template inline MKL_LONG fft_1d_inplace(const MKL_LONG n, Data x[], const Precision forward_scale, const Precision backward_scale, const DFTI_CONFIG_VALUE precision, const DFTI_CONFIG_VALUE domain, FFT fft) { - MKL_LONG status = 0; - DFTI_DESCRIPTOR_HANDLE descriptor = 0; + MKL_LONG status; + DFTI_DESCRIPTOR_HANDLE descriptor = nullptr; status = DftiCreateDescriptor(&descriptor, precision, domain, 1, n); if (0 != status) goto cleanup; @@ -25,7 +24,6 @@ inline MKL_LONG fft_1d_inplace(const MKL_LONG n, Data x[], const Precision forwa if (0 != status) goto cleanup; status = fft(descriptor, x); - if (0 != status) goto cleanup; cleanup: DftiFreeDescriptor(&descriptor); diff --git a/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj.filters b/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj.filters index 336ae134..bf583e27 100644 --- a/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj.filters +++ b/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj.filters @@ -33,7 +33,7 @@ Source Files - + Source Files diff --git a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs index 696f8f95..f20ece67 100644 --- a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs @@ -26,11 +26,15 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using System.Numerics; using MathNet.Numerics.IntegralTransforms; namespace MathNet.Numerics.Providers.FourierTransform { + +#if !NOSYSNUMERICS + using Complex = System.Numerics.Complex; +#endif + public class ManagedFourierTransformProvider : IFourierTransformProvider { public virtual void InitializeVerify() diff --git a/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs index c3c66484..f1fa08c9 100644 --- a/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs @@ -26,6 +26,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // +#if NATIVE + using System; using System.Numerics; @@ -88,3 +90,5 @@ namespace MathNet.Numerics.Providers.FourierTransform.Mkl } } } + +#endif From a4f5b618e4ebac37c2c5f99d04c785df2f619c1f Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 8 Oct 2016 09:34:11 +0200 Subject: [PATCH 09/17] FFT: MKL provider refactoring due to sharing between FFT and LA --- src/NativeProviders/MKL/capabilities.cpp | 5 +- src/Numerics/Numerics.csproj | 8 +- .../Providers/Common/Mkl/MklProvider.cs | 84 +++++++++++++++ .../Mkl/MklProviderCapabilities.cs | 6 +- .../Mkl/SafeNativeMethods.cs | 21 +++- .../{ => Common}/NativeProviderLoader.cs | 4 +- .../Mkl/MklFourierTransformProvider.cs | 15 +++ .../FourierTransform/Mkl/SafeNativeMethods.cs | 102 ------------------ .../Cuda/CudaLinearAlgebraProvider.cs | 1 + .../Mkl/MklLinearAlgebraProvider.Complex.cs | 1 + .../Mkl/MklLinearAlgebraProvider.Complex32.cs | 1 + .../Mkl/MklLinearAlgebraProvider.Double.cs | 1 + .../Mkl/MklLinearAlgebraProvider.Single.cs | 1 + .../Mkl/MklLinearAlgebraProvider.cs | 45 +------- .../OpenBlas/OpenBlasLinearAlgebraProvider.cs | 1 + 15 files changed, 142 insertions(+), 154 deletions(-) create mode 100644 src/Numerics/Providers/Common/Mkl/MklProvider.cs rename src/Numerics/Providers/{LinearAlgebra => Common}/Mkl/MklProviderCapabilities.cs (91%) rename src/Numerics/Providers/{LinearAlgebra => Common}/Mkl/SafeNativeMethods.cs (95%) rename src/Numerics/Providers/{ => Common}/NativeProviderLoader.cs (99%) delete mode 100644 src/Numerics/Providers/FourierTransform/Mkl/SafeNativeMethods.cs diff --git a/src/NativeProviders/MKL/capabilities.cpp b/src/NativeProviders/MKL/capabilities.cpp index 85a7ff64..0c422644 100644 --- a/src/NativeProviders/MKL/capabilities.cpp +++ b/src/NativeProviders/MKL/capabilities.cpp @@ -41,7 +41,7 @@ extern "C" { #endif // COMMON/SHARED - case 64: return 10; // revision + case 64: return 11; // revision case 65: return 1; // numerical consistency, precision and accuracy modes case 66: return 1; // threading control case 67: return 1; // memory management @@ -54,7 +54,8 @@ extern "C" { case 256: return 0; // basic optimization // FFT - case 384: return 0; // basic FFT + case 384: return 1; // basic FFT (major - breaking) + case 385: return 0; // basic FFT (minor - non-breaking) default: return 0; // unknown or not supported diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 601d96ac..06febca7 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -172,10 +172,10 @@ True Resources.resx + - @@ -188,7 +188,7 @@ - + @@ -206,7 +206,7 @@ - + @@ -225,7 +225,7 @@ - + diff --git a/src/Numerics/Providers/Common/Mkl/MklProvider.cs b/src/Numerics/Providers/Common/Mkl/MklProvider.cs new file mode 100644 index 00000000..c096513c --- /dev/null +++ b/src/Numerics/Providers/Common/Mkl/MklProvider.cs @@ -0,0 +1,84 @@ +// +// 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; + +namespace MathNet.Numerics.Providers.Common.Mkl +{ + internal static class MklProvider + { + static int _nativeRevision; + static bool _nativeX86; + static bool _nativeX64; + static bool _nativeIA64; + + public static void Load(int minRevision) + { + int a, b; + try + { + // Load the native library + NativeProviderLoader.TryLoad(SafeNativeMethods.DllName); + + a = SafeNativeMethods.query_capability(0); + b = SafeNativeMethods.query_capability(1); + + _nativeX86 = SafeNativeMethods.query_capability((int)ProviderPlatform.x86) > 0; + _nativeX64 = SafeNativeMethods.query_capability((int)ProviderPlatform.x64) > 0; + _nativeIA64 = SafeNativeMethods.query_capability((int)ProviderPlatform.ia64) > 0; + + _nativeRevision = SafeNativeMethods.query_capability((int)ProviderConfig.Revision); + } + catch (DllNotFoundException e) + { + throw new NotSupportedException("MKL Native Provider not found.", e); + } + catch (BadImageFormatException e) + { + throw new NotSupportedException("MKL Native Provider found but failed to load. Please verify that the platform matches (x64 vs x32, Windows vs Linux).", e); + } + catch (EntryPointNotFoundException e) + { + throw new NotSupportedException("MKL Native Provider does not support capability querying and is therefore not compatible. Consider upgrading to a newer version.", e); + } + + if (a != 0 || b != -1 || _nativeRevision < minRevision) + { + throw new NotSupportedException("MKL Native Provider too old. Consider upgrading to a newer version."); + } + } + + public static string Describe() + { + return string.Format("Intel MKL ({1}; revision {0})", + _nativeRevision, + _nativeX86 ? "x86" : _nativeX64 ? "x64" : _nativeIA64 ? "IA64" : "unknown"); + } + } +} diff --git a/src/Numerics/Providers/LinearAlgebra/Mkl/MklProviderCapabilities.cs b/src/Numerics/Providers/Common/Mkl/MklProviderCapabilities.cs similarity index 91% rename from src/Numerics/Providers/LinearAlgebra/Mkl/MklProviderCapabilities.cs rename to src/Numerics/Providers/Common/Mkl/MklProviderCapabilities.cs index 85aa8a43..b1676a65 100644 --- a/src/Numerics/Providers/LinearAlgebra/Mkl/MklProviderCapabilities.cs +++ b/src/Numerics/Providers/Common/Mkl/MklProviderCapabilities.cs @@ -27,7 +27,7 @@ // OTHER DEALINGS IN THE SOFTWARE. // -namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl +namespace MathNet.Numerics.Providers.Common.Mkl { internal enum ProviderPlatform : int { @@ -46,7 +46,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl internal enum ProviderCapability : int { - LinearAlgebra = 128, + LinearAlgebraMajor = 128, LinearAlgebraMinor = 129, + FourierTransformMajor = 384, + FourierTransformMinor = 385 } } diff --git a/src/Numerics/Providers/LinearAlgebra/Mkl/SafeNativeMethods.cs b/src/Numerics/Providers/Common/Mkl/SafeNativeMethods.cs similarity index 95% rename from src/Numerics/Providers/LinearAlgebra/Mkl/SafeNativeMethods.cs rename to src/Numerics/Providers/Common/Mkl/SafeNativeMethods.cs index 6fc00840..287b5c83 100644 --- a/src/Numerics/Providers/LinearAlgebra/Mkl/SafeNativeMethods.cs +++ b/src/Numerics/Providers/Common/Mkl/SafeNativeMethods.cs @@ -2,7 +2,7 @@ // Math.NET Numerics, part of the Math.NET Project // http://mathnet.opensourcedotnet.info // -// Copyright (c) 2009-2014 Math.NET +// 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 @@ -31,8 +31,9 @@ using System.Numerics; using System.Runtime.InteropServices; using System.Security; +using MathNet.Numerics.Providers.LinearAlgebra; -namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl +namespace MathNet.Numerics.Providers.Common.Mkl { /// /// P/Invoke methods to the native math libraries. @@ -367,6 +368,22 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl #endregion Vector Functions + #region FFT + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern long z_fft_forward_inplace(long n, double scaling, [In, Out] Complex[] x); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern long c_fft_forward_inplace(long n, float scaling, [In, Out] Complex32[] x); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern long z_fft_backward_inplace(long n, double scaling, [In, Out] Complex[] x); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern long c_fft_backward_inplace(long n, float scaling, [In, Out] Complex32[] x); + + #endregion FFT + // ReSharper restore InconsistentNaming } } diff --git a/src/Numerics/Providers/NativeProviderLoader.cs b/src/Numerics/Providers/Common/NativeProviderLoader.cs similarity index 99% rename from src/Numerics/Providers/NativeProviderLoader.cs rename to src/Numerics/Providers/Common/NativeProviderLoader.cs index 5f1e4193..0b4f19e4 100644 --- a/src/Numerics/Providers/NativeProviderLoader.cs +++ b/src/Numerics/Providers/Common/NativeProviderLoader.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2015 Math.NET +// 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 @@ -37,7 +37,7 @@ using System.Threading; #if NATIVE -namespace MathNet.Numerics.Providers +namespace MathNet.Numerics.Providers.Common { /// /// Helper class to load native libraries depending on the architecture of the OS and process. diff --git a/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs index f1fa08c9..50c08ab3 100644 --- a/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs @@ -30,6 +30,7 @@ using System; using System.Numerics; +using MathNet.Numerics.Providers.Common.Mkl; namespace MathNet.Numerics.Providers.FourierTransform.Mkl { @@ -37,6 +38,20 @@ namespace MathNet.Numerics.Providers.FourierTransform.Mkl { public void InitializeVerify() { + MklProvider.Load(minRevision: 11); + + // we only support exactly one major version, since major version changes imply a breaking change. + int fftMajor = SafeNativeMethods.query_capability((int)ProviderCapability.FourierTransformMajor); + int fftMinor = SafeNativeMethods.query_capability((int)ProviderCapability.FourierTransformMinor); + if (!(fftMajor == 1 && fftMinor >= 0)) + { + throw new NotSupportedException(string.Format("MKL Native Provider not compatible. Expecting fourier transform v1 but provider implements v{0}.", fftMajor)); + } + } + + public override string ToString() + { + return MklProvider.Describe(); } public void ForwardInplace(Complex[] complex, FourierTransformScaling scaling) diff --git a/src/Numerics/Providers/FourierTransform/Mkl/SafeNativeMethods.cs b/src/Numerics/Providers/FourierTransform/Mkl/SafeNativeMethods.cs deleted file mode 100644 index 02be1430..00000000 --- a/src/Numerics/Providers/FourierTransform/Mkl/SafeNativeMethods.cs +++ /dev/null @@ -1,102 +0,0 @@ -// -// 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, double scaling, [In, Out] Complex[] x); - - [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern long c_fft_forward_inplace(long n, float scaling, [In, Out] Complex32[] x); - - [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern long z_fft_backward_inplace(long n, double scaling, [In, Out] Complex[] x); - - [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern long c_fft_backward_inplace(long n, float scaling, [In, Out] Complex32[] x); - - #endregion FFT - - // ReSharper restore InconsistentNaming - } -} - -#endif diff --git a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.cs index ff93ac97..6c814870 100644 --- a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.cs +++ b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.cs @@ -28,6 +28,7 @@ // using System; +using MathNet.Numerics.Providers.Common; #if NATIVE diff --git a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex.cs b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex.cs index d698d188..880620c4 100644 --- a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex.cs +++ b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex.cs @@ -34,6 +34,7 @@ using System.Numerics; using System.Security; using MathNet.Numerics.LinearAlgebra.Factorization; using MathNet.Numerics.Properties; +using MathNet.Numerics.Providers.Common.Mkl; namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl { diff --git a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex32.cs b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex32.cs index be916438..901020c2 100644 --- a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex32.cs +++ b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Complex32.cs @@ -34,6 +34,7 @@ using System.Numerics; using System.Security; using MathNet.Numerics.LinearAlgebra.Factorization; using MathNet.Numerics.Properties; +using MathNet.Numerics.Providers.Common.Mkl; namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl { diff --git a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Double.cs b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Double.cs index 081ed4ff..f6712217 100644 --- a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Double.cs +++ b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Double.cs @@ -34,6 +34,7 @@ using System.Numerics; using System.Security; using MathNet.Numerics.LinearAlgebra.Factorization; using MathNet.Numerics.Properties; +using MathNet.Numerics.Providers.Common.Mkl; namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl { diff --git a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Single.cs b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Single.cs index c58d65a8..2a9a2c79 100644 --- a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Single.cs +++ b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.Single.cs @@ -34,6 +34,7 @@ using System.Numerics; using System.Security; using MathNet.Numerics.LinearAlgebra.Factorization; using MathNet.Numerics.Properties; +using MathNet.Numerics.Providers.Common.Mkl; namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl { diff --git a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs index 8c526e52..a107839a 100644 --- a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs +++ b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs @@ -28,6 +28,8 @@ // using System; +using MathNet.Numerics.Providers.Common; +using MathNet.Numerics.Providers.Common.Mkl; #if NATIVE @@ -105,11 +107,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl /// public partial class MklLinearAlgebraProvider : ManagedLinearAlgebraProvider { - int _nativeRevision; - bool _nativeIX86; - bool _nativeX64; - bool _nativeIA64; - readonly MklConsistency _consistency; readonly MklPrecision _precision; readonly MklAccuracy _accuracy; @@ -144,39 +141,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl /// public override void InitializeVerify() { - int a, b, linearAlgebra; - try - { - // Load the native library - NativeProviderLoader.TryLoad(SafeNativeMethods.DllName); - - a = SafeNativeMethods.query_capability(0); - b = SafeNativeMethods.query_capability(1); - - _nativeIX86 = SafeNativeMethods.query_capability((int)ProviderPlatform.x86) > 0; - _nativeX64 = SafeNativeMethods.query_capability((int)ProviderPlatform.x64) > 0; - _nativeIA64 = SafeNativeMethods.query_capability((int)ProviderPlatform.ia64) > 0; - - _nativeRevision = SafeNativeMethods.query_capability((int)ProviderConfig.Revision); - linearAlgebra = SafeNativeMethods.query_capability((int)ProviderCapability.LinearAlgebra); - } - catch (DllNotFoundException e) - { - throw new NotSupportedException("MKL Native Provider not found.", e); - } - catch (BadImageFormatException e) - { - throw new NotSupportedException("MKL Native Provider found but failed to load. Please verify that the platform matches (x64 vs x32, Windows vs Linux).", e); - } - catch (EntryPointNotFoundException e) - { - throw new NotSupportedException("MKL Native Provider does not support capability querying and is therefore not compatible. Consider upgrading to a newer version.", e); - } + MklProvider.Load(minRevision: 4); - if (a != 0 || b != -1 || _nativeRevision < 4) - { - throw new NotSupportedException("MKL Native Provider too old. Consider upgrading to a newer version."); - } + int linearAlgebra = SafeNativeMethods.query_capability((int)ProviderCapability.LinearAlgebraMajor); // we only support exactly one major version, since major version changes imply a breaking change. if (linearAlgebra != 2) @@ -295,9 +262,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl public override string ToString() { - return string.Format("Intel MKL ({1}; revision {0})", - _nativeRevision, - _nativeIX86 ? "x86" : _nativeX64 ? "x64" : _nativeIA64 ? "IA64" : "unknown"); + return MklProvider.Describe(); } } } diff --git a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.cs index 89e8bb5e..433243cf 100644 --- a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.cs +++ b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.cs @@ -30,6 +30,7 @@ #if NATIVE using System; +using MathNet.Numerics.Providers.Common; namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas { From d7440d51259bb03004b329051921b20938c0e55c Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 8 Oct 2016 10:02:30 +0200 Subject: [PATCH 10/17] FFT: MKL provider refactoring due to sharing between FFT and LA - 2 --- src/Numerics/Control.cs | 15 ++ src/Numerics/Numerics.csproj | 1 + .../Providers/Common/Mkl/MklProvider.cs | 139 ++++++++++++++++++ .../Common/Mkl/MklProviderCapabilities.cs | 2 +- .../Common/Mkl/MklProviderPrecision.cs | 66 +++++++++ .../Mkl/MklFourierTransformProvider.cs | 60 ++++++++ .../LinearAlgebra/ILinearAlgebraProvider.cs | 12 +- .../Mkl/MklLinearAlgebraProvider.cs | 119 +++++---------- 8 files changed, 323 insertions(+), 91 deletions(-) create mode 100644 src/Numerics/Providers/Common/Mkl/MklProviderPrecision.cs diff --git a/src/Numerics/Control.cs b/src/Numerics/Control.cs index dfd8fe32..4045a7d4 100644 --- a/src/Numerics/Control.cs +++ b/src/Numerics/Control.cs @@ -135,6 +135,7 @@ namespace MathNet.Numerics /// Throws if it is not available or failed to initialize, in which case the previous provider is still active. /// [CLSCompliant(false)] + [Obsolete("Will be removed in the next major version. Use the enums in the Common namespace instead.")] public static void UseNativeMKL( Providers.LinearAlgebra.Mkl.MklConsistency consistency = Providers.LinearAlgebra.Mkl.MklConsistency.Auto, Providers.LinearAlgebra.Mkl.MklPrecision precision = Providers.LinearAlgebra.Mkl.MklPrecision.Double, @@ -144,6 +145,20 @@ namespace MathNet.Numerics FourierTransformProvider = new Providers.FourierTransform.Mkl.MklFourierTransformProvider(); } + /// + /// Use the Intel MKL native provider for linear algebra, with the specified configuration parameters. + /// Throws if it is not available or failed to initialize, in which case the previous provider is still active. + /// + [CLSCompliant(false)] + public static void UseNativeMKL( + Providers.Common.Mkl.MklConsistency consistency = Providers.Common.Mkl.MklConsistency.Auto, + Providers.Common.Mkl.MklPrecision precision = Providers.Common.Mkl.MklPrecision.Double, + Providers.Common.Mkl.MklAccuracy accuracy = Providers.Common.Mkl.MklAccuracy.High) + { + LinearAlgebraProvider = new Providers.LinearAlgebra.Mkl.MklLinearAlgebraProvider(consistency, precision, accuracy); + FourierTransformProvider = new Providers.FourierTransform.Mkl.MklFourierTransformProvider(); + } + /// /// Try to use the Intel MKL native provider for linear algebra. /// diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 06febca7..842ae7c9 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -173,6 +173,7 @@ Resources.resx + diff --git a/src/Numerics/Providers/Common/Mkl/MklProvider.cs b/src/Numerics/Providers/Common/Mkl/MklProvider.cs index c096513c..86618320 100644 --- a/src/Numerics/Providers/Common/Mkl/MklProvider.cs +++ b/src/Numerics/Providers/Common/Mkl/MklProvider.cs @@ -72,6 +72,122 @@ namespace MathNet.Numerics.Providers.Common.Mkl { throw new NotSupportedException("MKL Native Provider too old. Consider upgrading to a newer version."); } + + ConfigureThreading(); + } + + static void ConfigureThreading() + { + // set threading settings, if supported + if (SafeNativeMethods.query_capability((int)ProviderConfig.Threading) > 0) + { + SafeNativeMethods.set_max_threads(Control.MaxDegreeOfParallelism); + } + } + + public static void ConfigurePrecision(MklConsistency consistency, MklPrecision precision, MklAccuracy accuracy) + { + // set numerical consistency, precision and accuracy modes, if supported + if (SafeNativeMethods.query_capability((int)ProviderConfig.Precision) > 0) + { + SafeNativeMethods.set_consistency_mode((int)consistency); + SafeNativeMethods.set_vml_mode((uint)precision | (uint)accuracy); + } + } + + /// + /// Frees the memory allocated to the MKL memory pool. + /// + public static void FreeBuffers() + { + if (SafeNativeMethods.query_capability((int)ProviderConfig.Memory) < 1) + { + throw new NotSupportedException("MKL Native Provider does not support memory management functions. Consider upgrading to a newer version."); + } + + SafeNativeMethods.free_buffers(); + } + + /// + /// Frees the memory allocated to the MKL memory pool on the current thread. + /// + public static void ThreadFreeBuffers() + { + if (SafeNativeMethods.query_capability((int)ProviderConfig.Memory) < 1) + { + throw new NotSupportedException("MKL Native Provider does not support memory management functions. Consider upgrading to a newer version."); + } + + SafeNativeMethods.thread_free_buffers(); + } + + /// + /// Disable the MKL memory pool. May impact performance. + /// + public static void DisableMemoryPool() + { + if (SafeNativeMethods.query_capability((int)ProviderConfig.Memory) < 1) + { + throw new NotSupportedException("MKL Native Provider does not support memory management functions. Consider upgrading to a newer version."); + } + + SafeNativeMethods.disable_fast_mm(); + } + + /// + /// Retrieves information about the MKL memory pool. + /// + /// On output, returns the number of memory buffers allocated. + /// Returns the number of bytes allocated to all memory buffers. + public static long MemoryStatistics(out int allocatedBuffers) + { + if (SafeNativeMethods.query_capability((int)ProviderConfig.Memory) < 1) + { + throw new NotSupportedException("MKL Native Provider does not support memory management functions. Consider upgrading to a newer version."); + } + + return SafeNativeMethods.mem_stat(out allocatedBuffers); + } + + /// + /// Enable gathering of peak memory statistics of the MKL memory pool. + /// + public static void EnablePeakMemoryStatistics() + { + if (SafeNativeMethods.query_capability((int)ProviderConfig.Memory) < 1) + { + throw new NotSupportedException("MKL Native Provider does not support memory management functions. Consider upgrading to a newer version."); + } + + SafeNativeMethods.peak_mem_usage((int)MklMemoryRequestMode.Enable); + } + + /// + /// Disable gathering of peak memory statistics of the MKL memory pool. + /// + public static void DisablePeakMemoryStatistics() + { + if (SafeNativeMethods.query_capability((int)ProviderConfig.Memory) < 1) + { + throw new NotSupportedException("MKL Native Provider does not support memory management functions. Consider upgrading to a newer version."); + } + + SafeNativeMethods.peak_mem_usage((int)MklMemoryRequestMode.Disable); + } + + /// + /// Measures peak memory usage of the MKL memory pool. + /// + /// Whether the usage counter should be reset. + /// The peak number of bytes allocated to all memory buffers. + public static long PeakMemoryStatistics(bool reset = true) + { + if (SafeNativeMethods.query_capability((int)ProviderConfig.Memory) < 1) + { + throw new NotSupportedException("MKL Native Provider does not support memory management functions. Consider upgrading to a newer version."); + } + + return SafeNativeMethods.peak_mem_usage((int)(reset ? MklMemoryRequestMode.PeakMemoryReset : MklMemoryRequestMode.PeakMemory)); } public static string Describe() @@ -80,5 +196,28 @@ namespace MathNet.Numerics.Providers.Common.Mkl _nativeRevision, _nativeX86 ? "x86" : _nativeX64 ? "x64" : _nativeIA64 ? "IA64" : "unknown"); } + + enum MklMemoryRequestMode : int + { + /// + /// Disable gathering memory usage + /// + Disable = 0, + + /// + /// Enable gathering memory usage + /// + Enable = 1, + + /// + /// Return peak memory usage + /// + PeakMemory = 2, + + /// + /// Return peak memory usage and reset counter + /// + PeakMemoryReset = -1 + } } } diff --git a/src/Numerics/Providers/Common/Mkl/MklProviderCapabilities.cs b/src/Numerics/Providers/Common/Mkl/MklProviderCapabilities.cs index b1676a65..908c04a2 100644 --- a/src/Numerics/Providers/Common/Mkl/MklProviderCapabilities.cs +++ b/src/Numerics/Providers/Common/Mkl/MklProviderCapabilities.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2015 Math.NET +// 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 diff --git a/src/Numerics/Providers/Common/Mkl/MklProviderPrecision.cs b/src/Numerics/Providers/Common/Mkl/MklProviderPrecision.cs new file mode 100644 index 00000000..45c5ed6c --- /dev/null +++ b/src/Numerics/Providers/Common/Mkl/MklProviderPrecision.cs @@ -0,0 +1,66 @@ +// +// 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; + +namespace MathNet.Numerics.Providers.Common.Mkl +{ + /// + /// Consistency vs. performance trade-off between runs on different machines. + /// + public enum MklConsistency : int + { + /// Consistent on the same CPU only (maximum performance) + Auto = 2, + /// Consistent on Intel and compatible CPUs with SSE2 support (maximum compatibility) + Compatible = 3, + /// Consistent on Intel CPUs supporting SSE2 or later + SSE2 = 4, + /// Consistent on Intel CPUs supporting SSE4.2 or later + SSE4_2 = 8, + /// Consistent on Intel CPUs supporting AVX or later + AVX = 9, + /// Consistent on Intel CPUs supporting AVX2 or later + AVX2 = 10 + } + + [CLSCompliant(false)] + public enum MklAccuracy : uint + { + Low = 0x1, + High = 0x2 + } + + [CLSCompliant(false)] + public enum MklPrecision : uint + { + Single = 0x10, + Double = 0x20 + } +} diff --git a/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs index 50c08ab3..afb9ac74 100644 --- a/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs @@ -49,6 +49,66 @@ namespace MathNet.Numerics.Providers.FourierTransform.Mkl } } + /// + /// Frees the memory allocated to the MKL memory pool. + /// + public void FreeBuffers() + { + MklProvider.FreeBuffers(); + } + + /// + /// Frees the memory allocated to the MKL memory pool on the current thread. + /// + public void ThreadFreeBuffers() + { + MklProvider.ThreadFreeBuffers(); + } + + /// + /// Disable the MKL memory pool. May impact performance. + /// + public void DisableMemoryPool() + { + MklProvider.DisableMemoryPool(); + } + + /// + /// Retrieves information about the MKL memory pool. + /// + /// On output, returns the number of memory buffers allocated. + /// Returns the number of bytes allocated to all memory buffers. + public long MemoryStatistics(out int allocatedBuffers) + { + return MklProvider.MemoryStatistics(out allocatedBuffers); + } + + /// + /// Enable gathering of peak memory statistics of the MKL memory pool. + /// + public void EnablePeakMemoryStatistics() + { + MklProvider.EnablePeakMemoryStatistics(); + } + + /// + /// Disable gathering of peak memory statistics of the MKL memory pool. + /// + public void DisablePeakMemoryStatistics() + { + MklProvider.DisablePeakMemoryStatistics(); + } + + /// + /// Measures peak memory usage of the MKL memory pool. + /// + /// Whether the usage counter should be reset. + /// The peak number of bytes allocated to all memory buffers. + public long PeakMemoryStatistics(bool reset = true) + { + return MklProvider.PeakMemoryStatistics(reset); + } + public override string ToString() { return MklProvider.Describe(); diff --git a/src/Numerics/Providers/LinearAlgebra/ILinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/ILinearAlgebraProvider.cs index 57ed0707..2298581a 100644 --- a/src/Numerics/Providers/LinearAlgebra/ILinearAlgebraProvider.cs +++ b/src/Numerics/Providers/LinearAlgebra/ILinearAlgebraProvider.cs @@ -2,9 +2,9 @@ // Math.NET Numerics, part of the Math.NET Project // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics -// +// // Copyright (c) 2009-2013 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 @@ -13,10 +13,10 @@ // 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 @@ -151,7 +151,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra T DotProduct(T[] x, T[] y); /// - /// Does a point wise add of two arrays z = x + y. This can be used + /// Does a point wise add of two arrays z = x + y. This can be used /// to add vectors or matrices. /// /// The array x. @@ -163,7 +163,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra void AddArrays(T[] x, T[] y, T[] result); /// - /// Does a point wise subtraction of two arrays z = x - y. This can be used + /// Does a point wise subtraction of two arrays z = x - y. This can be used /// to subtract vectors or matrices. /// /// The array x. diff --git a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs index a107839a..92f7a316 100644 --- a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs +++ b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs @@ -49,6 +49,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl /// /// Consistency vs. performance trade-off between runs on different machines. /// + [Obsolete("Will be removed in the next major version. Use the enums in the Common namespace instead.")] public enum MklConsistency : int { /// Consistent on the same CPU only (maximum performance) @@ -66,6 +67,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl } [CLSCompliant(false)] + [Obsolete("Will be removed in the next major version. Use the enums in the Common namespace instead.")] public enum MklAccuracy : uint { Low = 0x1, @@ -73,43 +75,21 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl } [CLSCompliant(false)] + [Obsolete("Will be removed in the next major version. Use the enums in the Common namespace instead.")] public enum MklPrecision : uint { Single = 0x10, Double = 0x20 } - internal enum MklMemoryRequestMode : int - { - /// - /// Disable gathering memory usage - /// - Disable = 0, - - /// - /// Enable gathering memory usage - /// - Enable = 1, - - /// - /// Return peak memory usage - /// - PeakMemory = 2, - - /// - /// Return peak memory usage and reset counter - /// - PeakMemoryReset = -1 - } - /// /// Intel's Math Kernel Library (MKL) linear algebra provider. /// public partial class MklLinearAlgebraProvider : ManagedLinearAlgebraProvider { - readonly MklConsistency _consistency; - readonly MklPrecision _precision; - readonly MklAccuracy _accuracy; + readonly Common.Mkl.MklConsistency _consistency; + readonly Common.Mkl.MklPrecision _precision; + readonly Common.Mkl.MklAccuracy _accuracy; /// /// Sets the desired bit consistency on repeated identical computations on varying CPU architectures, @@ -118,10 +98,28 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl /// VML optimal precision and rounding. /// VML accuracy mode. [CLSCompliant(false)] + [Obsolete("Will be removed in the next major version. Use the enums in the Common namespace instead.")] public MklLinearAlgebraProvider( MklConsistency consistency = MklConsistency.Auto, MklPrecision precision = MklPrecision.Double, MklAccuracy accuracy = MklAccuracy.High) + { + _consistency = (Common.Mkl.MklConsistency)consistency; + _precision = (Common.Mkl.MklPrecision)precision; + _accuracy = (Common.Mkl.MklAccuracy)accuracy; + } + + /// + /// Sets the desired bit consistency on repeated identical computations on varying CPU architectures, + /// as a trade-off with performance. + /// + /// VML optimal precision and rounding. + /// VML accuracy mode. + [CLSCompliant(false)] + public MklLinearAlgebraProvider( + Common.Mkl.MklConsistency consistency = Common.Mkl.MklConsistency.Auto, + Common.Mkl.MklPrecision precision = Common.Mkl.MklPrecision.Double, + Common.Mkl.MklAccuracy accuracy = Common.Mkl.MklAccuracy.High) { _consistency = consistency; _precision = precision; @@ -130,9 +128,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl public MklLinearAlgebraProvider() { - _consistency = MklConsistency.Auto; - _precision = MklPrecision.Double; - _accuracy = MklAccuracy.High; + _consistency = Common.Mkl.MklConsistency.Auto; + _precision = Common.Mkl.MklPrecision.Double; + _accuracy = Common.Mkl.MklAccuracy.High; } /// @@ -142,6 +140,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl public override void InitializeVerify() { MklProvider.Load(minRevision: 4); + MklProvider.ConfigurePrecision(_consistency, _precision, _accuracy); int linearAlgebra = SafeNativeMethods.query_capability((int)ProviderCapability.LinearAlgebraMajor); @@ -150,19 +149,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl { throw new NotSupportedException(string.Format("MKL Native Provider not compatible. Expecting linear algebra v2 but provider implements v{0}.", linearAlgebra)); } - - // set numerical consistency, precision and accuracy modes, if supported - if (SafeNativeMethods.query_capability((int)ProviderConfig.Precision) > 0) - { - SafeNativeMethods.set_consistency_mode((int)_consistency); - SafeNativeMethods.set_vml_mode((uint)_precision | (uint)_accuracy); - } - - // set threading settings, if supported - if (SafeNativeMethods.query_capability((int)ProviderConfig.Threading) > 0) - { - SafeNativeMethods.set_max_threads(Control.MaxDegreeOfParallelism); - } } /// @@ -170,12 +156,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl /// public void FreeBuffers() { - if (SafeNativeMethods.query_capability((int)ProviderConfig.Memory) < 1) - { - throw new NotSupportedException("MKL Native Provider does not support memory management functions. Consider upgrading to a newer version."); - } - - SafeNativeMethods.free_buffers(); + MklProvider.FreeBuffers(); } /// @@ -183,12 +164,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl /// public void ThreadFreeBuffers() { - if (SafeNativeMethods.query_capability((int)ProviderConfig.Memory) < 1) - { - throw new NotSupportedException("MKL Native Provider does not support memory management functions. Consider upgrading to a newer version."); - } - - SafeNativeMethods.thread_free_buffers(); + MklProvider.ThreadFreeBuffers(); } /// @@ -196,12 +172,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl /// public void DisableMemoryPool() { - if (SafeNativeMethods.query_capability((int)ProviderConfig.Memory) < 1) - { - throw new NotSupportedException("MKL Native Provider does not support memory management functions. Consider upgrading to a newer version."); - } - - SafeNativeMethods.disable_fast_mm(); + MklProvider.DisableMemoryPool(); } /// @@ -211,12 +182,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl /// Returns the number of bytes allocated to all memory buffers. public long MemoryStatistics(out int allocatedBuffers) { - if (SafeNativeMethods.query_capability((int)ProviderConfig.Memory) < 1) - { - throw new NotSupportedException("MKL Native Provider does not support memory management functions. Consider upgrading to a newer version."); - } - - return SafeNativeMethods.mem_stat(out allocatedBuffers); + return MklProvider.MemoryStatistics(out allocatedBuffers); } /// @@ -224,12 +190,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl /// public void EnablePeakMemoryStatistics() { - if (SafeNativeMethods.query_capability((int)ProviderConfig.Memory) < 1) - { - throw new NotSupportedException("MKL Native Provider does not support memory management functions. Consider upgrading to a newer version."); - } - - SafeNativeMethods.peak_mem_usage((int)MklMemoryRequestMode.Enable); + MklProvider.EnablePeakMemoryStatistics(); } /// @@ -237,12 +198,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl /// public void DisablePeakMemoryStatistics() { - if (SafeNativeMethods.query_capability((int)ProviderConfig.Memory) < 1) - { - throw new NotSupportedException("MKL Native Provider does not support memory management functions. Consider upgrading to a newer version."); - } - - SafeNativeMethods.peak_mem_usage((int)MklMemoryRequestMode.Disable); + MklProvider.DisablePeakMemoryStatistics(); } /// @@ -252,12 +208,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl /// The peak number of bytes allocated to all memory buffers. public long PeakMemoryStatistics(bool reset = true) { - if (SafeNativeMethods.query_capability((int)ProviderConfig.Memory) < 1) - { - throw new NotSupportedException("MKL Native Provider does not support memory management functions. Consider upgrading to a newer version."); - } - - return SafeNativeMethods.peak_mem_usage((int)(reset ? MklMemoryRequestMode.PeakMemoryReset : MklMemoryRequestMode.PeakMemory)); + return MklProvider.PeakMemoryStatistics(reset); } public override string ToString() From ea3b57c4b1ba2a65be17a446cd9e34ade76e43ad Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 8 Oct 2016 10:55:48 +0200 Subject: [PATCH 11/17] FFT: include MKL version in MKL provider description (ToString) --- src/NativeProviders/MKL/capabilities.cpp | 20 ++++++++++++++++ .../Providers/Common/Mkl/MklProvider.cs | 23 ++++++++++++++++--- .../Common/Mkl/MklProviderCapabilities.cs | 3 +++ 3 files changed, 43 insertions(+), 3 deletions(-) diff --git a/src/NativeProviders/MKL/capabilities.cpp b/src/NativeProviders/MKL/capabilities.cpp index 0c422644..0c86282e 100644 --- a/src/NativeProviders/MKL/capabilities.cpp +++ b/src/NativeProviders/MKL/capabilities.cpp @@ -40,6 +40,26 @@ extern "C" { return 0; #endif + // MKL VERSION + case 32: // major version + { + MKLVersion Version; + mkl_get_version(&Version); + return Version.MajorVersion; + } + case 33: // minor version + { + MKLVersion Version; + mkl_get_version(&Version); + return Version.MinorVersion; + } + case 34: // update version + { + MKLVersion Version; + mkl_get_version(&Version); + return Version.UpdateVersion; + } + // COMMON/SHARED case 64: return 11; // revision case 65: return 1; // numerical consistency, precision and accuracy modes diff --git a/src/Numerics/Providers/Common/Mkl/MklProvider.cs b/src/Numerics/Providers/Common/Mkl/MklProvider.cs index 86618320..9a0ee156 100644 --- a/src/Numerics/Providers/Common/Mkl/MklProvider.cs +++ b/src/Numerics/Providers/Common/Mkl/MklProvider.cs @@ -28,11 +28,13 @@ // using System; +using System.Collections.Generic; namespace MathNet.Numerics.Providers.Common.Mkl { internal static class MklProvider { + static Version _mklVersion; static int _nativeRevision; static bool _nativeX86; static bool _nativeX64; @@ -54,6 +56,11 @@ namespace MathNet.Numerics.Providers.Common.Mkl _nativeIA64 = SafeNativeMethods.query_capability((int)ProviderPlatform.ia64) > 0; _nativeRevision = SafeNativeMethods.query_capability((int)ProviderConfig.Revision); + + _mklVersion = new Version( + SafeNativeMethods.query_capability((int)ProviderConfig.MklMajorVersion), + SafeNativeMethods.query_capability((int)ProviderConfig.MklMinorVersion), + SafeNativeMethods.query_capability((int)ProviderConfig.MklUpdateVersion)); } catch (DllNotFoundException e) { @@ -192,9 +199,19 @@ namespace MathNet.Numerics.Providers.Common.Mkl public static string Describe() { - return string.Format("Intel MKL ({1}; revision {0})", - _nativeRevision, - _nativeX86 ? "x86" : _nativeX64 ? "x64" : _nativeIA64 ? "IA64" : "unknown"); + var parts = new List(); + if (_nativeX86) parts.Add("x86"); + if (_nativeX64) parts.Add("x64"); + if (_nativeIA64) parts.Add("IA64"); + parts.Add("revision " + _nativeRevision); + if (_mklVersion.Major > 0) + { + parts.Add(_mklVersion.Build == 0 + ? string.Concat("MKL ", _mklVersion.ToString(2)) + : string.Concat("MKL ", _mklVersion.ToString(2), " Update ", _mklVersion.Build)); + } + + return string.Concat("Intel MKL (", string.Join("; ", parts), ")"); } enum MklMemoryRequestMode : int diff --git a/src/Numerics/Providers/Common/Mkl/MklProviderCapabilities.cs b/src/Numerics/Providers/Common/Mkl/MklProviderCapabilities.cs index 908c04a2..af0f971d 100644 --- a/src/Numerics/Providers/Common/Mkl/MklProviderCapabilities.cs +++ b/src/Numerics/Providers/Common/Mkl/MklProviderCapabilities.cs @@ -38,6 +38,9 @@ namespace MathNet.Numerics.Providers.Common.Mkl internal enum ProviderConfig : int { + MklMajorVersion = 32, + MklMinorVersion = 33, + MklUpdateVersion = 34, Revision = 64, Precision = 65, Threading = 66, From 037c1043cb737dde0327f2849e6006d4a5c64a54 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 8 Oct 2016 14:24:25 +0200 Subject: [PATCH 12/17] FFT: rework provider discovery to throw less exceptions and work better when FFT and LA use different providers --- src/Numerics/Control.cs | 99 +++++---------- src/Numerics/Numerics.csproj | 2 + .../Providers/Common/Mkl/MklProvider.cs | 24 +++- .../FourierTransformControl.cs | 97 ++++++++++++++ .../IFourierTransformProvider.cs | 6 + .../ManagedFourierTransformProvider.cs | 12 ++ .../Mkl/MklFourierTransformProvider.cs | 12 ++ .../Cuda/CudaLinearAlgebraProvider.cs | 28 +++- .../LinearAlgebra/ILinearAlgebraProvider.cs | 6 + .../LinearAlgebra/LinearAlgebraControl.cs | 120 ++++++++++++++++++ .../ManagedLinearAlgebraProvider.cs | 9 ++ .../Mkl/MklLinearAlgebraProvider.cs | 11 +- .../OpenBlas/OpenBlasLinearAlgebraProvider.cs | 32 ++++- 13 files changed, 383 insertions(+), 75 deletions(-) create mode 100644 src/Numerics/Providers/FourierTransform/FourierTransformControl.cs create mode 100644 src/Numerics/Providers/LinearAlgebra/LinearAlgebraControl.cs diff --git a/src/Numerics/Control.cs b/src/Numerics/Control.cs index 4045a7d4..8cbae11c 100644 --- a/src/Numerics/Control.cs +++ b/src/Numerics/Control.cs @@ -27,10 +27,10 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using MathNet.Numerics.Providers.LinearAlgebra; using System; using System.Threading.Tasks; using MathNet.Numerics.Providers.FourierTransform; +using MathNet.Numerics.Providers.LinearAlgebra; namespace MathNet.Numerics { @@ -39,8 +39,6 @@ namespace MathNet.Numerics /// public static class Control { - const string EnvVarLAProvider = "MathNetNumericsLAProvider"; - static int _maxDegreeOfParallelism; static int _blockSize; static int _parallelizeOrder; @@ -68,54 +66,10 @@ namespace MathNet.Numerics TaskScheduler = TaskScheduler.Default; } - private static void InitializeDefaultProviders() - { - lock (_staticLock) - { - if (_linearAlgebraProvider == null || _fourierTransformProvider == null) - { -#if NATIVE - try - { - var value = Environment.GetEnvironmentVariable(EnvVarLAProvider); - switch (value != null ? value.ToUpperInvariant() : string.Empty) - { - case "MKL": - UseNativeMKL(); - break; - - case "CUDA": - UseNativeCUDA(); - break; - - case "OPENBLAS": - UseNativeOpenBLAS(); - break; - - default: - if (!TryUseNative()) - { - UseManaged(); - } - break; - } - } - catch - { - // We don't care about any failures here at all (because "auto") - UseManaged(); - } -#else - UseManaged(); -#endif - } - } - } - public static void UseManaged() { - LinearAlgebraProvider = new ManagedLinearAlgebraProvider(); - FourierTransformProvider = new ManagedFourierTransformProvider(); + LinearAlgebraControl.UseManaged(); + FourierTransformControl.UseManaged(); } #if NATIVE @@ -126,8 +80,8 @@ namespace MathNet.Numerics /// public static void UseNativeMKL() { - LinearAlgebraProvider = new Providers.LinearAlgebra.Mkl.MklLinearAlgebraProvider(); - FourierTransformProvider = new Providers.FourierTransform.Mkl.MklFourierTransformProvider(); + LinearAlgebraControl.UseNativeMKL(); + FourierTransformControl.UseNativeMKL(); } /// @@ -141,8 +95,11 @@ namespace MathNet.Numerics Providers.LinearAlgebra.Mkl.MklPrecision precision = Providers.LinearAlgebra.Mkl.MklPrecision.Double, Providers.LinearAlgebra.Mkl.MklAccuracy accuracy = Providers.LinearAlgebra.Mkl.MklAccuracy.High) { - LinearAlgebraProvider = new Providers.LinearAlgebra.Mkl.MklLinearAlgebraProvider(consistency, precision, accuracy); - FourierTransformProvider = new Providers.FourierTransform.Mkl.MklFourierTransformProvider(); + LinearAlgebraControl.UseNativeMKL( + (Providers.Common.Mkl.MklConsistency)consistency, + (Providers.Common.Mkl.MklPrecision)precision, + (Providers.Common.Mkl.MklAccuracy)accuracy); + FourierTransformControl.UseNativeMKL(); } /// @@ -155,8 +112,8 @@ namespace MathNet.Numerics Providers.Common.Mkl.MklPrecision precision = Providers.Common.Mkl.MklPrecision.Double, Providers.Common.Mkl.MklAccuracy accuracy = Providers.Common.Mkl.MklAccuracy.High) { - LinearAlgebraProvider = new Providers.LinearAlgebra.Mkl.MklLinearAlgebraProvider(consistency, precision, accuracy); - FourierTransformProvider = new Providers.FourierTransform.Mkl.MklFourierTransformProvider(); + LinearAlgebraControl.UseNativeMKL(consistency, precision, accuracy); + FourierTransformControl.UseNativeMKL(); } /// @@ -177,11 +134,7 @@ namespace MathNet.Numerics /// public static void UseNativeCUDA() { - LinearAlgebraProvider = new Providers.LinearAlgebra.Cuda.CudaLinearAlgebraProvider(); - if (_fourierTransformProvider == null) - { - FourierTransformProvider = new ManagedFourierTransformProvider(); - } + LinearAlgebraControl.UseNativeCUDA(); } /// @@ -202,11 +155,7 @@ namespace MathNet.Numerics /// public static void UseNativeOpenBLAS() { - LinearAlgebraProvider = new Providers.LinearAlgebra.OpenBlas.OpenBlasLinearAlgebraProvider(); - if (_fourierTransformProvider == null) - { - FourierTransformProvider = new ManagedFourierTransformProvider(); - } + LinearAlgebraControl.UseNativeOpenBLAS(); } /// @@ -296,7 +245,15 @@ namespace MathNet.Numerics get { if (_linearAlgebraProvider == null) - InitializeDefaultProviders(); + { + lock (_staticLock) + { + if (_linearAlgebraProvider == null) + { + LinearAlgebraControl.UseDefault(); + } + } + } return _linearAlgebraProvider; } @@ -318,7 +275,15 @@ namespace MathNet.Numerics get { if (_fourierTransformProvider == null) - InitializeDefaultProviders(); + { + lock (_staticLock) + { + if (_fourierTransformProvider == null) + { + FourierTransformControl.UseDefault(); + } + } + } return _fourierTransformProvider; } diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 842ae7c9..1c546f4c 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -175,6 +175,7 @@ + @@ -190,6 +191,7 @@ + diff --git a/src/Numerics/Providers/Common/Mkl/MklProvider.cs b/src/Numerics/Providers/Common/Mkl/MklProvider.cs index 9a0ee156..f9879b83 100644 --- a/src/Numerics/Providers/Common/Mkl/MklProvider.cs +++ b/src/Numerics/Providers/Common/Mkl/MklProvider.cs @@ -40,23 +40,41 @@ namespace MathNet.Numerics.Providers.Common.Mkl static bool _nativeX64; static bool _nativeIA64; + public static bool IsAvailable(int minRevision) + { + try + { + if (!NativeProviderLoader.TryLoad(SafeNativeMethods.DllName)) + { + return false; + } + + int a = SafeNativeMethods.query_capability(0); + int b = SafeNativeMethods.query_capability(1); + int nativeRevision = SafeNativeMethods.query_capability((int)ProviderConfig.Revision); + return a == 0 && b == -1 && nativeRevision >= minRevision; + } + catch + { + return false; + } + } + public static void Load(int minRevision) { int a, b; try { - // Load the native library NativeProviderLoader.TryLoad(SafeNativeMethods.DllName); a = SafeNativeMethods.query_capability(0); b = SafeNativeMethods.query_capability(1); + _nativeRevision = SafeNativeMethods.query_capability((int)ProviderConfig.Revision); _nativeX86 = SafeNativeMethods.query_capability((int)ProviderPlatform.x86) > 0; _nativeX64 = SafeNativeMethods.query_capability((int)ProviderPlatform.x64) > 0; _nativeIA64 = SafeNativeMethods.query_capability((int)ProviderPlatform.ia64) > 0; - _nativeRevision = SafeNativeMethods.query_capability((int)ProviderConfig.Revision); - _mklVersion = new Version( SafeNativeMethods.query_capability((int)ProviderConfig.MklMajorVersion), SafeNativeMethods.query_capability((int)ProviderConfig.MklMinorVersion), diff --git a/src/Numerics/Providers/FourierTransform/FourierTransformControl.cs b/src/Numerics/Providers/FourierTransform/FourierTransformControl.cs new file mode 100644 index 00000000..d989432a --- /dev/null +++ b/src/Numerics/Providers/FourierTransform/FourierTransformControl.cs @@ -0,0 +1,97 @@ +// +// 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; + +namespace MathNet.Numerics.Providers.FourierTransform +{ + internal static class FourierTransformControl + { + const string EnvVarFFTProvider = "MathNetNumericsFFTProvider"; + + public static void UseManaged() + { + Control.FourierTransformProvider = new ManagedFourierTransformProvider(); + } + +#if NATIVE + public static void UseNativeMKL() + { + Control.FourierTransformProvider = new Mkl.MklFourierTransformProvider(); + } +#endif + + public static bool TryUse(IFourierTransformProvider provider) + { + try + { + if (!provider.IsAvailable()) + { + return false; + } + + Control.FourierTransformProvider = provider; + return true; + } + catch + { + // intentionally swallow exceptions here - use the explicit variants if you're interested in why + return false; + } + } + + public static void UseBest() + { +#if NATIVE + if (!TryUse(new Mkl.MklFourierTransformProvider())) + { + UseManaged(); + } +#else + UseManaged(); +#endif + } + + public static void UseDefault() + { + var value = Environment.GetEnvironmentVariable(EnvVarFFTProvider); + switch (value != null ? value.ToUpperInvariant() : string.Empty) + { +#if NATIVE + case "MKL": + UseNativeMKL(); + break; +#endif + default: + UseBest(); + break; + } + } + } +} diff --git a/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs index 2e4dea88..ec44d664 100644 --- a/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs @@ -43,6 +43,12 @@ namespace MathNet.Numerics.Providers.FourierTransform public interface IFourierTransformProvider { + /// + /// Try to find out whether the provider is available, at least in principle. + /// Verification may still fail if available, but it will certainly fail if unavailable. + /// + bool IsAvailable(); + /// /// Initialize and verify that the provided is indeed available. If not, fall back to alternatives like the managed provider /// diff --git a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs index f20ece67..eb2bb5c3 100644 --- a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs @@ -37,6 +37,18 @@ namespace MathNet.Numerics.Providers.FourierTransform public class ManagedFourierTransformProvider : IFourierTransformProvider { + /// + /// Try to find out whether the provider is available, at least in principle. + /// Verification may still fail if available, but it will certainly fail if unavailable. + /// + public virtual bool IsAvailable() + { + return true; + } + + /// + /// Initialize and verify that the provided is indeed available. If not, fall back to alternatives like the managed provider + /// public virtual void InitializeVerify() { } diff --git a/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs index afb9ac74..cf033d09 100644 --- a/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs @@ -36,6 +36,18 @@ namespace MathNet.Numerics.Providers.FourierTransform.Mkl { public class MklFourierTransformProvider : IFourierTransformProvider { + /// + /// Try to find out whether the provider is available, at least in principle. + /// Verification may still fail if available, but it will certainly fail if unavailable. + /// + public bool IsAvailable() + { + return MklProvider.IsAvailable(minRevision: 11); + } + + /// + /// Initialize and verify that the provided is indeed available. If not, fall back to alternatives like the managed provider + /// public void InitializeVerify() { MklProvider.Load(minRevision: 11); diff --git a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.cs index 6c814870..b1c3de28 100644 --- a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.cs +++ b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2015 Math.NET +// 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 @@ -46,6 +46,30 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda IntPtr _blasHandle; IntPtr _solverHandle; + /// + /// Try to find out whether the provider is available, at least in principle. + /// Verification may still fail if available, but it will certainly fail if unavailable. + /// + public override bool IsAvailable() + { + try + { + if (!NativeProviderLoader.TryLoad(SafeNativeMethods.DllName)) + { + return false; + } + + int a = SafeNativeMethods.query_capability(0); + int b = SafeNativeMethods.query_capability(1); + int nativeRevision = SafeNativeMethods.query_capability((int)ProviderConfig.Revision); + return a == 0 && b == -1 && nativeRevision >= 1; + } + catch + { + return false; + } + } + /// /// Initialize and verify that the provided is indeed available. /// If calling this method fails, consider to fall back to alternatives like the managed provider. @@ -60,12 +84,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda a = SafeNativeMethods.query_capability(0); b = SafeNativeMethods.query_capability(1); + _nativeRevision = SafeNativeMethods.query_capability((int)ProviderConfig.Revision); _nativeIX86 = SafeNativeMethods.query_capability((int)ProviderPlatform.x86) > 0; _nativeX64 = SafeNativeMethods.query_capability((int)ProviderPlatform.x64) > 0; _nativeIA64 = SafeNativeMethods.query_capability((int)ProviderPlatform.ia64) > 0; - _nativeRevision = SafeNativeMethods.query_capability((int)ProviderConfig.Revision); linearAlgebra = SafeNativeMethods.query_capability((int)ProviderCapability.LinearAlgebra); } catch (DllNotFoundException e) diff --git a/src/Numerics/Providers/LinearAlgebra/ILinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/ILinearAlgebraProvider.cs index 2298581a..1717dbee 100644 --- a/src/Numerics/Providers/LinearAlgebra/ILinearAlgebraProvider.cs +++ b/src/Numerics/Providers/LinearAlgebra/ILinearAlgebraProvider.cs @@ -93,6 +93,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra ILinearAlgebraProvider, ILinearAlgebraProvider { + /// + /// Try to find out whether the provider is available, at least in principle. + /// Verification may still fail if available, but it will certainly fail if unavailable. + /// + bool IsAvailable(); + /// /// Initialize and verify that the provided is indeed available. If not, fall back to alternatives like the managed provider /// diff --git a/src/Numerics/Providers/LinearAlgebra/LinearAlgebraControl.cs b/src/Numerics/Providers/LinearAlgebra/LinearAlgebraControl.cs new file mode 100644 index 00000000..a4cf8882 --- /dev/null +++ b/src/Numerics/Providers/LinearAlgebra/LinearAlgebraControl.cs @@ -0,0 +1,120 @@ +// +// 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; + +namespace MathNet.Numerics.Providers.LinearAlgebra +{ + internal static class LinearAlgebraControl + { + const string EnvVarLAProvider = "MathNetNumericsLAProvider"; + + public static void UseManaged() + { + Control.LinearAlgebraProvider = new ManagedLinearAlgebraProvider(); + } + +#if NATIVE + public static void UseNativeMKL( + Common.Mkl.MklConsistency consistency = Common.Mkl.MklConsistency.Auto, + Common.Mkl.MklPrecision precision = Common.Mkl.MklPrecision.Double, + Common.Mkl.MklAccuracy accuracy = Common.Mkl.MklAccuracy.High) + { + Control.LinearAlgebraProvider = new Mkl.MklLinearAlgebraProvider(consistency, precision, accuracy); + } + + public static void UseNativeCUDA() + { + Control.LinearAlgebraProvider = new Cuda.CudaLinearAlgebraProvider(); + } + + public static void UseNativeOpenBLAS() + { + Control.LinearAlgebraProvider = new OpenBlas.OpenBlasLinearAlgebraProvider(); + } +#endif + + public static bool TryUse(ILinearAlgebraProvider provider) + { + try + { + if (!provider.IsAvailable()) + { + return false; + } + + Control.LinearAlgebraProvider = provider; + return true; + } + catch + { + // intentionally swallow exceptions here - use the explicit variants if you're interested in why + return false; + } + } + + public static void UseBest() + { +#if NATIVE + if (!(TryUse(new Cuda.CudaLinearAlgebraProvider()) + || TryUse(new Mkl.MklLinearAlgebraProvider()) + || TryUse(new OpenBlas.OpenBlasLinearAlgebraProvider()))) + { + UseManaged(); + } +#else + UseManaged(); +#endif + } + + public static void UseDefault() + { + var value = Environment.GetEnvironmentVariable(EnvVarLAProvider); + switch (value != null ? value.ToUpperInvariant() : string.Empty) + { +#if NATIVE + case "MKL": + UseNativeMKL(); + break; + + case "CUDA": + UseNativeCUDA(); + break; + + case "OPENBLAS": + UseNativeOpenBLAS(); + break; +#endif + default: + UseBest(); + break; + } + } + } +} diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.cs index ac7a3f88..6fbc98a8 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.cs @@ -34,6 +34,15 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// public partial class ManagedLinearAlgebraProvider : ILinearAlgebraProvider { + /// + /// Try to find out whether the provider is available, at least in principle. + /// Verification may still fail if available, but it will certainly fail if unavailable. + /// + public virtual bool IsAvailable() + { + return true; + } + /// /// Initialize and verify that the provided is indeed available. If not, fall back to alternatives like the managed provider /// diff --git a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs index 92f7a316..f6d1c34a 100644 --- a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs +++ b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2015 Math.NET +// 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 @@ -133,6 +133,15 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl _accuracy = Common.Mkl.MklAccuracy.High; } + /// + /// Try to find out whether the provider is available, at least in principle. + /// Verification may still fail if available, but it will certainly fail if unavailable. + /// + public override bool IsAvailable() + { + return MklProvider.IsAvailable(minRevision: 4); + } + /// /// Initialize and verify that the provided is indeed available. /// If calling this method fails, consider to fall back to alternatives like the managed provider. diff --git a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.cs index 433243cf..18b4cc52 100644 --- a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.cs +++ b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2015 Math.NET +// 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 @@ -63,6 +63,34 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas bool _nativeIA64; bool _nativeARM; + /// + /// Try to find out whether the provider is available, at least in principle. + /// Verification may still fail if available, but it will certainly fail if unavailable. + /// + public override bool IsAvailable() + { + try + { + if (!NativeProviderLoader.TryLoad(SafeNativeMethods.DllName)) + { + return false; + } + + int a = SafeNativeMethods.query_capability(0); + int b = SafeNativeMethods.query_capability(1); + int nativeRevision = SafeNativeMethods.query_capability((int)ProviderConfig.Revision); + return a == 0 && b == -1 && nativeRevision >= 1; + } + catch + { + return false; + } + } + + /// + /// Initialize and verify that the provided is indeed available. + /// If not, fall back to alternatives like the managed provider + /// public override void InitializeVerify() { int a, b, linearAlgebra; @@ -73,13 +101,13 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas a = SafeNativeMethods.query_capability(0); b = SafeNativeMethods.query_capability(1); + _nativeRevision = SafeNativeMethods.query_capability((int)ProviderConfig.Revision); _nativeIX86 = SafeNativeMethods.query_capability((int)ProviderPlatform.x86) > 0; _nativeX64 = SafeNativeMethods.query_capability((int)ProviderPlatform.x64) > 0; _nativeIA64 = SafeNativeMethods.query_capability((int)ProviderPlatform.ia64) > 0; _nativeARM = SafeNativeMethods.query_capability((int)ProviderPlatform.arm) > 0; - _nativeRevision = SafeNativeMethods.query_capability((int)ProviderConfig.Revision); linearAlgebra = SafeNativeMethods.query_capability((int)ProviderCapability.LinearAlgebra); } catch (DllNotFoundException e) From 25282489f357e037f80188718945efed5ce6f18b Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 8 Oct 2016 15:06:04 +0200 Subject: [PATCH 13/17] FFT: cleanup Cuda and OpenBlas providers similar to MKL refactoring --- src/NativeProviders/CUDA/lapack.cpp | 8 +- src/NativeProviders/Common/WindowsDLL.cpp | 2 +- src/NativeProviders/Common/lapack.cpp | 14 +- src/Numerics/Numerics.csproj | 11 +- .../Providers/Common/Cuda/CudaProvider.cs | 107 ++++++++++++++++ .../Cuda/CudaProviderCapabilities.cs | 9 +- .../Cuda/SafeNativeMethods.cs | 5 +- .../Common/OpenBlas/OpenBlasProvider.cs | 121 ++++++++++++++++++ .../OpenBlas/OpenBlasProviderCapabilities.cs | 9 +- .../OpenBlas/SafeNativeMethods.cs | 4 +- .../ManagedFourierTransformProvider.cs | 5 + .../Cuda/CudaLinearAlgebraProvider.Complex.cs | 3 +- .../CudaLinearAlgebraProvider.Complex32.cs | 5 +- .../Cuda/CudaLinearAlgebraProvider.Double.cs | 5 +- .../Cuda/CudaLinearAlgebraProvider.Single.cs | 5 +- .../Cuda/CudaLinearAlgebraProvider.cs | 62 ++------- .../Mkl/MklLinearAlgebraProvider.cs | 1 - .../OpenBlasLinearAlgebraProvider.Complex.cs | 3 +- ...OpenBlasLinearAlgebraProvider.Complex32.cs | 3 +- .../OpenBlasLinearAlgebraProvider.Double.cs | 3 +- .../OpenBlasLinearAlgebraProvider.Single.cs | 5 +- .../OpenBlas/OpenBlasLinearAlgebraProvider.cs | 70 ++-------- 22 files changed, 304 insertions(+), 156 deletions(-) create mode 100644 src/Numerics/Providers/Common/Cuda/CudaProvider.cs rename src/Numerics/Providers/{LinearAlgebra => Common}/Cuda/CudaProviderCapabilities.cs (90%) rename src/Numerics/Providers/{LinearAlgebra => Common}/Cuda/SafeNativeMethods.cs (99%) create mode 100644 src/Numerics/Providers/Common/OpenBlas/OpenBlasProvider.cs rename src/Numerics/Providers/{LinearAlgebra => Common}/OpenBlas/OpenBlasProviderCapabilities.cs (90%) rename src/Numerics/Providers/{LinearAlgebra => Common}/OpenBlas/SafeNativeMethods.cs (99%) diff --git a/src/NativeProviders/CUDA/lapack.cpp b/src/NativeProviders/CUDA/lapack.cpp index a7bbbfbd..a4f27923 100644 --- a/src/NativeProviders/CUDA/lapack.cpp +++ b/src/NativeProviders/CUDA/lapack.cpp @@ -65,7 +65,7 @@ inline int lu_inverse(cusolverDnHandle_t solverHandle, cublasHandle_t blasHandle getrf(solverHandle, n, n, d_A, n, work, d_I, d_info); cudaMemcpy(&info, d_info, sizeof(int), cudaMemcpyDeviceToHost); - + cudaFree(work); if (info != 0) @@ -134,7 +134,7 @@ inline int lu_inverse_factored(cublasHandle_t blasHandle, int n, T a[], int ipiv getri(blasHandle, n, d_Aarray, n, d_I, d_Carray, n, d_info, 1); cudaMemcpy(&info, d_info, sizeof(int), cudaMemcpyDeviceToHost); - cublasGetMatrix(n, n, sizeof(T), d_C, n, a, n); + cublasGetMatrix(n, n, sizeof(T), d_C, n, a, n); cublasGetVector(n, sizeof(int), d_I, 1, ipiv, 1); shift_ipiv_down(n, ipiv); @@ -299,7 +299,7 @@ inline int cholesky_solve(cusolverDnHandle_t solverHandle, int n, int nrhs, T a[ cudaMemcpy(&info, d_info, sizeof(int), cudaMemcpyDeviceToHost); cudaFree(work); - + if (info != 0) { cudaFree(d_A); @@ -425,7 +425,7 @@ inline int cholesky_solve_factored(cusolverDnHandle_t solverHandle, int n, int n // int info = 0; // ormqr(solverHandle, &side, &tran, &m, &bn, &n, r, &m, tau, clone_b, &m, work, &len, &info); // trsm(blasHandle, CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, 1.0, r, m, clone_b, m); -// +// // copyBtoX(m, n, bn, clone_b, x); // delete[] clone_b; // return info; diff --git a/src/NativeProviders/Common/WindowsDLL.cpp b/src/NativeProviders/Common/WindowsDLL.cpp index 83c4c1e1..7bafd83b 100644 --- a/src/NativeProviders/Common/WindowsDLL.cpp +++ b/src/NativeProviders/Common/WindowsDLL.cpp @@ -2,7 +2,7 @@ #define VC_EXTRALEAN #endif #include "windows.h" - + BOOL APIENTRY DllMain( HANDLE, DWORD, LPVOID ){ return TRUE; } diff --git a/src/NativeProviders/Common/lapack.cpp b/src/NativeProviders/Common/lapack.cpp index 4375f1ef..de6abbe6 100644 --- a/src/NativeProviders/Common/lapack.cpp +++ b/src/NativeProviders/Common/lapack.cpp @@ -16,12 +16,12 @@ inline lapack_int lu_factor(lapack_int m, T a[], lapack_int ipiv[], GETRF getrf) template inline lapack_int lu_inverse(lapack_int n, T a[], GETRF getrf, GETRI getri) { - try + try { auto ipiv = array_new(n); auto info = getrf(LAPACK_COL_MAJOR, n, n, a, n, ipiv.get()); - if (info != 0) + if (info != 0) { return info; } @@ -29,7 +29,7 @@ inline lapack_int lu_inverse(lapack_int n, T a[], GETRF getrf, GETRI getri) info = getri(LAPACK_COL_MAJOR, n, a, n, ipiv.get()); return info; } - catch (std::bad_alloc&) + catch (std::bad_alloc&) { return INSUFFICIENT_MEMORY; } @@ -56,20 +56,20 @@ inline lapack_int lu_solve_factored(lapack_int n, lapack_int nrhs, T a[], lapack template inline lapack_int lu_solve(lapack_int n, lapack_int nrhs, T a[], T b[], GETRF getrf, GETRS getrs) { - try + try { auto clone = array_clone(n * n, a); auto ipiv = array_new(n); auto info = getrf(LAPACK_COL_MAJOR, n, n, clone.get(), n, ipiv.get()); - if (info != 0) + if (info != 0) { return info; } return getrs(LAPACK_COL_MAJOR, 'N', n, nrhs, clone.get(), n, ipiv.get(), b, n); } - catch (std::bad_alloc&) + catch (std::bad_alloc&) { return INSUFFICIENT_MEMORY; } @@ -264,7 +264,7 @@ inline lapack_int svd_factor(bool compute_vectors, lapack_int m, lapack_int n, T template inline lapack_int complex_svd_factor(bool compute_vectors, lapack_int m, lapack_int n, T a[], T s[], T u[], T v[], GESVD gesvd) { - try + try { auto dim_s = std::min(m, n); auto s_local = array_new(dim_s); diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 1c546f4c..8a54e767 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -172,8 +172,10 @@ True Resources.resx + + @@ -188,8 +190,8 @@ - - + + @@ -197,7 +199,7 @@ - + @@ -227,7 +229,7 @@ - + @@ -483,5 +485,6 @@ Resources.Designer.cs + \ No newline at end of file diff --git a/src/Numerics/Providers/Common/Cuda/CudaProvider.cs b/src/Numerics/Providers/Common/Cuda/CudaProvider.cs new file mode 100644 index 00000000..53d2ded5 --- /dev/null +++ b/src/Numerics/Providers/Common/Cuda/CudaProvider.cs @@ -0,0 +1,107 @@ +// +// 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 System.Collections.Generic; + +namespace MathNet.Numerics.Providers.Common.Cuda +{ + internal static class CudaProvider + { + static int _nativeRevision; + static bool _nativeX86; + static bool _nativeX64; + static bool _nativeIA64; + + public static bool IsAvailable(int minRevision) + { + try + { + if (!NativeProviderLoader.TryLoad(SafeNativeMethods.DllName)) + { + return false; + } + + int a = SafeNativeMethods.query_capability(0); + int b = SafeNativeMethods.query_capability(1); + int nativeRevision = SafeNativeMethods.query_capability((int)ProviderConfig.Revision); + return a == 0 && b == -1 && nativeRevision >= minRevision; + } + catch + { + return false; + } + } + + public static void Load(int minRevision) + { + int a, b; + try + { + NativeProviderLoader.TryLoad(SafeNativeMethods.DllName); + + a = SafeNativeMethods.query_capability(0); + b = SafeNativeMethods.query_capability(1); + _nativeRevision = SafeNativeMethods.query_capability((int)ProviderConfig.Revision); + + _nativeX86 = SafeNativeMethods.query_capability((int)ProviderPlatform.x86) > 0; + _nativeX64 = SafeNativeMethods.query_capability((int)ProviderPlatform.x64) > 0; + _nativeIA64 = SafeNativeMethods.query_capability((int)ProviderPlatform.ia64) > 0; + } + catch (DllNotFoundException e) + { + throw new NotSupportedException("Cuda Native Provider not found.", e); + } + catch (BadImageFormatException e) + { + throw new NotSupportedException("Cuda Native Provider found but failed to load. Please verify that the platform matches (x64 vs x32, Windows vs Linux).", e); + } + catch (EntryPointNotFoundException e) + { + throw new NotSupportedException("Cuda Native Provider does not support capability querying and is therefore not compatible. Consider upgrading to a newer version.", e); + } + + if (a != 0 || b != -1 || _nativeRevision < minRevision) + { + throw new NotSupportedException("Cuda Native Provider too old. Consider upgrading to a newer version."); + } + } + + public static string Describe() + { + var parts = new List(); + if (_nativeX86) parts.Add("x86"); + if (_nativeX64) parts.Add("x64"); + if (_nativeIA64) parts.Add("IA64"); + parts.Add("revision " + _nativeRevision); + + return string.Concat("Nvidia CUDA (", string.Join("; ", parts), ")"); + } + } +} diff --git a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaProviderCapabilities.cs b/src/Numerics/Providers/Common/Cuda/CudaProviderCapabilities.cs similarity index 90% rename from src/Numerics/Providers/LinearAlgebra/Cuda/CudaProviderCapabilities.cs rename to src/Numerics/Providers/Common/Cuda/CudaProviderCapabilities.cs index f25d161d..90002973 100644 --- a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaProviderCapabilities.cs +++ b/src/Numerics/Providers/Common/Cuda/CudaProviderCapabilities.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2015 Math.NET +// 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 @@ -27,7 +27,7 @@ // OTHER DEALINGS IN THE SOFTWARE. // -namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda +namespace MathNet.Numerics.Providers.Common.Cuda { internal enum ProviderPlatform : int { @@ -43,6 +43,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda internal enum ProviderCapability : int { - LinearAlgebra = 128, + LinearAlgebraMajor = 128, + LinearAlgebraMinor = 129, } -} \ No newline at end of file +} diff --git a/src/Numerics/Providers/LinearAlgebra/Cuda/SafeNativeMethods.cs b/src/Numerics/Providers/Common/Cuda/SafeNativeMethods.cs similarity index 99% rename from src/Numerics/Providers/LinearAlgebra/Cuda/SafeNativeMethods.cs rename to src/Numerics/Providers/Common/Cuda/SafeNativeMethods.cs index ba58a483..d767e7b0 100644 --- a/src/Numerics/Providers/LinearAlgebra/Cuda/SafeNativeMethods.cs +++ b/src/Numerics/Providers/Common/Cuda/SafeNativeMethods.cs @@ -2,7 +2,7 @@ // Math.NET Numerics, part of the Math.NET Project // http://mathnet.opensourcedotnet.info // -// Copyright (c) 2009-2014 Math.NET +// 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 @@ -32,8 +32,9 @@ using System; using System.Numerics; using System.Runtime.InteropServices; using System.Security; +using MathNet.Numerics.Providers.LinearAlgebra; -namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda +namespace MathNet.Numerics.Providers.Common.Cuda { /// /// P/Invoke methods to the native math libraries. diff --git a/src/Numerics/Providers/Common/OpenBlas/OpenBlasProvider.cs b/src/Numerics/Providers/Common/OpenBlas/OpenBlasProvider.cs new file mode 100644 index 00000000..473ee263 --- /dev/null +++ b/src/Numerics/Providers/Common/OpenBlas/OpenBlasProvider.cs @@ -0,0 +1,121 @@ +// +// 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 System.Collections.Generic; + +namespace MathNet.Numerics.Providers.Common.OpenBlas +{ + internal static class OpenBlasProvider + { + static int _nativeRevision; + static bool _nativeX86; + static bool _nativeX64; + static bool _nativeIA64; + static bool _nativeARM; + + public static bool IsAvailable(int minRevision) + { + try + { + if (!NativeProviderLoader.TryLoad(SafeNativeMethods.DllName)) + { + return false; + } + + int a = SafeNativeMethods.query_capability(0); + int b = SafeNativeMethods.query_capability(1); + int nativeRevision = SafeNativeMethods.query_capability((int)ProviderConfig.Revision); + return a == 0 && b == -1 && nativeRevision >= minRevision; + } + catch + { + return false; + } + } + + public static void Load(int minRevision) + { + int a, b; + try + { + NativeProviderLoader.TryLoad(SafeNativeMethods.DllName); + + a = SafeNativeMethods.query_capability(0); + b = SafeNativeMethods.query_capability(1); + _nativeRevision = SafeNativeMethods.query_capability((int)ProviderConfig.Revision); + + _nativeX86 = SafeNativeMethods.query_capability((int)ProviderPlatform.x86) > 0; + _nativeX64 = SafeNativeMethods.query_capability((int)ProviderPlatform.x64) > 0; + _nativeIA64 = SafeNativeMethods.query_capability((int)ProviderPlatform.ia64) > 0; + _nativeARM = SafeNativeMethods.query_capability((int)ProviderPlatform.arm) > 0; + } + catch (DllNotFoundException e) + { + throw new NotSupportedException("OpenBLAS Native Provider not found.", e); + } + catch (BadImageFormatException e) + { + throw new NotSupportedException("OpenBLAS Native Provider found but failed to load. Please verify that the platform matches (x64 vs x32, Windows vs Linux).", e); + } + catch (EntryPointNotFoundException e) + { + throw new NotSupportedException("OpenBLAS Native Provider does not support capability querying and is therefore not compatible. Consider upgrading to a newer version.", e); + } + + if (a != 0 || b != -1 || _nativeRevision < minRevision) + { + throw new NotSupportedException("OpenBLAS Native Provider too old. Consider upgrading to a newer version."); + } + + ConfigureThreading(); + } + + static void ConfigureThreading() + { + // set threading settings, if supported + if (SafeNativeMethods.query_capability((int)ProviderConfig.Threading) > 0) + { + SafeNativeMethods.set_max_threads(Control.MaxDegreeOfParallelism); + } + } + + public static string Describe() + { + var parts = new List(); + if (_nativeX86) parts.Add("x86"); + if (_nativeX64) parts.Add("x64"); + if (_nativeIA64) parts.Add("IA64"); + if (_nativeARM) parts.Add("ARM"); + parts.Add("revision " + _nativeRevision); + + return string.Concat("OpenBLAS (", string.Join("; ", parts), ")"); + } + } +} diff --git a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasProviderCapabilities.cs b/src/Numerics/Providers/Common/OpenBlas/OpenBlasProviderCapabilities.cs similarity index 90% rename from src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasProviderCapabilities.cs rename to src/Numerics/Providers/Common/OpenBlas/OpenBlasProviderCapabilities.cs index b86c52db..858c710e 100644 --- a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasProviderCapabilities.cs +++ b/src/Numerics/Providers/Common/OpenBlas/OpenBlasProviderCapabilities.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2015 Math.NET +// 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 @@ -27,7 +27,7 @@ // OTHER DEALINGS IN THE SOFTWARE. // -namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas +namespace MathNet.Numerics.Providers.Common.OpenBlas { internal enum ProviderPlatform : int { @@ -45,6 +45,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas internal enum ProviderCapability : int { - LinearAlgebra = 128, + LinearAlgebraMajor = 128, + LinearAlgebraMinor = 129 } -} \ No newline at end of file +} diff --git a/src/Numerics/Providers/LinearAlgebra/OpenBlas/SafeNativeMethods.cs b/src/Numerics/Providers/Common/OpenBlas/SafeNativeMethods.cs similarity index 99% rename from src/Numerics/Providers/LinearAlgebra/OpenBlas/SafeNativeMethods.cs rename to src/Numerics/Providers/Common/OpenBlas/SafeNativeMethods.cs index fc024d90..aa62c2cc 100644 --- a/src/Numerics/Providers/LinearAlgebra/OpenBlas/SafeNativeMethods.cs +++ b/src/Numerics/Providers/Common/OpenBlas/SafeNativeMethods.cs @@ -31,8 +31,10 @@ using System.Numerics; using System.Runtime.InteropServices; using System.Security; +using MathNet.Numerics.Providers.LinearAlgebra; +using MathNet.Numerics.Providers.LinearAlgebra.OpenBlas; -namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas +namespace MathNet.Numerics.Providers.Common.OpenBlas { /// /// P/Invoke methods to the native math libraries. diff --git a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs index eb2bb5c3..58d42e3d 100644 --- a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs @@ -53,6 +53,11 @@ namespace MathNet.Numerics.Providers.FourierTransform { } + public override string ToString() + { + return "Managed"; + } + public virtual void ForwardInplace(Complex[] complex, FourierTransformScaling scaling) { Fourier.BluesteinForward(complex, Options(scaling)); diff --git a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex.cs b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex.cs index 65a4689b..28ae3741 100644 --- a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex.cs +++ b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex.cs @@ -33,6 +33,7 @@ using System; using System.Numerics; using System.Security; using MathNet.Numerics.Properties; +using MathNet.Numerics.Providers.Common.Cuda; namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda { @@ -569,7 +570,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda if (columnsA > rowsA || !computeVectors) // see remarks http://docs.nvidia.com/cuda/cusolver/index.html#cuds-lt-t-gt-gesvd base.SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt); else Solver(SafeNativeMethods.z_svd_factor(_solverHandle, computeVectors, rowsA, columnsA, a, s, u, vt)); - } + } } } diff --git a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex32.cs b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex32.cs index 5a2f8806..4e5674ff 100644 --- a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex32.cs +++ b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex32.cs @@ -32,6 +32,7 @@ using System; using System.Security; using MathNet.Numerics.Properties; +using MathNet.Numerics.Providers.Common.Cuda; namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda { @@ -466,7 +467,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda } Solver(SafeNativeMethods.c_cholesky_solve_factored(_solverHandle, orderA, columnsB, a, b)); - } + } /// /// Solves A*X=B for X using the singular value decomposition of A. @@ -568,7 +569,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda if (columnsA > rowsA || !computeVectors) // see remarks http://docs.nvidia.com/cuda/cusolver/index.html#cuds-lt-t-gt-gesvd base.SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt); else Solver(SafeNativeMethods.c_svd_factor(_solverHandle, computeVectors, rowsA, columnsA, a, s, u, vt)); - } + } } } diff --git a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Double.cs b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Double.cs index 5b3cc78c..1a60985b 100644 --- a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Double.cs +++ b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Double.cs @@ -32,6 +32,7 @@ using System; using System.Security; using MathNet.Numerics.Properties; +using MathNet.Numerics.Providers.Common.Cuda; namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda { @@ -466,7 +467,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda } Solver(SafeNativeMethods.d_cholesky_solve_factored(_solverHandle, orderA, columnsB, a, b)); - } + } /// /// Solves A*X=B for X using the singular value decomposition of A. @@ -568,7 +569,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda if (columnsA > rowsA || !computeVectors) // see remarks http://docs.nvidia.com/cuda/cusolver/index.html#cuds-lt-t-gt-gesvd base.SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt); else Solver (SafeNativeMethods.d_svd_factor(_solverHandle, computeVectors, rowsA, columnsA, a, s, u, vt)); - } + } } } diff --git a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Single.cs b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Single.cs index bfcbf73e..fa109558 100644 --- a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Single.cs +++ b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Single.cs @@ -32,6 +32,7 @@ using System; using System.Security; using MathNet.Numerics.Properties; +using MathNet.Numerics.Providers.Common.Cuda; namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda { @@ -466,7 +467,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda } Solver(SafeNativeMethods.s_cholesky_solve_factored(_solverHandle, orderA, columnsB, a, b)); - } + } /// /// Solves A*X=B for X using the singular value decomposition of A. @@ -568,7 +569,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda if (columnsA > rowsA || !computeVectors) // see remarks http://docs.nvidia.com/cuda/cusolver/index.html#cuds-lt-t-gt-gesvd base.SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt); else Solver(SafeNativeMethods.s_svd_factor(_solverHandle, computeVectors, rowsA, columnsA, a, s, u, vt)); - } + } } } diff --git a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.cs index b1c3de28..6ae60291 100644 --- a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.cs +++ b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.cs @@ -28,7 +28,7 @@ // using System; -using MathNet.Numerics.Providers.Common; +using MathNet.Numerics.Providers.Common.Cuda; #if NATIVE @@ -39,10 +39,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda /// public partial class CudaLinearAlgebraProvider : ManagedLinearAlgebraProvider, IDisposable { - int _nativeRevision; - bool _nativeIX86; - bool _nativeX64; - bool _nativeIA64; IntPtr _blasHandle; IntPtr _solverHandle; @@ -52,22 +48,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda /// public override bool IsAvailable() { - try - { - if (!NativeProviderLoader.TryLoad(SafeNativeMethods.DllName)) - { - return false; - } - - int a = SafeNativeMethods.query_capability(0); - int b = SafeNativeMethods.query_capability(1); - int nativeRevision = SafeNativeMethods.query_capability((int)ProviderConfig.Revision); - return a == 0 && b == -1 && nativeRevision >= 1; - } - catch - { - return false; - } + return CudaProvider.IsAvailable(minRevision: 1); } /// @@ -76,38 +57,14 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda /// public override void InitializeVerify() { - int a, b, linearAlgebra; - try - { - // Load the native library - NativeProviderLoader.TryLoad(SafeNativeMethods.DllName); - - a = SafeNativeMethods.query_capability(0); - b = SafeNativeMethods.query_capability(1); - _nativeRevision = SafeNativeMethods.query_capability((int)ProviderConfig.Revision); + CudaProvider.Load(minRevision: 1); - _nativeIX86 = SafeNativeMethods.query_capability((int)ProviderPlatform.x86) > 0; - _nativeX64 = SafeNativeMethods.query_capability((int)ProviderPlatform.x64) > 0; - _nativeIA64 = SafeNativeMethods.query_capability((int)ProviderPlatform.ia64) > 0; + int linearAlgebra = SafeNativeMethods.query_capability((int)ProviderCapability.LinearAlgebraMajor); - linearAlgebra = SafeNativeMethods.query_capability((int)ProviderCapability.LinearAlgebra); - } - catch (DllNotFoundException e) + // we only support exactly one major version, since major version changes imply a breaking change. + if (linearAlgebra != 1) { - throw new NotSupportedException("Cuda Native Provider not found.", e); - } - catch (BadImageFormatException e) - { - throw new NotSupportedException("Cuda Native Provider found but failed to load. Please verify that the platform matches (x64 vs x32, Windows vs Linux).", e); - } - catch (EntryPointNotFoundException e) - { - throw new NotSupportedException("Cuda Native Provider does not support capability querying and is therefore not compatible. Consider upgrading to a newer version.", e); - } - - if (a != 0 || b != -1 || linearAlgebra <=0 || _nativeRevision < 1) - { - throw new NotSupportedException("Cuda Native Provider too old or not compatible. Consider upgrading to a newer version."); + throw new NotSupportedException(string.Format("Cuda Native Provider not compatible. Expecting linear algebra v1 but provider implements v{0}.", linearAlgebra)); } BLAS(SafeNativeMethods.createBLASHandle(ref _blasHandle)); @@ -200,12 +157,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda public override string ToString() { - return string.Format("Nvidia CUDA ({1}; revision {0})", - _nativeRevision, - _nativeIX86 ? "x86" : _nativeX64 ? "x64" : _nativeIA64 ? "IA64" : "unknown"); + return CudaProvider.Describe(); } - public void Dispose() { BLAS(SafeNativeMethods.destroyBLASHandle(_blasHandle)); diff --git a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs index f6d1c34a..1de85759 100644 --- a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs +++ b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs @@ -28,7 +28,6 @@ // using System; -using MathNet.Numerics.Providers.Common; using MathNet.Numerics.Providers.Common.Mkl; #if NATIVE diff --git a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex.cs b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex.cs index 7a87dac9..88c15144 100644 --- a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex.cs +++ b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex.cs @@ -34,6 +34,7 @@ using MathNet.Numerics.Properties; using System; using System.Numerics; using System.Security; +using MathNet.Numerics.Providers.Common.OpenBlas; namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas { @@ -357,7 +358,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas if (info > 0) { throw new SingularUMatrixException(info); - } + } } /// diff --git a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex32.cs b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex32.cs index d431c8c9..d49562cf 100644 --- a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex32.cs +++ b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex32.cs @@ -34,6 +34,7 @@ using MathNet.Numerics.Properties; using System; using System.Numerics; using System.Security; +using MathNet.Numerics.Providers.Common.OpenBlas; namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas { @@ -357,7 +358,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas if (info > 0) { throw new SingularUMatrixException(info); - } + } } /// diff --git a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Double.cs b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Double.cs index d4e6b117..4cd9888c 100644 --- a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Double.cs +++ b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Double.cs @@ -34,6 +34,7 @@ using MathNet.Numerics.Properties; using System; using System.Numerics; using System.Security; +using MathNet.Numerics.Providers.Common.OpenBlas; namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas { @@ -357,7 +358,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas if (info > 0) { throw new SingularUMatrixException(info); - } + } } /// diff --git a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Single.cs b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Single.cs index 7de51a33..2752606f 100644 --- a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Single.cs +++ b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Single.cs @@ -34,6 +34,7 @@ using MathNet.Numerics.Properties; using System; using System.Numerics; using System.Security; +using MathNet.Numerics.Providers.Common.OpenBlas; namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas { @@ -357,7 +358,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas if (info > 0) { throw new SingularUMatrixException(info); - } + } } /// @@ -811,7 +812,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas if (method == QRMethod.Full) { var info = SafeNativeMethods.s_qr_solve_factored(rowsA, columnsA, columnsB, r, b, tau, x); - + if (info == (int)NativeError.MemoryAllocation) { throw new MemoryAllocationException(); diff --git a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.cs index 18b4cc52..c7306b06 100644 --- a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.cs +++ b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.cs @@ -30,7 +30,7 @@ #if NATIVE using System; -using MathNet.Numerics.Providers.Common; +using MathNet.Numerics.Providers.Common.OpenBlas; namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas { @@ -57,34 +57,13 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas /// public partial class OpenBlasLinearAlgebraProvider : ManagedLinearAlgebraProvider { - int _nativeRevision; - bool _nativeIX86; - bool _nativeX64; - bool _nativeIA64; - bool _nativeARM; - /// /// Try to find out whether the provider is available, at least in principle. /// Verification may still fail if available, but it will certainly fail if unavailable. /// public override bool IsAvailable() { - try - { - if (!NativeProviderLoader.TryLoad(SafeNativeMethods.DllName)) - { - return false; - } - - int a = SafeNativeMethods.query_capability(0); - int b = SafeNativeMethods.query_capability(1); - int nativeRevision = SafeNativeMethods.query_capability((int)ProviderConfig.Revision); - return a == 0 && b == -1 && nativeRevision >= 1; - } - catch - { - return false; - } + return OpenBlasProvider.IsAvailable(minRevision: 1); } /// @@ -93,53 +72,20 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas /// public override void InitializeVerify() { - int a, b, linearAlgebra; - try - { - // Load the native library - NativeProviderLoader.TryLoad(SafeNativeMethods.DllName); - - a = SafeNativeMethods.query_capability(0); - b = SafeNativeMethods.query_capability(1); - _nativeRevision = SafeNativeMethods.query_capability((int)ProviderConfig.Revision); + OpenBlasProvider.Load(minRevision: 1); - _nativeIX86 = SafeNativeMethods.query_capability((int)ProviderPlatform.x86) > 0; - _nativeX64 = SafeNativeMethods.query_capability((int)ProviderPlatform.x64) > 0; - _nativeIA64 = SafeNativeMethods.query_capability((int)ProviderPlatform.ia64) > 0; - _nativeARM = SafeNativeMethods.query_capability((int)ProviderPlatform.arm) > 0; - - linearAlgebra = SafeNativeMethods.query_capability((int)ProviderCapability.LinearAlgebra); - } - catch (DllNotFoundException e) - { - throw new NotSupportedException("OpenBLAS Native Provider not found.", e); - } - catch (BadImageFormatException e) - { - throw new NotSupportedException("OpenBLAS Native Provider found but failed to load. Please verify that the platform matches (x64 vs x32, Windows vs Linux).", e); - } - catch (EntryPointNotFoundException e) - { - throw new NotSupportedException("OpenBLAS Native Provider does not support capability querying and is therefore not compatible. Consider upgrading to a newer version.", e); - } - - if (a != 0 || b != -1 || linearAlgebra <=0 || _nativeRevision < 1) - { - throw new NotSupportedException("OpenBLAS Native Provider too old or not compatible. Consider upgrading to a newer version."); - } + int linearAlgebra = SafeNativeMethods.query_capability((int)ProviderCapability.LinearAlgebraMajor); - // set threading settings, if supported - if (SafeNativeMethods.query_capability((int)ProviderConfig.Threading) > 0) + // we only support exactly one major version, since major version changes imply a breaking change. + if (linearAlgebra != 1) { - SafeNativeMethods.set_max_threads(Control.MaxDegreeOfParallelism); + throw new NotSupportedException(string.Format("OpenBLAS Native Provider not compatible. Expecting linear algebra v1 but provider implements v{0}.", linearAlgebra)); } } public override string ToString() { - return string.Format("OpenBLAS ({1}; revision {0})", - _nativeRevision, - _nativeIX86 ? "x86" : _nativeX64 ? "x64" : _nativeIA64 ? "IA64" : _nativeARM ? "ARM" : "unknown"); + return OpenBlasProvider.Describe(); } } } From c472e8d4c8740d6606a869b81053d6f33215ed93 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 8 Oct 2016 15:20:48 +0200 Subject: [PATCH 14/17] FFT: fix PCL builds --- src/Numerics/Providers/Common/Cuda/CudaProvider.cs | 4 ++++ .../Providers/Common/Cuda/CudaProviderCapabilities.cs | 4 ++++ src/Numerics/Providers/Common/Mkl/MklProvider.cs | 4 ++++ .../Providers/Common/Mkl/MklProviderCapabilities.cs | 4 ++++ src/Numerics/Providers/Common/NativeProviderLoader.cs | 6 +++--- .../Providers/Common/OpenBlas/OpenBlasProvider.cs | 4 ++++ .../Common/OpenBlas/OpenBlasProviderCapabilities.cs | 4 ++++ .../Providers/FourierTransform/FourierTransformControl.cs | 8 ++++++-- .../LinearAlgebra/Cuda/CudaLinearAlgebraProvider.cs | 4 ++-- .../Providers/LinearAlgebra/LinearAlgebraControl.cs | 6 ++++-- .../LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs | 4 ++-- 11 files changed, 41 insertions(+), 11 deletions(-) diff --git a/src/Numerics/Providers/Common/Cuda/CudaProvider.cs b/src/Numerics/Providers/Common/Cuda/CudaProvider.cs index 53d2ded5..14c21127 100644 --- a/src/Numerics/Providers/Common/Cuda/CudaProvider.cs +++ b/src/Numerics/Providers/Common/Cuda/CudaProvider.cs @@ -27,6 +27,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // +#if NATIVE + using System; using System.Collections.Generic; @@ -105,3 +107,5 @@ namespace MathNet.Numerics.Providers.Common.Cuda } } } + +#endif diff --git a/src/Numerics/Providers/Common/Cuda/CudaProviderCapabilities.cs b/src/Numerics/Providers/Common/Cuda/CudaProviderCapabilities.cs index 90002973..6856ab49 100644 --- a/src/Numerics/Providers/Common/Cuda/CudaProviderCapabilities.cs +++ b/src/Numerics/Providers/Common/Cuda/CudaProviderCapabilities.cs @@ -27,6 +27,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // +#if NATIVE + namespace MathNet.Numerics.Providers.Common.Cuda { internal enum ProviderPlatform : int @@ -47,3 +49,5 @@ namespace MathNet.Numerics.Providers.Common.Cuda LinearAlgebraMinor = 129, } } + +#endif diff --git a/src/Numerics/Providers/Common/Mkl/MklProvider.cs b/src/Numerics/Providers/Common/Mkl/MklProvider.cs index f9879b83..b2007a8d 100644 --- a/src/Numerics/Providers/Common/Mkl/MklProvider.cs +++ b/src/Numerics/Providers/Common/Mkl/MklProvider.cs @@ -27,6 +27,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // +#if NATIVE + using System; using System.Collections.Generic; @@ -256,3 +258,5 @@ namespace MathNet.Numerics.Providers.Common.Mkl } } } + +#endif diff --git a/src/Numerics/Providers/Common/Mkl/MklProviderCapabilities.cs b/src/Numerics/Providers/Common/Mkl/MklProviderCapabilities.cs index af0f971d..4e692bab 100644 --- a/src/Numerics/Providers/Common/Mkl/MklProviderCapabilities.cs +++ b/src/Numerics/Providers/Common/Mkl/MklProviderCapabilities.cs @@ -27,6 +27,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // +#if NATIVE + namespace MathNet.Numerics.Providers.Common.Mkl { internal enum ProviderPlatform : int @@ -55,3 +57,5 @@ namespace MathNet.Numerics.Providers.Common.Mkl FourierTransformMinor = 385 } } + +#endif diff --git a/src/Numerics/Providers/Common/NativeProviderLoader.cs b/src/Numerics/Providers/Common/NativeProviderLoader.cs index 0b4f19e4..51c446bf 100644 --- a/src/Numerics/Providers/Common/NativeProviderLoader.cs +++ b/src/Numerics/Providers/Common/NativeProviderLoader.cs @@ -27,6 +27,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // +#if NATIVE + using System; using System.Collections.Generic; using System.IO; @@ -35,8 +37,6 @@ using System.Runtime.InteropServices; using System.Security; using System.Threading; -#if NATIVE - namespace MathNet.Numerics.Providers.Common { /// @@ -243,4 +243,4 @@ namespace MathNet.Numerics.Providers.Common } } -#endif \ No newline at end of file +#endif diff --git a/src/Numerics/Providers/Common/OpenBlas/OpenBlasProvider.cs b/src/Numerics/Providers/Common/OpenBlas/OpenBlasProvider.cs index 473ee263..0261ee1b 100644 --- a/src/Numerics/Providers/Common/OpenBlas/OpenBlasProvider.cs +++ b/src/Numerics/Providers/Common/OpenBlas/OpenBlasProvider.cs @@ -27,6 +27,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // +#if NATIVE + using System; using System.Collections.Generic; @@ -119,3 +121,5 @@ namespace MathNet.Numerics.Providers.Common.OpenBlas } } } + +#endif diff --git a/src/Numerics/Providers/Common/OpenBlas/OpenBlasProviderCapabilities.cs b/src/Numerics/Providers/Common/OpenBlas/OpenBlasProviderCapabilities.cs index 858c710e..fdabbf03 100644 --- a/src/Numerics/Providers/Common/OpenBlas/OpenBlasProviderCapabilities.cs +++ b/src/Numerics/Providers/Common/OpenBlas/OpenBlasProviderCapabilities.cs @@ -27,6 +27,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // +#if NATIVE + namespace MathNet.Numerics.Providers.Common.OpenBlas { internal enum ProviderPlatform : int @@ -49,3 +51,5 @@ namespace MathNet.Numerics.Providers.Common.OpenBlas LinearAlgebraMinor = 129 } } + +#endif diff --git a/src/Numerics/Providers/FourierTransform/FourierTransformControl.cs b/src/Numerics/Providers/FourierTransform/FourierTransformControl.cs index d989432a..d553cd27 100644 --- a/src/Numerics/Providers/FourierTransform/FourierTransformControl.cs +++ b/src/Numerics/Providers/FourierTransform/FourierTransformControl.cs @@ -80,18 +80,22 @@ namespace MathNet.Numerics.Providers.FourierTransform public static void UseDefault() { +#if NATIVE var value = Environment.GetEnvironmentVariable(EnvVarFFTProvider); switch (value != null ? value.ToUpperInvariant() : string.Empty) { -#if NATIVE + case "MKL": UseNativeMKL(); break; -#endif + default: UseBest(); break; } +#else + UseManaged(); +#endif } } } diff --git a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.cs index 6ae60291..bfd0cd99 100644 --- a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.cs +++ b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.cs @@ -27,11 +27,11 @@ // OTHER DEALINGS IN THE SOFTWARE. // +#if NATIVE + using System; using MathNet.Numerics.Providers.Common.Cuda; -#if NATIVE - namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda { /// diff --git a/src/Numerics/Providers/LinearAlgebra/LinearAlgebraControl.cs b/src/Numerics/Providers/LinearAlgebra/LinearAlgebraControl.cs index a4cf8882..e3712c46 100644 --- a/src/Numerics/Providers/LinearAlgebra/LinearAlgebraControl.cs +++ b/src/Numerics/Providers/LinearAlgebra/LinearAlgebraControl.cs @@ -95,10 +95,10 @@ namespace MathNet.Numerics.Providers.LinearAlgebra public static void UseDefault() { +#if NATIVE var value = Environment.GetEnvironmentVariable(EnvVarLAProvider); switch (value != null ? value.ToUpperInvariant() : string.Empty) { -#if NATIVE case "MKL": UseNativeMKL(); break; @@ -110,11 +110,13 @@ namespace MathNet.Numerics.Providers.LinearAlgebra case "OPENBLAS": UseNativeOpenBLAS(); break; -#endif default: UseBest(); break; } +#else + UseManaged(); +#endif } } } diff --git a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs index 1de85759..46a48110 100644 --- a/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs +++ b/src/Numerics/Providers/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs @@ -27,11 +27,11 @@ // OTHER DEALINGS IN THE SOFTWARE. // +#if NATIVE + using System; using MathNet.Numerics.Providers.Common.Mkl; -#if NATIVE - namespace MathNet.Numerics.Providers.LinearAlgebra.Mkl { /// From 432bf493a5aae37888bf83beb3419d1652b6a764 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 8 Oct 2016 16:35:49 +0200 Subject: [PATCH 15/17] FFT: force 64bit integers with MKL_INT64 --- src/NativeProviders/MKL/fft.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/NativeProviders/MKL/fft.cpp b/src/NativeProviders/MKL/fft.cpp index b4603d71..a3598c89 100644 --- a/src/NativeProviders/MKL/fft.cpp +++ b/src/NativeProviders/MKL/fft.cpp @@ -7,11 +7,11 @@ #include "mkl_dfti.h" template -inline MKL_LONG fft_1d_inplace(const MKL_LONG n, Data x[], const Precision forward_scale, const Precision backward_scale, const DFTI_CONFIG_VALUE precision, const DFTI_CONFIG_VALUE domain, FFT fft) +inline MKL_INT64 fft_1d_inplace(const MKL_INT64 n, Data x[], const Precision forward_scale, const Precision backward_scale, const DFTI_CONFIG_VALUE precision, const DFTI_CONFIG_VALUE domain, FFT fft) { MKL_LONG status; DFTI_DESCRIPTOR_HANDLE descriptor = nullptr; - status = DftiCreateDescriptor(&descriptor, precision, domain, 1, n); + status = DftiCreateDescriptor(&descriptor, precision, domain, 1, static_cast(n)); if (0 != status) goto cleanup; status = DftiSetValue(descriptor, DFTI_FORWARD_SCALE, forward_scale); @@ -27,27 +27,27 @@ inline MKL_LONG fft_1d_inplace(const MKL_LONG n, Data x[], const Precision forwa cleanup: DftiFreeDescriptor(&descriptor); - return status; + return static_cast(status); } extern "C" { - DLLEXPORT MKL_LONG z_fft_forward_inplace(const MKL_LONG n, const double scaling, MKL_Complex16 x[]) + DLLEXPORT MKL_INT64 z_fft_forward_inplace(const MKL_INT64 n, const double scaling, MKL_Complex16 x[]) { return fft_1d_inplace(n, x, scaling, 1.0, DFTI_DOUBLE, DFTI_COMPLEX, DftiComputeForward); } - DLLEXPORT MKL_LONG c_fft_forward_inplace(const MKL_LONG n, const float scaling, MKL_Complex8 x[]) + DLLEXPORT MKL_INT64 c_fft_forward_inplace(const MKL_INT64 n, const float scaling, MKL_Complex8 x[]) { return fft_1d_inplace(n, x, scaling, 1.0f, DFTI_SINGLE, DFTI_COMPLEX, DftiComputeForward); } - DLLEXPORT MKL_LONG z_fft_backward_inplace(const MKL_LONG n, const double scaling, MKL_Complex16 x[]) + DLLEXPORT MKL_INT64 z_fft_backward_inplace(const MKL_INT64 n, const double scaling, MKL_Complex16 x[]) { return fft_1d_inplace(n, x, 1.0, scaling, DFTI_DOUBLE, DFTI_COMPLEX, DftiComputeBackward); } - DLLEXPORT MKL_LONG c_fft_backward_inplace(const MKL_LONG n, const float scaling, MKL_Complex8 x[]) + DLLEXPORT MKL_INT64 c_fft_backward_inplace(const MKL_INT64 n, const float scaling, MKL_Complex8 x[]) { return fft_1d_inplace(n, x, 1.0f, scaling, DFTI_SINGLE, DFTI_COMPLEX, DftiComputeBackward); } From 4bad397c44f3d7a63fa592b39349541610907b58 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 8 Oct 2016 18:34:48 +0200 Subject: [PATCH 16/17] FFT: activate in Fourier.Forward/Backward, with full option mapping --- src/Numerics/IntegralTransforms/Fourier.cs | 47 ++++++++++++++-- .../IFourierTransformProvider.cs | 3 +- .../ManagedFourierTransformProvider.cs | 56 +++++++++++-------- .../Mkl/MklFourierTransformProvider.cs | 8 ++- .../FourierTransformProviderTests.cs | 2 +- 5 files changed, 83 insertions(+), 33 deletions(-) diff --git a/src/Numerics/IntegralTransforms/Fourier.cs b/src/Numerics/IntegralTransforms/Fourier.cs index 74a96752..c9703e4b 100644 --- a/src/Numerics/IntegralTransforms/Fourier.cs +++ b/src/Numerics/IntegralTransforms/Fourier.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2015 Math.NET +// 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 @@ -28,6 +28,7 @@ // using System; +using MathNet.Numerics.Providers.FourierTransform; namespace MathNet.Numerics.IntegralTransforms { @@ -46,7 +47,7 @@ namespace MathNet.Numerics.IntegralTransforms /// Sample vector, where the FFT is evaluated in place. public static void Forward(Complex[] samples) { - BluesteinForward(samples, FourierOptions.Default); + Control.FourierTransformProvider.ForwardInplace(samples, FourierTransformScaling.SymmetricScaling); } /// @@ -56,7 +57,23 @@ namespace MathNet.Numerics.IntegralTransforms /// Fourier Transform Convention Options. public static void Forward(Complex[] samples, FourierOptions options) { - BluesteinForward(samples, options); + switch (options) + { + case FourierOptions.NoScaling: + case FourierOptions.AsymmetricScaling: + Control.FourierTransformProvider.ForwardInplace(samples, FourierTransformScaling.NoScaling); + break; + case FourierOptions.InverseExponent: + Control.FourierTransformProvider.BackwardInplace(samples, FourierTransformScaling.SymmetricScaling); + break; + case FourierOptions.InverseExponent | FourierOptions.NoScaling: + case FourierOptions.InverseExponent | FourierOptions.AsymmetricScaling: + Control.FourierTransformProvider.BackwardInplace(samples, FourierTransformScaling.NoScaling); + break; + default: + Control.FourierTransformProvider.ForwardInplace(samples, FourierTransformScaling.SymmetricScaling); + break; + } } /// @@ -65,7 +82,7 @@ namespace MathNet.Numerics.IntegralTransforms /// Sample vector, where the FFT is evaluated in place. public static void Inverse(Complex[] samples) { - BluesteinInverse(samples, FourierOptions.Default); + Control.FourierTransformProvider.BackwardInplace(samples, FourierTransformScaling.SymmetricScaling); } /// @@ -75,7 +92,27 @@ namespace MathNet.Numerics.IntegralTransforms /// Fourier Transform Convention Options. public static void Inverse(Complex[] samples, FourierOptions options) { - BluesteinInverse(samples, options); + switch (options) + { + case FourierOptions.NoScaling: + Control.FourierTransformProvider.BackwardInplace(samples, FourierTransformScaling.NoScaling); + break; + case FourierOptions.AsymmetricScaling: + Control.FourierTransformProvider.BackwardInplace(samples, FourierTransformScaling.BackwardScaling); + break; + case FourierOptions.InverseExponent: + Control.FourierTransformProvider.ForwardInplace(samples, FourierTransformScaling.SymmetricScaling); + break; + case FourierOptions.InverseExponent | FourierOptions.NoScaling: + Control.FourierTransformProvider.ForwardInplace(samples, FourierTransformScaling.NoScaling); + break; + case FourierOptions.InverseExponent | FourierOptions.AsymmetricScaling: + Control.FourierTransformProvider.ForwardInplace(samples, FourierTransformScaling.ForwardScaling); + break; + default: + Control.FourierTransformProvider.BackwardInplace(samples, FourierTransformScaling.SymmetricScaling); + break; + } } /// diff --git a/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs index ec44d664..f5fea08e 100644 --- a/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs @@ -38,7 +38,8 @@ namespace MathNet.Numerics.Providers.FourierTransform { NoScaling = 0, SymmetricScaling = 1, - AsymmetricScaling = 2 + BackwardScaling = 2, + ForwardScaling = 3 } public interface IFourierTransformProvider diff --git a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs index 58d42e3d..b5c89c3e 100644 --- a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs @@ -26,6 +26,7 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System.Collections; using MathNet.Numerics.IntegralTransforms; namespace MathNet.Numerics.Providers.FourierTransform @@ -41,7 +42,7 @@ namespace MathNet.Numerics.Providers.FourierTransform /// Try to find out whether the provider is available, at least in principle. /// Verification may still fail if available, but it will certainly fail if unavailable. /// - public virtual bool IsAvailable() + public bool IsAvailable() { return true; } @@ -49,26 +50,49 @@ namespace MathNet.Numerics.Providers.FourierTransform /// /// Initialize and verify that the provided is indeed available. If not, fall back to alternatives like the managed provider /// - public virtual void InitializeVerify() + public void InitializeVerify() { } - public override string ToString() + public string ToString() { return "Managed"; } - public virtual void ForwardInplace(Complex[] complex, FourierTransformScaling scaling) + public void ForwardInplace(Complex[] complex, FourierTransformScaling scaling) { - Fourier.BluesteinForward(complex, Options(scaling)); + switch (scaling) + { + case FourierTransformScaling.SymmetricScaling: + Fourier.BluesteinForward(complex, FourierOptions.Default); + break; + case FourierTransformScaling.ForwardScaling: + // Only backward scaling can be expressed with options, hence the double-inverse + Fourier.BluesteinInverse(complex, FourierOptions.AsymmetricScaling | FourierOptions.InverseExponent); + break; + default: + Fourier.BluesteinForward(complex, FourierOptions.NoScaling); + break; + } } - public virtual void BackwardInplace(Complex[] complex, FourierTransformScaling scaling) + public void BackwardInplace(Complex[] complex, FourierTransformScaling scaling) { - Fourier.BluesteinInverse(complex, Options(scaling)); + switch (scaling) + { + case FourierTransformScaling.SymmetricScaling: + Fourier.BluesteinInverse(complex, FourierOptions.Default); + break; + case FourierTransformScaling.BackwardScaling: + Fourier.BluesteinInverse(complex, FourierOptions.AsymmetricScaling); + break; + default: + Fourier.BluesteinInverse(complex, FourierOptions.NoScaling); + break; + } } - public virtual Complex[] Forward(Complex[] complexTimeSpace, FourierTransformScaling scaling) + public Complex[] Forward(Complex[] complexTimeSpace, FourierTransformScaling scaling) { Complex[] work = new Complex[complexTimeSpace.Length]; complexTimeSpace.Copy(work); @@ -76,26 +100,12 @@ namespace MathNet.Numerics.Providers.FourierTransform return work; } - public virtual Complex[] Backward(Complex[] complexFrequenceSpace, FourierTransformScaling scaling) + public Complex[] Backward(Complex[] complexFrequenceSpace, FourierTransformScaling scaling) { Complex[] work = new Complex[complexFrequenceSpace.Length]; complexFrequenceSpace.Copy(work); BackwardInplace(work, scaling); return work; } - - private FourierOptions Options(FourierTransformScaling scaling) - { - switch (scaling) - { - case FourierTransformScaling.NoScaling: - return FourierOptions.NoScaling; - case FourierTransformScaling.AsymmetricScaling: - return FourierOptions.AsymmetricScaling; - case FourierTransformScaling.SymmetricScaling: - default: - return FourierOptions.Default; - } - } } } diff --git a/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs index cf033d09..eb7b413d 100644 --- a/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs @@ -152,24 +152,26 @@ namespace MathNet.Numerics.Providers.FourierTransform.Mkl return work; } - private double ForwardScaling(FourierTransformScaling scaling, int length) + static double ForwardScaling(FourierTransformScaling scaling, int length) { switch (scaling) { case FourierTransformScaling.SymmetricScaling: return Math.Sqrt(1.0/length); + case FourierTransformScaling.ForwardScaling: + return 1.0/length; default: return 1.0; } } - private double BackwardScaling(FourierTransformScaling scaling, int length) + static double BackwardScaling(FourierTransformScaling scaling, int length) { switch (scaling) { case FourierTransformScaling.SymmetricScaling: return Math.Sqrt(1.0/length); - case FourierTransformScaling.AsymmetricScaling: + case FourierTransformScaling.BackwardScaling: return 1.0/length; default: return 1.0; diff --git a/src/UnitTests/FourierTransformProviderTests/FourierTransformProviderTests.cs b/src/UnitTests/FourierTransformProviderTests/FourierTransformProviderTests.cs index 91e67507..c29e2689 100644 --- a/src/UnitTests/FourierTransformProviderTests/FourierTransformProviderTests.cs +++ b/src/UnitTests/FourierTransformProviderTests/FourierTransformProviderTests.cs @@ -55,7 +55,7 @@ namespace MathNet.Numerics.UnitTests.FourierTransformProviderTests // real-odd transforms to imaginary odd samples.Copy(spectrum); - Control.FourierTransformProvider.ForwardInplace(spectrum, FourierTransformScaling.AsymmetricScaling); + Control.FourierTransformProvider.ForwardInplace(spectrum, FourierTransformScaling.BackwardScaling); // all real components must be zero foreach (var c in spectrum) From efcb5d5b5b7fb0c60254ddac47a4d47ac917855d Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 8 Oct 2016 18:57:04 +0200 Subject: [PATCH 17/17] FFT: verify mapping between all fourier options --- .../MatchingNaiveTransformTest.cs | 70 +++++++++++++++++++ 1 file changed, 70 insertions(+) diff --git a/src/UnitTests/IntegralTransformsTests/MatchingNaiveTransformTest.cs b/src/UnitTests/IntegralTransformsTests/MatchingNaiveTransformTest.cs index 145eeeb4..85c75e9b 100644 --- a/src/UnitTests/IntegralTransformsTests/MatchingNaiveTransformTest.cs +++ b/src/UnitTests/IntegralTransformsTests/MatchingNaiveTransformTest.cs @@ -30,6 +30,7 @@ using System; using MathNet.Numerics.Distributions; using MathNet.Numerics.IntegralTransforms; +using MathNet.Numerics.Providers.FourierTransform; using NUnit.Framework; namespace MathNet.Numerics.UnitTests.IntegralTransformsTests @@ -69,6 +70,24 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests AssertHelpers.AlmostEqual(spectrumNaive, spectrumFast, maximumErrorDecimalPlaces); } + static void VerifyInplace( + Complex[] samples, + int maximumErrorDecimalPlaces, + FourierOptions options, + Action expected, + Action actual) + { + var spectrumExpected = new Complex[samples.Length]; + samples.CopyTo(spectrumExpected, 0); + expected(spectrumExpected, options); + + var spectrumActual = new Complex[samples.Length]; + samples.CopyTo(spectrumActual, 0); + actual(spectrumActual, options); + + AssertHelpers.AlmostEqual(spectrumExpected, spectrumActual, maximumErrorDecimalPlaces); + } + /// /// Fourier Radix2XX matches naive on real sine. /// @@ -144,6 +163,57 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests Verify(samples, 10, options, Fourier.NaiveInverse, Fourier.BluesteinInverse); } + /// + /// Fourier bluestein matches providers on random power of two. + /// + /// Fourier options. + [TestCase(FourierOptions.Default)] + [TestCase(FourierOptions.NoScaling)] + [TestCase(FourierOptions.AsymmetricScaling)] + [TestCase(FourierOptions.InverseExponent)] + [TestCase(FourierOptions.InverseExponent | FourierOptions.NoScaling)] + [TestCase(FourierOptions.InverseExponent | FourierOptions.AsymmetricScaling)] + public void FourierBluesteinMatchesProvider_Random_Arbitrary(FourierOptions options) + { + var samples = Generate.RandomComplex(0x7F, GetUniform(1)); + + VerifyInplace(samples, 10, options, Fourier.Forward, Fourier.BluesteinForward); + VerifyInplace(samples, 10, options, Fourier.Inverse, Fourier.BluesteinInverse); + } + + [Test] + public void AlgorithmsMatchProvider_PowerOfTwo_Large() + { + // 65536 = 2^16 + const FourierOptions options = FourierOptions.NoScaling; + var samples = Generate.RandomComplex(65536, GetUniform(1)); + var provider = Control.FourierTransformProvider.Forward(samples, FourierTransformScaling.NoScaling); + + Verify(samples, 10, options, (a, b) => provider, Fourier.Radix2Forward); + Verify(samples, 10, options, (a, b) => provider, Fourier.BluesteinForward); + } + + [Test] + public void AlgorithmsMatchProvider_Arbitrary_Large() + { + // 30870 = 2*3*3*5*7*7*7 + const FourierOptions options = FourierOptions.NoScaling; + var samples = Generate.RandomComplex(30870, GetUniform(1)); + var provider = Control.FourierTransformProvider.Forward(samples, FourierTransformScaling.NoScaling); + + Verify(samples, 10, options, (a, b) => provider, Fourier.BluesteinForward); + } + + [Test] + public void AlgorithmsMatchProvider_Arbitrary_Large_GH286() + { + const FourierOptions options = FourierOptions.NoScaling; + var samples = Generate.RandomComplex(46500, GetUniform(1)); + var provider = Control.FourierTransformProvider.Forward(samples, FourierTransformScaling.NoScaling); + + Verify(samples, 10, options, (a, b) => provider, Fourier.BluesteinForward); + } + [Test, Explicit("Long-Running")] public void AlgorithmsMatchNaive_PowerOfTwo_Large() {