Browse Source

OpenBlasLinearAlgerbraProvider implementation. It is basically just a copy of the MklLinearAlgebraProvider minus the vector methods not included in OpenBLAS. The OpenBlasLinearAlgerbraProvider replaces the GotoBlasLinearAlgebraProvider.

Added a new unit test project which is just a copy of the UnitTests-MKL project but with different compilation symbols (NATIVEMKL vs NATIVEOPENBLAS) to control which native provider is loaded. LinearAlgebraProviderTests should probably be in its own project.
pull/312/head
Kuan Bartel 11 years ago
committed by Christoph Ruegg
parent
commit
16851ac48d
  1. 20
      MathNet.Numerics.NativeProviders.sln
  2. 77
      src/NativeProviders/OpenBLAS/capabilities.cpp
  3. 39
      src/NativeProviders/OpenBLAS/complex.h
  4. 44
      src/NativeProviders/OpenBLAS/lapack.cpp
  5. 22
      src/NativeProviders/Windows/OpenBLAS/OpenBLASWrapper.vcxproj
  6. 8
      src/NativeProviders/Windows/OpenBLAS/OpenBLASWrapper.vcxproj.filters
  7. 6
      src/Numerics/Control.cs
  8. 13
      src/Numerics/Numerics.csproj
  9. 1
      src/Numerics/Properties/AssemblyInfo.cs
  10. 367
      src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.cs
  11. 263
      src/Numerics/Providers/LinearAlgebra/GotoBlas/SafeNativeMethods.cs
  12. 436
      src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex.cs
  13. 437
      src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex32.cs
  14. 437
      src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Double.cs
  15. 437
      src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Single.cs
  16. 107
      src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.cs
  17. 303
      src/Numerics/Providers/LinearAlgebra/OpenBlas/SafeNativeMethods.cs
  18. 55
      src/Numerics/Providers/LinearAlgebra/ProviderCapabilities.cs
  19. 7
      src/UnitTests/UnitTests-MKL.csproj
  20. 351
      src/UnitTests/UnitTests-OpenBLAS.csproj
  21. 6
      src/UnitTests/UseLinearAlgebraProvider.cs

20
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

77
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 */

39
src/NativeProviders/OpenBLAS/complex.h

@ -0,0 +1,39 @@
template <typename _T>
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<typename _Other> inline
complex& operator=(const complex<_Other>& right)
{
real = (_T)right.real;
imag = (_T)right.imag;
return *this;
}
};

44
src/NativeProviders/OpenBLAS/lapack.cpp

@ -1,47 +1,7 @@
#include "cblas.h"
#include "complex.h"
#define LAPACK_COMPLEX_CUSTOM
template <typename T>
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<typename _Other> inline
complex& operator=(const complex<_Other>& right)
{
real = (T)right.real;
imag = (T)right.imag;
return *this;
}
};
#define lapack_complex_float complex<float>
#define lapack_complex_double complex<double>
@ -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;

22
src/NativeProviders/Windows/OpenBLAS/OpenBLASWrapper.vcxproj

@ -62,9 +62,9 @@
<Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
</ImportGroup>
<PropertyGroup Label="UserMacros">
<OpenBLASIncludeDir>$(ProjectDir)..\..\..\..\..\libs\OpenBLAS\include</OpenBLASIncludeDir>
<OpenBLASLibDir>$(ProjectDir)..\..\..\..\..\libs\OpenBLAS\x86</OpenBLASLibDir>
<OpenBLASLibDir_x64>$(ProjectDir)..\..\..\..\..\libs\OpenBLAS\x64</OpenBLASLibDir_x64>
<OpenBLASIncludeDir>$(ProjectDir)..\..\..\..\..\libs\OpenBLAS\include\</OpenBLASIncludeDir>
<OpenBLASLibDir>$(ProjectDir)..\..\..\..\..\libs\OpenBLAS\x86\</OpenBLASLibDir>
<OpenBLASLibDir_x64>$(ProjectDir)..\..\..\..\..\libs\OpenBLAS\x64\</OpenBLASLibDir_x64>
</PropertyGroup>
<PropertyGroup>
<_ProjectFileVersion>11.0.50727.1</_ProjectFileVersion>
@ -114,8 +114,7 @@
<AdditionalLibraryDirectories>$(OpenBLASLibDir);%(AdditionalLibraryDirectories)</AdditionalLibraryDirectories>
</Link>
<PostBuildEvent>
<Command>
</Command>
<Command>copy "$(OpenBLASLibDir)*.dll" $(OutputPath)</Command>
</PostBuildEvent>
</ItemDefinitionGroup>
<ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">
@ -142,8 +141,7 @@
<AdditionalLibraryDirectories>$(OpenBLASLibDir_x64);%(AdditionalLibraryDirectories)</AdditionalLibraryDirectories>
</Link>
<PostBuildEvent>
<Command>
</Command>
<Command>copy "$(OpenBLASLibDir_x64)*.dll" $(OutputPath)</Command>
</PostBuildEvent>
</ItemDefinitionGroup>
<ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">
@ -170,8 +168,7 @@
<AdditionalLibraryDirectories>$(OpenBLASLibDir);%(AdditionalLibraryDirectories)</AdditionalLibraryDirectories>
</Link>
<PostBuildEvent>
<Command>
</Command>
<Command>copy "$(OpenBLASLibDir)*.dll" $(OutputPath)</Command>
</PostBuildEvent>
</ItemDefinitionGroup>
<ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
@ -201,18 +198,21 @@
<AdditionalLibraryDirectories>$(OpenBLASLibDir_x64);%(AdditionalLibraryDirectories)</AdditionalLibraryDirectories>
</Link>
<PostBuildEvent>
<Command>
</Command>
<Command>copy "$(OpenBLASLibDir_x64)*.dll" $(OutputPath)</Command>
</PostBuildEvent>
</ItemDefinitionGroup>
<ItemGroup>
<ClCompile Include="..\..\Common\WindowsDLL.cpp" />
<ClCompile Include="..\..\OpenBLAS\blas.c" />
<ClCompile Include="..\..\OpenBLAS\capabilities.cpp" />
<ClCompile Include="..\..\OpenBLAS\lapack.cpp" />
</ItemGroup>
<ItemGroup>
<ResourceCompile Include="..\..\Common\resource.rc" />
</ItemGroup>
<ItemGroup>
<ClInclude Include="..\..\OpenBLAS\complex.h" />
</ItemGroup>
<Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />
<ImportGroup Label="ExtensionTargets">
</ImportGroup>

