diff --git a/MathNet.Numerics.NativeProviders.sln b/MathNet.Numerics.NativeProviders.sln index 25b52c6c..4e790037 100644 --- a/MathNet.Numerics.NativeProviders.sln +++ b/MathNet.Numerics.NativeProviders.sln @@ -22,6 +22,10 @@ Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "CUDA", "src\NativeProviders EndProject Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "UnitTests-CUDA", "src\UnitTests\UnitTests-CUDA.csproj", "{E79C0395-01DC-4BC9-B86C-ED45790892C5}" 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 @@ -40,6 +44,10 @@ Global Release-MKL|Mixed Platforms = Release-MKL|Mixed Platforms Release-MKL|Win32 = Release-MKL|Win32 Release-MKL|x64 = Release-MKL|x64 + Release-OpenBLAS|Any CPU = Release-OpenBLAS|Any CPU + Release-OpenBLAS|Mixed Platforms = Release-OpenBLAS|Mixed Platforms + Release-OpenBLAS|Win32 = Release-OpenBLAS|Win32 + Release-OpenBLAS|x64 = Release-OpenBLAS|x64 Release-Signed|Any CPU = Release-Signed|Any CPU Release-Signed|Mixed Platforms = Release-Signed|Mixed Platforms Release-Signed|Win32 = Release-Signed|Win32 @@ -72,6 +80,11 @@ Global {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release-MKL|Win32.Build.0 = Release|Win32 {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release-MKL|x64.ActiveCfg = Release|x64 {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release-MKL|x64.Build.0 = Release|x64 + {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release-OpenBLAS|Any CPU.ActiveCfg = Release|Win32 + {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release-OpenBLAS|Mixed Platforms.ActiveCfg = Release|Win32 + {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release-OpenBLAS|Mixed Platforms.Build.0 = Release|Win32 + {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release-OpenBLAS|Win32.ActiveCfg = Release|Win32 + {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release-OpenBLAS|x64.ActiveCfg = Release|x64 {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release-Signed|Any CPU.ActiveCfg = Release|Win32 {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release-Signed|Mixed Platforms.ActiveCfg = Release|Win32 {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release-Signed|Mixed Platforms.Build.0 = Release|Win32 @@ -95,6 +108,10 @@ Global {2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release-MKL|Mixed Platforms.ActiveCfg = Release|Win32 {2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release-MKL|Win32.ActiveCfg = Release|Win32 {2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release-MKL|x64.ActiveCfg = Release|x64 + {2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release-OpenBLAS|Any CPU.ActiveCfg = Release|Win32 + {2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release-OpenBLAS|Mixed Platforms.ActiveCfg = Release|Win32 + {2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release-OpenBLAS|Win32.ActiveCfg = Release|Win32 + {2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release-OpenBLAS|x64.ActiveCfg = Release|x64 {2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release-Signed|Any CPU.ActiveCfg = Release|Win32 {2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release-Signed|Mixed Platforms.ActiveCfg = Release|Win32 {2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release-Signed|Win32.ActiveCfg = Release|Win32 @@ -133,6 +150,14 @@ Global {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-MKL|Win32.Build.0 = Release|Any CPU {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-MKL|x64.ActiveCfg = Release|Any CPU {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-MKL|x64.Build.0 = Release|Any CPU + {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-OpenBLAS|Any CPU.ActiveCfg = Release|Any CPU + {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-OpenBLAS|Any CPU.Build.0 = Release|Any CPU + {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-OpenBLAS|Mixed Platforms.ActiveCfg = Release|Any CPU + {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-OpenBLAS|Mixed Platforms.Build.0 = Release|Any CPU + {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-OpenBLAS|Win32.ActiveCfg = Release|Any CPU + {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-OpenBLAS|Win32.Build.0 = Release|Any CPU + {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-OpenBLAS|x64.ActiveCfg = Release|Any CPU + {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-OpenBLAS|x64.Build.0 = Release|Any CPU {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-Signed|Any CPU.ActiveCfg = Release-Signed|Any CPU {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-Signed|Any CPU.Build.0 = Release-Signed|Any CPU {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-Signed|Mixed Platforms.ActiveCfg = Release-Signed|Any CPU @@ -165,6 +190,12 @@ Global {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-MKL|Win32.Build.0 = Release|Any CPU {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-MKL|x64.ActiveCfg = Release|Any CPU {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-MKL|x64.Build.0 = Release|Any CPU + {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-OpenBLAS|Any CPU.ActiveCfg = Release|Any CPU + {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-OpenBLAS|Any CPU.Build.0 = Release|Any CPU + {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-OpenBLAS|Mixed Platforms.ActiveCfg = Release|Any CPU + {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-OpenBLAS|Mixed Platforms.Build.0 = Release|Any CPU + {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-OpenBLAS|Win32.ActiveCfg = Release|Any CPU + {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-OpenBLAS|x64.ActiveCfg = Release|Any CPU {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-Signed|Any CPU.ActiveCfg = Release|Any CPU {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-Signed|Any CPU.Build.0 = Release|Any CPU {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-Signed|Mixed Platforms.ActiveCfg = Release|Any CPU @@ -197,6 +228,11 @@ Global {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-MKL|Mixed Platforms.Build.0 = Release|Win32 {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-MKL|Win32.ActiveCfg = Release|Win32 {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-MKL|x64.ActiveCfg = Release|x64 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-OpenBLAS|Any CPU.ActiveCfg = Release|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-OpenBLAS|Mixed Platforms.ActiveCfg = Release|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-OpenBLAS|Mixed Platforms.Build.0 = Release|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-OpenBLAS|Win32.ActiveCfg = Release|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-OpenBLAS|x64.ActiveCfg = Release|x64 {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-Signed|Any CPU.ActiveCfg = Release|Win32 {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-Signed|Mixed Platforms.ActiveCfg = Release|Win32 {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-Signed|Mixed Platforms.Build.0 = Release|Win32 @@ -229,12 +265,85 @@ Global {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-MKL|Mixed Platforms.Build.0 = Release|Any CPU {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-MKL|Win32.ActiveCfg = Release|Any CPU {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-MKL|x64.ActiveCfg = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-OpenBLAS|Any CPU.ActiveCfg = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-OpenBLAS|Any CPU.Build.0 = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-OpenBLAS|Mixed Platforms.ActiveCfg = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-OpenBLAS|Mixed Platforms.Build.0 = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-OpenBLAS|Win32.ActiveCfg = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-OpenBLAS|x64.ActiveCfg = Release|Any CPU {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-Signed|Any CPU.ActiveCfg = Release|Any CPU {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-Signed|Any CPU.Build.0 = Release|Any CPU {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-Signed|Mixed Platforms.ActiveCfg = Release|Any CPU {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-Signed|Mixed Platforms.Build.0 = Release|Any CPU {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-Signed|Win32.ActiveCfg = Release|Any CPU {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-Signed|x64.ActiveCfg = Release|Any CPU + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Debug|Any CPU.ActiveCfg = Debug|Win32 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Debug|Mixed Platforms.ActiveCfg = Debug|Win32 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Debug|Mixed Platforms.Build.0 = Debug|Win32 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Debug|Win32.ActiveCfg = Debug|Win32 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Debug|Win32.Build.0 = Debug|Win32 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Debug|x64.ActiveCfg = Debug|x64 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Debug|x64.Build.0 = Debug|x64 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release|Any CPU.ActiveCfg = Release|Win32 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release|Mixed Platforms.ActiveCfg = Release|Win32 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release|Mixed Platforms.Build.0 = Release|Win32 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release|Win32.ActiveCfg = Release|Win32 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release|Win32.Build.0 = Release|Win32 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release|x64.ActiveCfg = Release|x64 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release|x64.Build.0 = Release|x64 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release-CUDA|Any CPU.ActiveCfg = Release|Win32 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release-CUDA|Mixed Platforms.ActiveCfg = Release|Win32 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release-CUDA|Win32.ActiveCfg = Release|Win32 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release-CUDA|x64.ActiveCfg = Release|x64 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release-MKL|Any CPU.ActiveCfg = Release|Win32 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release-MKL|Mixed Platforms.ActiveCfg = Release|Win32 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release-MKL|Win32.ActiveCfg = Release|Win32 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release-MKL|x64.ActiveCfg = Release|x64 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release-OpenBLAS|Any CPU.ActiveCfg = Release|Win32 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release-OpenBLAS|Mixed Platforms.ActiveCfg = Release|Win32 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release-OpenBLAS|Win32.ActiveCfg = Release|Win32 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release-OpenBLAS|Win32.Build.0 = Release|Win32 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release-OpenBLAS|x64.ActiveCfg = Release|x64 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release-OpenBLAS|x64.Build.0 = Release|x64 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release-Signed|Any CPU.ActiveCfg = Release|Win32 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release-Signed|Mixed Platforms.ActiveCfg = Release|Win32 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release-Signed|Mixed Platforms.Build.0 = Release|Win32 + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1}.Release-Signed|Win32.ActiveCfg = Release|Win32 + {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-CUDA|Any CPU.ActiveCfg = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-CUDA|Mixed Platforms.ActiveCfg = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-CUDA|Win32.ActiveCfg = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-CUDA|x64.ActiveCfg = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-MKL|Any CPU.ActiveCfg = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-MKL|Mixed Platforms.ActiveCfg = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-MKL|Win32.ActiveCfg = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-MKL|x64.ActiveCfg = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-OpenBLAS|Any CPU.ActiveCfg = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-OpenBLAS|Mixed Platforms.ActiveCfg = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-OpenBLAS|Win32.ActiveCfg = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-OpenBLAS|Win32.Build.0 = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-OpenBLAS|x64.ActiveCfg = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-OpenBLAS|x64.Build.0 = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-Signed|Any CPU.ActiveCfg = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-Signed|Any CPU.Build.0 = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-Signed|Mixed Platforms.ActiveCfg = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-Signed|Mixed Platforms.Build.0 = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-Signed|Win32.ActiveCfg = Release|Any CPU + {96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-Signed|x64.ActiveCfg = Release|Any CPU EndGlobalSection GlobalSection(SolutionProperties) = preSolution HideSolutionNode = FALSE diff --git a/src/NativeProviders/OpenBLAS/blas.c b/src/NativeProviders/OpenBLAS/blas.c new file mode 100644 index 00000000..a961425c --- /dev/null +++ b/src/NativeProviders/OpenBLAS/blas.c @@ -0,0 +1,89 @@ +#include "cblas.h" +#include "wrapper_common.h" + +#if GCC +extern "C" { +#endif +DLLEXPORT void s_axpy(const blasint n, const float alpha, const float x[], float y[]){ + cblas_saxpy(n, alpha, x, 1, y, 1); +} + +DLLEXPORT void d_axpy(const blasint n, const double alpha, const double x[], double y[]){ + cblas_daxpy(n, alpha, x, 1, y, 1); +} + +DLLEXPORT void c_axpy(const blasint n, const openblas_complex_float alpha, const openblas_complex_float x[], openblas_complex_float y[]){ + cblas_caxpy(n, (float*)&alpha, (float*)x, 1, (float*)y, 1); +} + +DLLEXPORT void z_axpy(const blasint n, const openblas_complex_double alpha, const openblas_complex_double x[], openblas_complex_double y[]){ + cblas_zaxpy(n, (double*)&alpha, (double*)x, 1, (double*)y, 1); +} + +DLLEXPORT void s_scale(const blasint n, const float alpha, float x[]){ + cblas_sscal(n, alpha, x, 1); +} + +DLLEXPORT void d_scale(const blasint n, const double alpha, double x[]){ + cblas_dscal(n, alpha, x, 1); +} + +DLLEXPORT void c_scale(const blasint n, const openblas_complex_float alpha, openblas_complex_float x[]){ + cblas_cscal(n, (float*)&alpha, (float*)x, 1); +} + +DLLEXPORT void z_scale(const blasint n, const openblas_complex_double alpha, openblas_complex_double x[]){ + cblas_zscal(n, (double*)&alpha, (double*)x, 1); +} + +DLLEXPORT float s_dot_product(const blasint n, const float x[], const float y[]){ + return cblas_sdot(n, x, 1, y, 1); +} + +DLLEXPORT double d_dot_product(const blasint n, const double x[], const double y[]){ + return cblas_ddot(n, x, 1, y, 1); +} + +DLLEXPORT openblas_complex_float c_dot_product(const blasint n, const openblas_complex_float x[], const openblas_complex_float y[]){ + openblas_complex_float ret; + cblas_cdotu_sub(n, (float*)x, 1, (float*)y, 1, &ret); + return ret; +} + +DLLEXPORT openblas_complex_double z_dot_product(const blasint n, const openblas_complex_double x[], const openblas_complex_double y[]){ + openblas_complex_double ret; + cblas_zdotu_sub(n, (double*)x, 1, (double*)y, 1, &ret); + return ret; +} + +DLLEXPORT void s_matrix_multiply(CBLAS_TRANSPOSE transA, CBLAS_TRANSPOSE transB, const blasint m, const blasint n, const blasint k, const float alpha, const float x[], const float y[], const float beta, float c[]){ + blasint lda = transA == CblasNoTrans ? m : k; + blasint ldb = transB == CblasNoTrans ? k : n; + + cblas_sgemm(CblasColMajor, transA, transB, m, n, k, alpha, x, lda, y, ldb, beta, c, m); +} + +DLLEXPORT void d_matrix_multiply(CBLAS_TRANSPOSE transA, CBLAS_TRANSPOSE transB, const blasint m, const blasint n, const blasint k, const double alpha, const double x[], const double y[], const double beta, double c[]){ + blasint lda = transA == CblasNoTrans ? m : k; + blasint ldb = transB == CblasNoTrans ? k : n; + + cblas_dgemm(CblasColMajor, transA, transB, m, n, k, alpha, x, lda, y, ldb, beta, c, m); +} + +DLLEXPORT void c_matrix_multiply(CBLAS_TRANSPOSE transA, CBLAS_TRANSPOSE transB, const blasint m, const blasint n, const blasint k, const openblas_complex_float alpha, const openblas_complex_float x[], const openblas_complex_float y[], const openblas_complex_float beta, openblas_complex_float c[]){ + blasint lda = transA == CblasNoTrans ? m : k; + blasint ldb = transB == CblasNoTrans ? k : n; + + cblas_cgemm(CblasColMajor, transA, transB, m, n, k, (float*)&alpha, (float*)x, lda, (float*)y, ldb, (float*)&beta, (float*)c, m); +} + +DLLEXPORT void z_matrix_multiply(CBLAS_TRANSPOSE transA, CBLAS_TRANSPOSE transB, const blasint m, const blasint n, const blasint k, const openblas_complex_double alpha, const openblas_complex_double x[], const openblas_complex_double y[], const openblas_complex_double beta, openblas_complex_double c[]){ + blasint lda = transA == CblasNoTrans ? m : k; + blasint ldb = transB == CblasNoTrans ? k : n; + + cblas_zgemm(CblasColMajor, transA, transB, m, n, k, (double*)&alpha, (double*)x, lda, (double*)y, ldb, (double*)&beta, (double*)c, m); +} + +#if GCC +} +#endif diff --git a/src/NativeProviders/OpenBLAS/capabilities.cpp b/src/NativeProviders/OpenBLAS/capabilities.cpp new file mode 100644 index 00000000..90dc57f2 --- /dev/null +++ b/src/NativeProviders/OpenBLAS/capabilities.cpp @@ -0,0 +1,77 @@ +#include "wrapper_common.h" +#include "cblas.h" + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ + + /* + Capability is supported if >0 + + Actual number can be increased over time to indicate + extensions/revisions (that do not break compatibility) + */ + DLLEXPORT int query_capability(const int capability) + { + switch (capability) + { + + // SANITY CHECKS + case 0: return 0; + case 1: return -1; + + // PLATFORM + case 8: +#ifdef _M_IX86 + return 1; +#else + return 0; +#endif + case 9: +#ifdef _M_X64 + return 1; +#else + return 0; +#endif + case 10: +#ifdef _M_IA64 + return 1; +#else + return 0; +#endif + + // COMMON/SHARED + case 64: return 7; // revision + case 66: return 1; // threading control + + // LINEAR ALGEBRA + case 128: return 1; // basic dense linear algebra + + default: return 0; // unknown or not supported + + } + } + + DLLEXPORT void set_max_threads(const blasint num_threads) + { + openblas_set_num_threads(num_threads); + } + + DLLEXPORT char* get_build_config() + { + return openblas_get_config(); + } + + DLLEXPORT char* get_cpu_core() + { + return openblas_get_corename(); + } + + DLLEXPORT int get_parallel_type() + { + return openblas_get_parallel(); + } + +#ifdef __cplusplus +} +#endif /* __cplusplus */ diff --git a/src/NativeProviders/OpenBLAS/complex.h b/src/NativeProviders/OpenBLAS/complex.h new file mode 100644 index 00000000..d107e400 --- /dev/null +++ b/src/NativeProviders/OpenBLAS/complex.h @@ -0,0 +1,39 @@ +template +struct complex +{ + _T real, imag; + + complex(_T _real = 0, _T _imag = 0) + { + real = _real; + imag = _imag; + } + + complex(const complex<_T>& right) + { + real = right.real; + imag = right.imag; + } + + complex& operator=(const complex& right) + { + real = right.real; + imag = right.imag; + return *this; + } + + complex& operator=(const _T& right) + { + real = right; + imag = 0; + return *this; + } + + template inline + complex& operator=(const complex<_Other>& right) + { + real = (_T)right.real; + imag = (_T)right.imag; + return *this; + } +}; \ No newline at end of file diff --git a/src/NativeProviders/OpenBLAS/lapack.cpp b/src/NativeProviders/OpenBLAS/lapack.cpp new file mode 100644 index 00000000..fee9abec --- /dev/null +++ b/src/NativeProviders/OpenBLAS/lapack.cpp @@ -0,0 +1,712 @@ +#include "cblas.h" + +#include "complex.h" +#define LAPACK_COMPLEX_CUSTOM +#define lapack_complex_float complex +#define lapack_complex_double complex + +#include "lapacke.h" +#include "lapack_common.h" +#include "wrapper_common.h" +#include + +template +inline lapack_int lu_factor(lapack_int m, T a[], lapack_int ipiv[], GETRF getrf) +{ + lapack_int info = 0; + getrf(&m, &m, a, &m, ipiv, &info); + shift_ipiv_down(m, ipiv); + return info; +}; + +template +inline lapack_int lu_inverse(lapack_int n, T a[], T work[], lapack_int lwork, GETRF getrf, GETRI getri) +{ + lapack_int* ipiv = new lapack_int[n]; + lapack_int info = 0; + getrf(&n, &n, a, &n, ipiv, &info); + + if (info != 0) + { + delete[] ipiv; + return info; + } + + getri(&n, a, &n, ipiv, work, &lwork, &info); + delete[] ipiv; + return info; +}; + +template +inline lapack_int lu_inverse_factored(lapack_int n, T a[], lapack_int ipiv[], T work[], lapack_int lwork, GETRI getri) +{ + shift_ipiv_up(n, ipiv); + lapack_int info = 0; + getri(&n, a, &n, ipiv, work, &lwork, &info); + shift_ipiv_down(n, ipiv); + return info; +} + +template +inline lapack_int lu_solve_factored(lapack_int n, lapack_int nrhs, T a[], lapack_int ipiv[], T b[], GETRS getrs) +{ + shift_ipiv_up(n, ipiv); + lapack_int info = 0; + char trans ='N'; + getrs(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info); + shift_ipiv_down(n, ipiv); + return info; +} + +template +inline lapack_int lu_solve(lapack_int n, lapack_int nrhs, T a[], T b[], GETRF getrf, GETRS getrs) +{ + T* clone = Clone(n, n, a); + lapack_int* ipiv = new lapack_int[n]; + lapack_int info = 0; + getrf(&n, &n, clone, &n, ipiv, &info); + + if (info != 0) + { + delete[] ipiv; + delete[] clone; + return info; + } + + char trans ='N'; + getrs(&trans, &n, &nrhs, clone, &n, ipiv, b, &n, &info); + delete[] ipiv; + delete[] clone; + return info; +} + + +template +inline lapack_int cholesky_factor(lapack_int n, T* a, POTRF potrf) +{ + char uplo = 'L'; + lapack_int info = 0; + potrf(&uplo, &n, a, &n, &info); + T zero = T(); + + for (lapack_int i = 0; i < n; ++i) + { + lapack_int index = i * n; + + for (lapack_int j = 0; j < n && i > j; ++j) + { + a[index + j] = zero; + } + } + + return info; +} + +template +inline lapack_int cholesky_solve(lapack_int n, lapack_int nrhs, T a[], T b[], POTRF potrf, POTRS potrs) +{ + T* clone = Clone(n, n, a); + char uplo = 'L'; + lapack_int info = 0; + potrf(&uplo, &n, clone, &n, &info); + + if (info != 0) + { + delete[] clone; + return info; + } + + potrs(&uplo, &n, &nrhs, clone, &n, b, &n, &info); + delete[] clone; + return info; +} + +template +inline lapack_int cholesky_solve_factored(lapack_int n, lapack_int nrhs, T a[], T b[], POTRS potrs) +{ + char uplo = 'L'; + lapack_int info = 0; + potrs(&uplo, &n, &nrhs, a, &n, b, &n, &info); + return info; +} + +template +inline lapack_int qr_factor(lapack_int m, lapack_int n, T r[], T tau[], T q[], T work[], lapack_int len, GEQRF geqrf, ORGQR orgqr) +{ + lapack_int info = 0; + geqrf(&m, &n, r, &m, tau, work, &len, &info); + + for (lapack_int i = 0; i < m; ++i) + { + for (lapack_int j = 0; j < m && j < n; ++j) + { + if (i > j) + { + q[j * m + i] = r[j * m + i]; + } + } + } + + //compute the q elements explicitly + if (m <= n) + { + orgqr(&m, &m, &m, q, &m, tau, work, &len, &info); + } + else + { + orgqr(&m, &m, &n, q, &m, tau, work, &len, &info); + } + + return info; +} + +template +inline lapack_int qr_thin_factor(lapack_int m, lapack_int n, T q[], T tau[], T r[], T work[], lapack_int len, GEQRF geqrf, ORGQR orgqr) +{ + lapack_int info = 0; + geqrf(&m, &n, q, &m, tau, work, &len, &info); + + for (lapack_int i = 0; i < n; ++i) + { + for (lapack_int j = 0; j < n; ++j) + { + if (i <= j) + { + r[j * n + i] = q[j * m + i]; + } + } + } + + orgqr(&m, &n, &n, q, &m, tau, work, &len, &info); + return info; +} + +template +inline lapack_int qr_solve(lapack_int m, lapack_int n, lapack_int bn, T a[], T b[], T x[], T work[], lapack_int len, GELS gels) +{ + T* clone_a = Clone(m, n, a); + T* clone_b = Clone(m, bn, b); + char N = 'N'; + lapack_int info = 0; + gels(&N, &m, &n, &bn, clone_a, &m, clone_b, &m, work, &len, &info); + copyBtoX(m, n, bn, clone_b, x); + delete[] clone_a; + delete[] clone_b; + return info; +} + +template +inline lapack_int qr_solve_factored(lapack_int m, lapack_int n, lapack_int bn, T r[], T b[], T tau[], T x[], T work[], lapack_int len, ORMQR ormqr, TRSM trsm) +{ + T* clone_b = Clone(m, bn, b); + char side ='L'; + char tran = 'T'; + lapack_int info = 0; + ormqr(&side, &tran, &m, &bn, &n, r, &m, tau, clone_b, &m, work, &len, &info); + trsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, 1.0, r, m, clone_b, m); + copyBtoX(m, n, bn, clone_b, x); + delete[] clone_b; + return info; +} + +template +inline lapack_int complex_qr_solve_factored(lapack_int m, lapack_int n, lapack_int bn, T r[], T b[], T tau[], T x[], T work[], lapack_int len, UNMQR unmqr, TRSM trsm) +{ + T* clone_b = Clone(m, bn, b); + char side ='L'; + char tran = 'C'; + 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, reinterpret_cast(&one), reinterpret_cast(r), m, reinterpret_cast(clone_b), m); + copyBtoX(m, n, bn, clone_b, x); + delete[] clone_b; + return info; +} + +template +inline lapack_int svd_factor(bool compute_vectors, lapack_int m, lapack_int n, T a[], T s[], T u[], T v[], T work[], lapack_int len, GESVD gesvd) +{ + lapack_int info = 0; + char job = compute_vectors ? 'A' : 'N'; + gesvd(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, &info); + return info; +} + +template +inline lapack_int complex_svd_factor(bool compute_vectors, lapack_int m, lapack_int n, T a[], T s[], T u[], T v[], T work[], lapack_int len, GESVD gesvd) +{ + lapack_int info = 0; + lapack_int dim_s = std::min(m,n); + R* rwork = new R[5 * dim_s]; + R* s_local = new R[dim_s]; + char job = compute_vectors ? 'A' : 'N'; + gesvd(&job, &job, &m, &n, a, &m, s_local, u, &m, v, &n, work, &len, rwork, &info); + + for (lapack_int index = 0; index < dim_s; ++index) + { + s[index] = s_local[index]; + } + + delete[] rwork; + delete[] s_local; + return info; +} + +template +inline lapack_int eigen_factor(lapack_int n, T a[], T vectors[], R values[], T d[], GEES gees, TREVC trevc) +{ + T* clone_a = Clone(n, n, a); + T* wr = new T[n]; + T* wi = new T[n]; + + lapack_int sdim; + lapack_int info = gees(LAPACK_COL_MAJOR, 'V', 'N', nullptr, n, clone_a, n, &sdim, wr, wi, vectors, n); + if (info != 0) + { + delete[] clone_a; + delete[] wr; + delete[] wi; + return info; + } + + lapack_int m; + info = trevc(LAPACK_COL_MAJOR, 'R', 'B', nullptr, n, clone_a, n, nullptr, n, vectors, n, n, &m); + if (info != 0) + { + delete[] clone_a; + delete[] wr; + delete[] wi; + return info; + } + + for (lapack_int index = 0; index < n; ++index) + { + values[index] = R(wr[index], wi[index]); + } + + for (lapack_int i = 0; i < n; ++i) + { + lapack_int in = i * n; + d[in + i] = wr[i]; + + if (wi[i] > 0) + { + d[in + n + i] = wi[i]; + } + else if (wi[i] < 0) + { + d[in - n + i] = wi[i]; + } + } + + delete[] clone_a; + delete[] wr; + delete[] wi; + return info; +} + +template +inline lapack_int eigen_complex_factor(lapack_int n, T a[], T vectors[], lapack_complex_double values[], T d[], GEES gees, TREVC trevc) +{ + T* clone_a = Clone(n, n, a); + T* w = new T[n]; + + lapack_int sdim; + lapack_int info = gees(LAPACK_COL_MAJOR, 'V', 'N', nullptr, n, clone_a, n, &sdim, w, vectors, n); + if (info != 0) + { + delete[] clone_a; + delete[] w; + return info; + } + + lapack_int m; + info = trevc(LAPACK_COL_MAJOR, 'R', 'B', nullptr, n, clone_a, n, nullptr, n, vectors, n, n, &m); + if (info != 0) + { + delete[] clone_a; + delete[] w; + return info; + } + + for (lapack_int i = 0; i < n; ++i) + { + values[i] = w[i]; + d[i * n + i] = w[i]; + } + + delete[] clone_a; + delete[] w; + return info; +} + +template +inline lapack_int sym_eigen_factor(lapack_int n, T a[], T vectors[], lapack_complex_double values[], T d[], SYEV syev) +{ + T* clone_a = Clone(n, n, a); + R* w = new R[n]; + + lapack_int info = syev(LAPACK_COL_MAJOR, 'V', 'U', n, clone_a, n, w); + if (info != 0) + { + delete[] clone_a; + delete[] w; + return info; + } + + memcpy(vectors, clone_a, n*n*sizeof(T)); + + for (lapack_int index = 0; index < n; ++index) + { + values[index] = lapack_complex_double(w[index]); + } + + for (lapack_int j = 0; j < n; ++j) + { + lapack_int jn = j*n; + + for (lapack_int i = 0; i < n; ++i) + { + if (i == j) + { + d[jn + i] = w[i]; + } + } + } + + delete[] clone_a; + delete[] w; + return info; +} + +extern "C" { + + DLLEXPORT float s_matrix_norm(char norm, lapack_int m, lapack_int n, float a[], float work[]) + { + return LAPACKE_slange_work(CblasColMajor, norm, m, n, a, m, work); + } + + DLLEXPORT double d_matrix_norm(char norm, lapack_int m, lapack_int n, double a[], double work[]) + { + return LAPACKE_dlange_work(CblasColMajor, norm, m, n, a, m, work); + } + + DLLEXPORT float c_matrix_norm(char norm, lapack_int m, lapack_int n, lapack_complex_float a[], float work[]) + { + return LAPACKE_clange_work(CblasColMajor, norm, m, n, a, m, work); + } + + DLLEXPORT double z_matrix_norm(char norm, lapack_int m, lapack_int n, lapack_complex_double a[], double work[]) + { + return LAPACKE_zlange_work(CblasColMajor, norm, m, n, a, m, work); + } + + DLLEXPORT lapack_int s_lu_factor(lapack_int m, float a[], lapack_int ipiv[]) + { + return lu_factor(m, a, ipiv, LAPACK_sgetrf); + } + + DLLEXPORT lapack_int d_lu_factor(lapack_int m, double a[], lapack_int ipiv[]) + { + return lu_factor(m, a, ipiv, LAPACK_dgetrf); + } + + DLLEXPORT lapack_int c_lu_factor(lapack_int m, lapack_complex_float a[], lapack_int ipiv[]) + { + return lu_factor(m, a, ipiv, LAPACK_cgetrf); + } + + DLLEXPORT lapack_int z_lu_factor(lapack_int m, lapack_complex_double a[], lapack_int ipiv[]) + { + return lu_factor(m, a, ipiv, LAPACK_zgetrf); + } + + DLLEXPORT lapack_int s_lu_inverse(lapack_int n, float a[], float work[], lapack_int lwork) + { + return lu_inverse(n, a, work, lwork, LAPACK_sgetrf, LAPACK_sgetri); + } + + DLLEXPORT lapack_int d_lu_inverse(lapack_int n, double a[], double work[], lapack_int lwork) + { + return lu_inverse(n, a, work, lwork, LAPACK_dgetrf, LAPACK_dgetri); + } + + DLLEXPORT lapack_int c_lu_inverse(lapack_int n, lapack_complex_float a[], lapack_complex_float work[], lapack_int lwork) + { + return lu_inverse(n, a, work, lwork, LAPACK_cgetrf, LAPACK_cgetri); + } + + DLLEXPORT lapack_int z_lu_inverse(lapack_int n, lapack_complex_double a[], lapack_complex_double work[], lapack_int lwork) + { + return lu_inverse(n, a, work, lwork, LAPACK_zgetrf, LAPACK_zgetri); + } + + DLLEXPORT lapack_int s_lu_inverse_factored(lapack_int n, float a[], lapack_int ipiv[], float work[], lapack_int lwork) + { + return lu_inverse_factored(n, a, ipiv, work, lwork, LAPACK_sgetri); + } + + DLLEXPORT lapack_int d_lu_inverse_factored(lapack_int n, double a[], lapack_int ipiv[], double work[], lapack_int lwork) + { + return lu_inverse_factored(n, a, ipiv, work, lwork, LAPACK_dgetri); + } + + DLLEXPORT lapack_int c_lu_inverse_factored(lapack_int n, lapack_complex_float a[], lapack_int ipiv[], lapack_complex_float work[], lapack_int lwork) + { + return lu_inverse_factored(n, a, ipiv, work, lwork, LAPACK_cgetri); + } + + DLLEXPORT lapack_int z_lu_inverse_factored(lapack_int n, lapack_complex_double a[], lapack_int ipiv[], lapack_complex_double work[], lapack_int lwork) + { + return lu_inverse_factored(n, a, ipiv, work, lwork, LAPACK_zgetri); + } + + DLLEXPORT lapack_int s_lu_solve_factored(lapack_int n, lapack_int nrhs, float a[], lapack_int ipiv[], float b[]) + { + return lu_solve_factored(n, nrhs, a, ipiv, b, LAPACK_sgetrs); + } + + DLLEXPORT lapack_int d_lu_solve_factored(lapack_int n, lapack_int nrhs, double a[], lapack_int ipiv[], double b[]) + { + return lu_solve_factored(n, nrhs, a, ipiv, b, LAPACK_dgetrs); + } + + DLLEXPORT lapack_int c_lu_solve_factored(lapack_int n, lapack_int nrhs, lapack_complex_float a[], lapack_int ipiv[], lapack_complex_float b[]) + { + return lu_solve_factored(n, nrhs, a, ipiv, b, LAPACK_cgetrs); + } + + DLLEXPORT lapack_int z_lu_solve_factored(lapack_int n, lapack_int nrhs, lapack_complex_double a[], lapack_int ipiv[], lapack_complex_double b[]) + { + return lu_solve_factored(n, nrhs, a, ipiv, b, LAPACK_zgetrs); + } + + DLLEXPORT lapack_int s_lu_solve(lapack_int n, lapack_int nrhs, float a[], float b[]) + { + return lu_solve(n, nrhs, a, b, LAPACK_sgetrf, LAPACK_sgetrs); + } + + DLLEXPORT lapack_int d_lu_solve(lapack_int n, lapack_int nrhs, double a[], double b[]) + { + return lu_solve(n, nrhs, a, b, LAPACK_dgetrf, LAPACK_dgetrs); + } + + DLLEXPORT lapack_int c_lu_solve(lapack_int n, lapack_int nrhs, lapack_complex_float a[], lapack_complex_float b[]) + { + return lu_solve(n, nrhs, a, b, LAPACK_cgetrf, LAPACK_cgetrs); + } + + DLLEXPORT lapack_int z_lu_solve(lapack_int n, lapack_int nrhs, lapack_complex_double a[], lapack_complex_double b[]) + { + return lu_solve(n, nrhs, a, b, LAPACK_zgetrf, LAPACK_zgetrs); + } + + DLLEXPORT lapack_int s_cholesky_factor(lapack_int n, float a[]) + { + return cholesky_factor(n, a, LAPACK_spotrf); + } + + DLLEXPORT lapack_int d_cholesky_factor(lapack_int n, double* a) + { + return cholesky_factor(n, a, LAPACK_dpotrf); + } + + DLLEXPORT lapack_int c_cholesky_factor(lapack_int n, lapack_complex_float a[]) + { + return cholesky_factor(n, a, LAPACK_cpotrf); + } + + DLLEXPORT lapack_int z_cholesky_factor(lapack_int n, lapack_complex_double a[]) + { + return cholesky_factor(n, a, LAPACK_zpotrf); + } + + DLLEXPORT lapack_int s_cholesky_solve(lapack_int n, lapack_int nrhs, float a[], float b[]) + { + return cholesky_solve(n, nrhs, a, b, LAPACK_spotrf, LAPACK_spotrs); + } + + DLLEXPORT lapack_int d_cholesky_solve(lapack_int n, lapack_int nrhs, double a[], double b[]) + { + return cholesky_solve(n, nrhs, a, b, LAPACK_dpotrf, LAPACK_dpotrs); + } + + DLLEXPORT lapack_int c_cholesky_solve(lapack_int n, lapack_int nrhs, lapack_complex_float a[], lapack_complex_float b[]) + { + return cholesky_solve(n, nrhs, a, b, LAPACK_cpotrf, LAPACK_cpotrs); + } + + DLLEXPORT lapack_int z_cholesky_solve(lapack_int n, lapack_int nrhs, lapack_complex_double a[], lapack_complex_double b[]) + { + return cholesky_solve(n, nrhs, a, b, LAPACK_zpotrf, LAPACK_zpotrs); + } + + DLLEXPORT lapack_int s_cholesky_solve_factored(lapack_int n, lapack_int nrhs, float a[], float b[]) + { + return cholesky_solve_factored(n, nrhs, a, b, LAPACK_spotrs); + } + + DLLEXPORT lapack_int d_cholesky_solve_factored(lapack_int n, lapack_int nrhs, double a[], double b[]) + { + return cholesky_solve_factored(n, nrhs, a, b, LAPACK_dpotrs); + } + + DLLEXPORT lapack_int c_cholesky_solve_factored(lapack_int n, lapack_int nrhs, lapack_complex_float a[], lapack_complex_float b[]) + { + return cholesky_solve_factored(n, nrhs, a, b, LAPACK_cpotrs); + } + + DLLEXPORT lapack_int z_cholesky_solve_factored(lapack_int n, lapack_int nrhs, lapack_complex_double a[], lapack_complex_double b[]) + { + return cholesky_solve_factored(n, nrhs, a, b, LAPACK_zpotrs); + } + + DLLEXPORT lapack_int s_qr_factor(lapack_int m, lapack_int n, float r[], float tau[], float q[], float work[], lapack_int len) + { + return qr_factor(m, n, r, tau, q, work, len, LAPACK_sgeqrf, LAPACK_sorgqr); + } + + DLLEXPORT lapack_int s_qr_thin_factor(lapack_int m, lapack_int n, float q[], float tau[], float r[], float work[], lapack_int len) + { + return qr_thin_factor(m, n, q, tau, r, work, len, LAPACK_sgeqrf, LAPACK_sorgqr); + } + + DLLEXPORT lapack_int d_qr_factor(lapack_int m, lapack_int n, double r[], double tau[], double q[], double work[], lapack_int len) + { + return qr_factor(m, n, r, tau, q, work, len, LAPACK_dgeqrf, LAPACK_dorgqr); + } + + DLLEXPORT lapack_int d_qr_thin_factor(lapack_int m, lapack_int n, double q[], double tau[], double r[], double work[], lapack_int len) + { + return qr_thin_factor(m, n, q, tau, r, work, len, LAPACK_dgeqrf, LAPACK_dorgqr); + } + + DLLEXPORT lapack_int c_qr_factor(lapack_int m, lapack_int n, lapack_complex_float r[], lapack_complex_float tau[], lapack_complex_float q[], lapack_complex_float work[], lapack_int len) + { + return qr_factor(m, n, r, tau, q, work, len, LAPACK_cgeqrf, LAPACK_cungqr); + } + + DLLEXPORT lapack_int c_qr_thin_factor(lapack_int m, lapack_int n, lapack_complex_float q[], lapack_complex_float tau[], lapack_complex_float r[], lapack_complex_float work[], lapack_int len) + { + return qr_thin_factor(m, n, q, tau, r, work, len, LAPACK_cgeqrf, LAPACK_cungqr); + } + + DLLEXPORT lapack_int z_qr_factor(lapack_int m, lapack_int n, lapack_complex_double r[], lapack_complex_double tau[], lapack_complex_double q[], lapack_complex_double work[], lapack_int len) + { + return qr_factor(m, n, r, tau, q, work, len, LAPACK_zgeqrf, LAPACK_zungqr); + } + + DLLEXPORT lapack_int z_qr_thin_factor(lapack_int m, lapack_int n, lapack_complex_double q[], lapack_complex_double tau[], lapack_complex_double r[], lapack_complex_double work[], lapack_int len) + { + return qr_thin_factor(m, n, q, tau, r, work, len, LAPACK_zgeqrf, LAPACK_zungqr); + } + + DLLEXPORT lapack_int s_qr_solve(lapack_int m, lapack_int n, lapack_int bn, float a[], float b[], float x[], float work[], lapack_int len) + { + return qr_solve(m, n, bn, a, b, x, work, len, LAPACK_sgels); + } + + DLLEXPORT lapack_int d_qr_solve(lapack_int m, lapack_int n, lapack_int bn, double a[], double b[], double x[], double work[], lapack_int len) + { + return qr_solve(m, n, bn, a, b, x, work, len, LAPACK_dgels); + } + + DLLEXPORT lapack_int c_qr_solve(lapack_int m, lapack_int n, lapack_int bn, lapack_complex_float a[], lapack_complex_float b[], lapack_complex_float x[], lapack_complex_float work[], lapack_int len) + { + return qr_solve(m, n, bn, a, b, x, work, len, LAPACK_cgels); + } + + DLLEXPORT lapack_int z_qr_solve(lapack_int m, lapack_int n, lapack_int bn, lapack_complex_double a[], lapack_complex_double b[], lapack_complex_double x[], lapack_complex_double work[], lapack_int len) + { + return qr_solve(m, n, bn, a, b, x, work, len, LAPACK_zgels); + } + + DLLEXPORT lapack_int s_qr_solve_factored(lapack_int m, lapack_int n, lapack_int bn, float r[], float b[], float tau[], float x[], float work[], lapack_int len) + { + return qr_solve_factored(m, n, bn, r, b, tau, x, work, len, LAPACK_sormqr, cblas_strsm); + } + + DLLEXPORT lapack_int d_qr_solve_factored(lapack_int m, lapack_int n, lapack_int bn, double r[], double b[], double tau[], double x[], double work[], lapack_int len) + { + return qr_solve_factored(m, n, bn, r, b, tau, x, work, len, LAPACK_dormqr, cblas_dtrsm); + } + + DLLEXPORT lapack_int c_qr_solve_factored(lapack_int m, lapack_int n, lapack_int bn, lapack_complex_float r[], lapack_complex_float b[], lapack_complex_float tau[], lapack_complex_float x[], lapack_complex_float work[], lapack_int len) + { + return complex_qr_solve_factored(m, n, bn, r, b, tau, x, work, len, LAPACK_cunmqr, cblas_ctrsm); + } + + DLLEXPORT lapack_int z_qr_solve_factored(lapack_int m, lapack_int n, lapack_int bn, lapack_complex_double r[], lapack_complex_double b[], lapack_complex_double tau[], lapack_complex_double x[], lapack_complex_double work[], lapack_int len) + { + return complex_qr_solve_factored(m, n, bn, r, b, tau, x, work, len, LAPACK_zunmqr, cblas_ztrsm); + } + + DLLEXPORT lapack_int s_svd_factor(bool compute_vectors, lapack_int m, lapack_int n, float a[], float s[], float u[], float v[], float work[], lapack_int len) + { + return svd_factor(compute_vectors, m, n, a, s, u, v, work, len, LAPACK_sgesvd); + } + + DLLEXPORT lapack_int d_svd_factor(bool compute_vectors, lapack_int m, lapack_int n, double a[], double s[], double u[], double v[], double work[], lapack_int len) + { + return svd_factor(compute_vectors, m, n, a, s, u, v, work, len, LAPACK_dgesvd); + } + + DLLEXPORT lapack_int c_svd_factor(bool compute_vectors, lapack_int m, lapack_int n, lapack_complex_float a[], lapack_complex_float s[], lapack_complex_float u[], lapack_complex_float v[], lapack_complex_float work[], lapack_int len) + { + return complex_svd_factor(compute_vectors, m, n, a, s, u, v, work, len, LAPACK_cgesvd); + } + + DLLEXPORT lapack_int z_svd_factor(bool compute_vectors, lapack_int m, lapack_int n, lapack_complex_double a[], lapack_complex_double s[], lapack_complex_double u[], lapack_complex_double v[], lapack_complex_double work[], lapack_int len) + { + return complex_svd_factor(compute_vectors, m, n, a, s, u, v, work, len, LAPACK_zgesvd); + } + + DLLEXPORT lapack_int s_eigen(bool isSymmetric, lapack_int n, float a[], float vectors[], lapack_complex_double values[], float d[]) + { + if (isSymmetric) + { + return sym_eigen_factor(n, a, vectors, values, d, LAPACKE_ssyev); + } + else + { + return eigen_factor(n, a, vectors, values, d, LAPACKE_sgees, LAPACKE_strevc); + } + } + + DLLEXPORT lapack_int d_eigen(bool isSymmetric, lapack_int n, double a[], double vectors[], lapack_complex_double values[], double d[]) + { + if (isSymmetric) + { + return sym_eigen_factor(n, a, vectors, values, d, LAPACKE_dsyev); + } + else + { + return eigen_factor(n, a, vectors, values, d, LAPACKE_dgees, LAPACKE_dtrevc); + } + } + + DLLEXPORT lapack_int c_eigen(bool isSymmetric, lapack_int n, lapack_complex_float a[], lapack_complex_float vectors[], lapack_complex_double values[], lapack_complex_float d[]) + { + if (isSymmetric) + { + return sym_eigen_factor(n, a, vectors, values, d, LAPACKE_cheev); + } + else + { + return eigen_complex_factor(n, a, vectors, values, d, LAPACKE_cgees, LAPACKE_ctrevc); + } + } + + DLLEXPORT lapack_int z_eigen(bool isSymmetric, lapack_int n, lapack_complex_double a[], lapack_complex_double vectors[], lapack_complex_double values[], lapack_complex_double d[]) + { + if (isSymmetric) + { + return sym_eigen_factor(n, a, vectors, values, d, LAPACKE_zheev); + } + else + { + return eigen_complex_factor(n, a, vectors, values, d, LAPACKE_zgees, LAPACKE_ztrevc); + } + } +} diff --git a/src/NativeProviders/Windows/OpenBLAS/OpenBLASWrapper.vcxproj b/src/NativeProviders/Windows/OpenBLAS/OpenBLASWrapper.vcxproj new file mode 100644 index 00000000..46523e0f --- /dev/null +++ b/src/NativeProviders/Windows/OpenBLAS/OpenBLASWrapper.vcxproj @@ -0,0 +1,220 @@ + + + + + Debug + Win32 + + + Debug + x64 + + + Release + Win32 + + + Release + x64 + + + + {CB4011B6-E9A7-480B-A7B1-8492039DAAD1} + OpenBLASWrapper + OpenBLAS + + + + DynamicLibrary + v120 + MultiByte + true + + + DynamicLibrary + v120 + MultiByte + + + DynamicLibrary + v120 + MultiByte + true + + + DynamicLibrary + v120 + MultiByte + + + + + + + + + + + + + + + + + + $(ProjectDir)..\..\..\..\..\libs\OpenBLAS + $(OpenBlasPath)\include\ + $(OpenBlasPath)\x86\ + $(OpenBlasPath)\x64\ + + + <_ProjectFileVersion>11.0.50727.1 + + + $(ProjectDir)..\..\..\..\out\OpenBLAS\Windows\x86\ + $(Platform)\$(Configuration)\ + MathNet.Numerics.OpenBLAS + $(OpenBLASIncludeDir);$(IncludePath) + + + $(ProjectDir)..\..\..\..\out\OpenBLAS\Windows\x64\ + $(Platform)\$(Configuration)\ + MathNet.Numerics.OpenBLAS + $(OpenBLASIncludeDir);$(IncludePath) + + + $(ProjectDir)..\..\..\..\out\OpenBLAS\Windows\x86\ + $(Platform)\$(Configuration)\ + MathNet.Numerics.OpenBLAS + $(OpenBLASIncludeDir);$(IncludePath) + + + $(ProjectDir)..\..\..\..\out\OpenBLAS\Windows\x64\ + $(Platform)\$(Configuration)\ + MathNet.Numerics.OpenBLAS + $(OpenBLASIncludeDir);$(IncludePath) + + + + Disabled + $(ProjectDir)..\..\Common;$(ProjectDir)..\..\OpenBLAS;$(OpenBLASIncludeDir);%(AdditionalIncludeDirectories) + _WINDOWS;%(PreprocessorDefinitions) + true + EnableFastChecks + MultiThreadedDebug + Level3 + EditAndContinue + Default + + + libopenblas.dll.a;%(AdditionalDependencies) + $(OutDir)$(TargetName)$(TargetExt) + true + MachineX86 + $(OutDir)$(TargetName).lib + $(OpenBLASLibDir);%(AdditionalLibraryDirectories) + + + copy "$(OpenBLASLibDir)*.dll" $(OutputPath) + + + + + X64 + + + Disabled + $(ProjectDir)..\..\Common;$(ProjectDir)..\..\OpenBLAS;$(OpenBLASIncludeDir);%(AdditionalIncludeDirectories) + _WINDOWS;%(PreprocessorDefinitions) + true + EnableFastChecks + MultiThreadedDebug + Level3 + ProgramDatabase + Default + + + libopenblas.dll.a;%(AdditionalDependencies) + $(OutDir)$(TargetName)$(TargetExt) + true + MachineX64 + $(OutDir)$(TargetName).lib + $(OpenBLASLibDir_x64);%(AdditionalLibraryDirectories) + + + copy "$(OpenBLASLibDir_x64)*.dll" $(OutputPath) + + + + + MaxSpeed + true + $(ProjectDir)..\..\Common;$(ProjectDir)..\..\OpenBLAS;$(OpenBLASIncludeDir);%(AdditionalIncludeDirectories) + _WINDOWS;%(PreprocessorDefinitions) + MultiThreaded + true + Level3 + ProgramDatabase + Default + /Qvec-report:1 %(AdditionalOptions) + + + libopenblas.dll.a;%(AdditionalDependencies) + $(OutDir)$(TargetName)$(TargetExt) + true + true + true + MachineX86 + $(OutDir)$(TargetName).lib + $(OpenBLASLibDir);%(AdditionalLibraryDirectories) + + + copy "$(OpenBLASLibDir)*.dll" $(OutputPath) + + + + + X64 + + + MaxSpeed + true + $(ProjectDir)..\..\Common;$(ProjectDir)..\..\OpenBLAS;$(OpenBLASIncludeDir);%(AdditionalIncludeDirectories) + _WINDOWS;%(PreprocessorDefinitions) + MultiThreaded + true + Level3 + ProgramDatabase + Default + /Qvec-report:1 %(AdditionalOptions) + + + libopenblas.dll.a;%(AdditionalDependencies) + $(OutDir)$(TargetName)$(TargetExt) + true + true + true + MachineX64 + $(OutDir)$(TargetName).lib + $(OpenBLASLibDir_x64);%(AdditionalLibraryDirectories) + + + copy "$(OpenBLASLibDir_x64)*.dll" $(OutputPath) + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/src/NativeProviders/Windows/OpenBLAS/OpenBLASWrapper.vcxproj.filters b/src/NativeProviders/Windows/OpenBLAS/OpenBLASWrapper.vcxproj.filters new file mode 100644 index 00000000..f3bab0ba --- /dev/null +++ b/src/NativeProviders/Windows/OpenBLAS/OpenBLASWrapper.vcxproj.filters @@ -0,0 +1,41 @@ + + + + + {4FC737F1-C7A5-4376-A066-2A32D752A2FF} + cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx + + + {93995380-89BD-4b04-88EB-625FBE52EBFB} + h;hpp;hxx;hm;inl;inc;xsd + + + {67DA6AB6-F800-4c08-8B7A-83BB121AAD01} + rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx;tiff;tif;png;wav + + + + + Source Files + + + Source Files + + + Source Files + + + Source Files + + + + + Resource Files + + + + + Header Files + + + \ No newline at end of file diff --git a/src/Numerics/Control.cs b/src/Numerics/Control.cs index 0e9570fd..804721ca 100644 --- a/src/Numerics/Control.cs +++ b/src/Numerics/Control.cs @@ -136,6 +136,11 @@ namespace MathNet.Numerics { LinearAlgebraProvider = new Providers.LinearAlgebra.Cuda.CudaLinearAlgebraProvider(); } + + public static void UseNativeOpenBLAS() + { + LinearAlgebraProvider = new Providers.LinearAlgebra.OpenBlas.OpenBlasLinearAlgebraProvider(); + } #endif /// diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 7d4fccd6..8255c31f 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -163,12 +163,12 @@ - - - - - - + + + + + + @@ -198,6 +198,7 @@ + diff --git a/src/Numerics/Properties/AssemblyInfo.cs b/src/Numerics/Properties/AssemblyInfo.cs index 930f5179..61e86993 100644 --- a/src/Numerics/Properties/AssemblyInfo.cs +++ b/src/Numerics/Properties/AssemblyInfo.cs @@ -77,6 +77,7 @@ using System.Runtime.InteropServices; [assembly: InternalsVisibleTo("MathNet.Numerics.UnitTests")] [assembly: InternalsVisibleTo("MathNet.Numerics.UnitTestsMKL")] [assembly: InternalsVisibleTo("MathNet.Numerics.UnitTestsCUDA")] +[assembly: InternalsVisibleTo("MathNet.Numerics.UnitTestsOpenBLAS")] [assembly: InternalsVisibleTo("Performance")] #endif diff --git a/src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.cs deleted file mode 100644 index 5644076b..00000000 --- a/src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.cs +++ /dev/null @@ -1,367 +0,0 @@ -// -// Math.NET Numerics, part of the Math.NET Project -// http://numerics.mathdotnet.com -// http://github.com/mathnet/mathnet-numerics -// http://mathnetnumerics.codeplex.com -// -// Copyright (c) 2009-2011 Math.NET -// -// Permission is hereby granted, free of charge, to any person -// obtaining a copy of this software and associated documentation -// files (the "Software"), to deal in the Software without -// restriction, including without limitation the rights to use, -// copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the -// Software is furnished to do so, subject to the following -// conditions: -// -// The above copyright notice and this permission notice shall be -// included in all copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, -// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES -// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND -// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT -// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, -// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING -// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR -// OTHER DEALINGS IN THE SOFTWARE. -// - -#if NATIVEGOTO - -using MathNet.Numerics.Properties; -using System; -using System.Numerics; -using System.Security; - -namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas -{ - /// - /// GotoBLAS2 linear algebra provider. - /// - public partial class GotoBlasLinearAlgebraProvider : ManagedLinearAlgebraProvider - { - /// - /// Computes the requested of the matrix. - /// - /// The type of norm to compute. - /// The number of rows in the matrix. - /// The number of columns in the matrix. - /// The matrix to compute the norm from. - /// - /// The requested of the matrix. - /// - [SecuritySafeCritical] - public override float MatrixNorm(Norm norm, int rows, int columns, float[] matrix) - { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - - if (rows <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "rows"); - } - - if (columns <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "columns"); - } - - if (matrix.Length < rows*columns) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows*columns), "matrix"); - } - - var work = new float[rows]; - return MatrixNorm(norm, rows, columns, matrix, work); - } - - /// - /// Computes the requested of the matrix. - /// - /// The type of norm to compute. - /// The number of rows in the matrix. - /// The number of columns in the matrix. - /// The matrix to compute the norm from. - /// The work array. Only used when - /// and needs to be have a length of at least M (number of rows of . - /// - /// The requested of the matrix. - /// - [SecuritySafeCritical] - public override float MatrixNorm(Norm norm, int rows, int columns, float[] matrix, float[] work) - { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - - if (rows <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "rows"); - } - - if (columns <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "columns"); - } - - if (matrix.Length < rows*columns) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows*columns), "matrix"); - } - - if (work.Length < rows) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows), "work"); - } - - return SafeNativeMethods.s_matrix_norm((byte) norm, rows, columns, matrix, work); - } - - /// - /// Computes the requested of the matrix. - /// - /// The type of norm to compute. - /// The number of rows in the matrix. - /// The number of columns in the matrix. - /// The matrix to compute the norm from. - /// - /// The requested of the matrix. - /// - [SecuritySafeCritical] - public override double MatrixNorm(Norm norm, int rows, int columns, double[] matrix) - { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - - if (rows <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "rows"); - } - - if (columns <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "columns"); - } - - if (matrix.Length < rows*columns) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows*columns), "matrix"); - } - - var work = new double[rows]; - return MatrixNorm(norm, rows, columns, matrix, work); - } - - /// - /// Computes the requested of the matrix. - /// - /// The type of norm to compute. - /// The number of rows in the matrix. - /// The number of columns in the matrix. - /// The matrix to compute the norm from. - /// The work array. Only used when - /// and needs to be have a length of at least M (number of rows of . - /// - /// The requested of the matrix. - /// - [SecuritySafeCritical] - public override double MatrixNorm(Norm norm, int rows, int columns, double[] matrix, double[] work) - { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - - if (rows <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "rows"); - } - - if (columns <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "columns"); - } - - if (matrix.Length < rows*columns) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows*columns), "matrix"); - } - - if (work.Length < rows) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows), "work"); - } - - return SafeNativeMethods.d_matrix_norm((byte) norm, rows, columns, matrix, work); - } - - /// - /// Computes the requested of the matrix. - /// - /// The type of norm to compute. - /// The number of rows in the matrix. - /// The number of columns in the matrix. - /// The matrix to compute the norm from. - /// - /// The requested of the matrix. - /// - [SecuritySafeCritical] - public override Complex32 MatrixNorm(Norm norm, int rows, int columns, Complex32[] matrix) - { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - - if (rows <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "rows"); - } - - if (columns <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "columns"); - } - - if (matrix.Length < rows*columns) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows*columns), "matrix"); - } - - var work = new float[rows]; - return MatrixNorm(norm, rows, columns, matrix, work); - } - - /// - /// Computes the requested of the matrix. - /// - /// The type of norm to compute. - /// The number of rows in the matrix. - /// The number of columns in the matrix. - /// The matrix to compute the norm from. - /// The work array. Only used when - /// and needs to be have a length of at least M (number of rows of . - /// - /// The requested of the matrix. - /// - [SecuritySafeCritical] - public override Complex32 MatrixNorm(Norm norm, int rows, int columns, Complex32[] matrix, float[] work) - { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - - if (rows <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "rows"); - } - - if (columns <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "columns"); - } - - if (matrix.Length < rows*columns) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows*columns), "matrix"); - } - - if (work.Length < rows) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows), "work"); - } - - return SafeNativeMethods.c_matrix_norm((byte) norm, rows, columns, matrix, work); - } - - /// - /// Computes the requested of the matrix. - /// - /// The type of norm to compute. - /// The number of rows in the matrix. - /// The number of columns in the matrix. - /// The matrix to compute the norm from. - /// - /// The requested of the matrix. - /// - [SecuritySafeCritical] - public override Complex MatrixNorm(Norm norm, int rows, int columns, Complex[] matrix) - { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - - if (rows <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "rows"); - } - - if (columns <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "columns"); - } - - if (matrix.Length < rows*columns) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows*columns), "matrix"); - } - - var work = new double[rows]; - return MatrixNorm(norm, rows, columns, matrix, work); - } - - /// - /// Computes the requested of the matrix. - /// - /// The type of norm to compute. - /// The number of rows in the matrix. - /// The number of columns in the matrix. - /// The matrix to compute the norm from. - /// The work array. Only used when - /// and needs to be have a length of at least M (number of rows of . - /// - /// The requested of the matrix. - /// - [SecuritySafeCritical] - public override Complex MatrixNorm(Norm norm, int rows, int columns, Complex[] matrix, double[] work) - { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - - if (rows <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "rows"); - } - - if (columns <= 0) - { - throw new ArgumentException(Resources.ArgumentMustBePositive, "columns"); - } - - if (matrix.Length < rows*columns) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows*columns), "matrix"); - } - - if (work.Length < rows) - { - throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows), "work"); - } - - return SafeNativeMethods.z_matrix_norm((byte) norm, rows, columns, matrix, work); - } - } -} - -#endif diff --git a/src/Numerics/Providers/LinearAlgebra/GotoBlas/SafeNativeMethods.cs b/src/Numerics/Providers/LinearAlgebra/GotoBlas/SafeNativeMethods.cs deleted file mode 100644 index 6ed47e53..00000000 --- a/src/Numerics/Providers/LinearAlgebra/GotoBlas/SafeNativeMethods.cs +++ /dev/null @@ -1,263 +0,0 @@ -// -// Math.NET Numerics, part of the Math.NET Project -// http://mathnet.opensourcedotnet.info -// -// Copyright (c) 2009-2010 Math.NET -// -// Permission is hereby granted, free of charge, to any person -// obtaining a copy of this software and associated documentation -// files (the "Software"), to deal in the Software without -// restriction, including without limitation the rights to use, -// copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the -// Software is furnished to do so, subject to the following -// conditions: -// -// The above copyright notice and this permission notice shall be -// included in all copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, -// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES -// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND -// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT -// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, -// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING -// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR -// OTHER DEALINGS IN THE SOFTWARE. -// - -#if NATIVEGOTO - -using System.Numerics; -using System.Runtime.InteropServices; -using System.Security; - -namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas -{ - /// - /// P/Invoke methods to the native math libraries. - /// - [SuppressUnmanagedCodeSecurity] - [SecurityCritical] - internal static class SafeNativeMethods - { - /// - /// Name of the native DLL. - /// - const string DllName = "MathNET.Numerics.GotoBLAS2.dll"; - - #region BLAS - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void s_axpy(int n, float alpha, float[] x, [In, Out] float[] y); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void d_axpy(int n, double alpha, double[] x, [In, Out] double[] y); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void c_axpy(int n, Complex32 alpha, Complex32[] x, [In, Out] Complex32[] y); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void z_axpy(int n, Complex alpha, Complex[] x, [In, Out] Complex[] y); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void s_scale(int n, float alpha, [Out] float[] x); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void d_scale(int n, double alpha, [Out] double[] x); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void c_scale(int n, Complex32 alpha, [In, Out] Complex32[] x); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void z_scale(int n, Complex alpha, [In, Out] Complex[] x); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern float s_dot_product(int n, float[] x, float[] y); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern double d_dot_product(int n, double[] x, double[] y); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern Complex32 c_dot_product(int n, Complex32[] x, Complex32[] y); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern Complex z_dot_product(int n, Complex[] x, Complex[] y); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void s_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, float alpha, float[] x, float[] y, float beta, [In, Out] float[] c); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void d_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, double alpha, double[] x, double[] y, double beta, [In, Out] double[] c); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void c_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, Complex32 alpha, Complex32[] x, Complex32[] y, Complex32 beta, [In, Out] Complex32[] c); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void z_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, Complex alpha, Complex[] x, Complex[] y, Complex beta, [In, Out] Complex[] c); - - #endregion BLAS - - #region LAPACK - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern float s_matrix_norm(byte norm, int rows, int columns, [In] float[] a, [In, Out] float[] work); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern float d_matrix_norm(byte norm, int rows, int columns, [In] double[] a, [In, Out] double[] work); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern float c_matrix_norm(byte norm, int rows, int columns, [In] Complex32[] a, [In, Out] float[] work); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern double z_matrix_norm(byte norm, int rows, int columns, [In] Complex[] a, [In, Out] double[] work); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_cholesky_factor(int n, [In, Out] float[] a); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_cholesky_factor(int n, [In, Out] double[] a); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_cholesky_factor(int n, [In, Out] Complex32[] a); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_cholesky_factor(int n, [In, Out] Complex[] a); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_lu_factor(int n, [In, Out] float[] a, [In, Out] int[] ipiv); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_lu_factor(int n, [In, Out] double[] a, [In, Out] int[] ipiv); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_lu_factor(int n, [In, Out] Complex32[] a, [In, Out] int[] ipiv); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_lu_factor(int n, [In, Out] Complex[] a, [In, Out] int[] ipiv); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_lu_inverse(int n, [In, Out] float[] a, [In, Out] float[] work, int lwork); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_lu_inverse(int n, [In, Out] double[] a, [In, Out] double[] work, int lwork); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_lu_inverse(int n, [In, Out] Complex32[] a, [In, Out] Complex32[] work, int lwork); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_lu_inverse(int n, [In, Out] Complex[] a, [In, Out] Complex[] work, int lwork); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_lu_inverse_factored(int n, [In, Out] float[] a, [In, Out] int[] ipiv, [In, Out] float[] work, int lwork); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_lu_inverse_factored(int n, [In, Out] double[] a, [In, Out] int[] ipiv, [In, Out] double[] work, int lwork); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_lu_inverse_factored(int n, [In, Out] Complex32[] a, [In, Out] int[] ipiv, [In, Out] Complex32[] work, int lwork); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_lu_inverse_factored(int n, [In, Out] Complex[] a, [In, Out] int[] ipiv, [In, Out] Complex[] work, int lwork); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_lu_solve_factored(int n, int nrhs, float[] a, [In, Out] int[] ipiv, [In, Out] float[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_lu_solve_factored(int n, int nrhs, double[] a, [In, Out] int[] ipiv, [In, Out] double[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_lu_solve_factored(int n, int nrhs, Complex32[] a, [In, Out] int[] ipiv, [In, Out] Complex32[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_lu_solve_factored(int n, int nrhs, Complex[] a, [In, Out] int[] ipiv, [In, Out] Complex[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_lu_solve(int n, int nrhs, float[] a, [In, Out] float[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_lu_solve(int n, int nrhs, double[] a, [In, Out] double[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_lu_solve(int n, int nrhs, Complex32[] a, [In, Out] Complex32[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_lu_solve(int n, int nrhs, Complex[] a, [In, Out] Complex[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_cholesky_solve(int n, int nrhs, float[] a, [In, Out] float[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_cholesky_solve(int n, int nrhs, double[] a, [In, Out] double[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_cholesky_solve(int n, int nrhs, Complex32[] a, [In, Out] Complex32[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_cholesky_solve(int n, int nrhs, Complex[] a, [In, Out] Complex[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_cholesky_solve_factored(int n, int nrhs, float[] a, [In, Out] float[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_cholesky_solve_factored(int n, int nrhs, double[] a, [In, Out] double[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_cholesky_solve_factored(int n, int nrhs, Complex32[] a, [In, Out] Complex32[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_cholesky_solve_factored(int n, int nrhs, Complex[] a, [In, Out] Complex[] b); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_qr_factor(int m, int n, [In, Out] float[] r, [In, Out] float[] tau, [In, Out] float[] q, [In, Out] float[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_qr_factor(int m, int n, [In, Out] double[] r, [In, Out] double[] tau, [In, Out] double[] q, [In, Out] double[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_qr_factor(int m, int n, [In, Out] Complex32[] r, [In, Out] Complex32[] tau, [In, Out] Complex32[] q, [In, Out] Complex32[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_qr_factor(int m, int n, [In, Out] Complex[] r, [In, Out] Complex[] tau, [In, Out] Complex[] q, [In, Out] Complex[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_qr_solve(int m, int n, int bn, float[] r, float[] b, [In, Out] float[] x, [In, Out] float[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_qr_solve(int m, int n, int bn, double[] r, double[] b, [In, Out] double[] x, [In, Out] double[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_qr_solve(int m, int n, int bn, Complex32[] r, Complex32[] b, [In, Out] Complex32[] x, [In, Out] Complex32[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_qr_solve(int m, int n, int bn, Complex[] r, Complex[] b, [In, Out] Complex[] x, [In, Out] Complex[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_qr_solve_factored(int m, int n, int bn, float[] r, float[] b, float[] tau, [In, Out] float[] x, [In, Out] float[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_qr_solve_factored(int m, int n, int bn, double[] r, double[] b, double[] tau, [In, Out] double[] x, [In, Out] double[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_qr_solve_factored(int m, int n, int bn, Complex32[] r, Complex32[] b, Complex32[] tau, [In, Out] Complex32[] x, [In, Out] Complex32[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_qr_solve_factored(int m, int n, int bn, Complex[] r, Complex[] b, Complex[] tau, [In, Out] Complex[] x, [In, Out] Complex[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int s_svd_factor(bool computeVectors, int m, int n, [In, Out] float[] a, [In, Out] float[] s, [In, Out] float[] u, [In, Out] float[] v, [In, Out] float[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int d_svd_factor(bool computeVectors, int m, int n, [In, Out] double[] a, [In, Out] double[] s, [In, Out] double[] u, [In, Out] double[] v, [In, Out] double[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int c_svd_factor(bool computeVectors, int m, int n, [In, Out] Complex32[] a, [In, Out] Complex32[] s, [In, Out] Complex32[] u, [In, Out] Complex32[] v, [In, Out] Complex32[] work, int len); - - [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern int z_svd_factor(bool computeVectors, int m, int n, [In, Out] Complex[] a, [In, Out] Complex[] s, [In, Out] Complex[] u, [In, Out] Complex[] v, [In, Out] Complex[] work, int len); - - #endregion LAPACK - } -} - -#endif diff --git a/src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Complex.cs b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex.cs similarity index 73% rename from src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Complex.cs rename to src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex.cs index efa55565..a683b7bf 100644 --- a/src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Complex.cs +++ b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex.cs @@ -1,10 +1,10 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2011 Math.NET +// Copyright (c) 2009-2015 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -28,7 +28,7 @@ // OTHER DEALINGS IN THE SOFTWARE. // -#if NATIVEGOTO +#if NATIVE using MathNet.Numerics.LinearAlgebra.Factorization; using MathNet.Numerics.Properties; @@ -36,13 +36,78 @@ using System; using System.Numerics; using System.Security; -namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas +namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas { /// - /// GotoBLAS2 linear algebra provider. + /// OpenBLAS linear algebra provider. /// - public partial class GotoBlasLinearAlgebraProvider + public partial class OpenBlasLinearAlgebraProvider { + /// + /// Computes the requested of the matrix. + /// + /// The type of norm to compute. + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The matrix to compute the norm from. + /// + /// The requested of the matrix. + /// + [SecuritySafeCritical] + public override double MatrixNorm(Norm norm, int rows, int columns, Complex[] matrix) + { + if (matrix == null) + { + throw new ArgumentNullException("matrix"); + } + + if (rows <= 0) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "rows"); + } + + if (columns <= 0) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "columns"); + } + + if (matrix.Length < rows * columns) + { + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows * columns), "matrix"); + } + + var work = new double[rows]; + return SafeNativeMethods.z_matrix_norm((byte)norm, rows, columns, matrix, work); + } + + /// + /// Computes the dot product of x and y. + /// + /// The vector x. + /// The vector y. + /// The dot product of x and y. + /// This is equivalent to the DOT BLAS routine. + [SecuritySafeCritical] + public override Complex DotProduct(Complex[] x, Complex[] y) + { + if (y == null) + { + throw new ArgumentNullException("y"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (x.Length != y.Length) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength); + } + + return SafeNativeMethods.z_dot_product(x.Length, x, y); + } + /// /// Adds a scaled vector to another: result = y + alpha*x. /// @@ -164,7 +229,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas var k = transposeA == Transpose.DontTranspose ? columnsA : rowsA; var l = transposeB == Transpose.DontTranspose ? rowsB : columnsB; - if (c.Length != m*n) + if (c.Length != m * n) { throw new ArgumentException(Resources.ArgumentMatrixDimensions); } @@ -199,7 +264,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (data.Length != order*order) + if (data.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "data"); } @@ -226,7 +291,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("a"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -255,7 +320,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -286,7 +351,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("a"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -327,7 +392,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -366,12 +431,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("a"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (b.Length != columnsOfB*order) + if (b.Length != columnsOfB * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -406,7 +471,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -416,7 +481,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); } - if (b.Length != columnsOfB*order) + if (b.Length != columnsOfB * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -449,7 +514,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -484,7 +549,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("b"); } - if (b.Length != orderA*columnsB) + if (b.Length != orderA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -518,7 +583,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("b"); } - if (b.Length != orderA*columnsB) + if (b.Length != orderA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -556,7 +621,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("q"); } - if (r.Length != rowsR*columnsR) + if (r.Length != rowsR * columnsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); } @@ -566,12 +631,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (q.Length != rowsR*rowsR) + if (q.Length != rowsR * rowsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } - var work = new Complex[columnsR*Control.BlockSize]; + var work = new Complex[columnsR * Control.BlockSize]; SafeNativeMethods.z_qr_factor(rowsR, columnsR, r, tau, q, work, work.Length); } @@ -608,7 +673,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (r.Length != rowsR*columnsR) + if (r.Length != rowsR * columnsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); } @@ -618,14 +683,14 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (q.Length != rowsR*rowsR) + if (q.Length != rowsR * rowsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } - if (work.Length < columnsR*Control.BlockSize) + if (work.Length < columnsR * Control.BlockSize) { - work[0] = columnsR*Control.BlockSize; + work[0] = columnsR * Control.BlockSize; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } @@ -633,54 +698,123 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas } /// - /// Solves A*X=B for X using QR factorization of A. + /// Computes the thin QR factorization of A where M > N. /// - /// The A matrix. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The B matrix. - /// The number of columns of B. - /// On exit, the solution matrix. - /// Rows must be greater or equal to columns. - public override void QRSolve(Complex[] a, int rows, int columns, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full) + /// On entry, it is the M by N A matrix to factor. On exit, + /// it is overwritten with the Q matrix of the QR factorization. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// On exit, A N by N matrix that holds the R matrix of the + /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + [SecuritySafeCritical] + public override void ThinQRFactor(Complex[] q, int rowsA, int columnsA, Complex[] r, Complex[] tau) { - if (a == null) + if (r == null) { - throw new ArgumentNullException("a"); + throw new ArgumentNullException("r"); } - if (b == null) + if (q == null) { - throw new ArgumentNullException("b"); + throw new ArgumentNullException("q"); } - if (x == null) + if (q.Length != rowsA * columnsA) { - throw new ArgumentNullException("x"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "q"); } - if (a.Length != rows*columns) + if (tau.Length < Math.Min(rowsA, columnsA)) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (b.Length != rows*columnsB) + if (r.Length != columnsA * columnsA) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r"); + } + + var work = new Complex[columnsA * Control.BlockSize]; + SafeNativeMethods.z_qr_thin_factor(rowsA, columnsA, q, tau, r, work, work.Length); + } + + /// + /// Computes the thin QR factorization of A where M > N. + /// + /// On entry, it is the M by N A matrix to factor. On exit, + /// it is overwritten with the Q matrix of the QR factorization. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// On exit, A N by N matrix that holds the R matrix of the + /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. + /// The work array. The array must have a length of at least N, + /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal + /// work size value. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + [SecuritySafeCritical] + public override void ThinQRFactor(Complex[] q, int rowsA, int columnsA, Complex[] r, Complex[] tau, Complex[] work) + { + if (r == null) + { + throw new ArgumentNullException("r"); } - if (x.Length != columns*columnsB) + if (q == null) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + throw new ArgumentNullException("q"); } - if (rows < columns) + if (work == null) { - throw new ArgumentException(Resources.RowsLessThanColumns); + throw new ArgumentNullException("q"); + } + + if (q.Length != rowsA * columnsA) + { + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "q"); + } + + if (tau.Length < Math.Min(rowsA, columnsA)) + { + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); + } + + if (r.Length != columnsA * columnsA) + { + throw new ArgumentException( + string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r"); + } + + if (work.Length < columnsA * Control.BlockSize) + { + work[0] = columnsA * Control.BlockSize; + throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - var work = new Complex[columns*Control.BlockSize]; - QRSolve(a, rows, columns, b, columnsB, x, work); + SafeNativeMethods.z_qr_thin_factor(rowsA, columnsA, q, tau, r, work, work.Length); + } + + /// + /// Solves A*X=B for X using QR factorization of A. + /// + /// The A matrix. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + /// The type of QR factorization to perform. + /// Rows must be greater or equal to columns. + [SecuritySafeCritical] + public override void QRSolve(Complex[] a, int rows, int columns, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full) + { + var work = new Complex[columns * Control.BlockSize]; + QRSolve(a, rows, columns, b, columnsB, x, work, method); } /// @@ -695,7 +829,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. + /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. + [SecuritySafeCritical] public override void QRSolve(Complex[] a, int rows, int columns, Complex[] b, int columnsB, Complex[] x, Complex[] work, QRMethod method = QRMethod.Full) { if (a == null) @@ -718,17 +854,17 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (a.Length != rows*columns) + if (a.Length != rows * columns) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (b.Length != rows*columnsB) + if (b.Length != rows * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (x.Length != columns*columnsB) + if (x.Length != columns * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); } @@ -740,7 +876,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas if (work.Length < 1) { - work[0] = rows*Control.BlockSize; + work[0] = rows * Control.BlockSize; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } @@ -750,8 +886,8 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// /// Solves A*X=B for X using a previously QR factored matrix. /// - /// The Q matrix obtained by calling . - /// The R matrix obtained by calling . + /// The Q matrix obtained by calling . + /// The R matrix obtained by calling . /// The number of rows in the A matrix. /// The number of columns in the A matrix. /// Contains additional information on Q. Only used for the native solver @@ -759,57 +895,13 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. + /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] public override void QRSolveFactored(Complex[] q, Complex[] r, int rowsR, int columnsR, Complex[] tau, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full) { - if (r == null) - { - throw new ArgumentNullException("r"); - } - - if (q == null) - { - throw new ArgumentNullException("q"); - } - - if (b == null) - { - throw new ArgumentNullException("q"); - } - - if (x == null) - { - throw new ArgumentNullException("q"); - } - - if (r.Length != rowsR*columnsR) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); - } - - if (q.Length != rowsR*rowsR) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); - } - - if (b.Length != rowsR*columnsB) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); - } - - if (x.Length != columnsR*columnsB) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); - } - - if (rowsR < columnsR) - { - throw new ArgumentException(Resources.RowsLessThanColumns); - } - - var work = new Complex[columnsR*Control.BlockSize]; - QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, work); + var work = new Complex[columnsR * Control.BlockSize]; + QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, work, method); } /// @@ -817,9 +909,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// /// The Q matrix obtained by QR factor. This is only used for the managed provider and can be /// null for the native provider. The native provider uses the Q portion stored in the R matrix. - /// The R matrix obtained by calling . - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. + /// The R matrix obtained by calling . + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. /// Contains additional information on Q. Only used for the native solver /// and can be null for the managed provider. /// On entry the B matrix; on exit the X matrix. @@ -828,8 +920,10 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// The work array - only used in the native provider. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. + /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. - public override void QRSolveFactored(Complex[] q, Complex[] r, int rowsR, int columnsR, Complex[] tau, Complex[] b, int columnsB, Complex[] x, Complex[] work, QRMethod method = QRMethod.Full) + [SecuritySafeCritical] + public override void QRSolveFactored(Complex[] q, Complex[] r, int rowsA, int columnsA, Complex[] tau, Complex[] b, int columnsB, Complex[] x, Complex[] work, QRMethod method = QRMethod.Full) { if (r == null) { @@ -856,38 +950,54 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (r.Length != rowsR*columnsR) + int rowsQ, columnsQ, rowsR, columnsR; + if (method == QRMethod.Full) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + rowsQ = columnsQ = rowsR = rowsA; + columnsR = columnsA; + } + else + { + rowsQ = rowsA; + columnsQ = rowsR = columnsR = columnsA; } - if (q.Length != rowsR*rowsR) + if (r.Length != rowsR * columnsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsR * columnsR), "r"); } - if (b.Length != rowsR*columnsB) + if (q.Length != rowsQ * columnsQ) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsQ * columnsQ), "q"); } - if (x.Length != columnsR*columnsB) + if (b.Length != rowsA * columnsB) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsA * columnsB), "b"); } - if (rowsR < columnsR) + if (x.Length != columnsA * columnsB) { - throw new ArgumentException(Resources.RowsLessThanColumns); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, columnsA * columnsB), "x"); } if (work.Length < 1) { - work[0] = rowsR*Control.BlockSize; + work[0] = rowsA * Control.BlockSize; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - SafeNativeMethods.z_qr_solve_factored(rowsR, columnsR, columnsB, r, b, tau, x, work, work.Length); + if (method == QRMethod.Full) + { + SafeNativeMethods.z_qr_solve_factored(rowsA, columnsA, columnsB, r, b, tau, x, work, work.Length); + } + else + { + // we don't have access to the raw Q matrix any more(it is stored in R in the full QR), need to think about this. + // let just call the managed version in the meantime. The heavy lifting has already been done. -marcus + base.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, QRMethod.Thin); + } } /// @@ -926,12 +1036,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("vt"); } - if (u.Length != rowsA*rowsA) + if (u.Length != rowsA * rowsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); } - if (vt.Length != columnsA*columnsA) + if (vt.Length != columnsA * columnsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); } @@ -941,7 +1051,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); } - var work = new Complex[(2*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)]; + var work = new Complex[(2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)]; SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, work); } @@ -971,20 +1081,20 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("x"); } - if (b.Length != rowsA*columnsB) + if (b.Length != rowsA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (x.Length != columnsA*columnsB) + if (x.Length != columnsA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - var work = new Complex[(2*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)]; + var work = new Complex[(2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)]; var s = new Complex[Math.Min(rowsA, columnsA)]; - var u = new Complex[rowsA*rowsA]; - var vt = new Complex[columnsA*columnsA]; + var u = new Complex[rowsA * rowsA]; + var vt = new Complex[columnsA * columnsA]; var clone = new Complex[a.Length]; a.Copy(clone); @@ -1036,12 +1146,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (u.Length != rowsA*rowsA) + if (u.Length != rowsA * rowsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); } - if (vt.Length != columnsA*columnsA) + if (vt.Length != columnsA * columnsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); } @@ -1056,13 +1166,73 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentSingleDimensionArray, "work"); } - if (work.Length < (2*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)) + if (work.Length < (2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)) { - work[0] = (2*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA); + work[0] = (2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA); throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - SafeNativeMethods.z_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length); + if (SafeNativeMethods.z_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length) > 0) + { + throw new NonConvergenceException(); + } + } + + /// + /// Computes the eigenvalues and eigenvectors of a matrix. + /// + /// Whether the matrix is symmetric or not. + /// The order of the matrix. + /// The matrix to decompose. The lenth of the array must be order * order. + /// On output, the matrix contains the eigen vectors. The lenth of the array must be order * order. + /// On output, the eigen values (λ) of matrix in ascending value. The length of the arry must . + /// On output, the block diagonal eigenvalue matrix. The lenth of the array must be order * order. + public override void EigenDecomp(bool isSymmetric, int order, Complex[] matrix, Complex[] matrixEv, Complex[] vectorEv, Complex[] matrixD) + { + if (matrix == null) + { + throw new ArgumentNullException("matrix"); + } + + if (matrix.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrix"); + } + + if (matrixEv == null) + { + throw new ArgumentNullException("matrixEv"); + } + + if (matrixEv.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixEv"); + } + + if (vectorEv == null) + { + throw new ArgumentNullException("vectorEv"); + } + + if (vectorEv.Length != order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order), "vectorEv"); + } + + if (matrixD == null) + { + throw new ArgumentNullException("matrixD"); + } + + if (matrixD.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixD"); + } + + if (SafeNativeMethods.z_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD) > 0) + { + throw new NonConvergenceException(); + } } } } diff --git a/src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Complex32.cs b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex32.cs similarity index 73% rename from src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Complex32.cs rename to src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex32.cs index e5de77d7..9f118068 100644 --- a/src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Complex32.cs +++ b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Complex32.cs @@ -1,10 +1,10 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2011 Math.NET +// Copyright (c) 2009-2015 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -28,20 +28,86 @@ // OTHER DEALINGS IN THE SOFTWARE. // -#if NATIVEGOTO +#if NATIVE using MathNet.Numerics.LinearAlgebra.Factorization; using MathNet.Numerics.Properties; using System; +using System.Numerics; using System.Security; -namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas +namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas { /// - /// GotoBLAS2 linear algebra provider. + /// OpenBLAS linear algebra provider. /// - public partial class GotoBlasLinearAlgebraProvider + public partial class OpenBlasLinearAlgebraProvider { + /// + /// Computes the requested of the matrix. + /// + /// The type of norm to compute. + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The matrix to compute the norm from. + /// + /// The requested of the matrix. + /// + [SecuritySafeCritical] + public override double MatrixNorm(Norm norm, int rows, int columns, Complex32[] matrix) + { + if (matrix == null) + { + throw new ArgumentNullException("matrix"); + } + + if (rows <= 0) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "rows"); + } + + if (columns <= 0) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "columns"); + } + + if (matrix.Length < rows * columns) + { + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows * columns), "matrix"); + } + + var work = new float[rows]; + return SafeNativeMethods.c_matrix_norm((byte)norm, rows, columns, matrix, work); + } + + /// + /// Computes the dot product of x and y. + /// + /// The vector x. + /// The vector y. + /// The dot product of x and y. + /// This is equivalent to the DOT BLAS routine. + [SecuritySafeCritical] + public override Complex32 DotProduct(Complex32[] x, Complex32[] y) + { + if (y == null) + { + throw new ArgumentNullException("y"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (x.Length != y.Length) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength); + } + + return SafeNativeMethods.c_dot_product(x.Length, x, y); + } + /// /// Adds a scaled vector to another: result = y + alpha*x. /// @@ -163,7 +229,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas var k = transposeA == Transpose.DontTranspose ? columnsA : rowsA; var l = transposeB == Transpose.DontTranspose ? rowsB : columnsB; - if (c.Length != m*n) + if (c.Length != m * n) { throw new ArgumentException(Resources.ArgumentMatrixDimensions); } @@ -198,7 +264,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (data.Length != order*order) + if (data.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "data"); } @@ -225,7 +291,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("a"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -254,7 +320,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -285,7 +351,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("a"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -326,7 +392,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -365,12 +431,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("a"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (b.Length != columnsOfB*order) + if (b.Length != columnsOfB * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -405,7 +471,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -415,7 +481,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); } - if (b.Length != columnsOfB*order) + if (b.Length != columnsOfB * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -448,7 +514,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -483,7 +549,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("b"); } - if (b.Length != orderA*columnsB) + if (b.Length != orderA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -517,7 +583,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("b"); } - if (b.Length != orderA*columnsB) + if (b.Length != orderA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -555,7 +621,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("q"); } - if (r.Length != rowsR*columnsR) + if (r.Length != rowsR * columnsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); } @@ -565,12 +631,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (q.Length != rowsR*rowsR) + if (q.Length != rowsR * rowsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } - var work = new Complex32[columnsR*Control.BlockSize]; + var work = new Complex32[columnsR * Control.BlockSize]; SafeNativeMethods.c_qr_factor(rowsR, columnsR, r, tau, q, work, work.Length); } @@ -607,7 +673,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (r.Length != rowsR*columnsR) + if (r.Length != rowsR * columnsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); } @@ -617,14 +683,14 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (q.Length != rowsR*rowsR) + if (q.Length != rowsR * rowsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } - if (work.Length < columnsR*Control.BlockSize) + if (work.Length < columnsR * Control.BlockSize) { - work[0] = columnsR*Control.BlockSize; + work[0] = columnsR * Control.BlockSize; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } @@ -632,54 +698,123 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas } /// - /// Solves A*X=B for X using QR factorization of A. + /// Computes the thin QR factorization of A where M > N. /// - /// The A matrix. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The B matrix. - /// The number of columns of B. - /// On exit, the solution matrix. - /// Rows must be greater or equal to columns. - public override void QRSolve(Complex32[] a, int rows, int columns, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full) + /// On entry, it is the M by N A matrix to factor. On exit, + /// it is overwritten with the Q matrix of the QR factorization. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// On exit, A N by N matrix that holds the R matrix of the + /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + [SecuritySafeCritical] + public override void ThinQRFactor(Complex32[] q, int rowsA, int columnsA, Complex32[] r, Complex32[] tau) { - if (a == null) + if (r == null) { - throw new ArgumentNullException("a"); + throw new ArgumentNullException("r"); } - if (b == null) + if (q == null) { - throw new ArgumentNullException("b"); + throw new ArgumentNullException("q"); } - if (x == null) + if (q.Length != rowsA * columnsA) { - throw new ArgumentNullException("x"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "q"); } - if (a.Length != rows*columns) + if (tau.Length < Math.Min(rowsA, columnsA)) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (b.Length != rows*columnsB) + if (r.Length != columnsA * columnsA) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r"); + } + + var work = new Complex32[columnsA * Control.BlockSize]; + SafeNativeMethods.c_qr_thin_factor(rowsA, columnsA, q, tau, r, work, work.Length); + } + + /// + /// Computes the thin QR factorization of A where M > N. + /// + /// On entry, it is the M by N A matrix to factor. On exit, + /// it is overwritten with the Q matrix of the QR factorization. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// On exit, A N by N matrix that holds the R matrix of the + /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. + /// The work array. The array must have a length of at least N, + /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal + /// work size value. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + [SecuritySafeCritical] + public override void ThinQRFactor(Complex32[] q, int rowsA, int columnsA, Complex32[] r, Complex32[] tau, Complex32[] work) + { + if (r == null) + { + throw new ArgumentNullException("r"); } - if (x.Length != columns*columnsB) + if (q == null) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + throw new ArgumentNullException("q"); } - if (rows < columns) + if (work == null) { - throw new ArgumentException(Resources.RowsLessThanColumns); + throw new ArgumentNullException("q"); + } + + if (q.Length != rowsA * columnsA) + { + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "q"); + } + + if (tau.Length < Math.Min(rowsA, columnsA)) + { + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); + } + + if (r.Length != columnsA * columnsA) + { + throw new ArgumentException( + string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r"); + } + + if (work.Length < columnsA * Control.BlockSize) + { + work[0] = columnsA * Control.BlockSize; + throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - var work = new Complex32[columns*Control.BlockSize]; - QRSolve(a, rows, columns, b, columnsB, x, work); + SafeNativeMethods.c_qr_thin_factor(rowsA, columnsA, q, tau, r, work, work.Length); + } + + /// + /// Solves A*X=B for X using QR factorization of A. + /// + /// The A matrix. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + /// The type of QR factorization to perform. + /// Rows must be greater or equal to columns. + [SecuritySafeCritical] + public override void QRSolve(Complex32[] a, int rows, int columns, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full) + { + var work = new Complex32[columns * Control.BlockSize]; + QRSolve(a, rows, columns, b, columnsB, x, work, method); } /// @@ -694,7 +829,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. + /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. + [SecuritySafeCritical] public override void QRSolve(Complex32[] a, int rows, int columns, Complex32[] b, int columnsB, Complex32[] x, Complex32[] work, QRMethod method = QRMethod.Full) { if (a == null) @@ -717,17 +854,17 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (a.Length != rows*columns) + if (a.Length != rows * columns) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (b.Length != rows*columnsB) + if (b.Length != rows * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (x.Length != columns*columnsB) + if (x.Length != columns * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); } @@ -739,7 +876,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas if (work.Length < 1) { - work[0] = rows*Control.BlockSize; + work[0] = rows * Control.BlockSize; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } @@ -749,8 +886,8 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// /// Solves A*X=B for X using a previously QR factored matrix. /// - /// The Q matrix obtained by calling . - /// The R matrix obtained by calling . + /// The Q matrix obtained by calling . + /// The R matrix obtained by calling . /// The number of rows in the A matrix. /// The number of columns in the A matrix. /// Contains additional information on Q. Only used for the native solver @@ -758,57 +895,13 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. + /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] public override void QRSolveFactored(Complex32[] q, Complex32[] r, int rowsR, int columnsR, Complex32[] tau, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full) { - if (r == null) - { - throw new ArgumentNullException("r"); - } - - if (q == null) - { - throw new ArgumentNullException("q"); - } - - if (b == null) - { - throw new ArgumentNullException("q"); - } - - if (x == null) - { - throw new ArgumentNullException("q"); - } - - if (r.Length != rowsR*columnsR) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); - } - - if (q.Length != rowsR*rowsR) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); - } - - if (b.Length != rowsR*columnsB) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); - } - - if (x.Length != columnsR*columnsB) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); - } - - if (rowsR < columnsR) - { - throw new ArgumentException(Resources.RowsLessThanColumns); - } - - var work = new Complex32[columnsR*Control.BlockSize]; - QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, work); + var work = new Complex32[columnsR * Control.BlockSize]; + QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, work, method); } /// @@ -816,9 +909,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// /// The Q matrix obtained by QR factor. This is only used for the managed provider and can be /// null for the native provider. The native provider uses the Q portion stored in the R matrix. - /// The R matrix obtained by calling . - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. + /// The R matrix obtained by calling . + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. /// Contains additional information on Q. Only used for the native solver /// and can be null for the managed provider. /// On entry the B matrix; on exit the X matrix. @@ -827,8 +920,10 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// The work array - only used in the native provider. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. + /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. - public override void QRSolveFactored(Complex32[] q, Complex32[] r, int rowsR, int columnsR, Complex32[] tau, Complex32[] b, int columnsB, Complex32[] x, Complex32[] work, QRMethod method = QRMethod.Full) + [SecuritySafeCritical] + public override void QRSolveFactored(Complex32[] q, Complex32[] r, int rowsA, int columnsA, Complex32[] tau, Complex32[] b, int columnsB, Complex32[] x, Complex32[] work, QRMethod method = QRMethod.Full) { if (r == null) { @@ -855,38 +950,54 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (r.Length != rowsR*columnsR) + int rowsQ, columnsQ, rowsR, columnsR; + if (method == QRMethod.Full) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + rowsQ = columnsQ = rowsR = rowsA; + columnsR = columnsA; + } + else + { + rowsQ = rowsA; + columnsQ = rowsR = columnsR = columnsA; } - if (q.Length != rowsR*rowsR) + if (r.Length != rowsR * columnsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsR * columnsR), "r"); } - if (b.Length != rowsR*columnsB) + if (q.Length != rowsQ * columnsQ) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsQ * columnsQ), "q"); } - if (x.Length != columnsR*columnsB) + if (b.Length != rowsA * columnsB) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsA * columnsB), "b"); } - if (rowsR < columnsR) + if (x.Length != columnsA * columnsB) { - throw new ArgumentException(Resources.RowsLessThanColumns); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, columnsA * columnsB), "x"); } if (work.Length < 1) { - work[0] = rowsR*Control.BlockSize; + work[0] = rowsA * Control.BlockSize; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - SafeNativeMethods.c_qr_solve_factored(rowsR, columnsR, columnsB, r, b, tau, x, work, work.Length); + if (method == QRMethod.Full) + { + SafeNativeMethods.c_qr_solve_factored(rowsA, columnsA, columnsB, r, b, tau, x, work, work.Length); + } + else + { + // we don't have access to the raw Q matrix any more(it is stored in R in the full QR), need to think about this. + // let just call the managed version in the meantime. The heavy lifting has already been done. -marcus + base.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, QRMethod.Thin); + } } /// @@ -925,12 +1036,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("vt"); } - if (u.Length != rowsA*rowsA) + if (u.Length != rowsA * rowsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); } - if (vt.Length != columnsA*columnsA) + if (vt.Length != columnsA * columnsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); } @@ -940,7 +1051,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); } - var work = new Complex32[(2*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)]; + var work = new Complex32[(2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)]; SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, work); } @@ -970,20 +1081,20 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("x"); } - if (b.Length != rowsA*columnsB) + if (b.Length != rowsA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (x.Length != columnsA*columnsB) + if (x.Length != columnsA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - var work = new Complex32[(2*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)]; + var work = new Complex32[(2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)]; var s = new Complex32[Math.Min(rowsA, columnsA)]; - var u = new Complex32[rowsA*rowsA]; - var vt = new Complex32[columnsA*columnsA]; + var u = new Complex32[rowsA * rowsA]; + var vt = new Complex32[columnsA * columnsA]; var clone = new Complex32[a.Length]; a.Copy(clone); @@ -1035,12 +1146,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (u.Length != rowsA*rowsA) + if (u.Length != rowsA * rowsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); } - if (vt.Length != columnsA*columnsA) + if (vt.Length != columnsA * columnsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); } @@ -1055,13 +1166,73 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentSingleDimensionArray, "work"); } - if (work.Length < (2*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)) + if (work.Length < (2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)) { - work[0] = (2*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA); + work[0] = (2 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA); throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - SafeNativeMethods.c_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length); + if (SafeNativeMethods.c_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length) > 0) + { + throw new NonConvergenceException(); + } + } + + /// + /// Computes the eigenvalues and eigenvectors of a matrix. + /// + /// Whether the matrix is symmetric or not. + /// The order of the matrix. + /// The matrix to decompose. The lenth of the array must be order * order. + /// On output, the matrix contains the eigen vectors. The lenth of the array must be order * order. + /// On output, the eigen values (λ) of matrix in ascending value. The length of the arry must . + /// On output, the block diagonal eigenvalue matrix. The lenth of the array must be order * order. + public override void EigenDecomp(bool isSymmetric, int order, Complex32[] matrix, Complex32[] matrixEv, Complex[] vectorEv, Complex32[] matrixD) + { + if (matrix == null) + { + throw new ArgumentNullException("matrix"); + } + + if (matrix.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrix"); + } + + if (matrixEv == null) + { + throw new ArgumentNullException("matrixEv"); + } + + if (matrixEv.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixEv"); + } + + if (vectorEv == null) + { + throw new ArgumentNullException("vectorEv"); + } + + if (vectorEv.Length != order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order), "vectorEv"); + } + + if (matrixD == null) + { + throw new ArgumentNullException("matrixD"); + } + + if (matrixD.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixD"); + } + + if (SafeNativeMethods.c_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD) > 0) + { + throw new NonConvergenceException(); + } } } } diff --git a/src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Double.cs b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Double.cs similarity index 73% rename from src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Double.cs rename to src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Double.cs index 8dfc7a2b..0143e0be 100644 --- a/src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Double.cs +++ b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Double.cs @@ -1,10 +1,10 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2011 Math.NET +// Copyright (c) 2009-2015 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -28,20 +28,86 @@ // OTHER DEALINGS IN THE SOFTWARE. // -#if NATIVEGOTO +#if NATIVE using MathNet.Numerics.LinearAlgebra.Factorization; using MathNet.Numerics.Properties; using System; +using System.Numerics; using System.Security; -namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas +namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas { /// - /// GotoBLAS2 linear algebra provider. + /// OpenBLAS linear algebra provider. /// - public partial class GotoBlasLinearAlgebraProvider + public partial class OpenBlasLinearAlgebraProvider { + /// + /// Computes the requested of the matrix. + /// + /// The type of norm to compute. + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The matrix to compute the norm from. + /// + /// The requested of the matrix. + /// + [SecuritySafeCritical] + public override double MatrixNorm(Norm norm, int rows, int columns, double[] matrix) + { + if (matrix == null) + { + throw new ArgumentNullException("matrix"); + } + + if (rows <= 0) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "rows"); + } + + if (columns <= 0) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "columns"); + } + + if (matrix.Length < rows * columns) + { + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows * columns), "matrix"); + } + + var work = new double[rows]; + return SafeNativeMethods.d_matrix_norm((byte)norm, rows, columns, matrix, work); + } + + /// + /// Computes the dot product of x and y. + /// + /// The vector x. + /// The vector y. + /// The dot product of x and y. + /// This is equivalent to the DOT BLAS routine. + [SecuritySafeCritical] + public override double DotProduct(double[] x, double[] y) + { + if (y == null) + { + throw new ArgumentNullException("y"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (x.Length != y.Length) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength); + } + + return SafeNativeMethods.d_dot_product(x.Length, x, y); + } + /// /// Adds a scaled vector to another: result = y + alpha*x. /// @@ -163,7 +229,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas var k = transposeA == Transpose.DontTranspose ? columnsA : rowsA; var l = transposeB == Transpose.DontTranspose ? rowsB : columnsB; - if (c.Length != m*n) + if (c.Length != m * n) { throw new ArgumentException(Resources.ArgumentMatrixDimensions); } @@ -198,7 +264,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (data.Length != order*order) + if (data.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "data"); } @@ -225,7 +291,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("a"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -254,7 +320,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -285,7 +351,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("a"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -326,7 +392,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -365,12 +431,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("a"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (b.Length != columnsOfB*order) + if (b.Length != columnsOfB * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -405,7 +471,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -415,7 +481,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); } - if (b.Length != columnsOfB*order) + if (b.Length != columnsOfB * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -448,7 +514,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -483,7 +549,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("b"); } - if (b.Length != orderA*columnsB) + if (b.Length != orderA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -517,7 +583,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("b"); } - if (b.Length != orderA*columnsB) + if (b.Length != orderA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -555,7 +621,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("q"); } - if (r.Length != rowsR*columnsR) + if (r.Length != rowsR * columnsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); } @@ -565,12 +631,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (q.Length != rowsR*rowsR) + if (q.Length != rowsR * rowsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } - var work = new double[columnsR*Control.BlockSize]; + var work = new double[columnsR * Control.BlockSize]; SafeNativeMethods.d_qr_factor(rowsR, columnsR, r, tau, q, work, work.Length); } @@ -607,7 +673,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (r.Length != rowsR*columnsR) + if (r.Length != rowsR * columnsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); } @@ -617,14 +683,14 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (q.Length != rowsR*rowsR) + if (q.Length != rowsR * rowsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } - if (work.Length < columnsR*Control.BlockSize) + if (work.Length < columnsR * Control.BlockSize) { - work[0] = columnsR*Control.BlockSize; + work[0] = columnsR * Control.BlockSize; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } @@ -632,54 +698,123 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas } /// - /// Solves A*X=B for X using QR factorization of A. + /// Computes the thin QR factorization of A where M > N. /// - /// The A matrix. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The B matrix. - /// The number of columns of B. - /// On exit, the solution matrix. - /// Rows must be greater or equal to columns. - public override void QRSolve(double[] a, int rows, int columns, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full) + /// On entry, it is the M by N A matrix to factor. On exit, + /// it is overwritten with the Q matrix of the QR factorization. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// On exit, A N by N matrix that holds the R matrix of the + /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + [SecuritySafeCritical] + public override void ThinQRFactor(double[] q, int rowsA, int columnsA, double[] r, double[] tau) { - if (a == null) + if (r == null) { - throw new ArgumentNullException("a"); + throw new ArgumentNullException("r"); } - if (b == null) + if (q == null) { - throw new ArgumentNullException("b"); + throw new ArgumentNullException("q"); } - if (x == null) + if (q.Length != rowsA * columnsA) { - throw new ArgumentNullException("x"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "q"); } - if (a.Length != rows*columns) + if (tau.Length < Math.Min(rowsA, columnsA)) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (b.Length != rows*columnsB) + if (r.Length != columnsA * columnsA) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r"); + } + + var work = new double[columnsA * Control.BlockSize]; + SafeNativeMethods.d_qr_thin_factor(rowsA, columnsA, q, tau, r, work, work.Length); + } + + /// + /// Computes the thin QR factorization of A where M > N. + /// + /// On entry, it is the M by N A matrix to factor. On exit, + /// it is overwritten with the Q matrix of the QR factorization. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// On exit, A N by N matrix that holds the R matrix of the + /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. + /// The work array. The array must have a length of at least N, + /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal + /// work size value. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + [SecuritySafeCritical] + public override void ThinQRFactor(double[] q, int rowsA, int columnsA, double[] r, double[] tau, double[] work) + { + if (r == null) + { + throw new ArgumentNullException("r"); } - if (x.Length != columns*columnsB) + if (q == null) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + throw new ArgumentNullException("q"); } - if (rows < columns) + if (work == null) { - throw new ArgumentException(Resources.RowsLessThanColumns); + throw new ArgumentNullException("q"); + } + + if (q.Length != rowsA * columnsA) + { + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "q"); + } + + if (tau.Length < Math.Min(rowsA, columnsA)) + { + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); + } + + if (r.Length != columnsA * columnsA) + { + throw new ArgumentException( + string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r"); + } + + if (work.Length < columnsA * Control.BlockSize) + { + work[0] = columnsA * Control.BlockSize; + throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - var work = new double[columns*Control.BlockSize]; - QRSolve(a, rows, columns, b, columnsB, x, work); + SafeNativeMethods.d_qr_thin_factor(rowsA, columnsA, q, tau, r, work, work.Length); + } + + /// + /// Solves A*X=B for X using QR factorization of A. + /// + /// The A matrix. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + /// The type of QR factorization to perform. + /// Rows must be greater or equal to columns. + [SecuritySafeCritical] + public override void QRSolve(double[] a, int rows, int columns, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full) + { + var work = new double[columns * Control.BlockSize]; + QRSolve(a, rows, columns, b, columnsB, x, work, method); } /// @@ -694,7 +829,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. + /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. + [SecuritySafeCritical] public override void QRSolve(double[] a, int rows, int columns, double[] b, int columnsB, double[] x, double[] work, QRMethod method = QRMethod.Full) { if (a == null) @@ -717,17 +854,17 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (a.Length != rows*columns) + if (a.Length != rows * columns) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (b.Length != rows*columnsB) + if (b.Length != rows * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (x.Length != columns*columnsB) + if (x.Length != columns * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); } @@ -739,7 +876,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas if (work.Length < 1) { - work[0] = rows*Control.BlockSize; + work[0] = rows * Control.BlockSize; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } @@ -749,8 +886,8 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// /// Solves A*X=B for X using a previously QR factored matrix. /// - /// The Q matrix obtained by calling . - /// The R matrix obtained by calling . + /// The Q matrix obtained by calling . + /// The R matrix obtained by calling . /// The number of rows in the A matrix. /// The number of columns in the A matrix. /// Contains additional information on Q. Only used for the native solver @@ -758,57 +895,13 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. + /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] public override void QRSolveFactored(double[] q, double[] r, int rowsR, int columnsR, double[] tau, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full) { - if (r == null) - { - throw new ArgumentNullException("r"); - } - - if (q == null) - { - throw new ArgumentNullException("q"); - } - - if (b == null) - { - throw new ArgumentNullException("q"); - } - - if (x == null) - { - throw new ArgumentNullException("q"); - } - - if (r.Length != rowsR*columnsR) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); - } - - if (q.Length != rowsR*rowsR) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); - } - - if (b.Length != rowsR*columnsB) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); - } - - if (x.Length != columnsR*columnsB) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); - } - - if (rowsR < columnsR) - { - throw new ArgumentException(Resources.RowsLessThanColumns); - } - - var work = new double[columnsR*Control.BlockSize]; - QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, work); + var work = new double[columnsR * Control.BlockSize]; + QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, work, method); } /// @@ -816,9 +909,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// /// The Q matrix obtained by QR factor. This is only used for the managed provider and can be /// null for the native provider. The native provider uses the Q portion stored in the R matrix. - /// The R matrix obtained by calling . - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. + /// The R matrix obtained by calling . + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. /// Contains additional information on Q. Only used for the native solver /// and can be null for the managed provider. /// On entry the B matrix; on exit the X matrix. @@ -827,8 +920,10 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// The work array - only used in the native provider. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. + /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. - public override void QRSolveFactored(double[] q, double[] r, int rowsR, int columnsR, double[] tau, double[] b, int columnsB, double[] x, double[] work, QRMethod method = QRMethod.Full) + [SecuritySafeCritical] + public override void QRSolveFactored(double[] q, double[] r, int rowsA, int columnsA, double[] tau, double[] b, int columnsB, double[] x, double[] work, QRMethod method = QRMethod.Full) { if (r == null) { @@ -855,38 +950,54 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (r.Length != rowsR*columnsR) + int rowsQ, columnsQ, rowsR, columnsR; + if (method == QRMethod.Full) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + rowsQ = columnsQ = rowsR = rowsA; + columnsR = columnsA; + } + else + { + rowsQ = rowsA; + columnsQ = rowsR = columnsR = columnsA; } - if (q.Length != rowsR*rowsR) + if (r.Length != rowsR * columnsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsR * columnsR), "r"); } - if (b.Length != rowsR*columnsB) + if (q.Length != rowsQ * columnsQ) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsQ * columnsQ), "q"); } - if (x.Length != columnsR*columnsB) + if (b.Length != rowsA * columnsB) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsA * columnsB), "b"); } - if (rowsR < columnsR) + if (x.Length != columnsA * columnsB) { - throw new ArgumentException(Resources.RowsLessThanColumns); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, columnsA * columnsB), "x"); } if (work.Length < 1) { - work[0] = rowsR*Control.BlockSize; + work[0] = rowsA * Control.BlockSize; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - SafeNativeMethods.d_qr_solve_factored(rowsR, columnsR, columnsB, r, b, tau, x, work, work.Length); + if (method == QRMethod.Full) + { + SafeNativeMethods.d_qr_solve_factored(rowsA, columnsA, columnsB, r, b, tau, x, work, work.Length); + } + else + { + // we don't have access to the raw Q matrix any more(it is stored in R in the full QR), need to think about this. + // let just call the managed version in the meantime. The heavy lifting has already been done. -marcus + base.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, QRMethod.Thin); + } } /// @@ -925,12 +1036,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("vt"); } - if (u.Length != rowsA*rowsA) + if (u.Length != rowsA * rowsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); } - if (vt.Length != columnsA*columnsA) + if (vt.Length != columnsA * columnsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); } @@ -940,7 +1051,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); } - var work = new double[Math.Max((3*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA), 5*Math.Min(rowsA, columnsA))]; + var work = new double[Math.Max((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA), 5 * Math.Min(rowsA, columnsA))]; SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, work); } @@ -970,20 +1081,20 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("x"); } - if (b.Length != rowsA*columnsB) + if (b.Length != rowsA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (x.Length != columnsA*columnsB) + if (x.Length != columnsA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - var work = new double[Math.Max((3*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA), 5*Math.Min(rowsA, columnsA))]; + var work = new double[Math.Max((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA), 5 * Math.Min(rowsA, columnsA))]; var s = new double[Math.Min(rowsA, columnsA)]; - var u = new double[rowsA*rowsA]; - var vt = new double[columnsA*columnsA]; + var u = new double[rowsA * rowsA]; + var vt = new double[columnsA * columnsA]; var clone = new double[a.Length]; a.Copy(clone); @@ -1035,12 +1146,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (u.Length != rowsA*rowsA) + if (u.Length != rowsA * rowsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); } - if (vt.Length != columnsA*columnsA) + if (vt.Length != columnsA * columnsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); } @@ -1055,13 +1166,73 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentSingleDimensionArray, "work"); } - if (work.Length < Math.Max((3*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA), 5*Math.Min(rowsA, columnsA))) + if (work.Length < Math.Max((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA), 5 * Math.Min(rowsA, columnsA))) { - work[0] = Math.Max((3*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA), 5*Math.Min(rowsA, columnsA)); + work[0] = Math.Max((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA), 5 * Math.Min(rowsA, columnsA)); throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - SafeNativeMethods.d_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length); + if (SafeNativeMethods.d_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length) > 0) + { + throw new NonConvergenceException(); + } + } + + /// + /// Computes the eigenvalues and eigenvectors of a matrix. + /// + /// Whether the matrix is symmetric or not. + /// The order of the matrix. + /// The matrix to decompose. The lenth of the array must be order * order. + /// On output, the matrix contains the eigen vectors. The lenth of the array must be order * order. + /// On output, the eigen values (λ) of matrix in ascending value. The length of the arry must . + /// On output, the block diagonal eigenvalue matrix. The lenth of the array must be order * order. + public override void EigenDecomp(bool isSymmetric, int order, double[] matrix, double[] matrixEv, Complex[] vectorEv, double[] matrixD) + { + if (matrix == null) + { + throw new ArgumentNullException("matrix"); + } + + if (matrix.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrix"); + } + + if (matrixEv == null) + { + throw new ArgumentNullException("matrixEv"); + } + + if (matrixEv.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixEv"); + } + + if (vectorEv == null) + { + throw new ArgumentNullException("vectorEv"); + } + + if (vectorEv.Length != order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order), "vectorEv"); + } + + if (matrixD == null) + { + throw new ArgumentNullException("matrixD"); + } + + if (matrixD.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixD"); + } + + if (SafeNativeMethods.d_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD) > 0) + { + throw new NonConvergenceException(); + } } } } diff --git a/src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Single.cs b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Single.cs similarity index 73% rename from src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Single.cs rename to src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Single.cs index 6c4af16e..61c934c1 100644 --- a/src/Numerics/Providers/LinearAlgebra/GotoBlas/GotoBlasLinearAlgebraProvider.Single.cs +++ b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.Single.cs @@ -1,10 +1,10 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2011 Math.NET +// Copyright (c) 2009-2015 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -28,20 +28,86 @@ // OTHER DEALINGS IN THE SOFTWARE. // -#if NATIVEGOTO +#if NATIVE using MathNet.Numerics.LinearAlgebra.Factorization; using MathNet.Numerics.Properties; using System; +using System.Numerics; using System.Security; -namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas +namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas { /// - /// GotoBLAS2 linear algebra provider. + /// OpenBLAS linear algebra provider. /// - public partial class GotoBlasLinearAlgebraProvider + public partial class OpenBlasLinearAlgebraProvider { + /// + /// Computes the requested of the matrix. + /// + /// The type of norm to compute. + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The matrix to compute the norm from. + /// + /// The requested of the matrix. + /// + [SecuritySafeCritical] + public override double MatrixNorm(Norm norm, int rows, int columns, float[] matrix) + { + if (matrix == null) + { + throw new ArgumentNullException("matrix"); + } + + if (rows <= 0) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "rows"); + } + + if (columns <= 0) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "columns"); + } + + if (matrix.Length < rows * columns) + { + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows * columns), "matrix"); + } + + var work = new float[rows]; + return SafeNativeMethods.s_matrix_norm((byte)norm, rows, columns, matrix, work); + } + + /// + /// Computes the dot product of x and y. + /// + /// The vector x. + /// The vector y. + /// The dot product of x and y. + /// This is equivalent to the DOT BLAS routine. + [SecuritySafeCritical] + public override float DotProduct(float[] x, float[] y) + { + if (y == null) + { + throw new ArgumentNullException("y"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (x.Length != y.Length) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength); + } + + return SafeNativeMethods.s_dot_product(x.Length, x, y); + } + /// /// Adds a scaled vector to another: result = y + alpha*x. /// @@ -163,7 +229,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas var k = transposeA == Transpose.DontTranspose ? columnsA : rowsA; var l = transposeB == Transpose.DontTranspose ? rowsB : columnsB; - if (c.Length != m*n) + if (c.Length != m * n) { throw new ArgumentException(Resources.ArgumentMatrixDimensions); } @@ -198,7 +264,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (data.Length != order*order) + if (data.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "data"); } @@ -225,7 +291,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("a"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -254,7 +320,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -285,7 +351,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("a"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -326,7 +392,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -365,12 +431,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("a"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (b.Length != columnsOfB*order) + if (b.Length != columnsOfB * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -405,7 +471,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("ipiv"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -415,7 +481,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); } - if (b.Length != columnsOfB*order) + if (b.Length != columnsOfB * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -448,7 +514,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); } - if (a.Length != order*order) + if (a.Length != order * order) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } @@ -483,7 +549,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("b"); } - if (b.Length != orderA*columnsB) + if (b.Length != orderA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -517,7 +583,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("b"); } - if (b.Length != orderA*columnsB) + if (b.Length != orderA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } @@ -555,7 +621,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("q"); } - if (r.Length != rowsR*columnsR) + if (r.Length != rowsR * columnsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); } @@ -565,12 +631,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (q.Length != rowsR*rowsR) + if (q.Length != rowsR * rowsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } - var work = new float[columnsR*Control.BlockSize]; + var work = new float[columnsR * Control.BlockSize]; SafeNativeMethods.s_qr_factor(rowsR, columnsR, r, tau, q, work, work.Length); } @@ -607,7 +673,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (r.Length != rowsR*columnsR) + if (r.Length != rowsR * columnsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r"); } @@ -617,14 +683,14 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (q.Length != rowsR*rowsR) + if (q.Length != rowsR * rowsR) { throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q"); } - if (work.Length < columnsR*Control.BlockSize) + if (work.Length < columnsR * Control.BlockSize) { - work[0] = columnsR*Control.BlockSize; + work[0] = columnsR * Control.BlockSize; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } @@ -632,54 +698,123 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas } /// - /// Solves A*X=B for X using QR factorization of A. + /// Computes the thin QR factorization of A where M > N. /// - /// The A matrix. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The B matrix. - /// The number of columns of B. - /// On exit, the solution matrix. - /// Rows must be greater or equal to columns. - public override void QRSolve(float[] a, int rows, int columns, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full) + /// On entry, it is the M by N A matrix to factor. On exit, + /// it is overwritten with the Q matrix of the QR factorization. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// On exit, A N by N matrix that holds the R matrix of the + /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + [SecuritySafeCritical] + public override void ThinQRFactor(float[] q, int rowsA, int columnsA, float[] r, float[] tau) { - if (a == null) + if (r == null) { - throw new ArgumentNullException("a"); + throw new ArgumentNullException("r"); } - if (b == null) + if (q == null) { - throw new ArgumentNullException("b"); + throw new ArgumentNullException("q"); } - if (x == null) + if (q.Length != rowsA * columnsA) { - throw new ArgumentNullException("x"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "q"); } - if (a.Length != rows*columns) + if (tau.Length < Math.Min(rowsA, columnsA)) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); } - if (b.Length != rows*columnsB) + if (r.Length != columnsA * columnsA) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r"); + } + + var work = new float[columnsA * Control.BlockSize]; + SafeNativeMethods.s_qr_thin_factor(rowsA, columnsA, q, tau, r, work, work.Length); + } + + /// + /// Computes the thin QR factorization of A where M > N. + /// + /// On entry, it is the M by N A matrix to factor. On exit, + /// it is overwritten with the Q matrix of the QR factorization. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// On exit, A N by N matrix that holds the R matrix of the + /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. + /// The work array. The array must have a length of at least N, + /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal + /// work size value. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + [SecuritySafeCritical] + public override void ThinQRFactor(float[] q, int rowsA, int columnsA, float[] r, float[] tau, float[] work) + { + if (r == null) + { + throw new ArgumentNullException("r"); } - if (x.Length != columns*columnsB) + if (q == null) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + throw new ArgumentNullException("q"); } - if (rows < columns) + if (work == null) { - throw new ArgumentException(Resources.RowsLessThanColumns); + throw new ArgumentNullException("q"); + } + + if (q.Length != rowsA * columnsA) + { + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "q"); + } + + if (tau.Length < Math.Min(rowsA, columnsA)) + { + throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau"); + } + + if (r.Length != columnsA * columnsA) + { + throw new ArgumentException( + string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r"); + } + + if (work.Length < columnsA * Control.BlockSize) + { + work[0] = columnsA * Control.BlockSize; + throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - var work = new float[columns*Control.BlockSize]; - QRSolve(a, rows, columns, b, columnsB, x, work); + SafeNativeMethods.s_qr_thin_factor(rowsA, columnsA, q, tau, r, work, work.Length); + } + + /// + /// Solves A*X=B for X using QR factorization of A. + /// + /// The A matrix. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + /// The type of QR factorization to perform. + /// Rows must be greater or equal to columns. + [SecuritySafeCritical] + public override void QRSolve(float[] a, int rows, int columns, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full) + { + var work = new float[columns * Control.BlockSize]; + QRSolve(a, rows, columns, b, columnsB, x, work, method); } /// @@ -694,7 +829,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. + /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. + [SecuritySafeCritical] public override void QRSolve(float[] a, int rows, int columns, float[] b, int columnsB, float[] x, float[] work, QRMethod method = QRMethod.Full) { if (a == null) @@ -717,17 +854,17 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (a.Length != rows*columns) + if (a.Length != rows * columns) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); } - if (b.Length != rows*columnsB) + if (b.Length != rows * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (x.Length != columns*columnsB) + if (x.Length != columns * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); } @@ -739,7 +876,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas if (work.Length < 1) { - work[0] = rows*Control.BlockSize; + work[0] = rows * Control.BlockSize; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } @@ -749,8 +886,8 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// /// Solves A*X=B for X using a previously QR factored matrix. /// - /// The Q matrix obtained by calling . - /// The R matrix obtained by calling . + /// The Q matrix obtained by calling . + /// The R matrix obtained by calling . /// The number of rows in the A matrix. /// The number of columns in the A matrix. /// Contains additional information on Q. Only used for the native solver @@ -758,57 +895,13 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. + /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] public override void QRSolveFactored(float[] q, float[] r, int rowsR, int columnsR, float[] tau, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full) { - if (r == null) - { - throw new ArgumentNullException("r"); - } - - if (q == null) - { - throw new ArgumentNullException("q"); - } - - if (b == null) - { - throw new ArgumentNullException("q"); - } - - if (x == null) - { - throw new ArgumentNullException("q"); - } - - if (r.Length != rowsR*columnsR) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); - } - - if (q.Length != rowsR*rowsR) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); - } - - if (b.Length != rowsR*columnsB) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); - } - - if (x.Length != columnsR*columnsB) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); - } - - if (rowsR < columnsR) - { - throw new ArgumentException(Resources.RowsLessThanColumns); - } - - var work = new float[columnsR*Control.BlockSize]; - QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, work); + var work = new float[columnsR * Control.BlockSize]; + QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, work, method); } /// @@ -816,9 +909,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// /// The Q matrix obtained by QR factor. This is only used for the managed provider and can be /// null for the native provider. The native provider uses the Q portion stored in the R matrix. - /// The R matrix obtained by calling . - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. + /// The R matrix obtained by calling . + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. /// Contains additional information on Q. Only used for the native solver /// and can be null for the managed provider. /// On entry the B matrix; on exit the X matrix. @@ -827,8 +920,10 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas /// The work array - only used in the native provider. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. + /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. - public override void QRSolveFactored(float[] q, float[] r, int rowsR, int columnsR, float[] tau, float[] b, int columnsB, float[] x, float[] work, QRMethod method = QRMethod.Full) + [SecuritySafeCritical] + public override void QRSolveFactored(float[] q, float[] r, int rowsA, int columnsA, float[] tau, float[] b, int columnsB, float[] x, float[] work, QRMethod method = QRMethod.Full) { if (r == null) { @@ -855,38 +950,54 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (r.Length != rowsR*columnsR) + int rowsQ, columnsQ, rowsR, columnsR; + if (method == QRMethod.Full) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + rowsQ = columnsQ = rowsR = rowsA; + columnsR = columnsA; + } + else + { + rowsQ = rowsA; + columnsQ = rowsR = columnsR = columnsA; } - if (q.Length != rowsR*rowsR) + if (r.Length != rowsR * columnsR) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsR * columnsR), "r"); } - if (b.Length != rowsR*columnsB) + if (q.Length != rowsQ * columnsQ) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsQ * columnsQ), "q"); } - if (x.Length != columnsR*columnsB) + if (b.Length != rowsA * columnsB) { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsA * columnsB), "b"); } - if (rowsR < columnsR) + if (x.Length != columnsA * columnsB) { - throw new ArgumentException(Resources.RowsLessThanColumns); + throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, columnsA * columnsB), "x"); } if (work.Length < 1) { - work[0] = rowsR*Control.BlockSize; + work[0] = rowsA * Control.BlockSize; throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - SafeNativeMethods.s_qr_solve_factored(rowsR, columnsR, columnsB, r, b, tau, x, work, work.Length); + if (method == QRMethod.Full) + { + SafeNativeMethods.s_qr_solve_factored(rowsA, columnsA, columnsB, r, b, tau, x, work, work.Length); + } + else + { + // we don't have access to the raw Q matrix any more(it is stored in R in the full QR), need to think about this. + // let just call the managed version in the meantime. The heavy lifting has already been done. -marcus + base.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, QRMethod.Thin); + } } /// @@ -925,12 +1036,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("vt"); } - if (u.Length != rowsA*rowsA) + if (u.Length != rowsA * rowsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); } - if (vt.Length != columnsA*columnsA) + if (vt.Length != columnsA * columnsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); } @@ -940,7 +1051,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); } - var work = new float[Math.Max(((3*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)), 5*Math.Min(rowsA, columnsA))]; + var work = new float[Math.Max(((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)), 5 * Math.Min(rowsA, columnsA))]; SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, work); } @@ -970,20 +1081,20 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("x"); } - if (b.Length != rowsA*columnsB) + if (b.Length != rowsA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - if (x.Length != columnsA*columnsB) + if (x.Length != columnsA * columnsB) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - var work = new float[Math.Max(((3*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)), 5*Math.Min(rowsA, columnsA))]; + var work = new float[Math.Max(((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)), 5 * Math.Min(rowsA, columnsA))]; var s = new float[Math.Min(rowsA, columnsA)]; - var u = new float[rowsA*rowsA]; - var vt = new float[columnsA*columnsA]; + var u = new float[rowsA * rowsA]; + var vt = new float[columnsA * columnsA]; var clone = new float[a.Length]; a.Copy(clone); @@ -1035,12 +1146,12 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentNullException("work"); } - if (u.Length != rowsA*rowsA) + if (u.Length != rowsA * rowsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); } - if (vt.Length != columnsA*columnsA) + if (vt.Length != columnsA * columnsA) { throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); } @@ -1055,13 +1166,73 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.GotoBlas throw new ArgumentException(Resources.ArgumentSingleDimensionArray, "work"); } - if (work.Length < Math.Max(((3*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)), 5*Math.Min(rowsA, columnsA))) + if (work.Length < Math.Max(((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)), 5 * Math.Min(rowsA, columnsA))) { - work[0] = Math.Max(((3*Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)), 5*Math.Min(rowsA, columnsA)); + work[0] = Math.Max(((3 * Math.Min(rowsA, columnsA)) + Math.Max(rowsA, columnsA)), 5 * Math.Min(rowsA, columnsA)); throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); } - SafeNativeMethods.s_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length); + if (SafeNativeMethods.s_svd_factor(computeVectors, rowsA, columnsA, a, s, u, vt, work, work.Length) > 0) + { + throw new NonConvergenceException(); + } + } + + /// + /// Computes the eigenvalues and eigenvectors of a matrix. + /// + /// Whether the matrix is symmetric or not. + /// The order of the matrix. + /// The matrix to decompose. The lenth of the array must be order * order. + /// On output, the matrix contains the eigen vectors. The lenth of the array must be order * order. + /// On output, the eigen values (λ) of matrix in ascending value. The length of the arry must . + /// On output, the block diagonal eigenvalue matrix. The lenth of the array must be order * order. + public override void EigenDecomp(bool isSymmetric, int order, float[] matrix, float[] matrixEv, Complex[] vectorEv, float[] matrixD) + { + if (matrix == null) + { + throw new ArgumentNullException("matrix"); + } + + if (matrix.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrix"); + } + + if (matrixEv == null) + { + throw new ArgumentNullException("matrixEv"); + } + + if (matrixEv.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixEv"); + } + + if (vectorEv == null) + { + throw new ArgumentNullException("vectorEv"); + } + + if (vectorEv.Length != order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order), "vectorEv"); + } + + if (matrixD == null) + { + throw new ArgumentNullException("matrixD"); + } + + if (matrixD.Length != order * order) + { + throw new ArgumentException(String.Format(Resources.ArgumentArrayWrongLength, order * order), "matrixD"); + } + + if (SafeNativeMethods.s_eigen(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD) > 0) + { + throw new NonConvergenceException(); + } } } } diff --git a/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.cs new file mode 100644 index 00000000..6a70a194 --- /dev/null +++ b/src/Numerics/Providers/LinearAlgebra/OpenBlas/OpenBlasLinearAlgebraProvider.cs @@ -0,0 +1,107 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2015 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +#if NATIVE + +using MathNet.Numerics.Properties; +using System; +using System.Numerics; +using System.Security; + +namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas +{ + public enum ParallelType : int + { + Sequential = 0, + Thread = 1, + OpenMP = 2 + } + + /// + /// OpenBLAS linear algebra provider. + /// + public partial class OpenBlasLinearAlgebraProvider : ManagedLinearAlgebraProvider + { + bool _nativeIX86; + bool _nativeX64; + bool _nativeIA64; + bool _nativeARM; + + public OpenBlasLinearAlgebraProvider() + { + + } + + public override void InitializeVerify() + { + int linearAlgebra; + try + { + // Load the native library + NativeProviderLoader.TryLoad(SafeNativeMethods.DllName); + + _nativeIX86 = SafeNativeMethods.query_capability((int)ProviderPlatform.x86) > 0; + _nativeX64 = SafeNativeMethods.query_capability((int)ProviderPlatform.x64) > 0; + _nativeIA64 = SafeNativeMethods.query_capability((int)ProviderPlatform.ia64) > 0; + _nativeARM = SafeNativeMethods.query_capability((int)ProviderPlatform.arm) > 0; + + linearAlgebra = SafeNativeMethods.query_capability((int)ProviderCapability.LinearAlgebra); + } + catch (DllNotFoundException e) + { + throw new NotSupportedException("OpenBLAS Native Provider not found.", e); + } + catch (BadImageFormatException e) + { + throw new NotSupportedException("OpenBLAS Native Provider found but failed to load. Please verify that the platform matches (x64 vs x32, Windows vs Linux).", e); + } + catch (EntryPointNotFoundException e) + { + throw new NotSupportedException("OpenBLAS Native Provider does not support capability querying and is therefore not compatible. Consider upgrading to a newer version.", e); + } + + // set threading settings, if supported + if (SafeNativeMethods.query_capability((int)ProviderConfig.Threading) > 0) + { + SafeNativeMethods.set_max_threads(Control.MaxDegreeOfParallelism); + } + } + + public override string ToString() + { + return string.Format("OpenBLAS\r\nProvider revision: {0}\r\nCPU core name: {1}\r\nLibrary config: {1})", + SafeNativeMethods.query_capability((int)ProviderConfig.Revision), + SafeNativeMethods.get_cpu_core(), + SafeNativeMethods.get_build_config()); + } + } +} + +#endif diff --git a/src/Numerics/Providers/LinearAlgebra/OpenBlas/SafeNativeMethods.cs b/src/Numerics/Providers/LinearAlgebra/OpenBlas/SafeNativeMethods.cs new file mode 100644 index 00000000..fa59a5d7 --- /dev/null +++ b/src/Numerics/Providers/LinearAlgebra/OpenBlas/SafeNativeMethods.cs @@ -0,0 +1,303 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009-2010 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +#if NATIVE + +using System.Numerics; +using System.Runtime.InteropServices; +using System.Security; + +namespace MathNet.Numerics.Providers.LinearAlgebra.OpenBlas +{ + /// + /// P/Invoke methods to the native math libraries. + /// + [SuppressUnmanagedCodeSecurity] + [SecurityCritical] + internal static class SafeNativeMethods + { + /// + /// Name of the native DLL. + /// + const string _DllName = "MathNET.Numerics.OpenBLAS.dll"; + internal static string DllName { get { return _DllName; } } + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int query_capability(int capability); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern string get_build_config(); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern string get_cpu_core(); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern ParallelType get_parallel_type(); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void set_max_threads(int num_threads); + + #region BLAS + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void s_axpy(int n, float alpha, float[] x, [In, Out] float[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void d_axpy(int n, double alpha, double[] x, [In, Out] double[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void c_axpy(int n, Complex32 alpha, Complex32[] x, [In, Out] Complex32[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void z_axpy(int n, Complex alpha, Complex[] x, [In, Out] Complex[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void s_scale(int n, float alpha, [Out] float[] x); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void d_scale(int n, double alpha, [Out] double[] x); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void c_scale(int n, Complex32 alpha, [In, Out] Complex32[] x); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void z_scale(int n, Complex alpha, [In, Out] Complex[] x); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern float s_dot_product(int n, float[] x, float[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern double d_dot_product(int n, double[] x, double[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern Complex32 c_dot_product(int n, Complex32[] x, Complex32[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern Complex z_dot_product(int n, Complex[] x, Complex[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void s_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, float alpha, float[] x, float[] y, float beta, [In, Out] float[] c); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void d_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, double alpha, double[] x, double[] y, double beta, [In, Out] double[] c); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void c_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, Complex32 alpha, Complex32[] x, Complex32[] y, Complex32 beta, [In, Out] Complex32[] c); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void z_matrix_multiply(Transpose transA, Transpose transB, int m, int n, int k, Complex alpha, Complex[] x, Complex[] y, Complex beta, [In, Out] Complex[] c); + + #endregion BLAS + + #region LAPACK + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern float s_matrix_norm(byte norm, int rows, int columns, [In] float[] a, [In, Out] float[] work); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern float d_matrix_norm(byte norm, int rows, int columns, [In] double[] a, [In, Out] double[] work); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern float c_matrix_norm(byte norm, int rows, int columns, [In] Complex32[] a, [In, Out] float[] work); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern double z_matrix_norm(byte norm, int rows, int columns, [In] Complex[] a, [In, Out] double[] work); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_lu_factor(int n, [In, Out] float[] a, [In, Out] int[] ipiv); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_lu_factor(int n, [In, Out] double[] a, [In, Out] int[] ipiv); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_lu_factor(int n, [In, Out] Complex32[] a, [In, Out] int[] ipiv); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_lu_factor(int n, [In, Out] Complex[] a, [In, Out] int[] ipiv); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_lu_inverse(int n, [In, Out] float[] a, [In, Out] float[] work, int lwork); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_lu_inverse(int n, [In, Out] double[] a, [In, Out] double[] work, int lwork); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_lu_inverse(int n, [In, Out] Complex32[] a, [In, Out] Complex32[] work, int lwork); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_lu_inverse(int n, [In, Out] Complex[] a, [In, Out] Complex[] work, int lwork); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_lu_inverse_factored(int n, [In, Out] float[] a, [In, Out] int[] ipiv, [In, Out] float[] work, int lwork); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_lu_inverse_factored(int n, [In, Out] double[] a, [In, Out] int[] ipiv, [In, Out] double[] work, int lwork); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_lu_inverse_factored(int n, [In, Out] Complex32[] a, [In, Out] int[] ipiv, [In, Out] Complex32[] work, int lwork); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_lu_inverse_factored(int n, [In, Out] Complex[] a, [In, Out] int[] ipiv, [In, Out] Complex[] work, int lwork); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_lu_solve_factored(int n, int nrhs, float[] a, [In, Out] int[] ipiv, [In, Out] float[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_lu_solve_factored(int n, int nrhs, double[] a, [In, Out] int[] ipiv, [In, Out] double[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_lu_solve_factored(int n, int nrhs, Complex32[] a, [In, Out] int[] ipiv, [In, Out] Complex32[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_lu_solve_factored(int n, int nrhs, Complex[] a, [In, Out] int[] ipiv, [In, Out] Complex[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_lu_solve(int n, int nrhs, float[] a, [In, Out] float[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_lu_solve(int n, int nrhs, double[] a, [In, Out] double[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_lu_solve(int n, int nrhs, Complex32[] a, [In, Out] Complex32[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_lu_solve(int n, int nrhs, Complex[] a, [In, Out] Complex[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_cholesky_factor(int n, [In, Out] float[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_cholesky_factor(int n, [In, Out] double[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_cholesky_factor(int n, [In, Out] Complex32[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_cholesky_factor(int n, [In, Out] Complex[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_cholesky_solve(int n, int nrhs, float[] a, [In, Out] float[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_cholesky_solve(int n, int nrhs, double[] a, [In, Out] double[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_cholesky_solve(int n, int nrhs, Complex32[] a, [In, Out] Complex32[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_cholesky_solve(int n, int nrhs, Complex[] a, [In, Out] Complex[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_cholesky_solve_factored(int n, int nrhs, float[] a, [In, Out] float[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_cholesky_solve_factored(int n, int nrhs, double[] a, [In, Out] double[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_cholesky_solve_factored(int n, int nrhs, Complex32[] a, [In, Out] Complex32[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_cholesky_solve_factored(int n, int nrhs, Complex[] a, [In, Out] Complex[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_qr_factor(int m, int n, [In, Out] float[] r, [In, Out] float[] tau, [In, Out] float[] q, [In, Out] float[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_qr_thin_factor(int m, int n, [In, Out] float[] q, [In, Out] float[] tau, [In, Out] float[] r, [In, Out] float[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_qr_factor(int m, int n, [In, Out] double[] r, [In, Out] double[] tau, [In, Out] double[] q, [In, Out] double[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_qr_thin_factor(int m, int n, [In, Out] double[] q, [In, Out] double[] tau, [In, Out] double[] r, [In, Out] double[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_qr_factor(int m, int n, [In, Out] Complex32[] r, [In, Out] Complex32[] tau, [In, Out] Complex32[] q, [In, Out] Complex32[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_qr_thin_factor(int m, int n, [In, Out] Complex32[] q, [In, Out] Complex32[] tau, [In, Out] Complex32[] r, [In, Out] Complex32[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_qr_factor(int m, int n, [In, Out] Complex[] r, [In, Out] Complex[] tau, [In, Out] Complex[] q, [In, Out] Complex[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_qr_thin_factor(int m, int n, [In, Out] Complex[] q, [In, Out] Complex[] tau, [In, Out] Complex[] r, [In, Out] Complex[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_qr_solve(int m, int n, int bn, float[] r, float[] b, [In, Out] float[] x, [In, Out] float[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_qr_solve(int m, int n, int bn, double[] r, double[] b, [In, Out] double[] x, [In, Out] double[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_qr_solve(int m, int n, int bn, Complex32[] r, Complex32[] b, [In, Out] Complex32[] x, [In, Out] Complex32[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_qr_solve(int m, int n, int bn, Complex[] r, Complex[] b, [In, Out] Complex[] x, [In, Out] Complex[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_qr_solve_factored(int m, int n, int bn, float[] r, float[] b, float[] tau, [In, Out] float[] x, [In, Out] float[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_qr_solve_factored(int m, int n, int bn, double[] r, double[] b, double[] tau, [In, Out] double[] x, [In, Out] double[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_qr_solve_factored(int m, int n, int bn, Complex32[] r, Complex32[] b, Complex32[] tau, [In, Out] Complex32[] x, [In, Out] Complex32[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_qr_solve_factored(int m, int n, int bn, Complex[] r, Complex[] b, Complex[] tau, [In, Out] Complex[] x, [In, Out] Complex[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_svd_factor(bool computeVectors, int m, int n, [In, Out] float[] a, [In, Out] float[] s, [In, Out] float[] u, [In, Out] float[] v, [In, Out] float[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_svd_factor(bool computeVectors, int m, int n, [In, Out] double[] a, [In, Out] double[] s, [In, Out] double[] u, [In, Out] double[] v, [In, Out] double[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_svd_factor(bool computeVectors, int m, int n, [In, Out] Complex32[] a, [In, Out] Complex32[] s, [In, Out] Complex32[] u, [In, Out] Complex32[] v, [In, Out] Complex32[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_svd_factor(bool computeVectors, int m, int n, [In, Out] Complex[] a, [In, Out] Complex[] s, [In, Out] Complex[] u, [In, Out] Complex[] v, [In, Out] Complex[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_eigen(bool isSymmetric, int n, [In, Out] float[] a, [In, Out] float[] vectors, [In, Out] Complex[] values, [In, Out] float[] d); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_eigen(bool isSymmetric, int n, [In, Out] double[] a, [In, Out] double[] vectors, [In, Out] Complex[] values, [In, Out] double[] d); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_eigen(bool isSymmetric, int n, [In, Out] Complex32[] a, [In, Out] Complex32[] vectors, [In, Out] Complex[] values, [In, Out] Complex32[] d); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_eigen(bool isSymmetric, int n, [In, Out] Complex[] a, [In, Out] Complex[] vectors, [In, Out] Complex[] values, [In, Out] Complex[] d); + + #endregion LAPACK + } +} + +#endif diff --git a/src/Numerics/Providers/LinearAlgebra/ProviderCapabilities.cs b/src/Numerics/Providers/LinearAlgebra/ProviderCapabilities.cs new file mode 100644 index 00000000..b07d4e0f --- /dev/null +++ b/src/Numerics/Providers/LinearAlgebra/ProviderCapabilities.cs @@ -0,0 +1,55 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2015 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +namespace MathNet.Numerics.Providers.LinearAlgebra +{ + public enum ProviderPlatform : int + { + x86 = 8, + x64 = 9, + ia64 = 10, + arm = 11, + } + + public enum ProviderConfig : int + { + Revision = 64, + Precision = 65, + Threading = 66, + Memory = 67, + } + + public enum ProviderCapability : int + { + LinearAlgebra = 128, + Optimization = 256, + FFT = 384, + } +} \ No newline at end of file diff --git a/src/UnitTests/UnitTests-MKL.csproj b/src/UnitTests/UnitTests-MKL.csproj index 07bd43ec..df1bc55f 100644 --- a/src/UnitTests/UnitTests-MKL.csproj +++ b/src/UnitTests/UnitTests-MKL.csproj @@ -345,4 +345,7 @@ True + + + \ No newline at end of file diff --git a/src/UnitTests/UnitTests-OpenBLAS.csproj b/src/UnitTests/UnitTests-OpenBLAS.csproj new file mode 100644 index 00000000..94b0ac3e --- /dev/null +++ b/src/UnitTests/UnitTests-OpenBLAS.csproj @@ -0,0 +1,351 @@ + + + + 10.0 + Debug + AnyCPU + 8.0.30703 + 2.0 + {96B903EF-3EE1-4569-803C-0482D2F5ED37} + Library + Properties + MathNet.Numerics.UnitTests + MathNet.Numerics.UnitTestsOpenBLAS + v4.5 + 512 + ..\..\ + + + TRACE;NATIVE;OPENBLAS; + ..\..\out\OpenBLAS\Windows\ + ..\..\obj\OpenBLAS\Windows\AnyCPU\ + ..\..\obj\OpenBLAS\Windows\AnyCPU\ + true + pdbonly + prompt + MinimumRecommendedRules.ruleset + 1591 + AnyCPU + + + TRACE;DEBUG;NATIVE;OPENBLAS; + ..\..\out\OpenBLAS\Windows\ + ..\..\obj\OpenBLAS\Windows\AnyCPU\ + ..\..\obj\OpenBLAS\Windows\AnyCPU\ + false + full + true + prompt + 4 + 1591 + AnyCPU + + + + + + + + + + + + + + + + + + + data\Codeplex-5667.csv + Always + + + data\Github-Cureos-1.csv + Always + + + data\Matlab\A.mat + Always + + + data\Matlab\collection-nocompress.mat + Always + + + data\Matlab\collection.mat + Always + + + data\Matlab\complex.mat + Always + + + data\Matlab\sparse-large.mat + Always + + + data\Matlab\sparse-small.mat + Always + + + data\Matlab\sparse_complex.mat + Always + + + data\Matlab\v.mat + Always + + + data\NIST\AtmWtAgt.dat + Always + + + data\NIST\Bennett5.dat + Always + + + data\NIST\BoxBOD.dat + Always + + + data\NIST\Chwirut1.dat + Always + + + data\NIST\Chwirut2.dat + Always + + + data\NIST\DanWood.dat + Always + + + data\NIST\Eckerle4.dat + Always + + + data\NIST\ENSO.dat + Always + + + data\NIST\Filip.dat + Always + + + data\NIST\Gauss1.dat + Always + + + data\NIST\Gauss2.dat + Always + + + data\NIST\Gauss3.dat + Always + + + data\NIST\Hahn1.dat + Always + + + data\NIST\Kirby2.dat + Always + + + data\NIST\Lanczos1.dat + Always + + + data\NIST\Lanczos2.dat + Always + + + data\NIST\Lanczos3.dat + Always + + + data\NIST\Lew.dat + Always + + + data\NIST\Longley.dat + Always + + + data\NIST\Lottery.dat + Always + + + data\NIST\Mavro.dat + Always + + + data\NIST\MGH09.dat + Always + + + data\NIST\MGH10.dat + Always + + + data\NIST\MGH17.dat + Always + + + data\NIST\Michelso.dat + Always + + + data\NIST\Misra1a.dat + Always + + + data\NIST\Misra1b.dat + Always + + + data\NIST\Misra1c.dat + Always + + + data\NIST\Misra1d.dat + Always + + + data\NIST\Nelson.dat + Always + + + data\NIST\NoInt1.dat + Always + + + data\NIST\NoInt2.dat + Always + + + data\NIST\Norris.dat + Always + + + data\NIST\NumAcc1.dat + Always + + + data\NIST\NumAcc2.dat + Always + + + data\NIST\NumAcc3.dat + Always + + + data\NIST\NumAcc4.dat + Always + + + data\NIST\Pontius.dat + Always + + + data\NIST\Rat42.dat + Always + + + data\NIST\Rat43.dat + Always + + + data\NIST\Roszman1.dat + Always + + + data\NIST\SiRstvt.dat + Always + + + data\NIST\SmLs01t.dat + Always + + + data\NIST\SmLs02t.dat + Always + + + data\NIST\SmLs03t.dat + Always + + + data\NIST\SmLs04t.dat + Always + + + data\NIST\SmLs05t.dat + Always + + + data\NIST\SmLs06t.dat + Always + + + data\NIST\SmLs07t.dat + Always + + + data\NIST\SmLs08t.dat + Always + + + data\NIST\SmLs09t.dat + Always + + + data\NIST\Thurber.dat + Always + + + data\NIST\Wampler1.dat + Always + + + data\NIST\Wampler2.dat + Always + + + data\NIST\Wampler3.dat + Always + + + data\NIST\Wampler4.dat + Always + + + data\NIST\Wampler5.dat + Always + + + + data\NIST\Meixner.dat + Always + + + + + + {b7cae5f4-a23f-4438-b5be-41226618b695} + Numerics + + + + + + ..\..\packages\NUnit\lib\nunit.framework.dll + True + True + + + + + + \ No newline at end of file diff --git a/src/UnitTests/UseLinearAlgebraProvider.cs b/src/UnitTests/UseLinearAlgebraProvider.cs index 00266687..13975002 100644 --- a/src/UnitTests/UseLinearAlgebraProvider.cs +++ b/src/UnitTests/UseLinearAlgebraProvider.cs @@ -43,6 +43,8 @@ namespace MathNet.Numerics.UnitTests Control.UseNativeMKL(); #elif CUDA Control.UseNativeCUDA(); +#elif OPENBLAS + Control.UseNativeOpenBLAS(); #endif #endif }