diff --git a/MathNet.Numerics.NativeProviders.sln b/MathNet.Numerics.NativeProviders.sln index 25a8d6e1..c8e8144a 100644 --- a/MathNet.Numerics.NativeProviders.sln +++ b/MathNet.Numerics.NativeProviders.sln @@ -22,6 +22,8 @@ Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "UnitTests-MKL", "src\UnitTe EndProject Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "OpenBLAS", "src\NativeProviders\Windows\OpenBLAS\OpenBLASWrapper.vcxproj", "{CB4011B6-E9A7-480B-A7B1-8492039DAAD1}" EndProject +Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "UnitTests-OpenBLAS", "src\UnitTests\UnitTests-OpenBLAS.csproj", "{96B903EF-3EE1-4569-803C-0482D2F5ED37}" +EndProject Global GlobalSection(SolutionConfigurationPlatforms) = preSolution Debug|Any CPU = Debug|Any CPU @@ -134,6 +136,24 @@ Global {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release-Signed|Win32.Build.0 = Release|Win32 {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release-Signed|x64.ActiveCfg = Release|x64 {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release-Signed|x64.Build.0 = Release|x64 + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Debug|Any CPU.ActiveCfg = Debug|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Debug|Any CPU.Build.0 = Debug|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Debug|Mixed Platforms.ActiveCfg = Debug|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Debug|Mixed Platforms.Build.0 = Debug|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Debug|Win32.ActiveCfg = Debug|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Debug|x64.ActiveCfg = Debug|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release|Any CPU.ActiveCfg = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release|Any CPU.Build.0 = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release|Mixed Platforms.ActiveCfg = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release|Mixed Platforms.Build.0 = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release|Win32.ActiveCfg = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release|x64.ActiveCfg = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-Signed|Any CPU.ActiveCfg = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-Signed|Any CPU.Build.0 = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-Signed|Mixed Platforms.ActiveCfg = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-Signed|Mixed Platforms.Build.0 = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-Signed|Win32.ActiveCfg = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-Signed|x64.ActiveCfg = Release|Any CPU EndGlobalSection GlobalSection(SolutionProperties) = preSolution HideSolutionNode = FALSE diff --git a/src/NativeProviders/OpenBLAS/capabilities.cpp b/src/NativeProviders/OpenBLAS/capabilities.cpp new file mode 100644 index 00000000..90dc57f2 --- /dev/null +++ b/src/NativeProviders/OpenBLAS/capabilities.cpp @@ -0,0 +1,77 @@ +#include "wrapper_common.h" +#include "cblas.h" + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ + + /* + Capability is supported if >0 + + Actual number can be increased over time to indicate + extensions/revisions (that do not break compatibility) + */ + DLLEXPORT int query_capability(const int capability) + { + switch (capability) + { + + // SANITY CHECKS + case 0: return 0; + case 1: return -1; + + // PLATFORM + case 8: +#ifdef _M_IX86 + return 1; +#else + return 0; +#endif + case 9: +#ifdef _M_X64 + return 1; +#else + return 0; +#endif + case 10: +#ifdef _M_IA64 + return 1; +#else + return 0; +#endif + + // COMMON/SHARED + case 64: return 7; // revision + case 66: return 1; // threading control + + // LINEAR ALGEBRA + case 128: return 1; // basic dense linear algebra + + default: return 0; // unknown or not supported + + } + } + + DLLEXPORT void set_max_threads(const blasint num_threads) + { + openblas_set_num_threads(num_threads); + } + + DLLEXPORT char* get_build_config() + { + return openblas_get_config(); + } + + DLLEXPORT char* get_cpu_core() + { + return openblas_get_corename(); + } + + DLLEXPORT int get_parallel_type() + { + return openblas_get_parallel(); + } + +#ifdef __cplusplus +} +#endif /* __cplusplus */ diff --git a/src/NativeProviders/OpenBLAS/complex.h b/src/NativeProviders/OpenBLAS/complex.h new file mode 100644 index 00000000..d107e400 --- /dev/null +++ b/src/NativeProviders/OpenBLAS/complex.h @@ -0,0 +1,39 @@ +template +struct complex +{ + _T real, imag; + + complex(_T _real = 0, _T _imag = 0) + { + real = _real; + imag = _imag; + } + + complex(const complex<_T>& right) + { + real = right.real; + imag = right.imag; + } + + complex& operator=(const complex& right) + { + real = right.real; + imag = right.imag; + return *this; + } + + complex& operator=(const _T& right) + { + real = right; + imag = 0; + return *this; + } + + template inline + complex& operator=(const complex<_Other>& right) + { + real = (_T)right.real; + imag = (_T)right.imag; + return *this; + } +}; \ No newline at end of file diff --git a/src/NativeProviders/OpenBLAS/lapack.cpp b/src/NativeProviders/OpenBLAS/lapack.cpp index 0f641319..7d2287ff 100644 --- a/src/NativeProviders/OpenBLAS/lapack.cpp +++ b/src/NativeProviders/OpenBLAS/lapack.cpp @@ -1,47 +1,7 @@ #include "cblas.h" +#include "complex.h" #define LAPACK_COMPLEX_CUSTOM - -template -struct complex -{ - T real, imag; - - complex(T _real = 0, T _imag = 0) - { - real = _real; - imag = _imag; - } - - complex(const complex& right) - { - real = right.real; - imag = right.imag; - } - - complex& operator=(const complex& right) - { - real = right.real; - imag = right.imag; - return *this; - } - - complex& operator=(const T& right) - { - real = right; - imag = 0; - return *this; - } - - template inline - complex& operator=(const complex<_Other>& right) - { - real = (T)right.real; - imag = (T)right.imag; - return *this; - } -}; - #define lapack_complex_float complex #define lapack_complex_double complex @@ -258,7 +218,7 @@ inline lapack_int complex_qr_solve_factored(lapack_int m, lapack_int n, lapack_i lapack_int info = 0; unmqr(&side, &tran, &m, &bn, &n, r, &m, tau, clone_b, &m, work, &len, &info); T one = { 1.0f, 0.0f }; - trsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, &(one.real), &(r->real), m, &(clone_b->real), m); + trsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, &(one.real), &(r->real), m, &(clone_b->real), m); copyBtoX(m, n, bn, clone_b, x); delete[] clone_b; return info; diff --git a/src/NativeProviders/Windows/OpenBLAS/OpenBLASWrapper.vcxproj b/src/NativeProviders/Windows/OpenBLAS/OpenBLASWrapper.vcxproj index 540a55f0..5a24f9b1 100644 --- a/src/NativeProviders/Windows/OpenBLAS/OpenBLASWrapper.vcxproj +++ b/src/NativeProviders/Windows/OpenBLAS/OpenBLASWrapper.vcxproj @@ -62,9 +62,9 @@ - $(ProjectDir)..\..\..\..\..\libs\OpenBLAS\include - $(ProjectDir)..\..\..\..\..\libs\OpenBLAS\x86 - $(ProjectDir)..\..\..\..\..\libs\OpenBLAS\x64 + $(ProjectDir)..\..\..\..\..\libs\OpenBLAS\include\ + $(ProjectDir)..\..\..\..\..\libs\OpenBLAS\x86\ + $(ProjectDir)..\..\..\..\..\libs\OpenBLAS\x64\ <_ProjectFileVersion>11.0.50727.1 @@ -114,8 +114,7 @@ $(OpenBLASLibDir);%(AdditionalLibraryDirectories) - - + copy "$(OpenBLASLibDir)*.dll" $(OutputPath) @@ -142,8 +141,7 @@ $(OpenBLASLibDir_x64);%(AdditionalLibraryDirectories) - - + copy "$(OpenBLASLibDir_x64)*.dll" $(OutputPath) @@ -170,8 +168,7 @@ $(OpenBLASLibDir);%(AdditionalLibraryDirectories) - - + copy "$(OpenBLASLibDir)*.dll" $(OutputPath) @@ -201,18 +198,21 @@ $(OpenBLASLibDir_x64);%(AdditionalLibraryDirectories) - - + copy "$(OpenBLASLibDir_x64)*.dll" $(OutputPath) + + + + diff --git a/src/NativeProviders/Windows/OpenBLAS/OpenBLASWrapper.vcxproj.filters b/src/NativeProviders/Windows/OpenBLAS/OpenBLASWrapper.vcxproj.filters index 62e1f976..f3bab0ba 100644 --- a/src/NativeProviders/Windows/OpenBLAS/OpenBLASWrapper.vcxproj.filters +++ b/src/NativeProviders/Windows/OpenBLAS/OpenBLASWrapper.vcxproj.filters @@ -24,10 +24,18 @@ Source Files + + Source Files + Resource Files + + + Header Files + + \ No newline at end of file diff --git a/src/Numerics/Control.cs b/src/Numerics/Control.cs index fa1028fc..2ee43d3b 100644 --- a/src/Numerics/Control.cs +++ b/src/Numerics/Control.cs @@ -127,6 +127,12 @@ namespace MathNet.Numerics { LinearAlgebraProvider = new Providers.LinearAlgebra.Mkl.MklLinearAlgebraProvider(consistency, precision, accuracy); } + + public static void UseNativeOpenBLAS() + { + LinearAlgebraProvider = new Providers.LinearAlgebra.OpenBlas.OpenBlasLinearAlgebraProvider(); + } + #endif /// diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index b783234e..471fff0d 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -157,12 +157,12 @@ - - - - - - + + + + + + @@ -192,6 +192,7 @@ + diff --git a/src/Numerics/Properties/AssemblyInfo.cs b/src/Numerics/Properties/AssemblyInfo.cs index 66cf3bdf..4da0e42c 100644 --- a/src/Numerics/Properties/AssemblyInfo.cs +++ b/src/Numerics/Properties/AssemblyInfo.cs @@ -76,6 +76,7 @@ using System.Runtime.InteropServices; #else [assembly: InternalsVisibleTo("MathNet.Numerics.UnitTests")] [assembly: InternalsVisibleTo("MathNet.Numerics.UnitTestsMKL")] +[assembly: InternalsVisibleTo("MathNet.Numerics.UnitTestsOpenBLAS")] [assembly: InternalsVisibleTo("Performance")] #endif diff --git a/src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.cs deleted file mode 100644 index 5644076b..00000000 --- a/src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.cs +++ /dev/null @@ -1,367 +0,0 @@ -// -// Math.NET Numerics, part of the Math.NET Project -// http://numerics.mathdotnet.com -// http://github.com/mathnet/mathnet-numerics -// http://mathnetnumerics.codeplex.com -// -// Copyright (c) 2009-2011 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 NATIVEGOTO - -using MathNet.Numerics.Properties; -using System; -using System.Numerics; -using System.Security; - -namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas -{ - /// - /// GotoBLAS2 linear algebra provider. - /// - public partial class GotoBlasLinearAlgebraProvider : ManagedLinearAlgebraProvider - { - /// - /// Computes the requested of the matrix. - /// - /// The type of norm to compute. - /// The number of rows in the matrix. - /// The number of columns in the matrix. - /// The matrix to compute the norm from. - /// - /// The requested of the matrix. - /// - [SecuritySafeCritical] - public override float MatrixNorm(Norm norm, int rows, int columns, float[] matrix) - { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - - if (rows <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "rows"); - } - - if (columns <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "columns"); - } - - if (matrix.Length < rows*columns) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows*columns), "matrix"); - } - - var work = new float[rows]; - return MatrixNorm(norm, rows, columns, matrix, work); - } - - /// - /// Computes the requested of the matrix. - /// - /// The type of norm to compute. - /// The number of rows in the matrix. - /// The number of columns in the matrix. - /// The matrix to compute the norm from. - /// The work array. Only used when - /// and needs to be have a length of at least M (number of rows of . - /// - /// The requested of the matrix. - /// - [SecuritySafeCritical] - public override float MatrixNorm(Norm norm, int rows, int columns, float[] matrix, float[] work) - { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - - if (rows <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "rows"); - } - - if (columns <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "columns"); - } - - if (matrix.Length < rows*columns) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows*columns), "matrix"); - } - - if (work.Length < rows) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows), "work"); - } - - return SafeNativeMethods.s_matrix_norm((byte) norm, rows, columns, matrix, work); - } - - /// - /// Computes the requested of the matrix. - /// - /// The type of norm to compute. - /// The number of rows in the matrix. - /// The number of columns in the matrix. - /// The matrix to compute the norm from. - /// - /// The requested of the matrix. - /// - [SecuritySafeCritical] - public override double MatrixNorm(Norm norm, int rows, int columns, double[] matrix) - { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - - if (rows <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "rows"); - } - - if (columns <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "columns"); - } - - if (matrix.Length < rows*columns) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows*columns), "matrix"); - } - - var work = new double[rows]; - return MatrixNorm(norm, rows, columns, matrix, work); - } - - /// - /// Computes the requested of the matrix. - /// - /// The type of norm to compute. - /// The number of rows in the matrix. - /// The number of columns in the matrix. - /// The matrix to compute the norm from. - /// The work array. Only used when - /// and needs to be have a length of at least M (number of rows of . - /// - /// The requested of the matrix. - /// - [SecuritySafeCritical] - public override double MatrixNorm(Norm norm, int rows, int columns, double[] matrix, double[] work) - { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - - if (rows <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "rows"); - } - - if (columns <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "columns"); - } - - if (matrix.Length < rows*columns) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows*columns), "matrix"); - } - - if (work.Length < rows) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows), "work"); - } - - return SafeNativeMethods.d_matrix_norm((byte) norm, rows, columns, matrix, work); - } - - /// - /// Computes the requested of the matrix. - /// - /// The type of norm to compute. - /// The number of rows in the matrix. - /// The number of columns in the matrix. - /// The matrix to compute the norm from. - /// - /// The requested of the matrix. - /// - [SecuritySafeCritical] - public override Complex32 MatrixNorm(Norm norm, int rows, int columns, Complex32[] matrix) - { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - - if (rows <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "rows"); - } - - if (columns <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "columns"); - } - - if (matrix.Length < rows*columns) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows*columns), "matrix"); - } - - var work = new float[rows]; - return MatrixNorm(norm, rows, columns, matrix, work); - } - - /// - /// Computes the requested of the matrix. - /// - /// The type of norm to compute. - /// The number of rows in the matrix. - /// The number of columns in the matrix. - /// The matrix to compute the norm from. - /// The work array. Only used when - /// and needs to be have a length of at least M (number of rows of . - /// - /// The requested of the matrix. - /// - [SecuritySafeCritical] - public override Complex32 MatrixNorm(Norm norm, int rows, int columns, Complex32[] matrix, float[] work) - { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - - if (rows <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "rows"); - } - - if (columns <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "columns"); - } - - if (matrix.Length < rows*columns) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows*columns), "matrix"); - } - - if (work.Length < rows) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows), "work"); - } - - return SafeNativeMethods.c_matrix_norm((byte) norm, rows, columns, matrix, work); - } - - /// - /// Computes the requested of the matrix. - /// - /// The type of norm to compute. - /// The number of rows in the matrix. - /// The number of columns in the matrix. - /// The matrix to compute the norm from. - /// - /// The requested of the matrix. - /// - [SecuritySafeCritical] - public override Complex MatrixNorm(Norm norm, int rows, int columns, Complex[] matrix) - { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - - if (rows <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "rows"); - } - - if (columns <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "columns"); - } - - if (matrix.Length < rows*columns) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows*columns), "matrix"); - } - - var work = new double[rows]; - return MatrixNorm(norm, rows, columns, matrix, work); - } - - /// - /// Computes the requested of the matrix. - /// - /// The type of norm to compute. - /// The number of rows in the matrix. - /// The number of columns in the matrix. - /// The matrix to compute the norm from. - /// The work array. Only used when - /// and needs to be have a length of at least M (number of rows of . - /// - /// The requested of the matrix. - /// - [SecuritySafeCritical] - public override Complex MatrixNorm(Norm norm, int rows, int columns, Complex[] matrix, double[] work) - { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - - if (rows <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "rows"); - } - - if (columns <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "columns"); - } - - if (matrix.Length < rows*columns) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows*columns), "matrix"); - } - - if (work.Length < rows) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows), "work"); - } - - return SafeNativeMethods.z_matrix_norm((byte) norm, rows, columns, matrix, work); - } - } -} - -#endif diff --git a/src/Numerics/Providers/LinearAlgebra/GotoBlas/SafeNativeMethods.cs b/src/Numerics/Providers/LinearAlgebra/GotoBlas/SafeNativeMethods.cs deleted file mode 100644 index 6ed47e53..00000000 --- a/src/Numerics/Providers/LinearAlgebra/GotoBlas/SafeNativeMethods.cs +++ /dev/null @@ -1,263 +0,0 @@ -// -// Math.NET Numerics, part of the Math.NET Project -// http://mathnet.opensourcedotnet.info -// -// Copyright (c) 2009-2010 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 NATIVEGOTO - -using System.Numerics; -using System.Runtime.InteropServices; -using System.Security; - -namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas -{ - /// - /// P/Invoke methods to the native math libraries. - /// - [SuppressUnmanagedCodeSecurity] - [SecurityCritical] - internal static class SafeNativeMethods - { - /// - /// Name of the native DLL. - /// - const string DllName = "MathNET.Numerics.GotoBLAS2.dll"; - - #region BLAS - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void s_axpy(int n, float alpha, float[] x, [In, Out] float[] y); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void d_axpy(int n, double alpha, double[] x, [In, Out] double[] y); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void c_axpy(int n, Complex32 alpha, Complex32[] x, [In, Out] Complex32[] y); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void z_axpy(int n, Complex alpha, Complex[] x, [In, Out] Complex[] y); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void s_scale(int n, float alpha, [Out] float[] x); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void d_scale(int n, double alpha, [Out] double[] x); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void c_scale(int n, Complex32 alpha, [In, Out] Complex32[] x); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void z_scale(int n, Complex alpha, [In, Out] Complex[] x); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern float s_dot_product(int n, float[] x, float[] y); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern double d_dot_product(int n, double[] x, double[] y); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern Complex32 c_dot_product(int n, Complex32[] x, Complex32[] y); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern Complex z_dot_product(int n, Complex[] x, Complex[] y); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void s_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, float alpha, float[] x, float[] y, float beta, [In, Out] float[] c); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void d_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, double alpha, double[] x, double[] y, double beta, [In, Out] double[] c); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void c_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, Complex32 alpha, Complex32[] x, Complex32[] y, Complex32 beta, [In, Out] Complex32[] c); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void z_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, Complex alpha, Complex[] x, Complex[] y, Complex beta, [In, Out] Complex[] c); - - #endregion BLAS - - #region LAPACK - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern float s_matrix_norm(byte norm, int rows, int columns, [In] float[] a, [In, Out] float[] work); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern float d_matrix_norm(byte norm, int rows, int columns, [In] double[] a, [In, Out] double[] work); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern float c_matrix_norm(byte norm, int rows, int columns, [In] Complex32[] a, [In, Out] float[] work); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern double z_matrix_norm(byte norm, int rows, int columns, [In] Complex[] a, [In, Out] double[] work); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_cholesky_factor(int n, [In, Out] float[] a); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_cholesky_factor(int n, [In, Out] double[] a); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_cholesky_factor(int n, [In, Out] Complex32[] a); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_cholesky_factor(int n, [In, Out] Complex[] a); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_lu_factor(int n, [In, Out] float[] a, [In, Out] int[] ipiv); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_lu_factor(int n, [In, Out] double[] a, [In, Out] int[] ipiv); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_lu_factor(int n, [In, Out] Complex32[] a, [In, Out] int[] ipiv); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_lu_factor(int n, [In, Out] Complex[] a, [In, Out] int[] ipiv); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_lu_inverse(int n, [In, Out] float[] a, [In, Out] float[] work, int lwork); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_lu_inverse(int n, [In, Out] double[] a, [In, Out] double[] work, int lwork); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_lu_inverse(int n, [In, Out] Complex32[] a, [In, Out] Complex32[] work, int lwork); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_lu_inverse(int n, [In, Out] Complex[] a, [In, Out] Complex[] work, int lwork); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_lu_inverse_factored(int n, [In, Out] float[] a, [In, Out] int[] ipiv, [In, Out] float[] work, int lwork); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_lu_inverse_factored(int n, [In, Out] double[] a, [In, Out] int[] ipiv, [In, Out] double[] work, int lwork); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_lu_inverse_factored(int n, [In, Out] Complex32[] a, [In, Out] int[] ipiv, [In, Out] Complex32[] work, int lwork); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_lu_inverse_factored(int n, [In, Out] Complex[] a, [In, Out] int[] ipiv, [In, Out] Complex[] work, int lwork); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_lu_solve_factored(int n, int nrhs, float[] a, [In, Out] int[] ipiv, [In, Out] float[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_lu_solve_factored(int n, int nrhs, double[] a, [In, Out] int[] ipiv, [In, Out] double[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_lu_solve_factored(int n, int nrhs, Complex32[] a, [In, Out] int[] ipiv, [In, Out] Complex32[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_lu_solve_factored(int n, int nrhs, Complex[] a, [In, Out] int[] ipiv, [In, Out] Complex[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_lu_solve(int n, int nrhs, float[] a, [In, Out] float[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_lu_solve(int n, int nrhs, double[] a, [In, Out] double[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_lu_solve(int n, int nrhs, Complex32[] a, [In, Out] Complex32[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_lu_solve(int n, int nrhs, Complex[] a, [In, Out] Complex[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_cholesky_solve(int n, int nrhs, float[] a, [In, Out] float[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_cholesky_solve(int n, int nrhs, double[] a, [In, Out] double[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_cholesky_solve(int n, int nrhs, Complex32[] a, [In, Out] Complex32[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_cholesky_solve(int n, int nrhs, Complex[] a, [In, Out] Complex[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_cholesky_solve_factored(int n, int nrhs, float[] a, [In, Out] float[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_cholesky_solve_factored(int n, int nrhs, double[] a, [In, Out] double[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_cholesky_solve_factored(int n, int nrhs, Complex32[] a, [In, Out] Complex32[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_cholesky_solve_factored(int n, int nrhs, Complex[] a, [In, Out] Complex[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_qr_factor(int m, int n, [In, Out] float[] r, [In, Out] float[] tau, [In, Out] float[] q, [In, Out] float[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_qr_factor(int m, int n, [In, Out] double[] r, [In, Out] double[] tau, [In, Out] double[] q, [In, Out] double[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_qr_factor(int m, int n, [In, Out] Complex32[] r, [In, Out] Complex32[] tau, [In, Out] Complex32[] q, [In, Out] Complex32[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_qr_factor(int m, int n, [In, Out] Complex[] r, [In, Out] Complex[] tau, [In, Out] Complex[] q, [In, Out] Complex[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_qr_solve(int m, int n, int bn, float[] r, float[] b, [In, Out] float[] x, [In, Out] float[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_qr_solve(int m, int n, int bn, double[] r, double[] b, [In, Out] double[] x, [In, Out] double[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_qr_solve(int m, int n, int bn, Complex32[] r, Complex32[] b, [In, Out] Complex32[] x, [In, Out] Complex32[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_qr_solve(int m, int n, int bn, Complex[] r, Complex[] b, [In, Out] Complex[] x, [In, Out] Complex[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_qr_solve_factored(int m, int n, int bn, float[] r, float[] b, float[] tau, [In, Out] float[] x, [In, Out] float[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_qr_solve_factored(int m, int n, int bn, double[] r, double[] b, double[] tau, [In, Out] double[] x, [In, Out] double[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_qr_solve_factored(int m, int n, int bn, Complex32[] r, Complex32[] b, Complex32[] tau, [In, Out] Complex32[] x, [In, Out] Complex32[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_qr_solve_factored(int m, int n, int bn, Complex[] r, Complex[] b, Complex[] tau, [In, Out] Complex[] x, [In, Out] Complex[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_svd_factor(bool computeVectors, int m, int n, [In, Out] float[] a, [In, Out] float[] s, [In, Out] float[] u, [In, Out] float[] v, [In, Out] float[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_svd_factor(bool computeVectors, int m, int n, [In, Out] double[] a, [In, Out] double[] s, [In, Out] double[] u, [In, Out] double[] v, [In, Out] double[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_svd_factor(bool computeVectors, int m, int n, [In, Out] Complex32[] a, [In, Out] Complex32[] s, [In, Out] Complex32[] u, [In, Out] Complex32[] v, [In, Out] Complex32[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_svd_factor(bool computeVectors, int m, int n, [In, Out] Complex[] a, [In, Out] Complex[] s, [In, Out] Complex[] u, [In, Out] Complex[] v, [In, Out] Complex[] work, int len); - - #endregion LAPACK - } -} - -#endif diff --git a/src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Complex.cs b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex.cs similarity index 73% rename from src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Complex.cs rename to src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex.cs index efa55565..a683b7bf 100644 --- a/src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Complex.cs +++ b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex.cs @@ -1,10 +1,10 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2011 Math.NET +// Copyright (c) 2009-2015 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -28,7 +28,7 @@ // OTHER DEALINGS IN THE SOFTWARE. // -#if NATIVEGOTO +#if NATIVE using MathNet.Numerics.LinearAlgebra.Factorization; using MathNet.Numerics.Properties; @@ -36,13 +36,78 @@ using System; using System.Numerics; using System.Security; -namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas +namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas { /// - /// GotoBLAS2 linear algebra provider. + /// OpenBLAS linear algebra provider. /// - public partial class GotoBlasLinearAlgebraProvider + public partial class OpenBlasLinearAlgebraProvider { + /// + /// Computes the requested of the matrix. + /// + /// The type of norm to compute. + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The matrix to compute the norm from. + /// + /// The requested of the matrix. + /// + [SecuritySafeCritical] + public override double MatrixNorm(Norm norm, int rows, int columns, Complex[] matrix) + { + if (matrix == null) + { + throw new ArgumentNullException("matrix"); + } + + if (rows <= 0) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "rows"); + } + + if (columns <= 0) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "columns"); + } + + if (matrix.Length < rows * columns) + { + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows * columns), "matrix"); + } + + var work = new double[rows]; + return SafeNativeMethods.z_matrix_norm((byte)norm, rows, columns, matrix, work); + } + + /// + /// Computes the dot product of x and y. + /// + /// The vector x. + /// The vector y. + /// The dot product of x and y. + /// This is equivalent to the DOT BLAS routine. + [SecuritySafeCritical] + public override Complex DotProduct(Complex[] x, Complex[] y) + { + if (y == null) + { + throw new ArgumentNullException("y"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (x.Length != y.Length) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength); + } + + return SafeNativeMethods.z_dot_product(x.Length, x, y); + } + /// /// Adds a scaled vector to another: result = y + alpha*x. /// @@ -164,7 +229,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas var k = transposeA == Transpose.DontTranspose ? columnsA : rowsA; var l = transposeB == Transpose.DontTranspose ? rowsB : columnsB; - if (c.Length != m*n) + if (c.Length != m * n) { throw new ArgumentException(Resources.ArgumentMatrixDimensions); } @@ -199,7 +264,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (data.Length != order*order) + if (data.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "data"); } @@ -226,7 +291,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("a"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -255,7 +320,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -286,7 +351,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("a"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -327,7 +392,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -366,12 +431,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("a"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (b.Length != columnsOfB*order) + if (b.Length != columnsOfB * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -406,7 +471,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -416,7 +481,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); } - if (b.Length != columnsOfB*order) + if (b.Length != columnsOfB * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -449,7 +514,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -484,7 +549,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("b"); } - if (b.Length != orderA*columnsB) + if (b.Length != orderA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -518,7 +583,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("b"); } - if (b.Length != orderA*columnsB) + if (b.Length != orderA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -556,7 +621,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("q"); } - if (r.Length != rowsR*columnsR) + if (r.Length != rowsR * columnsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); } @@ -566,12 +631,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (q.Length != rowsR*rowsR) + if (q.Length != rowsR * rowsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } - var work = new Complex[columnsR*Control.BlockSize]; + var work = new Complex[columnsR * Control.BlockSize]; SafeNativeMethods.z_qr_factor(rowsR, columnsR, r, tau, q, work, work.Length); } @@ -608,7 +673,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (r.Length != rowsR*columnsR) + if (r.Length != rowsR * columnsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); } @@ -618,14 +683,14 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (q.Length != rowsR*rowsR) + if (q.Length != rowsR * rowsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } - if (work.Length < columnsR*Control.BlockSize) + if (work.Length < columnsR * Control.BlockSize) { - work[0] = columnsR*Control.BlockSize; + work[0] = columnsR * Control.BlockSize; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } @@ -633,54 +698,123 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas } /// - /// Solves A*X=B for X using QR factorization of A. + /// Computes the thin QR factorization of A where M > N. /// - /// The A matrix. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The B matrix. - /// The number of columns of B. - /// On exit, the solution matrix. - /// Rows must be greater or equal to columns. - public override void QRSolve(Complex[] a, int rows, int columns, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full) + /// On entry, it is the M by N A matrix to factor. On exit, + /// it is overwritten with the Q matrix of the QR factorization. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// On exit, A N by N matrix that holds the R matrix of the + /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + [SecuritySafeCritical] + public override void ThinQRFactor(Complex[] q, int rowsA, int columnsA, Complex[] r, Complex[] tau) { - if (a == null) + if (r == null) { - throw new ArgumentNullException("a"); + throw new ArgumentNullException("r"); } - if (b == null) + if (q == null) { - throw new ArgumentNullException("b"); + throw new ArgumentNullException("q"); } - if (x == null) + if (q.Length != rowsA * columnsA) { - throw new ArgumentNullException("x"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "q"); } - if (a.Length != rows*columns) + if (tau.Length < Math.Min(rowsA, columnsA)) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (b.Length != rows*columnsB) + if (r.Length != columnsA * columnsA) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r"); + } + + var work = new Complex[columnsA * Control.BlockSize]; + SafeNativeMethods.z_qr_thin_factor(rowsA, columnsA, q, tau, r, work, work.Length); + } + + /// + /// Computes the thin QR factorization of A where M > N. + /// + /// On entry, it is the M by N A matrix to factor. On exit, + /// it is overwritten with the Q matrix of the QR factorization. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// On exit, A N by N matrix that holds the R matrix of the + /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. + /// The work array. The array must have a length of at least N, + /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal + /// work size value. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + [SecuritySafeCritical] + public override void ThinQRFactor(Complex[] q, int rowsA, int columnsA, Complex[] r, Complex[] tau, Complex[] work) + { + if (r == null) + { + throw new ArgumentNullException("r"); } - if (x.Length != columns*columnsB) + if (q == null) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + throw new ArgumentNullException("q"); } - if (rows < columns) + if (work == null) { - throw new ArgumentException(Resources.RowsLessThanColumns); + throw new ArgumentNullException("q"); + } + + if (q.Length != rowsA * columnsA) + { + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "q"); + } + + if (tau.Length < Math.Min(rowsA, columnsA)) + { + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); + } + + if (r.Length != columnsA * columnsA) + { + throw new ArgumentException( + string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r"); + } + + if (work.Length < columnsA * Control.BlockSize) + { + work[0] = columnsA * Control.BlockSize; + throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - var work = new Complex[columns*Control.BlockSize]; - QRSolve(a, rows, columns, b, columnsB, x, work); + SafeNativeMethods.z_qr_thin_factor(rowsA, columnsA, q, tau, r, work, work.Length); + } + + /// + /// Solves A*X=B for X using QR factorization of A. + /// + /// The A matrix. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + /// The type of QR factorization to perform. + /// Rows must be greater or equal to columns. + [SecuritySafeCritical] + public override void QRSolve(Complex[] a, int rows, int columns, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full) + { + var work = new Complex[columns * Control.BlockSize]; + QRSolve(a, rows, columns, b, columnsB, x, work, method); } /// @@ -695,7 +829,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. + /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. + [SecuritySafeCritical] public override void QRSolve(Complex[] a, int rows, int columns, Complex[] b, int columnsB, Complex[] x, Complex[] work, QRMethod method = QRMethod.Full) { if (a == null) @@ -718,17 +854,17 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (a.Length != rows*columns) + if (a.Length != rows * columns) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (b.Length != rows*columnsB) + if (b.Length != rows * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (x.Length != columns*columnsB) + if (x.Length != columns * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); } @@ -740,7 +876,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas if (work.Length < 1) { - work[0] = rows*Control.BlockSize; + work[0] = rows * Control.BlockSize; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } @@ -750,8 +886,8 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// /// Solves A*X=B for X using a previously QR factored matrix. /// - /// The Q matrix obtained by calling . - /// The R matrix obtained by calling . + /// The Q matrix obtained by calling . + /// The R matrix obtained by calling . /// The number of rows in the A matrix. /// The number of columns in the A matrix. /// Contains additional information on Q. Only used for the native solver @@ -759,57 +895,13 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. + /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] public override void QRSolveFactored(Complex[] q, Complex[] r, int rowsR, int columnsR, Complex[] tau, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full) { - if (r == null) - { - throw new ArgumentNullException("r"); - } - - if (q == null) - { - throw new ArgumentNullException("q"); - } - - if (b == null) - { - throw new ArgumentNullException("q"); - } - - if (x == null) - { - throw new ArgumentNullException("q"); - } - - if (r.Length != rowsR*columnsR) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); - } - - if (q.Length != rowsR*rowsR) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); - } - - if (b.Length != rowsR*columnsB) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); - } - - if (x.Length != columnsR*columnsB) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); - } - - if (rowsR < columnsR) - { - throw new ArgumentException(Resources.RowsLessThanColumns); - } - - var work = new Complex[columnsR*Control.BlockSize]; - QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, work); + var work = new Complex[columnsR * Control.BlockSize]; + QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, work, method); } /// @@ -817,9 +909,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// /// The Q matrix obtained by QR factor. This is only used for the managed provider and can be /// null for the native provider. The native provider uses the Q portion stored in the R matrix. - /// The R matrix obtained by calling . - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. + /// The R matrix obtained by calling . + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. /// Contains additional information on Q. Only used for the native solver /// and can be null for the managed provider. /// On entry the B matrix; on exit the X matrix. @@ -828,8 +920,10 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// The work array - only used in the native provider. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. + /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. - public override void QRSolveFactored(Complex[] q, Complex[] r, int rowsR, int columnsR, Complex[] tau, Complex[] b, int columnsB, Complex[] x, Complex[] work, QRMethod method = QRMethod.Full) + [SecuritySafeCritical] + public override void QRSolveFactored(Complex[] q, Complex[] r, int rowsA, int columnsA, Complex[] tau, Complex[] b, int columnsB, Complex[] x, Complex[] work, QRMethod method = QRMethod.Full) { if (r == null) { @@ -856,38 +950,54 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (r.Length != rowsR*columnsR) + int rowsQ, columnsQ, rowsR, columnsR; + if (method == QRMethod.Full) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + rowsQ = columnsQ = rowsR = rowsA; + columnsR = columnsA; + } + else + { + rowsQ = rowsA; + columnsQ = rowsR = columnsR = columnsA; } - if (q.Length != rowsR*rowsR) + if (r.Length != rowsR * columnsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsR * columnsR), "r"); } - if (b.Length != rowsR*columnsB) + if (q.Length != rowsQ * columnsQ) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsQ * columnsQ), "q"); } - if (x.Length != columnsR*columnsB) + if (b.Length != rowsA * columnsB) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsA * columnsB), "b"); } - if (rowsR < columnsR) + if (x.Length != columnsA * columnsB) { - throw new ArgumentException(Resources.RowsLessThanColumns); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, columnsA * columnsB), "x"); } if (work.Length < 1) { - work[0] = rowsR*Control.BlockSize; + work[0] = rowsA * Control.BlockSize; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - SafeNativeMethods.z_qr_solve_factored(rowsR, columnsR, columnsB, r, b, tau, x, work, work.Length); + if (method == QRMethod.Full) + { + SafeNativeMethods.z_qr_solve_factored(rowsA, columnsA, columnsB, r, b, tau, x, work, work.Length); + } + else + { + // we don't have access to the raw Q matrix any more(it is stored in R in the full QR), need to think about this. + // let just call the managed version in the meantime. The heavy lifting has already been done. -marcus + base.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, QRMethod.Thin); + } } /// @@ -926,12 +1036,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("vt"); } - if (u.Length != rowsA*rowsA) + if (u.Length != rowsA * rowsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); } - if (vt.Length != columnsA*columnsA) + if (vt.Length != columnsA * columnsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); } @@ -941,7 +1051,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); } - var work = new Complex[(2*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)]; + var work = new Complex[(2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)]; SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, work); } @@ -971,20 +1081,20 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("x"); } - if (b.Length != rowsA*columnsB) + if (b.Length != rowsA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (x.Length != columnsA*columnsB) + if (x.Length != columnsA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - var work = new Complex[(2*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)]; + var work = new Complex[(2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)]; var s = new Complex[Math.Min(rowsA, columnsA)]; - var u = new Complex[rowsA*rowsA]; - var vt = new Complex[columnsA*columnsA]; + var u = new Complex[rowsA * rowsA]; + var vt = new Complex[columnsA * columnsA]; var clone = new Complex[a.Length]; a.Copy(clone); @@ -1036,12 +1146,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (u.Length != rowsA*rowsA) + if (u.Length != rowsA * rowsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); } - if (vt.Length != columnsA*columnsA) + if (vt.Length != columnsA * columnsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); } @@ -1056,13 +1166,73 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentSingleDimensionArray, "work"); } - if (work.Length < (2*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)) + if (work.Length < (2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)) { - work[0] = (2*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA); + work[0] = (2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA); throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - SafeNativeMethods.z_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length); + if (SafeNativeMethods.z_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length) > 0) + { + throw new NonConvergenceException(); + } + } + + /// + /// Computes the eigenvalues and eigenvectors of a matrix. + /// + /// Whether the matrix is symmetric or not. + /// The order of the matrix. + /// The matrix to decompose. The lenth of the array must be order * order. + /// On output, the matrix contains the eigen vectors. The lenth of the array must be order * order. + /// On output, the eigen values (λ) of matrix in ascending value. The length of the arry must . + /// On output, the block diagonal eigenvalue matrix. The lenth of the array must be order * order. + public override void EigenDecomp(bool isSymmetric, int order, Complex[] matrix, Complex[] matrixEv, Complex[] vectorEv, Complex[] matrixD) + { + if (matrix == null) + { + throw new ArgumentNullException("matrix"); + } + + if (matrix.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrix"); + } + + if (matrixEv == null) + { + throw new ArgumentNullException("matrixEv"); + } + + if (matrixEv.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixEv"); + } + + if (vectorEv == null) + { + throw new ArgumentNullException("vectorEv"); + } + + if (vectorEv.Length != order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order), "vectorEv"); + } + + if (matrixD == null) + { + throw new ArgumentNullException("matrixD"); + } + + if (matrixD.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixD"); + } + + if (SafeNativeMethods.z_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD) > 0) + { + throw new NonConvergenceException(); + } } } } diff --git a/src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Complex32.cs b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex32.cs similarity index 73% rename from src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Complex32.cs rename to src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex32.cs index e5de77d7..9f118068 100644 --- a/src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Complex32.cs +++ b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex32.cs @@ -1,10 +1,10 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2011 Math.NET +// Copyright (c) 2009-2015 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -28,20 +28,86 @@ // OTHER DEALINGS IN THE SOFTWARE. // -#if NATIVEGOTO +#if NATIVE using MathNet.Numerics.LinearAlgebra.Factorization; using MathNet.Numerics.Properties; using System; +using System.Numerics; using System.Security; -namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas +namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas { /// - /// GotoBLAS2 linear algebra provider. + /// OpenBLAS linear algebra provider. /// - public partial class GotoBlasLinearAlgebraProvider + public partial class OpenBlasLinearAlgebraProvider { + /// + /// Computes the requested of the matrix. + /// + /// The type of norm to compute. + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The matrix to compute the norm from. + /// + /// The requested of the matrix. + /// + [SecuritySafeCritical] + public override double MatrixNorm(Norm norm, int rows, int columns, Complex32[] matrix) + { + if (matrix == null) + { + throw new ArgumentNullException("matrix"); + } + + if (rows <= 0) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "rows"); + } + + if (columns <= 0) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "columns"); + } + + if (matrix.Length < rows * columns) + { + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows * columns), "matrix"); + } + + var work = new float[rows]; + return SafeNativeMethods.c_matrix_norm((byte)norm, rows, columns, matrix, work); + } + + /// + /// Computes the dot product of x and y. + /// + /// The vector x. + /// The vector y. + /// The dot product of x and y. + /// This is equivalent to the DOT BLAS routine. + [SecuritySafeCritical] + public override Complex32 DotProduct(Complex32[] x, Complex32[] y) + { + if (y == null) + { + throw new ArgumentNullException("y"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (x.Length != y.Length) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength); + } + + return SafeNativeMethods.c_dot_product(x.Length, x, y); + } + /// /// Adds a scaled vector to another: result = y + alpha*x. /// @@ -163,7 +229,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas var k = transposeA == Transpose.DontTranspose ? columnsA : rowsA; var l = transposeB == Transpose.DontTranspose ? rowsB : columnsB; - if (c.Length != m*n) + if (c.Length != m * n) { throw new ArgumentException(Resources.ArgumentMatrixDimensions); } @@ -198,7 +264,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (data.Length != order*order) + if (data.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "data"); } @@ -225,7 +291,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("a"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -254,7 +320,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -285,7 +351,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("a"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -326,7 +392,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -365,12 +431,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("a"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (b.Length != columnsOfB*order) + if (b.Length != columnsOfB * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -405,7 +471,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -415,7 +481,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); } - if (b.Length != columnsOfB*order) + if (b.Length != columnsOfB * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -448,7 +514,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -483,7 +549,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("b"); } - if (b.Length != orderA*columnsB) + if (b.Length != orderA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -517,7 +583,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("b"); } - if (b.Length != orderA*columnsB) + if (b.Length != orderA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -555,7 +621,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("q"); } - if (r.Length != rowsR*columnsR) + if (r.Length != rowsR * columnsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); } @@ -565,12 +631,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (q.Length != rowsR*rowsR) + if (q.Length != rowsR * rowsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } - var work = new Complex32[columnsR*Control.BlockSize]; + var work = new Complex32[columnsR * Control.BlockSize]; SafeNativeMethods.c_qr_factor(rowsR, columnsR, r, tau, q, work, work.Length); } @@ -607,7 +673,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (r.Length != rowsR*columnsR) + if (r.Length != rowsR * columnsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); } @@ -617,14 +683,14 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (q.Length != rowsR*rowsR) + if (q.Length != rowsR * rowsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } - if (work.Length < columnsR*Control.BlockSize) + if (work.Length < columnsR * Control.BlockSize) { - work[0] = columnsR*Control.BlockSize; + work[0] = columnsR * Control.BlockSize; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } @@ -632,54 +698,123 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas } /// - /// Solves A*X=B for X using QR factorization of A. + /// Computes the thin QR factorization of A where M > N. /// - /// The A matrix. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The B matrix. - /// The number of columns of B. - /// On exit, the solution matrix. - /// Rows must be greater or equal to columns. - public override void QRSolve(Complex32[] a, int rows, int columns, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full) + /// On entry, it is the M by N A matrix to factor. On exit, + /// it is overwritten with the Q matrix of the QR factorization. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// On exit, A N by N matrix that holds the R matrix of the + /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + [SecuritySafeCritical] + public override void ThinQRFactor(Complex32[] q, int rowsA, int columnsA, Complex32[] r, Complex32[] tau) { - if (a == null) + if (r == null) { - throw new ArgumentNullException("a"); + throw new ArgumentNullException("r"); } - if (b == null) + if (q == null) { - throw new ArgumentNullException("b"); + throw new ArgumentNullException("q"); } - if (x == null) + if (q.Length != rowsA * columnsA) { - throw new ArgumentNullException("x"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "q"); } - if (a.Length != rows*columns) + if (tau.Length < Math.Min(rowsA, columnsA)) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (b.Length != rows*columnsB) + if (r.Length != columnsA * columnsA) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r"); + } + + var work = new Complex32[columnsA * Control.BlockSize]; + SafeNativeMethods.c_qr_thin_factor(rowsA, columnsA, q, tau, r, work, work.Length); + } + + /// + /// Computes the thin QR factorization of A where M > N. + /// + /// On entry, it is the M by N A matrix to factor. On exit, + /// it is overwritten with the Q matrix of the QR factorization. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// On exit, A N by N matrix that holds the R matrix of the + /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. + /// The work array. The array must have a length of at least N, + /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal + /// work size value. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + [SecuritySafeCritical] + public override void ThinQRFactor(Complex32[] q, int rowsA, int columnsA, Complex32[] r, Complex32[] tau, Complex32[] work) + { + if (r == null) + { + throw new ArgumentNullException("r"); } - if (x.Length != columns*columnsB) + if (q == null) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + throw new ArgumentNullException("q"); } - if (rows < columns) + if (work == null) { - throw new ArgumentException(Resources.RowsLessThanColumns); + throw new ArgumentNullException("q"); + } + + if (q.Length != rowsA * columnsA) + { + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "q"); + } + + if (tau.Length < Math.Min(rowsA, columnsA)) + { + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); + } + + if (r.Length != columnsA * columnsA) + { + throw new ArgumentException( + string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r"); + } + + if (work.Length < columnsA * Control.BlockSize) + { + work[0] = columnsA * Control.BlockSize; + throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - var work = new Complex32[columns*Control.BlockSize]; - QRSolve(a, rows, columns, b, columnsB, x, work); + SafeNativeMethods.c_qr_thin_factor(rowsA, columnsA, q, tau, r, work, work.Length); + } + + /// + /// Solves A*X=B for X using QR factorization of A. + /// + /// The A matrix. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + /// The type of QR factorization to perform. + /// Rows must be greater or equal to columns. + [SecuritySafeCritical] + public override void QRSolve(Complex32[] a, int rows, int columns, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full) + { + var work = new Complex32[columns * Control.BlockSize]; + QRSolve(a, rows, columns, b, columnsB, x, work, method); } /// @@ -694,7 +829,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. + /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. + [SecuritySafeCritical] public override void QRSolve(Complex32[] a, int rows, int columns, Complex32[] b, int columnsB, Complex32[] x, Complex32[] work, QRMethod method = QRMethod.Full) { if (a == null) @@ -717,17 +854,17 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (a.Length != rows*columns) + if (a.Length != rows * columns) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (b.Length != rows*columnsB) + if (b.Length != rows * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (x.Length != columns*columnsB) + if (x.Length != columns * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); } @@ -739,7 +876,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas if (work.Length < 1) { - work[0] = rows*Control.BlockSize; + work[0] = rows * Control.BlockSize; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } @@ -749,8 +886,8 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// /// Solves A*X=B for X using a previously QR factored matrix. /// - /// The Q matrix obtained by calling . - /// The R matrix obtained by calling . + /// The Q matrix obtained by calling . + /// The R matrix obtained by calling . /// The number of rows in the A matrix. /// The number of columns in the A matrix. /// Contains additional information on Q. Only used for the native solver @@ -758,57 +895,13 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. + /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] public override void QRSolveFactored(Complex32[] q, Complex32[] r, int rowsR, int columnsR, Complex32[] tau, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full) { - if (r == null) - { - throw new ArgumentNullException("r"); - } - - if (q == null) - { - throw new ArgumentNullException("q"); - } - - if (b == null) - { - throw new ArgumentNullException("q"); - } - - if (x == null) - { - throw new ArgumentNullException("q"); - } - - if (r.Length != rowsR*columnsR) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); - } - - if (q.Length != rowsR*rowsR) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); - } - - if (b.Length != rowsR*columnsB) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); - } - - if (x.Length != columnsR*columnsB) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); - } - - if (rowsR < columnsR) - { - throw new ArgumentException(Resources.RowsLessThanColumns); - } - - var work = new Complex32[columnsR*Control.BlockSize]; - QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, work); + var work = new Complex32[columnsR * Control.BlockSize]; + QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, work, method); } /// @@ -816,9 +909,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// /// The Q matrix obtained by QR factor. This is only used for the managed provider and can be /// null for the native provider. The native provider uses the Q portion stored in the R matrix. - /// The R matrix obtained by calling . - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. + /// The R matrix obtained by calling . + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. /// Contains additional information on Q. Only used for the native solver /// and can be null for the managed provider. /// On entry the B matrix; on exit the X matrix. @@ -827,8 +920,10 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// The work array - only used in the native provider. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. + /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. - public override void QRSolveFactored(Complex32[] q, Complex32[] r, int rowsR, int columnsR, Complex32[] tau, Complex32[] b, int columnsB, Complex32[] x, Complex32[] work, QRMethod method = QRMethod.Full) + [SecuritySafeCritical] + public override void QRSolveFactored(Complex32[] q, Complex32[] r, int rowsA, int columnsA, Complex32[] tau, Complex32[] b, int columnsB, Complex32[] x, Complex32[] work, QRMethod method = QRMethod.Full) { if (r == null) { @@ -855,38 +950,54 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (r.Length != rowsR*columnsR) + int rowsQ, columnsQ, rowsR, columnsR; + if (method == QRMethod.Full) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + rowsQ = columnsQ = rowsR = rowsA; + columnsR = columnsA; + } + else + { + rowsQ = rowsA; + columnsQ = rowsR = columnsR = columnsA; } - if (q.Length != rowsR*rowsR) + if (r.Length != rowsR * columnsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsR * columnsR), "r"); } - if (b.Length != rowsR*columnsB) + if (q.Length != rowsQ * columnsQ) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsQ * columnsQ), "q"); } - if (x.Length != columnsR*columnsB) + if (b.Length != rowsA * columnsB) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsA * columnsB), "b"); } - if (rowsR < columnsR) + if (x.Length != columnsA * columnsB) { - throw new ArgumentException(Resources.RowsLessThanColumns); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, columnsA * columnsB), "x"); } if (work.Length < 1) { - work[0] = rowsR*Control.BlockSize; + work[0] = rowsA * Control.BlockSize; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - SafeNativeMethods.c_qr_solve_factored(rowsR, columnsR, columnsB, r, b, tau, x, work, work.Length); + if (method == QRMethod.Full) + { + SafeNativeMethods.c_qr_solve_factored(rowsA, columnsA, columnsB, r, b, tau, x, work, work.Length); + } + else + { + // we don't have access to the raw Q matrix any more(it is stored in R in the full QR), need to think about this. + // let just call the managed version in the meantime. The heavy lifting has already been done. -marcus + base.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, QRMethod.Thin); + } } /// @@ -925,12 +1036,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("vt"); } - if (u.Length != rowsA*rowsA) + if (u.Length != rowsA * rowsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); } - if (vt.Length != columnsA*columnsA) + if (vt.Length != columnsA * columnsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); } @@ -940,7 +1051,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); } - var work = new Complex32[(2*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)]; + var work = new Complex32[(2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)]; SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, work); } @@ -970,20 +1081,20 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("x"); } - if (b.Length != rowsA*columnsB) + if (b.Length != rowsA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (x.Length != columnsA*columnsB) + if (x.Length != columnsA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - var work = new Complex32[(2*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)]; + var work = new Complex32[(2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)]; var s = new Complex32[Math.Min(rowsA, columnsA)]; - var u = new Complex32[rowsA*rowsA]; - var vt = new Complex32[columnsA*columnsA]; + var u = new Complex32[rowsA * rowsA]; + var vt = new Complex32[columnsA * columnsA]; var clone = new Complex32[a.Length]; a.Copy(clone); @@ -1035,12 +1146,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (u.Length != rowsA*rowsA) + if (u.Length != rowsA * rowsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); } - if (vt.Length != columnsA*columnsA) + if (vt.Length != columnsA * columnsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); } @@ -1055,13 +1166,73 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentSingleDimensionArray, "work"); } - if (work.Length < (2*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)) + if (work.Length < (2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)) { - work[0] = (2*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA); + work[0] = (2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA); throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - SafeNativeMethods.c_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length); + if (SafeNativeMethods.c_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length) > 0) + { + throw new NonConvergenceException(); + } + } + + /// + /// Computes the eigenvalues and eigenvectors of a matrix. + /// + /// Whether the matrix is symmetric or not. + /// The order of the matrix. + /// The matrix to decompose. The lenth of the array must be order * order. + /// On output, the matrix contains the eigen vectors. The lenth of the array must be order * order. + /// On output, the eigen values (λ) of matrix in ascending value. The length of the arry must . + /// On output, the block diagonal eigenvalue matrix. The lenth of the array must be order * order. + public override void EigenDecomp(bool isSymmetric, int order, Complex32[] matrix, Complex32[] matrixEv, Complex[] vectorEv, Complex32[] matrixD) + { + if (matrix == null) + { + throw new ArgumentNullException("matrix"); + } + + if (matrix.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrix"); + } + + if (matrixEv == null) + { + throw new ArgumentNullException("matrixEv"); + } + + if (matrixEv.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixEv"); + } + + if (vectorEv == null) + { + throw new ArgumentNullException("vectorEv"); + } + + if (vectorEv.Length != order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order), "vectorEv"); + } + + if (matrixD == null) + { + throw new ArgumentNullException("matrixD"); + } + + if (matrixD.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixD"); + } + + if (SafeNativeMethods.c_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD) > 0) + { + throw new NonConvergenceException(); + } } } } diff --git a/src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Double.cs b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Double.cs similarity index 73% rename from src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Double.cs rename to src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Double.cs index 8dfc7a2b..0143e0be 100644 --- a/src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Double.cs +++ b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Double.cs @@ -1,10 +1,10 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2011 Math.NET +// Copyright (c) 2009-2015 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -28,20 +28,86 @@ // OTHER DEALINGS IN THE SOFTWARE. // -#if NATIVEGOTO +#if NATIVE using MathNet.Numerics.LinearAlgebra.Factorization; using MathNet.Numerics.Properties; using System; +using System.Numerics; using System.Security; -namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas +namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas { /// - /// GotoBLAS2 linear algebra provider. + /// OpenBLAS linear algebra provider. /// - public partial class GotoBlasLinearAlgebraProvider + public partial class OpenBlasLinearAlgebraProvider { + /// + /// Computes the requested of the matrix. + /// + /// The type of norm to compute. + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The matrix to compute the norm from. + /// + /// The requested of the matrix. + /// + [SecuritySafeCritical] + public override double MatrixNorm(Norm norm, int rows, int columns, double[] matrix) + { + if (matrix == null) + { + throw new ArgumentNullException("matrix"); + } + + if (rows <= 0) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "rows"); + } + + if (columns <= 0) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "columns"); + } + + if (matrix.Length < rows * columns) + { + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows * columns), "matrix"); + } + + var work = new double[rows]; + return SafeNativeMethods.d_matrix_norm((byte)norm, rows, columns, matrix, work); + } + + /// + /// Computes the dot product of x and y. + /// + /// The vector x. + /// The vector y. + /// The dot product of x and y. + /// This is equivalent to the DOT BLAS routine. + [SecuritySafeCritical] + public override double DotProduct(double[] x, double[] y) + { + if (y == null) + { + throw new ArgumentNullException("y"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (x.Length != y.Length) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength); + } + + return SafeNativeMethods.d_dot_product(x.Length, x, y); + } + /// /// Adds a scaled vector to another: result = y + alpha*x. /// @@ -163,7 +229,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas var k = transposeA == Transpose.DontTranspose ? columnsA : rowsA; var l = transposeB == Transpose.DontTranspose ? rowsB : columnsB; - if (c.Length != m*n) + if (c.Length != m * n) { throw new ArgumentException(Resources.ArgumentMatrixDimensions); } @@ -198,7 +264,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (data.Length != order*order) + if (data.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "data"); } @@ -225,7 +291,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("a"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -254,7 +320,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -285,7 +351,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("a"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -326,7 +392,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -365,12 +431,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("a"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (b.Length != columnsOfB*order) + if (b.Length != columnsOfB * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -405,7 +471,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -415,7 +481,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); } - if (b.Length != columnsOfB*order) + if (b.Length != columnsOfB * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -448,7 +514,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -483,7 +549,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("b"); } - if (b.Length != orderA*columnsB) + if (b.Length != orderA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -517,7 +583,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("b"); } - if (b.Length != orderA*columnsB) + if (b.Length != orderA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -555,7 +621,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("q"); } - if (r.Length != rowsR*columnsR) + if (r.Length != rowsR * columnsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); } @@ -565,12 +631,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (q.Length != rowsR*rowsR) + if (q.Length != rowsR * rowsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } - var work = new double[columnsR*Control.BlockSize]; + var work = new double[columnsR * Control.BlockSize]; SafeNativeMethods.d_qr_factor(rowsR, columnsR, r, tau, q, work, work.Length); } @@ -607,7 +673,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (r.Length != rowsR*columnsR) + if (r.Length != rowsR * columnsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); } @@ -617,14 +683,14 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (q.Length != rowsR*rowsR) + if (q.Length != rowsR * rowsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } - if (work.Length < columnsR*Control.BlockSize) + if (work.Length < columnsR * Control.BlockSize) { - work[0] = columnsR*Control.BlockSize; + work[0] = columnsR * Control.BlockSize; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } @@ -632,54 +698,123 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas } /// - /// Solves A*X=B for X using QR factorization of A. + /// Computes the thin QR factorization of A where M > N. /// - /// The A matrix. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The B matrix. - /// The number of columns of B. - /// On exit, the solution matrix. - /// Rows must be greater or equal to columns. - public override void QRSolve(double[] a, int rows, int columns, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full) + /// On entry, it is the M by N A matrix to factor. On exit, + /// it is overwritten with the Q matrix of the QR factorization. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// On exit, A N by N matrix that holds the R matrix of the + /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + [SecuritySafeCritical] + public override void ThinQRFactor(double[] q, int rowsA, int columnsA, double[] r, double[] tau) { - if (a == null) + if (r == null) { - throw new ArgumentNullException("a"); + throw new ArgumentNullException("r"); } - if (b == null) + if (q == null) { - throw new ArgumentNullException("b"); + throw new ArgumentNullException("q"); } - if (x == null) + if (q.Length != rowsA * columnsA) { - throw new ArgumentNullException("x"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "q"); } - if (a.Length != rows*columns) + if (tau.Length < Math.Min(rowsA, columnsA)) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (b.Length != rows*columnsB) + if (r.Length != columnsA * columnsA) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r"); + } + + var work = new double[columnsA * Control.BlockSize]; + SafeNativeMethods.d_qr_thin_factor(rowsA, columnsA, q, tau, r, work, work.Length); + } + + /// + /// Computes the thin QR factorization of A where M > N. + /// + /// On entry, it is the M by N A matrix to factor. On exit, + /// it is overwritten with the Q matrix of the QR factorization. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// On exit, A N by N matrix that holds the R matrix of the + /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. + /// The work array. The array must have a length of at least N, + /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal + /// work size value. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + [SecuritySafeCritical] + public override void ThinQRFactor(double[] q, int rowsA, int columnsA, double[] r, double[] tau, double[] work) + { + if (r == null) + { + throw new ArgumentNullException("r"); } - if (x.Length != columns*columnsB) + if (q == null) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + throw new ArgumentNullException("q"); } - if (rows < columns) + if (work == null) { - throw new ArgumentException(Resources.RowsLessThanColumns); + throw new ArgumentNullException("q"); + } + + if (q.Length != rowsA * columnsA) + { + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "q"); + } + + if (tau.Length < Math.Min(rowsA, columnsA)) + { + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); + } + + if (r.Length != columnsA * columnsA) + { + throw new ArgumentException( + string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r"); + } + + if (work.Length < columnsA * Control.BlockSize) + { + work[0] = columnsA * Control.BlockSize; + throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - var work = new double[columns*Control.BlockSize]; - QRSolve(a, rows, columns, b, columnsB, x, work); + SafeNativeMethods.d_qr_thin_factor(rowsA, columnsA, q, tau, r, work, work.Length); + } + + /// + /// Solves A*X=B for X using QR factorization of A. + /// + /// The A matrix. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + /// The type of QR factorization to perform. + /// Rows must be greater or equal to columns. + [SecuritySafeCritical] + public override void QRSolve(double[] a, int rows, int columns, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full) + { + var work = new double[columns * Control.BlockSize]; + QRSolve(a, rows, columns, b, columnsB, x, work, method); } /// @@ -694,7 +829,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. + /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. + [SecuritySafeCritical] public override void QRSolve(double[] a, int rows, int columns, double[] b, int columnsB, double[] x, double[] work, QRMethod method = QRMethod.Full) { if (a == null) @@ -717,17 +854,17 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (a.Length != rows*columns) + if (a.Length != rows * columns) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (b.Length != rows*columnsB) + if (b.Length != rows * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (x.Length != columns*columnsB) + if (x.Length != columns * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); } @@ -739,7 +876,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas if (work.Length < 1) { - work[0] = rows*Control.BlockSize; + work[0] = rows * Control.BlockSize; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } @@ -749,8 +886,8 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// /// Solves A*X=B for X using a previously QR factored matrix. /// - /// The Q matrix obtained by calling . - /// The R matrix obtained by calling . + /// The Q matrix obtained by calling . + /// The R matrix obtained by calling . /// The number of rows in the A matrix. /// The number of columns in the A matrix. /// Contains additional information on Q. Only used for the native solver @@ -758,57 +895,13 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. + /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] public override void QRSolveFactored(double[] q, double[] r, int rowsR, int columnsR, double[] tau, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full) { - if (r == null) - { - throw new ArgumentNullException("r"); - } - - if (q == null) - { - throw new ArgumentNullException("q"); - } - - if (b == null) - { - throw new ArgumentNullException("q"); - } - - if (x == null) - { - throw new ArgumentNullException("q"); - } - - if (r.Length != rowsR*columnsR) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); - } - - if (q.Length != rowsR*rowsR) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); - } - - if (b.Length != rowsR*columnsB) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); - } - - if (x.Length != columnsR*columnsB) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); - } - - if (rowsR < columnsR) - { - throw new ArgumentException(Resources.RowsLessThanColumns); - } - - var work = new double[columnsR*Control.BlockSize]; - QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, work); + var work = new double[columnsR * Control.BlockSize]; + QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, work, method); } /// @@ -816,9 +909,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// /// The Q matrix obtained by QR factor. This is only used for the managed provider and can be /// null for the native provider. The native provider uses the Q portion stored in the R matrix. - /// The R matrix obtained by calling . - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. + /// The R matrix obtained by calling . + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. /// Contains additional information on Q. Only used for the native solver /// and can be null for the managed provider. /// On entry the B matrix; on exit the X matrix. @@ -827,8 +920,10 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// The work array - only used in the native provider. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. + /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. - public override void QRSolveFactored(double[] q, double[] r, int rowsR, int columnsR, double[] tau, double[] b, int columnsB, double[] x, double[] work, QRMethod method = QRMethod.Full) + [SecuritySafeCritical] + public override void QRSolveFactored(double[] q, double[] r, int rowsA, int columnsA, double[] tau, double[] b, int columnsB, double[] x, double[] work, QRMethod method = QRMethod.Full) { if (r == null) { @@ -855,38 +950,54 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (r.Length != rowsR*columnsR) + int rowsQ, columnsQ, rowsR, columnsR; + if (method == QRMethod.Full) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + rowsQ = columnsQ = rowsR = rowsA; + columnsR = columnsA; + } + else + { + rowsQ = rowsA; + columnsQ = rowsR = columnsR = columnsA; } - if (q.Length != rowsR*rowsR) + if (r.Length != rowsR * columnsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsR * columnsR), "r"); } - if (b.Length != rowsR*columnsB) + if (q.Length != rowsQ * columnsQ) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsQ * columnsQ), "q"); } - if (x.Length != columnsR*columnsB) + if (b.Length != rowsA * columnsB) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsA * columnsB), "b"); } - if (rowsR < columnsR) + if (x.Length != columnsA * columnsB) { - throw new ArgumentException(Resources.RowsLessThanColumns); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, columnsA * columnsB), "x"); } if (work.Length < 1) { - work[0] = rowsR*Control.BlockSize; + work[0] = rowsA * Control.BlockSize; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - SafeNativeMethods.d_qr_solve_factored(rowsR, columnsR, columnsB, r, b, tau, x, work, work.Length); + if (method == QRMethod.Full) + { + SafeNativeMethods.d_qr_solve_factored(rowsA, columnsA, columnsB, r, b, tau, x, work, work.Length); + } + else + { + // we don't have access to the raw Q matrix any more(it is stored in R in the full QR), need to think about this. + // let just call the managed version in the meantime. The heavy lifting has already been done. -marcus + base.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, QRMethod.Thin); + } } /// @@ -925,12 +1036,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("vt"); } - if (u.Length != rowsA*rowsA) + if (u.Length != rowsA * rowsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); } - if (vt.Length != columnsA*columnsA) + if (vt.Length != columnsA * columnsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); } @@ -940,7 +1051,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); } - var work = new double[Math.Max((3*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA), 5*Math.Min(rowsA, columnsA))]; + var work = new double[Math.Max((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA), 5 * Math.Min(rowsA, columnsA))]; SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, work); } @@ -970,20 +1081,20 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("x"); } - if (b.Length != rowsA*columnsB) + if (b.Length != rowsA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (x.Length != columnsA*columnsB) + if (x.Length != columnsA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - var work = new double[Math.Max((3*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA), 5*Math.Min(rowsA, columnsA))]; + var work = new double[Math.Max((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA), 5 * Math.Min(rowsA, columnsA))]; var s = new double[Math.Min(rowsA, columnsA)]; - var u = new double[rowsA*rowsA]; - var vt = new double[columnsA*columnsA]; + var u = new double[rowsA * rowsA]; + var vt = new double[columnsA * columnsA]; var clone = new double[a.Length]; a.Copy(clone); @@ -1035,12 +1146,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (u.Length != rowsA*rowsA) + if (u.Length != rowsA * rowsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); } - if (vt.Length != columnsA*columnsA) + if (vt.Length != columnsA * columnsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); } @@ -1055,13 +1166,73 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentSingleDimensionArray, "work"); } - if (work.Length < Math.Max((3*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA), 5*Math.Min(rowsA, columnsA))) + if (work.Length < Math.Max((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA), 5 * Math.Min(rowsA, columnsA))) { - work[0] = Math.Max((3*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA), 5*Math.Min(rowsA, columnsA)); + work[0] = Math.Max((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA), 5 * Math.Min(rowsA, columnsA)); throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - SafeNativeMethods.d_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length); + if (SafeNativeMethods.d_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length) > 0) + { + throw new NonConvergenceException(); + } + } + + /// + /// Computes the eigenvalues and eigenvectors of a matrix. + /// + /// Whether the matrix is symmetric or not. + /// The order of the matrix. + /// The matrix to decompose. The lenth of the array must be order * order. + /// On output, the matrix contains the eigen vectors. The lenth of the array must be order * order. + /// On output, the eigen values (λ) of matrix in ascending value. The length of the arry must . + /// On output, the block diagonal eigenvalue matrix. The lenth of the array must be order * order. + public override void EigenDecomp(bool isSymmetric, int order, double[] matrix, double[] matrixEv, Complex[] vectorEv, double[] matrixD) + { + if (matrix == null) + { + throw new ArgumentNullException("matrix"); + } + + if (matrix.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrix"); + } + + if (matrixEv == null) + { + throw new ArgumentNullException("matrixEv"); + } + + if (matrixEv.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixEv"); + } + + if (vectorEv == null) + { + throw new ArgumentNullException("vectorEv"); + } + + if (vectorEv.Length != order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order), "vectorEv"); + } + + if (matrixD == null) + { + throw new ArgumentNullException("matrixD"); + } + + if (matrixD.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixD"); + } + + if (SafeNativeMethods.d_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD) > 0) + { + throw new NonConvergenceException(); + } } } } diff --git a/src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Single.cs b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Single.cs similarity index 73% rename from src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Single.cs rename to src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Single.cs index 6c4af16e..61c934c1 100644 --- a/src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Single.cs +++ b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Single.cs @@ -1,10 +1,10 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2011 Math.NET +// Copyright (c) 2009-2015 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -28,20 +28,86 @@ // OTHER DEALINGS IN THE SOFTWARE. // -#if NATIVEGOTO +#if NATIVE using MathNet.Numerics.LinearAlgebra.Factorization; using MathNet.Numerics.Properties; using System; +using System.Numerics; using System.Security; -namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas +namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas { /// - /// GotoBLAS2 linear algebra provider. + /// OpenBLAS linear algebra provider. /// - public partial class GotoBlasLinearAlgebraProvider + public partial class OpenBlasLinearAlgebraProvider { + /// + /// Computes the requested of the matrix. + /// + /// The type of norm to compute. + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The matrix to compute the norm from. + /// + /// The requested of the matrix. + /// + [SecuritySafeCritical] + public override double MatrixNorm(Norm norm, int rows, int columns, float[] matrix) + { + if (matrix == null) + { + throw new ArgumentNullException("matrix"); + } + + if (rows <= 0) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "rows"); + } + + if (columns <= 0) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "columns"); + } + + if (matrix.Length < rows * columns) + { + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows * columns), "matrix"); + } + + var work = new float[rows]; + return SafeNativeMethods.s_matrix_norm((byte)norm, rows, columns, matrix, work); + } + + /// + /// Computes the dot product of x and y. + /// + /// The vector x. + /// The vector y. + /// The dot product of x and y. + /// This is equivalent to the DOT BLAS routine. + [SecuritySafeCritical] + public override float DotProduct(float[] x, float[] y) + { + if (y == null) + { + throw new ArgumentNullException("y"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (x.Length != y.Length) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength); + } + + return SafeNativeMethods.s_dot_product(x.Length, x, y); + } + /// /// Adds a scaled vector to another: result = y + alpha*x. /// @@ -163,7 +229,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas var k = transposeA == Transpose.DontTranspose ? columnsA : rowsA; var l = transposeB == Transpose.DontTranspose ? rowsB : columnsB; - if (c.Length != m*n) + if (c.Length != m * n) { throw new ArgumentException(Resources.ArgumentMatrixDimensions); } @@ -198,7 +264,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (data.Length != order*order) + if (data.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "data"); } @@ -225,7 +291,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("a"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -254,7 +320,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -285,7 +351,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("a"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -326,7 +392,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -365,12 +431,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("a"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (b.Length != columnsOfB*order) + if (b.Length != columnsOfB * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -405,7 +471,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -415,7 +481,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); } - if (b.Length != columnsOfB*order) + if (b.Length != columnsOfB * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -448,7 +514,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -483,7 +549,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("b"); } - if (b.Length != orderA*columnsB) + if (b.Length != orderA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -517,7 +583,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("b"); } - if (b.Length != orderA*columnsB) + if (b.Length != orderA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -555,7 +621,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("q"); } - if (r.Length != rowsR*columnsR) + if (r.Length != rowsR * columnsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); } @@ -565,12 +631,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (q.Length != rowsR*rowsR) + if (q.Length != rowsR * rowsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } - var work = new float[columnsR*Control.BlockSize]; + var work = new float[columnsR * Control.BlockSize]; SafeNativeMethods.s_qr_factor(rowsR, columnsR, r, tau, q, work, work.Length); } @@ -607,7 +673,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (r.Length != rowsR*columnsR) + if (r.Length != rowsR * columnsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); } @@ -617,14 +683,14 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (q.Length != rowsR*rowsR) + if (q.Length != rowsR * rowsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } - if (work.Length < columnsR*Control.BlockSize) + if (work.Length < columnsR * Control.BlockSize) { - work[0] = columnsR*Control.BlockSize; + work[0] = columnsR * Control.BlockSize; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } @@ -632,54 +698,123 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas } /// - /// Solves A*X=B for X using QR factorization of A. + /// Computes the thin QR factorization of A where M > N. /// - /// The A matrix. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The B matrix. - /// The number of columns of B. - /// On exit, the solution matrix. - /// Rows must be greater or equal to columns. - public override void QRSolve(float[] a, int rows, int columns, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full) + /// On entry, it is the M by N A matrix to factor. On exit, + /// it is overwritten with the Q matrix of the QR factorization. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// On exit, A N by N matrix that holds the R matrix of the + /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + [SecuritySafeCritical] + public override void ThinQRFactor(float[] q, int rowsA, int columnsA, float[] r, float[] tau) { - if (a == null) + if (r == null) { - throw new ArgumentNullException("a"); + throw new ArgumentNullException("r"); } - if (b == null) + if (q == null) { - throw new ArgumentNullException("b"); + throw new ArgumentNullException("q"); } - if (x == null) + if (q.Length != rowsA * columnsA) { - throw new ArgumentNullException("x"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "q"); } - if (a.Length != rows*columns) + if (tau.Length < Math.Min(rowsA, columnsA)) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (b.Length != rows*columnsB) + if (r.Length != columnsA * columnsA) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r"); + } + + var work = new float[columnsA * Control.BlockSize]; + SafeNativeMethods.s_qr_thin_factor(rowsA, columnsA, q, tau, r, work, work.Length); + } + + /// + /// Computes the thin QR factorization of A where M > N. + /// + /// On entry, it is the M by N A matrix to factor. On exit, + /// it is overwritten with the Q matrix of the QR factorization. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// On exit, A N by N matrix that holds the R matrix of the + /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. + /// The work array. The array must have a length of at least N, + /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal + /// work size value. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + [SecuritySafeCritical] + public override void ThinQRFactor(float[] q, int rowsA, int columnsA, float[] r, float[] tau, float[] work) + { + if (r == null) + { + throw new ArgumentNullException("r"); } - if (x.Length != columns*columnsB) + if (q == null) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + throw new ArgumentNullException("q"); } - if (rows < columns) + if (work == null) { - throw new ArgumentException(Resources.RowsLessThanColumns); + throw new ArgumentNullException("q"); + } + + if (q.Length != rowsA * columnsA) + { + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "q"); + } + + if (tau.Length < Math.Min(rowsA, columnsA)) + { + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); + } + + if (r.Length != columnsA * columnsA) + { + throw new ArgumentException( + string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r"); + } + + if (work.Length < columnsA * Control.BlockSize) + { + work[0] = columnsA * Control.BlockSize; + throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - var work = new float[columns*Control.BlockSize]; - QRSolve(a, rows, columns, b, columnsB, x, work); + SafeNativeMethods.s_qr_thin_factor(rowsA, columnsA, q, tau, r, work, work.Length); + } + + /// + /// Solves A*X=B for X using QR factorization of A. + /// + /// The A matrix. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + /// The type of QR factorization to perform. + /// Rows must be greater or equal to columns. + [SecuritySafeCritical] + public override void QRSolve(float[] a, int rows, int columns, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full) + { + var work = new float[columns * Control.BlockSize]; + QRSolve(a, rows, columns, b, columnsB, x, work, method); } /// @@ -694,7 +829,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. + /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. + [SecuritySafeCritical] public override void QRSolve(float[] a, int rows, int columns, float[] b, int columnsB, float[] x, float[] work, QRMethod method = QRMethod.Full) { if (a == null) @@ -717,17 +854,17 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (a.Length != rows*columns) + if (a.Length != rows * columns) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (b.Length != rows*columnsB) + if (b.Length != rows * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (x.Length != columns*columnsB) + if (x.Length != columns * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); } @@ -739,7 +876,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas if (work.Length < 1) { - work[0] = rows*Control.BlockSize; + work[0] = rows * Control.BlockSize; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } @@ -749,8 +886,8 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// /// Solves A*X=B for X using a previously QR factored matrix. /// - /// The Q matrix obtained by calling . - /// The R matrix obtained by calling . + /// The Q matrix obtained by calling . + /// The R matrix obtained by calling . /// The number of rows in the A matrix. /// The number of columns in the A matrix. /// Contains additional information on Q. Only used for the native solver @@ -758,57 +895,13 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. + /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] public override void QRSolveFactored(float[] q, float[] r, int rowsR, int columnsR, float[] tau, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full) { - if (r == null) - { - throw new ArgumentNullException("r"); - } - - if (q == null) - { - throw new ArgumentNullException("q"); - } - - if (b == null) - { - throw new ArgumentNullException("q"); - } - - if (x == null) - { - throw new ArgumentNullException("q"); - } - - if (r.Length != rowsR*columnsR) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); - } - - if (q.Length != rowsR*rowsR) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); - } - - if (b.Length != rowsR*columnsB) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); - } - - if (x.Length != columnsR*columnsB) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); - } - - if (rowsR < columnsR) - { - throw new ArgumentException(Resources.RowsLessThanColumns); - } - - var work = new float[columnsR*Control.BlockSize]; - QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, work); + var work = new float[columnsR * Control.BlockSize]; + QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, work, method); } /// @@ -816,9 +909,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// /// The Q matrix obtained by QR factor. This is only used for the managed provider and can be /// null for the native provider. The native provider uses the Q portion stored in the R matrix. - /// The R matrix obtained by calling . - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. + /// The R matrix obtained by calling . + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. /// Contains additional information on Q. Only used for the native solver /// and can be null for the managed provider. /// On entry the B matrix; on exit the X matrix. @@ -827,8 +920,10 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// The work array - only used in the native provider. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. + /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. - public override void QRSolveFactored(float[] q, float[] r, int rowsR, int columnsR, float[] tau, float[] b, int columnsB, float[] x, float[] work, QRMethod method = QRMethod.Full) + [SecuritySafeCritical] + public override void QRSolveFactored(float[] q, float[] r, int rowsA, int columnsA, float[] tau, float[] b, int columnsB, float[] x, float[] work, QRMethod method = QRMethod.Full) { if (r == null) { @@ -855,38 +950,54 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (r.Length != rowsR*columnsR) + int rowsQ, columnsQ, rowsR, columnsR; + if (method == QRMethod.Full) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + rowsQ = columnsQ = rowsR = rowsA; + columnsR = columnsA; + } + else + { + rowsQ = rowsA; + columnsQ = rowsR = columnsR = columnsA; } - if (q.Length != rowsR*rowsR) + if (r.Length != rowsR * columnsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsR * columnsR), "r"); } - if (b.Length != rowsR*columnsB) + if (q.Length != rowsQ * columnsQ) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsQ * columnsQ), "q"); } - if (x.Length != columnsR*columnsB) + if (b.Length != rowsA * columnsB) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsA * columnsB), "b"); } - if (rowsR < columnsR) + if (x.Length != columnsA * columnsB) { - throw new ArgumentException(Resources.RowsLessThanColumns); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, columnsA * columnsB), "x"); } if (work.Length < 1) { - work[0] = rowsR*Control.BlockSize; + work[0] = rowsA * Control.BlockSize; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - SafeNativeMethods.s_qr_solve_factored(rowsR, columnsR, columnsB, r, b, tau, x, work, work.Length); + if (method == QRMethod.Full) + { + SafeNativeMethods.s_qr_solve_factored(rowsA, columnsA, columnsB, r, b, tau, x, work, work.Length); + } + else + { + // we don't have access to the raw Q matrix any more(it is stored in R in the full QR), need to think about this. + // let just call the managed version in the meantime. The heavy lifting has already been done. -marcus + base.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, QRMethod.Thin); + } } /// @@ -925,12 +1036,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("vt"); } - if (u.Length != rowsA*rowsA) + if (u.Length != rowsA * rowsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); } - if (vt.Length != columnsA*columnsA) + if (vt.Length != columnsA * columnsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); } @@ -940,7 +1051,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); } - var work = new float[Math.Max(((3*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)), 5*Math.Min(rowsA, columnsA))]; + var work = new float[Math.Max(((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)), 5 * Math.Min(rowsA, columnsA))]; SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, work); } @@ -970,20 +1081,20 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("x"); } - if (b.Length != rowsA*columnsB) + if (b.Length != rowsA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (x.Length != columnsA*columnsB) + if (x.Length != columnsA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - var work = new float[Math.Max(((3*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)), 5*Math.Min(rowsA, columnsA))]; + var work = new float[Math.Max(((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)), 5 * Math.Min(rowsA, columnsA))]; var s = new float[Math.Min(rowsA, columnsA)]; - var u = new float[rowsA*rowsA]; - var vt = new float[columnsA*columnsA]; + var u = new float[rowsA * rowsA]; + var vt = new float[columnsA * columnsA]; var clone = new float[a.Length]; a.Copy(clone); @@ -1035,12 +1146,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (u.Length != rowsA*rowsA) + if (u.Length != rowsA * rowsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); } - if (vt.Length != columnsA*columnsA) + if (vt.Length != columnsA * columnsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); } @@ -1055,13 +1166,73 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentSingleDimensionArray, "work"); } - if (work.Length < Math.Max(((3*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)), 5*Math.Min(rowsA, columnsA))) + if (work.Length < Math.Max(((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)), 5 * Math.Min(rowsA, columnsA))) { - work[0] = Math.Max(((3*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)), 5*Math.Min(rowsA, columnsA)); + work[0] = Math.Max(((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)), 5 * Math.Min(rowsA, columnsA)); throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - SafeNativeMethods.s_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length); + if (SafeNativeMethods.s_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length) > 0) + { + throw new NonConvergenceException(); + } + } + + /// + /// Computes the eigenvalues and eigenvectors of a matrix. + /// + /// Whether the matrix is symmetric or not. + /// The order of the matrix. + /// The matrix to decompose. The lenth of the array must be order * order. + /// On output, the matrix contains the eigen vectors. The lenth of the array must be order * order. + /// On output, the eigen values (λ) of matrix in ascending value. The length of the arry must . + /// On output, the block diagonal eigenvalue matrix. The lenth of the array must be order * order. + public override void EigenDecomp(bool isSymmetric, int order, float[] matrix, float[] matrixEv, Complex[] vectorEv, float[] matrixD) + { + if (matrix == null) + { + throw new ArgumentNullException("matrix"); + } + + if (matrix.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrix"); + } + + if (matrixEv == null) + { + throw new ArgumentNullException("matrixEv"); + } + + if (matrixEv.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixEv"); + } + + if (vectorEv == null) + { + throw new ArgumentNullException("vectorEv"); + } + + if (vectorEv.Length != order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order), "vectorEv"); + } + + if (matrixD == null) + { + throw new ArgumentNullException("matrixD"); + } + + if (matrixD.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixD"); + } + + if (SafeNativeMethods.s_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD) > 0) + { + throw new NonConvergenceException(); + } } } } diff --git a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.cs new file mode 100644 index 00000000..6a70a194 --- /dev/null +++ b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.cs @@ -0,0 +1,107 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2015 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 MathNet.Numerics.Properties; +using System; +using System.Numerics; +using System.Security; + +namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas +{ + public enum ParallelType : int + { + Sequential = 0, + Thread = 1, + OpenMP = 2 + } + + /// + /// OpenBLAS linear algebra provider. + /// + public partial class OpenBlasLinearAlgebraProvider : ManagedLinearAlgebraProvider + { + bool _nativeIX86; + bool _nativeX64; + bool _nativeIA64; + bool _nativeARM; + + public OpenBlasLinearAlgebraProvider() + { + + } + + public override void InitializeVerify() + { + int linearAlgebra; + try + { + // Load the native library + NativeProviderLoader.TryLoad(SafeNativeMethods.DllName); + + _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); + } + + // set threading settings, if supported + if (SafeNativeMethods.query_capability((int)ProviderConfig.Threading) > 0) + { + SafeNativeMethods.set_max_threads(Control.MaxDegreeOfParallelism); + } + } + + public override string ToString() + { + return string.Format("OpenBLAS\r\nProvider revision: {0}\r\nCPU core name: {1}\r\nLibrary config: {1})", + SafeNativeMethods.query_capability((int)ProviderConfig.Revision), + SafeNativeMethods.get_cpu_core(), + SafeNativeMethods.get_build_config()); + } + } +} + +#endif diff --git a/src/Numerics/Providers/LinearAlgebra/OpenBlas/SafeNativeMethods.cs b/src/Numerics/Providers/LinearAlgebra/OpenBlas/SafeNativeMethods.cs new file mode 100644 index 00000000..fa59a5d7 --- /dev/null +++ b/src/Numerics/Providers/LinearAlgebra/OpenBlas/SafeNativeMethods.cs @@ -0,0 +1,303 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009-2010 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.LinearAlgebra.OpenBlas +{ + /// + /// P/Invoke methods to the native math libraries. + /// + [SuppressUnmanagedCodeSecurity] + [SecurityCritical] + internal static class SafeNativeMethods + { + /// + /// Name of the native DLL. + /// + const string _DllName = "MathNET.Numerics.OpenBLAS.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 string get_build_config(); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern string get_cpu_core(); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern ParallelType get_parallel_type(); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void set_max_threads(int num_threads); + + #region BLAS + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void s_axpy(int n, float alpha, float[] x, [In, Out] float[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void d_axpy(int n, double alpha, double[] x, [In, Out] double[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void c_axpy(int n, Complex32 alpha, Complex32[] x, [In, Out] Complex32[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void z_axpy(int n, Complex alpha, Complex[] x, [In, Out] Complex[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void s_scale(int n, float alpha, [Out] float[] x); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void d_scale(int n, double alpha, [Out] double[] x); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void c_scale(int n, Complex32 alpha, [In, Out] Complex32[] x); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void z_scale(int n, Complex alpha, [In, Out] Complex[] x); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern float s_dot_product(int n, float[] x, float[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern double d_dot_product(int n, double[] x, double[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern Complex32 c_dot_product(int n, Complex32[] x, Complex32[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern Complex z_dot_product(int n, Complex[] x, Complex[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void s_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, float alpha, float[] x, float[] y, float beta, [In, Out] float[] c); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void d_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, double alpha, double[] x, double[] y, double beta, [In, Out] double[] c); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void c_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, Complex32 alpha, Complex32[] x, Complex32[] y, Complex32 beta, [In, Out] Complex32[] c); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void z_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, Complex alpha, Complex[] x, Complex[] y, Complex beta, [In, Out] Complex[] c); + + #endregion BLAS + + #region LAPACK + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern float s_matrix_norm(byte norm, int rows, int columns, [In] float[] a, [In, Out] float[] work); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern float d_matrix_norm(byte norm, int rows, int columns, [In] double[] a, [In, Out] double[] work); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern float c_matrix_norm(byte norm, int rows, int columns, [In] Complex32[] a, [In, Out] float[] work); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern double z_matrix_norm(byte norm, int rows, int columns, [In] Complex[] a, [In, Out] double[] work); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_lu_factor(int n, [In, Out] float[] a, [In, Out] int[] ipiv); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_lu_factor(int n, [In, Out] double[] a, [In, Out] int[] ipiv); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_lu_factor(int n, [In, Out] Complex32[] a, [In, Out] int[] ipiv); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_lu_factor(int n, [In, Out] Complex[] a, [In, Out] int[] ipiv); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_lu_inverse(int n, [In, Out] float[] a, [In, Out] float[] work, int lwork); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_lu_inverse(int n, [In, Out] double[] a, [In, Out] double[] work, int lwork); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_lu_inverse(int n, [In, Out] Complex32[] a, [In, Out] Complex32[] work, int lwork); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_lu_inverse(int n, [In, Out] Complex[] a, [In, Out] Complex[] work, int lwork); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_lu_inverse_factored(int n, [In, Out] float[] a, [In, Out] int[] ipiv, [In, Out] float[] work, int lwork); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_lu_inverse_factored(int n, [In, Out] double[] a, [In, Out] int[] ipiv, [In, Out] double[] work, int lwork); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_lu_inverse_factored(int n, [In, Out] Complex32[] a, [In, Out] int[] ipiv, [In, Out] Complex32[] work, int lwork); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_lu_inverse_factored(int n, [In, Out] Complex[] a, [In, Out] int[] ipiv, [In, Out] Complex[] work, int lwork); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_lu_solve_factored(int n, int nrhs, float[] a, [In, Out] int[] ipiv, [In, Out] float[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_lu_solve_factored(int n, int nrhs, double[] a, [In, Out] int[] ipiv, [In, Out] double[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_lu_solve_factored(int n, int nrhs, Complex32[] a, [In, Out] int[] ipiv, [In, Out] Complex32[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_lu_solve_factored(int n, int nrhs, Complex[] a, [In, Out] int[] ipiv, [In, Out] Complex[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_lu_solve(int n, int nrhs, float[] a, [In, Out] float[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_lu_solve(int n, int nrhs, double[] a, [In, Out] double[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_lu_solve(int n, int nrhs, Complex32[] a, [In, Out] Complex32[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_lu_solve(int n, int nrhs, Complex[] a, [In, Out] Complex[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_cholesky_factor(int n, [In, Out] float[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_cholesky_factor(int n, [In, Out] double[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_cholesky_factor(int n, [In, Out] Complex32[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_cholesky_factor(int n, [In, Out] Complex[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_cholesky_solve(int n, int nrhs, float[] a, [In, Out] float[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_cholesky_solve(int n, int nrhs, double[] a, [In, Out] double[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_cholesky_solve(int n, int nrhs, Complex32[] a, [In, Out] Complex32[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_cholesky_solve(int n, int nrhs, Complex[] a, [In, Out] Complex[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_cholesky_solve_factored(int n, int nrhs, float[] a, [In, Out] float[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_cholesky_solve_factored(int n, int nrhs, double[] a, [In, Out] double[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_cholesky_solve_factored(int n, int nrhs, Complex32[] a, [In, Out] Complex32[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_cholesky_solve_factored(int n, int nrhs, Complex[] a, [In, Out] Complex[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_qr_factor(int m, int n, [In, Out] float[] r, [In, Out] float[] tau, [In, Out] float[] q, [In, Out] float[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_qr_thin_factor(int m, int n, [In, Out] float[] q, [In, Out] float[] tau, [In, Out] float[] r, [In, Out] float[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_qr_factor(int m, int n, [In, Out] double[] r, [In, Out] double[] tau, [In, Out] double[] q, [In, Out] double[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_qr_thin_factor(int m, int n, [In, Out] double[] q, [In, Out] double[] tau, [In, Out] double[] r, [In, Out] double[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_qr_factor(int m, int n, [In, Out] Complex32[] r, [In, Out] Complex32[] tau, [In, Out] Complex32[] q, [In, Out] Complex32[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_qr_thin_factor(int m, int n, [In, Out] Complex32[] q, [In, Out] Complex32[] tau, [In, Out] Complex32[] r, [In, Out] Complex32[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_qr_factor(int m, int n, [In, Out] Complex[] r, [In, Out] Complex[] tau, [In, Out] Complex[] q, [In, Out] Complex[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_qr_thin_factor(int m, int n, [In, Out] Complex[] q, [In, Out] Complex[] tau, [In, Out] Complex[] r, [In, Out] Complex[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_qr_solve(int m, int n, int bn, float[] r, float[] b, [In, Out] float[] x, [In, Out] float[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_qr_solve(int m, int n, int bn, double[] r, double[] b, [In, Out] double[] x, [In, Out] double[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_qr_solve(int m, int n, int bn, Complex32[] r, Complex32[] b, [In, Out] Complex32[] x, [In, Out] Complex32[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_qr_solve(int m, int n, int bn, Complex[] r, Complex[] b, [In, Out] Complex[] x, [In, Out] Complex[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_qr_solve_factored(int m, int n, int bn, float[] r, float[] b, float[] tau, [In, Out] float[] x, [In, Out] float[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_qr_solve_factored(int m, int n, int bn, double[] r, double[] b, double[] tau, [In, Out] double[] x, [In, Out] double[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_qr_solve_factored(int m, int n, int bn, Complex32[] r, Complex32[] b, Complex32[] tau, [In, Out] Complex32[] x, [In, Out] Complex32[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_qr_solve_factored(int m, int n, int bn, Complex[] r, Complex[] b, Complex[] tau, [In, Out] Complex[] x, [In, Out] Complex[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_svd_factor(bool computeVectors, int m, int n, [In, Out] float[] a, [In, Out] float[] s, [In, Out] float[] u, [In, Out] float[] v, [In, Out] float[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_svd_factor(bool computeVectors, int m, int n, [In, Out] double[] a, [In, Out] double[] s, [In, Out] double[] u, [In, Out] double[] v, [In, Out] double[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_svd_factor(bool computeVectors, int m, int n, [In, Out] Complex32[] a, [In, Out] Complex32[] s, [In, Out] Complex32[] u, [In, Out] Complex32[] v, [In, Out] Complex32[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_svd_factor(bool computeVectors, int m, int n, [In, Out] Complex[] a, [In, Out] Complex[] s, [In, Out] Complex[] u, [In, Out] Complex[] v, [In, Out] Complex[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_eigen(bool isSymmetric, int n, [In, Out] float[] a, [In, Out] float[] vectors, [In, Out] Complex[] values, [In, Out] float[] d); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_eigen(bool isSymmetric, int n, [In, Out] double[] a, [In, Out] double[] vectors, [In, Out] Complex[] values, [In, Out] double[] d); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_eigen(bool isSymmetric, int n, [In, Out] Complex32[] a, [In, Out] Complex32[] vectors, [In, Out] Complex[] values, [In, Out] Complex32[] d); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_eigen(bool isSymmetric, int n, [In, Out] Complex[] a, [In, Out] Complex[] vectors, [In, Out] Complex[] values, [In, Out] Complex[] d); + + #endregion LAPACK + } +} + +#endif diff --git a/src/Numerics/Providers/LinearAlgebra/ProviderCapabilities.cs b/src/Numerics/Providers/LinearAlgebra/ProviderCapabilities.cs new file mode 100644 index 00000000..b07d4e0f --- /dev/null +++ b/src/Numerics/Providers/LinearAlgebra/ProviderCapabilities.cs @@ -0,0 +1,55 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2015 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.LinearAlgebra +{ + public enum ProviderPlatform : int + { + x86 = 8, + x64 = 9, + ia64 = 10, + arm = 11, + } + + public enum ProviderConfig : int + { + Revision = 64, + Precision = 65, + Threading = 66, + Memory = 67, + } + + public enum ProviderCapability : int + { + LinearAlgebra = 128, + Optimization = 256, + FFT = 384, + } +} \ No newline at end of file diff --git a/src/UnitTests/UnitTests-MKL.csproj b/src/UnitTests/UnitTests-MKL.csproj index bbd73bbf..3d8c040f 100644 --- a/src/UnitTests/UnitTests-MKL.csproj +++ b/src/UnitTests/UnitTests-MKL.csproj @@ -16,7 +16,7 @@ ..\..\ - TRACE;NATIVE + TRACE;NATIVE;NATIVEMKL; ..\..\out\MKL\Windows\ ..\..\obj\MKL\Windows\x86\ ..\..\obj\MKL\Windows\x86\ @@ -28,7 +28,7 @@ AnyCPU - TRACE;DEBUG;NATIVE + TRACE;DEBUG;NATIVE;NATIVEMKL; ..\..\out\MKL\Windows\ ..\..\obj\MKL\Windows\x86\ ..\..\obj\MKL\Windows\x86\ @@ -345,4 +345,7 @@ True + + + \ No newline at end of file diff --git a/src/UnitTests/UnitTests-OpenBLAS.csproj b/src/UnitTests/UnitTests-OpenBLAS.csproj new file mode 100644 index 00000000..e49fc5bf --- /dev/null +++ b/src/UnitTests/UnitTests-OpenBLAS.csproj @@ -0,0 +1,351 @@ + + + + 10.0 + Debug + AnyCPU + 8.0.30703 + 2.0 + {96B903EF-3EE1-4569-803C-0482D2F5ED37} + Library + Properties + MathNet.Numerics.UnitTests + MathNet.Numerics.UnitTestsOpenBLAS + v4.5 + 512 + ..\..\ + + + TRACE;NATIVE;NATIVEOPENBLAS; + ..\..\out\OpenBLAS\Windows\ + ..\..\obj\OpenBLAS\Windows\AnyCPU\ + ..\..\obj\OpenBLAS\Windows\AnyCPU\ + true + pdbonly + prompt + MinimumRecommendedRules.ruleset + 1591 + AnyCPU + + + TRACE;DEBUG;NATIVE;NATIVEOPENBLAS; + ..\..\out\OpenBLAS\Windows\ + ..\..\obj\OpenBLAS\Windows\AnyCPU\ + ..\..\obj\OpenBLAS\Windows\AnyCPU\ + false + full + true + prompt + 4 + 1591 + AnyCPU + + + + + + + + + + + + + + + + + + + data\Codeplex-5667.csv + Always + + + data\Github-Cureos-1.csv + Always + + + data\Matlab\A.mat + Always + + + data\Matlab\collection-nocompress.mat + Always + + + data\Matlab\collection.mat + Always + + + data\Matlab\complex.mat + Always + + + data\Matlab\sparse-large.mat + Always + + + data\Matlab\sparse-small.mat + Always + + + data\Matlab\sparse_complex.mat + Always + + + data\Matlab\v.mat + Always + + + data\NIST\AtmWtAgt.dat + Always + + + data\NIST\Bennett5.dat + Always + + + data\NIST\BoxBOD.dat + Always + + + data\NIST\Chwirut1.dat + Always + + + data\NIST\Chwirut2.dat + Always + + + data\NIST\DanWood.dat + Always + + + data\NIST\Eckerle4.dat + Always + + + data\NIST\ENSO.dat + Always + + + data\NIST\Filip.dat + Always + + + data\NIST\Gauss1.dat + Always + + + data\NIST\Gauss2.dat + Always + + + data\NIST\Gauss3.dat + Always + + + data\NIST\Hahn1.dat + Always + + + data\NIST\Kirby2.dat + Always + + + data\NIST\Lanczos1.dat + Always + + + data\NIST\Lanczos2.dat + Always + + + data\NIST\Lanczos3.dat + Always + + + data\NIST\Lew.dat + Always + + + data\NIST\Longley.dat + Always + + + data\NIST\Lottery.dat + Always + + + data\NIST\Mavro.dat + Always + + + data\NIST\MGH09.dat + Always + + + data\NIST\MGH10.dat + Always + + + data\NIST\MGH17.dat + Always + + + data\NIST\Michelso.dat + Always + + + data\NIST\Misra1a.dat + Always + + + data\NIST\Misra1b.dat + Always + + + data\NIST\Misra1c.dat + Always + + + data\NIST\Misra1d.dat + Always + + + data\NIST\Nelson.dat + Always + + + data\NIST\NoInt1.dat + Always + + + data\NIST\NoInt2.dat + Always + + + data\NIST\Norris.dat + Always + + + data\NIST\NumAcc1.dat + Always + + + data\NIST\NumAcc2.dat + Always + + + data\NIST\NumAcc3.dat + Always + + + data\NIST\NumAcc4.dat + Always + + + data\NIST\Pontius.dat + Always + + + data\NIST\Rat42.dat + Always + + + data\NIST\Rat43.dat + Always + + + data\NIST\Roszman1.dat + Always + + + data\NIST\SiRstvt.dat + Always + + + data\NIST\SmLs01t.dat + Always + + + data\NIST\SmLs02t.dat + Always + + + data\NIST\SmLs03t.dat + Always + + + data\NIST\SmLs04t.dat + Always + + + data\NIST\SmLs05t.dat + Always + + + data\NIST\SmLs06t.dat + Always + + + data\NIST\SmLs07t.dat + Always + + + data\NIST\SmLs08t.dat + Always + + + data\NIST\SmLs09t.dat + Always + + + data\NIST\Thurber.dat + Always + + + data\NIST\Wampler1.dat + Always + + + data\NIST\Wampler2.dat + Always + + + data\NIST\Wampler3.dat + Always + + + data\NIST\Wampler4.dat + Always + + + data\NIST\Wampler5.dat + Always + + + + data\NIST\Meixner.dat + Always + + + + + + {b7cae5f4-a23f-4438-b5be-41226618b695} + Numerics + + + + + + ..\..\packages\NUnit\lib\nunit.framework.dll + True + True + + + + + + \ No newline at end of file diff --git a/src/UnitTests/UseLinearAlgebraProvider.cs b/src/UnitTests/UseLinearAlgebraProvider.cs index fc4a42d6..9476d5bc 100644 --- a/src/UnitTests/UseLinearAlgebraProvider.cs +++ b/src/UnitTests/UseLinearAlgebraProvider.cs @@ -38,8 +38,12 @@ namespace MathNet.Numerics.UnitTests { public void BeforeTest(TestDetails testDetails) { -#if !NET35 && NATIVE +#if !NET35 +#if NATIVEMKL Control.UseNativeMKL(); +#elif NATIVEOPENBLAS + Control.UseNativeOpenBLAS(); +#endif #endif }