8
src/NativeProviders/Windows/OpenBLAS/OpenBLASWrapper.vcxproj.filters

@ -24,10 +24,18 @@
<ClCompile Include="..\..\Common\WindowsDLL.cpp">
<Filter>Source Files</Filter>
</ClCompile>
<ClCompile Include="..\..\OpenBLAS\capabilities.cpp">
<Filter>Source Files</Filter>
</ClCompile>
</ItemGroup>
<ItemGroup>
<ResourceCompile Include="..\..\Common\resource.rc">
<Filter>Resource Files</Filter>
</ResourceCompile>
</ItemGroup>
<ItemGroup>
<ClInclude Include="..\..\OpenBLAS\complex.h">
<Filter>Header Files</Filter>
</ClInclude>
</ItemGroup>
</Project>

6
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
/// <summary>

13
src/Numerics/Numerics.csproj

@ -157,12 +157,12 @@
<Compile Include="Providers\LinearAlgebra\Acml\AcmlLinearAlgebraProvider.Double.cs" />
<Compile Include="Providers\LinearAlgebra\Acml\AcmlLinearAlgebraProvider.Single.cs" />
<Compile Include="Providers\LinearAlgebra\Acml\SafeNativeMethods.cs" />
<Compile Include="Providers\LinearAlgebra\GotoBlas\GotoBlasLinearAlgebraProvider.cs" />
<Compile Include="Providers\LinearAlgebra\GotoBlas\GotoBlasLinearAlgebraProvider.Complex.cs" />
<Compile Include="Providers\LinearAlgebra\GotoBlas\GotoBlasLinearAlgebraProvider.Complex32.cs" />
<Compile Include="Providers\LinearAlgebra\GotoBlas\GotoBlasLinearAlgebraProvider.Double.cs" />
<Compile Include="Providers\LinearAlgebra\GotoBlas\GotoBlasLinearAlgebraProvider.Single.cs" />
<Compile Include="Providers\LinearAlgebra\GotoBlas\SafeNativeMethods.cs" />
<Compile Include="Providers\LinearAlgebra\OpenBlas\OpenBlasLinearAlgebraProvider.cs" />
<Compile Include="Providers\LinearAlgebra\OpenBlas\OpenBlasLinearAlgebraProvider.Complex.cs" />
<Compile Include="Providers\LinearAlgebra\OpenBlas\OpenBlasLinearAlgebraProvider.Complex32.cs" />
<Compile Include="Providers\LinearAlgebra\OpenBlas\OpenBlasLinearAlgebraProvider.Double.cs" />
<Compile Include="Providers\LinearAlgebra\OpenBlas\OpenBlasLinearAlgebraProvider.Single.cs" />
<Compile Include="Providers\LinearAlgebra\OpenBlas\SafeNativeMethods.cs" />
<Compile Include="Providers\LinearAlgebra\ManagedLinearAlgebraProvider.Complex32.cs" />
<Compile Include="Providers\LinearAlgebra\ManagedLinearAlgebraProvider.Complex.cs" />
<Compile Include="Providers\LinearAlgebra\ManagedLinearAlgebraProvider.Double.cs" />
@ -192,6 +192,7 @@
<Compile Include="LinearAlgebra\Vector.BCL.cs" />
<Compile Include="LinearAlgebra\Vector.Operators.cs" />
<Compile Include="Exceptions.cs" />
<Compile Include="Providers\LinearAlgebra\ProviderCapabilities.cs" />
<Compile Include="Providers\NativeProviderLoader.cs" />
<Compile Include="Random\SystemRandomSource.cs" />
<Compile Include="Random\RandomSeed.cs" />

1
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

367
src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.cs

@ -1,367 +0,0 @@
// <copyright file="GotoBlasLinearAlgebraProvider.cs" company="Math.NET">
// 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.
// </copyright>
#if NATIVEGOTO
using MathNet.Numerics.Properties;
using System;
using System.Numerics;
using System.Security;
namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas
{
/// <summary>
/// GotoBLAS2 linear algebra provider.
/// </summary>
public partial class GotoBlasLinearAlgebraProvider : ManagedLinearAlgebraProvider
{
/// <summary>
/// Computes the requested <see cref="Norm"/> of the matrix.
/// </summary>
/// <param name="norm">The type of norm to compute.</param>
/// <param name="rows">The number of rows in the matrix.</param>
/// <param name="columns">The number of columns in the matrix.</param>
/// <param name="matrix">The matrix to compute the norm from.</param>
/// <returns>
/// The requested <see cref="Norm"/> of the matrix.
/// </returns>
[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);
}
/// <summary>
/// Computes the requested <see cref="Norm"/> of the matrix.
/// </summary>
/// <param name="norm">The type of norm to compute.</param>
/// <param name="rows">The number of rows in the matrix.</param>
/// <param name="columns">The number of columns in the matrix.</param>
/// <param name="matrix">The matrix to compute the norm from.</param>
/// <param name="work">The work array. Only used when <see cref="Norm.InfinityNorm"/>
/// and needs to be have a length of at least M (number of rows of <paramref name="matrix"/>.</param>
/// <returns>
/// The requested <see cref="Norm"/> of the matrix.
/// </returns>
[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);
}
/// <summary>
/// Computes the requested <see cref="Norm"/> of the matrix.
/// </summary>
/// <param name="norm">The type of norm to compute.</param>
/// <param name="rows">The number of rows in the matrix.</param>
/// <param name="columns">The number of columns in the matrix.</param>
/// <param name="matrix">The matrix to compute the norm from.</param>
/// <returns>
/// The requested <see cref="Norm"/> of the matrix.
/// </returns>
[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);
}
/// <summary>
/// Computes the requested <see cref="Norm"/> of the matrix.
/// </summary>
/// <param name="norm">The type of norm to compute.</param>
/// <param name="rows">The number of rows in the matrix.</param>
/// <param name="columns">The number of columns in the matrix.</param>
/// <param name="matrix">The matrix to compute the norm from.</param>
/// <param name="work">The work array. Only used when <see cref="Norm.InfinityNorm"/>
/// and needs to be have a length of at least M (number of rows of <paramref name="matrix"/>.</param>
/// <returns>
/// The requested <see cref="Norm"/> of the matrix.
/// </returns>
[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);
}
/// <summary>
/// Computes the requested <see cref="Norm"/> of the matrix.
/// </summary>
/// <param name="norm">The type of norm to compute.</param>
/// <param name="rows">The number of rows in the matrix.</param>
/// <param name="columns">The number of columns in the matrix.</param>
/// <param name="matrix">The matrix to compute the norm from.</param>
/// <returns>
/// The requested <see cref="Norm"/> of the matrix.
/// </returns>
[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);
}
/// <summary>
/// Computes the requested <see cref="Norm"/> of the matrix.
/// </summary>
/// <param name="norm">The type of norm to compute.</param>
/// <param name="rows">The number of rows in the matrix.</param>
/// <param name="columns">The number of columns in the matrix.</param>
/// <param name="matrix">The matrix to compute the norm from.</param>
/// <param name="work">The work array. Only used when <see cref="Norm.InfinityNorm"/>
/// and needs to be have a length of at least M (number of rows of <paramref name="matrix"/>.</param>
/// <returns>
/// The requested <see cref="Norm"/> of the matrix.
/// </returns>
[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);
}
/// <summary>
/// Computes the requested <see cref="Norm"/> of the matrix.
/// </summary>
/// <param name="norm">The type of norm to compute.</param>
/// <param name="rows">The number of rows in the matrix.</param>
/// <param name="columns">The number of columns in the matrix.</param>
/// <param name="matrix">The matrix to compute the norm from.</param>
/// <returns>
/// The requested <see cref="Norm"/> of the matrix.
/// </returns>
[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);
}
/// <summary>
/// Computes the requested <see cref="Norm"/> of the matrix.
/// </summary>
/// <param name="norm">The type of norm to compute.</param>
/// <param name="rows">The number of rows in the matrix.</param>
/// <param name="columns">The number of columns in the matrix.</param>
/// <param name="matrix">The matrix to compute the norm from.</param>
/// <param name="work">The work array. Only used when <see cref="Norm.InfinityNorm"/>
/// and needs to be have a length of at least M (number of rows of <paramref name="matrix"/>.</param>
/// <returns>
/// The requested <see cref="Norm"/> of the matrix.
/// </returns>
[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

263
src/Numerics/Providers/LinearAlgebra/GotoBlas/SafeNativeMethods.cs

@ -1,263 +0,0 @@
// <copyright file="SafeNativeMethods.cs" company="Math.NET">
// 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.
// </copyright>
#if NATIVEGOTO
using System.Numerics;
using System.Runtime.InteropServices;
using System.Security;
namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas
{
/// <summary>
/// P/Invoke methods to the native math libraries.
/// </summary>
[SuppressUnmanagedCodeSecurity]
[SecurityCritical]
internal static class SafeNativeMethods
{
/// <summary>
/// Name of the native DLL.
/// </summary>
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

436
src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Complex.cs → src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex.cs

@ -1,10 +1,10 @@
// <copyright file="GotoBlasLinearAlgebraProvider.Complex.cs" company="Math.NET">
// <copyright file="OpenBlasLinearAlgebraProvider.Complex.cs" company="Math.NET">
// 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.
// </copyright>
#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
{
/// <summary>
/// GotoBLAS2 linear algebra provider.
/// OpenBLAS linear algebra provider.
/// </summary>
public partial class GotoBlasLinearAlgebraProvider
public partial class OpenBlasLinearAlgebraProvider
{
/// <summary>
/// Computes the requested <see cref="Norm"/> of the matrix.
/// </summary>
/// <param name="norm">The type of norm to compute.</param>
/// <param name="rows">The number of rows in the matrix.</param>
/// <param name="columns">The number of columns in the matrix.</param>
/// <param name="matrix">The matrix to compute the norm from.</param>
/// <returns>
/// The requested <see cref="Norm"/> of the matrix.
/// </returns>
[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);
}
/// <summary>
/// Computes the dot product of x and y.
/// </summary>
/// <param name="x">The vector x.</param>
/// <param name="y">The vector y.</param>
/// <returns>The dot product of x and y.</returns>
/// <remarks>This is equivalent to the DOT BLAS routine.</remarks>
[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);
}
/// <summary>
/// Adds a scaled vector to another: <c>result = y + alpha*x</c>.
/// </summary>
@ -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
}
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
/// Computes the thin QR factorization of A where M &gt; N.
/// </summary>
/// <param name="a">The A matrix.</param>
/// <param name="rows">The number of rows in the A matrix.</param>
/// <param name="columns">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
public override void QRSolve(Complex[] a, int rows, int columns, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full)
/// <param name="q">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.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="r">On exit, A N by N matrix that holds the R matrix of the
/// QR factorization.</param>
/// <param name="tau">A min(m,n) vector. On exit, contains additional information
/// to be used by the QR solve routine.</param>
/// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
[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);
}
/// <summary>
/// Computes the thin QR factorization of A where M &gt; N.
/// </summary>
/// <param name="q">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.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="r">On exit, A N by N matrix that holds the R matrix of the
/// QR factorization.</param>
/// <param name="tau">A min(m,n) vector. On exit, contains additional information
/// to be used by the QR solve routine.</param>
/// <param name="work">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.</param>
/// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
[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);
}
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
/// </summary>
/// <param name="a">The A matrix.</param>
/// <param name="rows">The number of rows in the A matrix.</param>
/// <param name="columns">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
[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);
}
/// <summary>
@ -695,7 +829,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas
/// <param name="work">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.</param>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
[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
/// <summary>
/// Solves A*X=B for X using a previously QR factored matrix.
/// </summary>
/// <param name="q">The Q matrix obtained by calling <see cref="QRFactor(Complex[],int,int,Complex[],Complex[],QRMethod)"/>.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(Complex[],int,int,Complex[],Complex[],QRMethod)"/>. </param>
/// <param name="q">The Q matrix obtained by calling <see cref="QRFactor(Complex[],int,int,Complex[],Complex[])"/>.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(Complex[],int,int,Complex[],Complex[])"/>. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="tau">Contains additional information on Q. Only used for the native solver
@ -759,57 +895,13 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
[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);
}
/// <summary>
@ -817,9 +909,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas
/// </summary>
/// <param name="q">The Q matrix obtained by QR factor. This is only used for the managed provider and can be
/// <c>null</c> for the native provider. The native provider uses the Q portion stored in the R matrix.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(Complex[],int,int,Complex[],Complex[],QRMethod)"/>. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(Complex[],int,int,Complex[],Complex[])"/>. </param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="tau">Contains additional information on Q. Only used for the native solver
/// and can be <c>null</c> for the managed provider.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
@ -828,8 +920,10 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas
/// <param name="work">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.</param>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
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);
}
}
/// <summary>
@ -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();
}
}
/// <summary>
/// Computes the eigenvalues and eigenvectors of a matrix.
/// </summary>
/// <param name="isSymmetric">Whether the matrix is symmetric or not.</param>
/// <param name="order">The order of the matrix.</param>
/// <param name="matrix">The matrix to decompose. The lenth of the array must be order * order.</param>
/// <param name="matrixEv">On output, the matrix contains the eigen vectors. The lenth of the array must be order * order.</param>
/// <param name="vectorEv">On output, the eigen values (λ) of matrix in ascending value. The length of the arry must <paramref name="order"/>.</param>
/// <param name="matrixD">On output, the block diagonal eigenvalue matrix. The lenth of the array must be order * order.</param>
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();
}
}
}
}

437
src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Complex32.cs → src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex32.cs

@ -1,10 +1,10 @@
// <copyright file="GotoBlasLinearAlgebraProvider.Complex32.cs" company="Math.NET">
// <copyright file="OpenBlasLinearAlgebraProvider.Complex32.cs" company="Math.NET">
// 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.
// </copyright>
#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
{
/// <summary>
/// GotoBLAS2 linear algebra provider.
/// OpenBLAS linear algebra provider.
/// </summary>
public partial class GotoBlasLinearAlgebraProvider
public partial class OpenBlasLinearAlgebraProvider
{
/// <summary>
/// Computes the requested <see cref="Norm"/> of the matrix.
/// </summary>
/// <param name="norm">The type of norm to compute.</param>
/// <param name="rows">The number of rows in the matrix.</param>
/// <param name="columns">The number of columns in the matrix.</param>
/// <param name="matrix">The matrix to compute the norm from.</param>
/// <returns>
/// The requested <see cref="Norm"/> of the matrix.
/// </returns>
[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);
}
/// <summary>
/// Computes the dot product of x and y.
/// </summary>
/// <param name="x">The vector x.</param>
/// <param name="y">The vector y.</param>
/// <returns>The dot product of x and y.</returns>
/// <remarks>This is equivalent to the DOT BLAS routine.</remarks>
[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);
}
/// <summary>
/// Adds a scaled vector to another: <c>result = y + alpha*x</c>.
/// </summary>
@ -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
}
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
/// Computes the thin QR factorization of A where M &gt; N.
/// </summary>
/// <param name="a">The A matrix.</param>
/// <param name="rows">The number of rows in the A matrix.</param>
/// <param name="columns">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
public override void QRSolve(Complex32[] a, int rows, int columns, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full)
/// <param name="q">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.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="r">On exit, A N by N matrix that holds the R matrix of the
/// QR factorization.</param>
/// <param name="tau">A min(m,n) vector. On exit, contains additional information
/// to be used by the QR solve routine.</param>
/// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
[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);
}
/// <summary>
/// Computes the thin QR factorization of A where M &gt; N.
/// </summary>
/// <param name="q">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.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="r">On exit, A N by N matrix that holds the R matrix of the
/// QR factorization.</param>
/// <param name="tau">A min(m,n) vector. On exit, contains additional information
/// to be used by the QR solve routine.</param>
/// <param name="work">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.</param>
/// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
[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);
}
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
/// </summary>
/// <param name="a">The A matrix.</param>
/// <param name="rows">The number of rows in the A matrix.</param>
/// <param name="columns">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
[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);
}
/// <summary>
@ -694,7 +829,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas
/// <param name="work">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.</param>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
[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
/// <summary>
/// Solves A*X=B for X using a previously QR factored matrix.
/// </summary>
/// <param name="q">The Q matrix obtained by calling <see cref="QRFactor(Complex32[],int,int,Complex32[],Complex32[],QRMethod)"/>.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(Complex32[],int,int,Complex32[],Complex32[],QRMethod)"/>. </param>
/// <param name="q">The Q matrix obtained by calling <see cref="QRFactor(Complex32[],int,int,Complex32[],Complex32[])"/>.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(Complex32[],int,int,Complex32[],Complex32[])"/>. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="tau">Contains additional information on Q. Only used for the native solver
@ -758,57 +895,13 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
[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);
}
/// <summary>
@ -816,9 +909,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas
/// </summary>
/// <param name="q">The Q matrix obtained by QR factor. This is only used for the managed provider and can be
/// <c>null</c> for the native provider. The native provider uses the Q portion stored in the R matrix.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(Complex32[],int,int,Complex32[],Complex32[],QRMethod)"/>. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(Complex32[],int,int,Complex32[],Complex32[])"/>. </param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="tau">Contains additional information on Q. Only used for the native solver
/// and can be <c>null</c> for the managed provider.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
@ -827,8 +920,10 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas
/// <param name="work">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.</param>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
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);
}
}
/// <summary>
@ -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();
}
}
/// <summary>
/// Computes the eigenvalues and eigenvectors of a matrix.
/// </summary>
/// <param name="isSymmetric">Whether the matrix is symmetric or not.</param>
/// <param name="order">The order of the matrix.</param>
/// <param name="matrix">The matrix to decompose. The lenth of the array must be order * order.</param>
/// <param name="matrixEv">On output, the matrix contains the eigen vectors. The lenth of the array must be order * order.</param>
/// <param name="vectorEv">On output, the eigen values (λ) of matrix in ascending value. The length of the arry must <paramref name="order"/>.</param>
/// <param name="matrixD">On output, the block diagonal eigenvalue matrix. The lenth of the array must be order * order.</param>
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();
}
}
}
}

437
src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Double.cs → src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Double.cs

@ -1,10 +1,10 @@
// <copyright file="GotoBlasLinearAlgebraProvider.Double.cs" company="Math.NET">
// <copyright file="OpenBlasLinearAlgebraProvider.Double.cs" company="Math.NET">
// 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.
// </copyright>
#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
{
/// <summary>
/// GotoBLAS2 linear algebra provider.
/// OpenBLAS linear algebra provider.
/// </summary>
public partial class GotoBlasLinearAlgebraProvider
public partial class OpenBlasLinearAlgebraProvider
{
/// <summary>
/// Computes the requested <see cref="Norm"/> of the matrix.
/// </summary>
/// <param name="norm">The type of norm to compute.</param>
/// <param name="rows">The number of rows in the matrix.</param>
/// <param name="columns">The number of columns in the matrix.</param>
/// <param name="matrix">The matrix to compute the norm from.</param>
/// <returns>
/// The requested <see cref="Norm"/> of the matrix.
/// </returns>
[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);
}
/// <summary>
/// Computes the dot product of x and y.
/// </summary>
/// <param name="x">The vector x.</param>
/// <param name="y">The vector y.</param>
/// <returns>The dot product of x and y.</returns>
/// <remarks>This is equivalent to the DOT BLAS routine.</remarks>
[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);
}
/// <summary>
/// Adds a scaled vector to another: <c>result = y + alpha*x</c>.
/// </summary>
@ -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
}
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
/// Computes the thin QR factorization of A where M &gt; N.
/// </summary>
/// <param name="a">The A matrix.</param>
/// <param name="rows">The number of rows in the A matrix.</param>
/// <param name="columns">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
public override void QRSolve(double[] a, int rows, int columns, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full)
/// <param name="q">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.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="r">On exit, A N by N matrix that holds the R matrix of the
/// QR factorization.</param>
/// <param name="tau">A min(m,n) vector. On exit, contains additional information
/// to be used by the QR solve routine.</param>
/// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
[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);
}
/// <summary>
/// Computes the thin QR factorization of A where M &gt; N.
/// </summary>
/// <param name="q">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.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="r">On exit, A N by N matrix that holds the R matrix of the
/// QR factorization.</param>
/// <param name="tau">A min(m,n) vector. On exit, contains additional information
/// to be used by the QR solve routine.</param>
/// <param name="work">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.</param>
/// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
[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);
}
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
/// </summary>
/// <param name="a">The A matrix.</param>
/// <param name="rows">The number of rows in the A matrix.</param>
/// <param name="columns">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
[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);
}
/// <summary>
@ -694,7 +829,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas
/// <param name="work">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.</param>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
[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
/// <summary>
/// Solves A*X=B for X using a previously QR factored matrix.
/// </summary>
/// <param name="q">The Q matrix obtained by calling <see cref="QRFactor(double[],int,int,double[],double[],QRMethod)"/>.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(double[],int,int,double[],double[],QRMethod)"/>. </param>
/// <param name="q">The Q matrix obtained by calling <see cref="QRFactor(double[],int,int,double[],double[])"/>.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(double[],int,int,double[],double[])"/>. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="tau">Contains additional information on Q. Only used for the native solver
@ -758,57 +895,13 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
[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);
}
/// <summary>
@ -816,9 +909,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas
/// </summary>
/// <param name="q">The Q matrix obtained by QR factor. This is only used for the managed provider and can be
/// <c>null</c> for the native provider. The native provider uses the Q portion stored in the R matrix.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(double[],int,int,double[],double[],QRMethod)"/>. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(double[],int,int,double[],double[])"/>. </param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="tau">Contains additional information on Q. Only used for the native solver
/// and can be <c>null</c> for the managed provider.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
@ -827,8 +920,10 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas
/// <param name="work">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.</param>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
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);
}
}
/// <summary>
@ -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();
}
}
/// <summary>
/// Computes the eigenvalues and eigenvectors of a matrix.
/// </summary>
/// <param name="isSymmetric">Whether the matrix is symmetric or not.</param>
/// <param name="order">The order of the matrix.</param>
/// <param name="matrix">The matrix to decompose. The lenth of the array must be order * order.</param>
/// <param name="matrixEv">On output, the matrix contains the eigen vectors. The lenth of the array must be order * order.</param>
/// <param name="vectorEv">On output, the eigen values (λ) of matrix in ascending value. The length of the arry must <paramref name="order"/>.</param>
/// <param name="matrixD">On output, the block diagonal eigenvalue matrix. The lenth of the array must be order * order.</param>
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();
}
}
}
}

437
src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Single.cs → src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Single.cs

@ -1,10 +1,10 @@
// <copyright file="GotoBlasLinearAlgebraProvider.Single.cs" company="Math.NET">
// <copyright file="OpenBlasLinearAlgebraProvider.Single.cs" company="Math.NET">
// 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.
// </copyright>
#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
{
/// <summary>
/// GotoBLAS2 linear algebra provider.
/// OpenBLAS linear algebra provider.
/// </summary>
public partial class GotoBlasLinearAlgebraProvider
public partial class OpenBlasLinearAlgebraProvider
{
/// <summary>
/// Computes the requested <see cref="Norm"/> of the matrix.
/// </summary>
/// <param name="norm">The type of norm to compute.</param>
/// <param name="rows">The number of rows in the matrix.</param>
/// <param name="columns">The number of columns in the matrix.</param>
/// <param name="matrix">The matrix to compute the norm from.</param>
/// <returns>
/// The requested <see cref="Norm"/> of the matrix.
/// </returns>
[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);
}
/// <summary>
/// Computes the dot product of x and y.
/// </summary>
/// <param name="x">The vector x.</param>
/// <param name="y">The vector y.</param>
/// <returns>The dot product of x and y.</returns>
/// <remarks>This is equivalent to the DOT BLAS routine.</remarks>
[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);
}
/// <summary>
/// Adds a scaled vector to another: <c>result = y + alpha*x</c>.
/// </summary>
@ -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
}
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
/// Computes the thin QR factorization of A where M &gt; N.
/// </summary>
/// <param name="a">The A matrix.</param>
/// <param name="rows">The number of rows in the A matrix.</param>
/// <param name="columns">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
public override void QRSolve(float[] a, int rows, int columns, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full)
/// <param name="q">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.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="r">On exit, A N by N matrix that holds the R matrix of the
/// QR factorization.</param>
/// <param name="tau">A min(m,n) vector. On exit, contains additional information
/// to be used by the QR solve routine.</param>
/// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
[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);
}
/// <summary>
/// Computes the thin QR factorization of A where M &gt; N.
/// </summary>
/// <param name="q">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.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="r">On exit, A N by N matrix that holds the R matrix of the
/// QR factorization.</param>
/// <param name="tau">A min(m,n) vector. On exit, contains additional information
/// to be used by the QR solve routine.</param>
/// <param name="work">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.</param>
/// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
[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);
}
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
/// </summary>
/// <param name="a">The A matrix.</param>
/// <param name="rows">The number of rows in the A matrix.</param>
/// <param name="columns">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
[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);
}
/// <summary>
@ -694,7 +829,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas
/// <param name="work">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.</param>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
[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
/// <summary>
/// Solves A*X=B for X using a previously QR factored matrix.
/// </summary>
/// <param name="q">The Q matrix obtained by calling <see cref="QRFactor(float[],int,int,float[],float[],QRMethod)"/>.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(float[],int,int,float[],float[],QRMethod)"/>. </param>
/// <param name="q">The Q matrix obtained by calling <see cref="QRFactor(float[],int,int,float[],float[])"/>.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(float[],int,int,float[],float[])"/>. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="tau">Contains additional information on Q. Only used for the native solver
@ -758,57 +895,13 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
[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);
}
/// <summary>
@ -816,9 +909,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas
/// </summary>
/// <param name="q">The Q matrix obtained by QR factor. This is only used for the managed provider and can be
/// <c>null</c> for the native provider. The native provider uses the Q portion stored in the R matrix.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(float[],int,int,float[],float[],QRMethod)"/>. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(float[],int,int,float[],float[])"/>. </param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="tau">Contains additional information on Q. Only used for the native solver
/// and can be <c>null</c> for the managed provider.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
@ -827,8 +920,10 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas
/// <param name="work">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.</param>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
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);
}
}
/// <summary>
@ -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();
}
}
/// <summary>
/// Computes the eigenvalues and eigenvectors of a matrix.
/// </summary>
/// <param name="isSymmetric">Whether the matrix is symmetric or not.</param>
/// <param name="order">The order of the matrix.</param>
/// <param name="matrix">The matrix to decompose. The lenth of the array must be order * order.</param>
/// <param name="matrixEv">On output, the matrix contains the eigen vectors. The lenth of the array must be order * order.</param>
/// <param name="vectorEv">On output, the eigen values (λ) of matrix in ascending value. The length of the arry must <paramref name="order"/>.</param>
/// <param name="matrixD">On output, the block diagonal eigenvalue matrix. The lenth of the array must be order * order.</param>
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();
}
}
}
}

107
src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.cs

@ -0,0 +1,107 @@
// <copyright file="OpenBlasLinearAlgebraProvider.cs" company="Math.NET">
// 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.
// </copyright>
#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
}
/// <summary>
/// OpenBLAS linear algebra provider.
/// </summary>
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

303
src/Numerics/Providers/LinearAlgebra/OpenBlas/SafeNativeMethods.cs

@ -0,0 +1,303 @@
// <copyright file="SafeNativeMethods.cs" company="Math.NET">
// 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.
// </copyright>
#if NATIVE
using System.Numerics;
using System.Runtime.InteropServices;
using System.Security;
namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas
{
/// <summary>
/// P/Invoke methods to the native math libraries.
/// </summary>
[SuppressUnmanagedCodeSecurity]
[SecurityCritical]
internal static class SafeNativeMethods
{
/// <summary>
/// Name of the native DLL.
/// </summary>
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

55
src/Numerics/Providers/LinearAlgebra/ProviderCapabilities.cs

@ -0,0 +1,55 @@
// <copyright file="ProviderCapabilities.cs" company="Math.NET">
// 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.
// </copyright>
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,
}
}

7
src/UnitTests/UnitTests-MKL.csproj

@ -16,7 +16,7 @@
<SolutionDir Condition="$(SolutionDir) == '' Or $(SolutionDir) == '*Undefined*'">..\..\</SolutionDir>
</PropertyGroup>
<PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Release|AnyCPU' ">
<DefineConstants>TRACE;NATIVE</DefineConstants>
<DefineConstants>TRACE;NATIVE;NATIVEMKL;</DefineConstants>
<OutputPath>..\..\out\MKL\Windows\</OutputPath>
<IntermediateOutputPath>..\..\obj\MKL\Windows\x86\</IntermediateOutputPath>
<BaseIntermediateOutputPath>..\..\obj\MKL\Windows\x86\</BaseIntermediateOutputPath>
@ -28,7 +28,7 @@
<PlatformTarget>AnyCPU</PlatformTarget>
</PropertyGroup>
<PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Debug|AnyCPU' ">
<DefineConstants>TRACE;DEBUG;NATIVE</DefineConstants>
<DefineConstants>TRACE;DEBUG;NATIVE;NATIVEMKL;</DefineConstants>
<OutputPath>..\..\out\MKL\Windows\</OutputPath>
<IntermediateOutputPath>..\..\obj\MKL\Windows\x86\</IntermediateOutputPath>
<BaseIntermediateOutputPath>..\..\obj\MKL\Windows\x86\</BaseIntermediateOutputPath>
@ -345,4 +345,7 @@
<Paket>True</Paket>
</Reference>
</ItemGroup>
<ItemGroup>
<Service Include="{82A7F48D-3B50-4B1E-B82E-3ADA8210C358}" />
</ItemGroup>
</Project>

351
src/UnitTests/UnitTests-OpenBLAS.csproj

@ -0,0 +1,351 @@
<?xml version="1.0" encoding="utf-8"?>
<Project ToolsVersion="4.0" DefaultTargets="Build" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
<PropertyGroup>
<MinimumVisualStudioVersion>10.0</MinimumVisualStudioVersion>
<Configuration Condition=" '$(Configuration)' == '' ">Debug</Configuration>
<Platform Condition=" '$(Platform)' == '' ">AnyCPU</Platform>
<ProductVersion>8.0.30703</ProductVersion>
<SchemaVersion>2.0</SchemaVersion>
<ProjectGuid>{96B903EF-3EE1-4569-803C-0482D2F5ED37}</ProjectGuid>
<OutputType>Library</OutputType>
<AppDesignerFolder>Properties</AppDesignerFolder>
<RootNamespace>MathNet.Numerics.UnitTests</RootNamespace>
<AssemblyName>MathNet.Numerics.UnitTestsOpenBLAS</AssemblyName>
<TargetFrameworkVersion>v4.5</TargetFrameworkVersion>
<FileAlignment>512</FileAlignment>
<SolutionDir Condition="$(SolutionDir) == '' Or $(SolutionDir) == '*Undefined*'">..\..\</SolutionDir>
</PropertyGroup>
<PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Release|AnyCPU' ">
<DefineConstants>TRACE;NATIVE;NATIVEOPENBLAS;</DefineConstants>
<OutputPath>..\..\out\OpenBLAS\Windows\</OutputPath>
<IntermediateOutputPath>..\..\obj\OpenBLAS\Windows\AnyCPU\</IntermediateOutputPath>
<BaseIntermediateOutputPath>..\..\obj\OpenBLAS\Windows\AnyCPU\</BaseIntermediateOutputPath>
<Optimize>true</Optimize>
<DebugType>pdbonly</DebugType>
<ErrorReport>prompt</ErrorReport>
<CodeAnalysisRuleSet>MinimumRecommendedRules.ruleset</CodeAnalysisRuleSet>
<NoWarn>1591</NoWarn>
<PlatformTarget>AnyCPU</PlatformTarget>
</PropertyGroup>
<PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Debug|AnyCPU' ">
<DefineConstants>TRACE;DEBUG;NATIVE;NATIVEOPENBLAS;</DefineConstants>
<OutputPath>..\..\out\OpenBLAS\Windows\</OutputPath>
<IntermediateOutputPath>..\..\obj\OpenBLAS\Windows\AnyCPU\</IntermediateOutputPath>
<BaseIntermediateOutputPath>..\..\obj\OpenBLAS\Windows\AnyCPU\</BaseIntermediateOutputPath>
<Optimize>false</Optimize>
<DebugType>full</DebugType>
<DebugSymbols>true</DebugSymbols>
<ErrorReport>prompt</ErrorReport>
<WarningLevel>4</WarningLevel>
<NoWarn>1591</NoWarn>
<PlatformTarget>AnyCPU</PlatformTarget>
</PropertyGroup>
<ItemGroup>
<Reference Include="System" />
<Reference Include="System" />
<Reference Include="System.Core" />
<Reference Include="System.Numerics" />
<Reference Include="System.Xml.Linq" />
<Reference Include="System.Data.DataSetExtensions" />
<Reference Include="Microsoft.CSharp" />
<Reference Include="System.Data" />
<Reference Include="System.Xml" />
</ItemGroup>
<ItemGroup>
<Compile Include="**\*.cs" Exclude="Properties\Settings.Designer.cs">
</Compile>
</ItemGroup>
<ItemGroup>
<None Include="..\..\data\Codeplex-5667.csv">
<Link>data\Codeplex-5667.csv</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\Github-Cureos-1.csv">
<Link>data\Github-Cureos-1.csv</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\Matlab\A.mat">
<Link>data\Matlab\A.mat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\Matlab\collection-nocompress.mat">
<Link>data\Matlab\collection-nocompress.mat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\Matlab\collection.mat">
<Link>data\Matlab\collection.mat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\Matlab\complex.mat">
<Link>data\Matlab\complex.mat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\Matlab\sparse-large.mat">
<Link>data\Matlab\sparse-large.mat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\Matlab\sparse-small.mat">
<Link>data\Matlab\sparse-small.mat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\Matlab\sparse_complex.mat">
<Link>data\Matlab\sparse_complex.mat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\Matlab\v.mat">
<Link>data\Matlab\v.mat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\AtmWtAgt.dat">
<Link>data\NIST\AtmWtAgt.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Bennett5.dat">
<Link>data\NIST\Bennett5.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\BoxBOD.dat">
<Link>data\NIST\BoxBOD.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Chwirut1.dat">
<Link>data\NIST\Chwirut1.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Chwirut2.dat">
<Link>data\NIST\Chwirut2.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\DanWood.dat">
<Link>data\NIST\DanWood.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Eckerle4.dat">
<Link>data\NIST\Eckerle4.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\ENSO.dat">
<Link>data\NIST\ENSO.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Filip.dat">
<Link>data\NIST\Filip.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Gauss1.dat">
<Link>data\NIST\Gauss1.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Gauss2.dat">
<Link>data\NIST\Gauss2.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Gauss3.dat">
<Link>data\NIST\Gauss3.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Hahn1.dat">
<Link>data\NIST\Hahn1.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Kirby2.dat">
<Link>data\NIST\Kirby2.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Lanczos1.dat">
<Link>data\NIST\Lanczos1.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Lanczos2.dat">
<Link>data\NIST\Lanczos2.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Lanczos3.dat">
<Link>data\NIST\Lanczos3.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Lew.dat">
<Link>data\NIST\Lew.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Longley.dat">
<Link>data\NIST\Longley.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Lottery.dat">
<Link>data\NIST\Lottery.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Mavro.dat">
<Link>data\NIST\Mavro.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\MGH09.dat">
<Link>data\NIST\MGH09.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\MGH10.dat">
<Link>data\NIST\MGH10.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\MGH17.dat">
<Link>data\NIST\MGH17.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Michelso.dat">
<Link>data\NIST\Michelso.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Misra1a.dat">
<Link>data\NIST\Misra1a.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Misra1b.dat">
<Link>data\NIST\Misra1b.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Misra1c.dat">
<Link>data\NIST\Misra1c.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Misra1d.dat">
<Link>data\NIST\Misra1d.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Nelson.dat">
<Link>data\NIST\Nelson.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\NoInt1.dat">
<Link>data\NIST\NoInt1.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\NoInt2.dat">
<Link>data\NIST\NoInt2.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Norris.dat">
<Link>data\NIST\Norris.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\NumAcc1.dat">
<Link>data\NIST\NumAcc1.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\NumAcc2.dat">
<Link>data\NIST\NumAcc2.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\NumAcc3.dat">
<Link>data\NIST\NumAcc3.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\NumAcc4.dat">
<Link>data\NIST\NumAcc4.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Pontius.dat">
<Link>data\NIST\Pontius.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Rat42.dat">
<Link>data\NIST\Rat42.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Rat43.dat">
<Link>data\NIST\Rat43.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Roszman1.dat">
<Link>data\NIST\Roszman1.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\SiRstvt.dat">
<Link>data\NIST\SiRstvt.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\SmLs01t.dat">
<Link>data\NIST\SmLs01t.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\SmLs02t.dat">
<Link>data\NIST\SmLs02t.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\SmLs03t.dat">
<Link>data\NIST\SmLs03t.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\SmLs04t.dat">
<Link>data\NIST\SmLs04t.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\SmLs05t.dat">
<Link>data\NIST\SmLs05t.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\SmLs06t.dat">
<Link>data\NIST\SmLs06t.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\SmLs07t.dat">
<Link>data\NIST\SmLs07t.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\SmLs08t.dat">
<Link>data\NIST\SmLs08t.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\SmLs09t.dat">
<Link>data\NIST\SmLs09t.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Thurber.dat">
<Link>data\NIST\Thurber.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Wampler1.dat">
<Link>data\NIST\Wampler1.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Wampler2.dat">
<Link>data\NIST\Wampler2.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Wampler3.dat">
<Link>data\NIST\Wampler3.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Wampler4.dat">
<Link>data\NIST\Wampler4.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="..\..\data\NIST\Wampler5.dat">
<Link>data\NIST\Wampler5.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="App.config" />
<None Include="..\..\data\NIST\Meixner.dat">
<Link>data\NIST\Meixner.dat</Link>
<CopyToOutputDirectory>Always</CopyToOutputDirectory>
</None>
<None Include="paket.references" />
</ItemGroup>
<ItemGroup>
<ProjectReference Include="..\Numerics\Numerics.csproj">
<Project>{b7cae5f4-a23f-4438-b5be-41226618b695}</Project>
<Name>Numerics</Name>
</ProjectReference>
</ItemGroup>
<Import Project="$(MSBuildToolsPath)\Microsoft.CSharp.targets" />
<ItemGroup>
<Reference Include="nunit.framework">
<HintPath>..\..\packages\NUnit\lib\nunit.framework.dll</HintPath>
<Private>True</Private>
<Paket>True</Paket>
</Reference>
</ItemGroup>
<ItemGroup>
<Service Include="{82A7F48D-3B50-4B1E-B82E-3ADA8210C358}" />
</ItemGroup>
</Project>

6
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
}

Loading…
Cancel
Save