From 8b7b617481103845ee7f29206c373884c435a503 Mon Sep 17 00:00:00 2001 From: Matthew Johnson Date: Sat, 25 Apr 2015 02:42:37 +0100 Subject: [PATCH] Some things are working, but most aren't. There appears to be something wrong with my usage of cublas?getriBatched, which isn't entirely surprising. I'll need to figure that out next. --- MathNet.Numerics.NativeProviders.sln | 44 +- src/NativeProviders/CUDA/blas.cpp | 6 +- src/NativeProviders/CUDA/capabilities.cpp | 79 ++ src/NativeProviders/CUDA/lapack.cpp | 125 +++- .../Windows/CUDA/CUDAWrapper.vcxproj | 18 +- .../Windows/CUDA/CUDAWrapper.vcxproj.filters | 3 + src/Numerics/Control.cs | 9 + src/Numerics/Numerics.csproj | 6 + src/Numerics/Properties/AssemblyInfo.cs | 1 + src/Numerics/Properties/Resources.Designer.cs | 21 +- src/Numerics/Properties/Resources.resx | 3 + .../Cuda/CudaLinearAlgebraProvider.Complex.cs | 703 +++++++++++++++++ .../CudaLinearAlgebraProvider.Complex32.cs | 703 +++++++++++++++++ .../Cuda/CudaLinearAlgebraProvider.Double.cs | 704 ++++++++++++++++++ .../Cuda/CudaLinearAlgebraProvider.Single.cs | 703 +++++++++++++++++ .../Cuda/CudaLinearAlgebraProvider.cs | 209 ++++++ .../LinearAlgebra/Cuda/SafeNativeMethods.cs | 378 ++++++++++ src/UnitTests/UnitTests-CUDA.csproj | 348 +++++++++ src/UnitTests/UseLinearAlgebraProvider.cs | 3 + 19 files changed, 4020 insertions(+), 46 deletions(-) create mode 100644 src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex.cs create mode 100644 src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex32.cs create mode 100644 src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Double.cs create mode 100644 src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Single.cs create mode 100644 src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.cs create mode 100644 src/Numerics/Providers/LinearAlgebra/Cuda/SafeNativeMethods.cs create mode 100644 src/UnitTests/UnitTests-CUDA.csproj diff --git a/MathNet.Numerics.NativeProviders.sln b/MathNet.Numerics.NativeProviders.sln index 1c7bf18d..0297e6f8 100644 --- a/MathNet.Numerics.NativeProviders.sln +++ b/MathNet.Numerics.NativeProviders.sln @@ -22,6 +22,10 @@ Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "UnitTests-MKL", "src\UnitTe EndProject Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "CUDA", "src\NativeProviders\Windows\CUDA\CUDAWrapper.vcxproj", "{5A52B796-7F41-4C90-8DE2-F3F391C4482C}" EndProject +Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "UnitTests-CUDA", "src\UnitTests\UnitTests-CUDA.csproj", "{E79C0395-01DC-4BC9-B86C-ED45790892C5}" +EndProject +Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "Scratch", "Scratch\Scratch.csproj", "{2386FAD1-BB99-4597-885C-8EF81D0637BA}" +EndProject Global GlobalSection(SolutionConfigurationPlatforms) = preSolution Debug|Any CPU = Debug|Any CPU @@ -125,13 +129,51 @@ Global {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release|Mixed Platforms.Build.0 = Release|Win32 {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release|Win32.ActiveCfg = Release|Win32 {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release|Win32.Build.0 = Release|Win32 - {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release|x64.ActiveCfg = Release|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release|x64.ActiveCfg = Release|x64 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release|x64.Build.0 = 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 {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-Signed|Win32.ActiveCfg = Release|Win32 {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-Signed|Win32.Build.0 = Release|Win32 {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-Signed|x64.ActiveCfg = Release|Win32 + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Debug|Any CPU.ActiveCfg = Debug|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Debug|Any CPU.Build.0 = Debug|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Debug|Mixed Platforms.ActiveCfg = Debug|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Debug|Mixed Platforms.Build.0 = Debug|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Debug|Win32.ActiveCfg = Debug|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Debug|x64.ActiveCfg = Debug|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release|Any CPU.ActiveCfg = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release|Any CPU.Build.0 = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release|Mixed Platforms.ActiveCfg = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release|Mixed Platforms.Build.0 = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release|Win32.ActiveCfg = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release|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 + {2386FAD1-BB99-4597-885C-8EF81D0637BA}.Debug|Any CPU.ActiveCfg = Debug|Any CPU + {2386FAD1-BB99-4597-885C-8EF81D0637BA}.Debug|Any CPU.Build.0 = Debug|Any CPU + {2386FAD1-BB99-4597-885C-8EF81D0637BA}.Debug|Mixed Platforms.ActiveCfg = Debug|Any CPU + {2386FAD1-BB99-4597-885C-8EF81D0637BA}.Debug|Mixed Platforms.Build.0 = Debug|Any CPU + {2386FAD1-BB99-4597-885C-8EF81D0637BA}.Debug|Win32.ActiveCfg = Debug|Any CPU + {2386FAD1-BB99-4597-885C-8EF81D0637BA}.Debug|x64.ActiveCfg = Debug|Any CPU + {2386FAD1-BB99-4597-885C-8EF81D0637BA}.Debug|x64.Build.0 = Debug|Any CPU + {2386FAD1-BB99-4597-885C-8EF81D0637BA}.Release|Any CPU.ActiveCfg = Release|Any CPU + {2386FAD1-BB99-4597-885C-8EF81D0637BA}.Release|Any CPU.Build.0 = Release|Any CPU + {2386FAD1-BB99-4597-885C-8EF81D0637BA}.Release|Mixed Platforms.ActiveCfg = Release|Any CPU + {2386FAD1-BB99-4597-885C-8EF81D0637BA}.Release|Mixed Platforms.Build.0 = Release|Any CPU + {2386FAD1-BB99-4597-885C-8EF81D0637BA}.Release|Win32.ActiveCfg = Release|Any CPU + {2386FAD1-BB99-4597-885C-8EF81D0637BA}.Release|x64.ActiveCfg = Release|Any CPU + {2386FAD1-BB99-4597-885C-8EF81D0637BA}.Release-Signed|Any CPU.ActiveCfg = Release|Any CPU + {2386FAD1-BB99-4597-885C-8EF81D0637BA}.Release-Signed|Any CPU.Build.0 = Release|Any CPU + {2386FAD1-BB99-4597-885C-8EF81D0637BA}.Release-Signed|Mixed Platforms.ActiveCfg = Release|Any CPU + {2386FAD1-BB99-4597-885C-8EF81D0637BA}.Release-Signed|Mixed Platforms.Build.0 = Release|Any CPU + {2386FAD1-BB99-4597-885C-8EF81D0637BA}.Release-Signed|Win32.ActiveCfg = Release|Any CPU + {2386FAD1-BB99-4597-885C-8EF81D0637BA}.Release-Signed|x64.ActiveCfg = Release|Any CPU EndGlobalSection GlobalSection(SolutionProperties) = preSolution HideSolutionNode = FALSE diff --git a/src/NativeProviders/CUDA/blas.cpp b/src/NativeProviders/CUDA/blas.cpp index 47233b72..c8348a32 100644 --- a/src/NativeProviders/CUDA/blas.cpp +++ b/src/NativeProviders/CUDA/blas.cpp @@ -75,9 +75,8 @@ void cuda_gemm(const cublasHandle_t handle, const cublasOperation_t transa, cons cudaFree(d_C); } -#if GCC extern "C" { -#endif + DLLEXPORT void s_axpy(const cublasHandle_t blasHandle, const int n, const float alpha, const float x[], float y[]){ cuda_axpy(blasHandle, n, &alpha, x, 1, y, 1, cublasSaxpy); } @@ -162,6 +161,5 @@ extern "C" { cuda_gemm(blasHandle, transA, transB, m, n, k, &alpha, x, lda, y, ldb, &beta, c, m, cublasZgemm); } -#if GCC } -#endif + diff --git a/src/NativeProviders/CUDA/capabilities.cpp b/src/NativeProviders/CUDA/capabilities.cpp index e69de29b..75fb1d59 100644 --- a/src/NativeProviders/CUDA/capabilities.cpp +++ b/src/NativeProviders/CUDA/capabilities.cpp @@ -0,0 +1,79 @@ +#include "wrapper_common.h" +#include "cublas_v2.h" +#include "cusolverDn.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 1; // revision + + // LINEAR ALGEBRA + case 128: return 1; // basic dense linear algebra + + // OPTIMIZATION + case 256: return 0; // basic optimization + + // FFT + case 384: return 0; // basic FFT + + default: return 0; // unknown or not supported + + } + } + + DLLEXPORT cublasStatus_t createBLASHandle(cublasHandle_t *blasHandle){ + return cublasCreate(blasHandle); + } + + DLLEXPORT cublasStatus_t destroyBLASHandle(cublasHandle_t blasHandle){ + return cublasDestroy(blasHandle); + } + + DLLEXPORT cusolverStatus_t createSolverHandle(cusolverDnHandle_t *solverHandle){ + return cusolverDnCreate(solverHandle); + } + + DLLEXPORT cusolverStatus_t destroySolverHandle(cusolverDnHandle_t solverHandle){ + return cusolverDnDestroy(solverHandle); + } + +#ifdef __cplusplus +} +#endif /* __cplusplus */ diff --git a/src/NativeProviders/CUDA/lapack.cpp b/src/NativeProviders/CUDA/lapack.cpp index 7798df12..718c4909 100644 --- a/src/NativeProviders/CUDA/lapack.cpp +++ b/src/NativeProviders/CUDA/lapack.cpp @@ -10,8 +10,6 @@ template inline int lu_factor(cusolverDnHandle_t solverHandle, int m, T a[], int ipiv[], GETRF getrf, GETRFBSIZE getrfbsize) { int info = 0; - T* work = NULL; - int lwork = 0; T* d_A = NULL; cudaMalloc((void**)&d_A, m*m*sizeof(T)); @@ -20,10 +18,17 @@ inline int lu_factor(cusolverDnHandle_t solverHandle, int m, T a[], int ipiv[], int* d_I = NULL; cudaMalloc((void**)&d_I, m*sizeof(int)); + T* work = NULL; + int lwork = 0; getrfbsize(solverHandle, m, m, a, m, &lwork); - cudaMalloc((void**)lwork, sizeof(T)*lwork); + cudaMalloc((void**)&work, sizeof(T)*lwork); + + int* d_info = NULL; + cudaMalloc((void**)&d_info, sizeof(int)); - getrf(solverHandle, m, m, d_A, m, work, d_I, &info); + getrf(solverHandle, m, m, d_A, m, work, d_I, d_info); + + cudaMemcpy(&info, d_info, 1, cudaMemcpyDeviceToHost); cublasGetMatrix(m, m, sizeof(T), d_A, m, a, m); cublasGetVector(m, sizeof(T), d_I, 1, ipiv, 1); @@ -32,6 +37,7 @@ inline int lu_factor(cusolverDnHandle_t solverHandle, int m, T a[], int ipiv[], cudaFree(d_A); cudaFree(d_I); + cudaFree(d_info); cudaFree(work); return info; @@ -41,8 +47,6 @@ template inline int lu_inverse(cusolverDnHandle_t solverHandle, cublasHandle_t blasHandle, int n, T a[], GETRF getrf, GETRI getri, GETRFBSIZE getrfbsize) { int info = 0; - T* work = NULL; - int lwork = 0; int* d_I = NULL; cudaMalloc((void**)&d_I, n*sizeof(T)); @@ -51,30 +55,48 @@ inline int lu_inverse(cusolverDnHandle_t solverHandle, cublasHandle_t blasHandle cudaMalloc((void**)&d_A, n*n*sizeof(T)); cublasSetMatrix(n, n, sizeof(T), a, n, d_A, n); + T* work = NULL; + int lwork = 0; getrfbsize(solverHandle, n, n, d_A, n, &lwork); - cudaMalloc((void**)lwork, sizeof(T)*lwork); + cudaMalloc((void**)&work, sizeof(T)*lwork); + + int* d_info = NULL; + cudaMalloc((void**)&d_info, sizeof(int)); + + printf("initial %f %f %f %f %f %f %f %f %f\r\n", a[0], a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8]); - getrf(solverHandle, n, n, d_A, n, work, d_I, &info); + getrf(solverHandle, n, n, d_A, n, work, d_I, d_info); + cudaMemcpy(&info, d_info, 1, cudaMemcpyDeviceToHost); + cublasGetMatrix(n, n, sizeof(T), d_A, n, a, n); + printf("after factor %f %f %f %f %f %f %f %f %f\r\n", a[0], a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8]); + cudaFree(work); if (info != 0) { cudaFree(d_A); cudaFree(d_I); + cudaFree(d_info); return info; } T* d_C = NULL; cudaMalloc((void**)&d_C, n*n*sizeof(T)); - getri(blasHandle, n, d_A, n, d_I, d_C, n, &info); + getri(blasHandle, n, d_A, n, d_I, d_C, n, d_info); + cudaMemcpy(&info, d_info, 1, cudaMemcpyDeviceToHost); + + cublasGetMatrix(n, n, sizeof(T), d_A, n, a, n); + printf("a inverse %f %f %f %f %f %f %f %f %f\r\n", a[0], a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8]); cublasGetMatrix(n, n, sizeof(T), d_C, n, a, n); + printf("c inverse %f %f %f %f %f %f %f %f %f\r\n", a[0], a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8]); cudaFree(d_A); cudaFree(d_I); cudaFree(d_C); + cudaFree(d_info); return info; }; @@ -82,9 +104,10 @@ inline int lu_inverse(cusolverDnHandle_t solverHandle, cublasHandle_t blasHandle template inline int lu_inverse_factored(cublasHandle_t blasHandle, int n, T a[], int ipiv[], GETRI getri) { - shift_ipiv_up(n, ipiv); int info = 0; + shift_ipiv_up(n, ipiv); + T* d_A = NULL; cudaMalloc((void**)&d_A, n*n*sizeof(T)); cublasSetMatrix(n, n, sizeof(T), a, n, d_A, n); @@ -96,7 +119,11 @@ inline int lu_inverse_factored(cublasHandle_t blasHandle, int n, T a[], int ipiv cudaMalloc((void**)&d_I, n*sizeof(int)); cublasSetVector(n, sizeof(int), ipiv, 1, d_I, 1); - getri(blasHandle, n, d_A, n, d_I, d_C, n, &info); + int* d_info = NULL; + cudaMalloc((void**)&d_info, sizeof(int)); + + getri(blasHandle, n, d_A, n, d_I, d_C, n, d_info); + cudaMemcpy(&info, d_info, 1, cudaMemcpyDeviceToHost); cublasGetMatrix(n, n, sizeof(T), d_C, n, a, n); cublasGetVector(n, sizeof(int), d_I, 1, ipiv, 1); @@ -106,6 +133,7 @@ inline int lu_inverse_factored(cublasHandle_t blasHandle, int n, T a[], int ipiv cudaFree(d_A); cudaFree(d_I); cudaFree(d_C); + cudaFree(d_info); return info; } @@ -113,9 +141,10 @@ inline int lu_inverse_factored(cublasHandle_t blasHandle, int n, T a[], int ipiv template inline int lu_solve_factored(cusolverDnHandle_t solverHandle, int n, int nrhs, T a[], int ipiv[], T b[], GETRS getrs) { - shift_ipiv_up(n, ipiv); int info = 0; + shift_ipiv_up(n, ipiv); + T* d_A = NULL; cudaMalloc((void**)&d_A, n*n*sizeof(T)); cublasSetMatrix(n, n, sizeof(T), a, n, d_A, n); @@ -128,7 +157,11 @@ inline int lu_solve_factored(cusolverDnHandle_t solverHandle, int n, int nrhs, T cudaMalloc((void**)&d_I, n*sizeof(int)); cublasSetVector(n, sizeof(int), ipiv, 1, d_I, 1); - getrs(solverHandle, CUBLAS_OP_N, n, nrhs, d_A, n, d_I, d_B, n, &info); + int* d_info = NULL; + cudaMalloc((void**)&d_info, sizeof(int)); + + getrs(solverHandle, CUBLAS_OP_N, n, nrhs, d_A, n, d_I, d_B, n, d_info); + cudaMemcpy(&info, d_info, 1, cudaMemcpyDeviceToHost); cublasGetMatrix(n, nrhs, sizeof(T), d_B, n, b, n); @@ -137,6 +170,7 @@ inline int lu_solve_factored(cusolverDnHandle_t solverHandle, int n, int nrhs, T cudaFree(d_A); cudaFree(d_B); cudaFree(d_I); + cudaFree(d_info); return info; } @@ -145,8 +179,6 @@ template inline int lu_solve(cusolverDnHandle_t solverHandle, int n, int nrhs, T a[], T b[], GETRF getrf, GETRS getrs, GETRFBSIZE getrfbsize) { int info = 0; - T* work = NULL; - int lwork = 0; int* d_I = NULL; cudaMalloc((void**)&d_I, n*sizeof(T)); @@ -155,15 +187,22 @@ inline int lu_solve(cusolverDnHandle_t solverHandle, int n, int nrhs, T a[], T b cudaMalloc((void**)&d_A, n*n*sizeof(T)); cublasSetMatrix(n, n, sizeof(T), a, n, d_A, n); + T* work = NULL; + int lwork = 0; getrfbsize(solverHandle, n, n, a, n, &lwork); - cudaMalloc((void**)lwork, sizeof(T)*lwork); + cudaMalloc((void**)&work, sizeof(T)*lwork); + + int* d_info = NULL; + cudaMalloc((void**)&d_info, sizeof(int)); - getrf(solverHandle, n, n, d_A, n, work, d_I, &info); + getrf(solverHandle, n, n, d_A, n, work, d_I, d_info); + cudaMemcpy(&info, d_info, 1, cudaMemcpyDeviceToHost); if (info != 0) { cudaFree(d_I); cudaFree(d_A); + cudaFree(d_info); return info; } @@ -171,13 +210,15 @@ inline int lu_solve(cusolverDnHandle_t solverHandle, int n, int nrhs, T a[], T b cudaMalloc((void**)&d_B, n*nrhs*sizeof(T)); cublasSetMatrix(n, nrhs, sizeof(T), b, n, d_B, n); - getrs(solverHandle, CUBLAS_OP_N, n, nrhs, d_A, n, d_I, d_B, n, &info); + getrs(solverHandle, CUBLAS_OP_N, n, nrhs, d_A, n, d_I, d_B, n, d_info); + cudaMemcpy(&info, d_info, 1, cudaMemcpyDeviceToHost); cublasGetMatrix(n, nrhs, sizeof(T), d_B, n, b, n); cudaFree(d_A); cudaFree(d_B); cudaFree(d_I); + cudaFree(d_info); return info; } @@ -197,7 +238,11 @@ inline int cholesky_factor(cusolverDnHandle_t solverHandle, int n, T* a, POTRF p potrfbsize(solverHandle, CUBLAS_FILL_MODE_LOWER, n, d_A, n, &lWork); cudaMalloc((void**)&work, sizeof(T)*lWork); - potrf(solverHandle, CUBLAS_FILL_MODE_LOWER, n, d_A, n, work, lWork, &info); + int* d_info = NULL; + cudaMalloc((void**)&d_info, sizeof(int)); + + potrf(solverHandle, CUBLAS_FILL_MODE_LOWER, n, d_A, n, work, lWork, d_info); + cudaMemcpy(&info, d_info, 1, cudaMemcpyDeviceToHost); cublasGetMatrix(n, n, sizeof(T), d_A, n, a, n); @@ -214,6 +259,7 @@ inline int cholesky_factor(cusolverDnHandle_t solverHandle, int n, T* a, POTRF p } cudaFree(d_A); + cudaFree(d_info); cudaFree(work); return info; @@ -233,13 +279,18 @@ inline int cholesky_solve(cusolverDnHandle_t solverHandle, int n, int nrhs, T a[ potrfbsize(solverHandle, CUBLAS_FILL_MODE_LOWER, n, d_A, n, &lWork); cudaMalloc((void**)&work, sizeof(T)*lWork); - potrf(solverHandle, CUBLAS_FILL_MODE_LOWER, n, d_A, n, work, lWork, &info); + int* d_info = NULL; + cudaMalloc((void**)&d_info, sizeof(int)); + + potrf(solverHandle, CUBLAS_FILL_MODE_LOWER, n, d_A, n, work, lWork, d_info); + cudaMemcpy(&info, d_info, 1, cudaMemcpyDeviceToHost); cudaFree(work); if (info != 0) { cudaFree(d_A); + cudaFree(d_info); return info; } @@ -247,12 +298,14 @@ inline int cholesky_solve(cusolverDnHandle_t solverHandle, int n, int nrhs, T a[ cudaMalloc((void**)d_B, n*nrhs*sizeof(T)); cublasSetMatrix(n, nrhs, sizeof(T), b, n, d_B, n); - potrs(solverHandle, CUBLAS_FILL_MODE_LOWER, n, nrhs, d_A, n, d_B, n, &info); + potrs(solverHandle, CUBLAS_FILL_MODE_LOWER, n, nrhs, d_A, n, d_B, n, d_info); + cudaMemcpy(&info, d_info, 1, cudaMemcpyDeviceToHost); cublasGetMatrix(n, nrhs, sizeof(T), d_B, n, b, n); cudaFree(d_A); cudaFree(d_B); + cudaFree(d_info); return info; } @@ -270,12 +323,17 @@ inline int cholesky_solve_factored(cusolverDnHandle_t solverHandle, int n, int n cudaMalloc((void**)d_B, n*nrhs*sizeof(T)); cublasSetMatrix(n, nrhs, sizeof(T), b, n, d_B, n); - potrs(solverHandle, CUBLAS_FILL_MODE_LOWER, n, nrhs, d_A, n, d_B, n, &info); + int* d_info = NULL; + cudaMalloc((void**)&d_info, sizeof(int)); + + potrs(solverHandle, CUBLAS_FILL_MODE_LOWER, n, nrhs, d_A, n, d_B, n, d_info); + cudaMemcpy(&info, d_info, 1, cudaMemcpyDeviceToHost); cublasGetMatrix(n, nrhs, sizeof(T), d_B, n, b, n); cudaFree(d_A); cudaFree(d_B); + cudaFree(d_info); return info; } @@ -402,8 +460,13 @@ inline int svd_factor(cusolverDnHandle_t solverHandle, bool compute_vectors, int T* rwork = NULL; cudaMalloc((void**)&rwork, 5 * dim_s * sizeof(T)); + int* d_info = NULL; + cudaMalloc((void**)&d_info, sizeof(int)); + + char job = compute_vectors ? 'A' : 'N'; - gesvd(solverHandle, job, job, m, n, d_A, m, d_S, d_U, m, d_V, n, work, lWork, rwork, &info); + gesvd(solverHandle, job, job, m, n, d_A, m, d_S, d_U, m, d_V, n, work, lWork, rwork, d_info); + cudaMemcpy(&info, d_info, 1, cudaMemcpyDeviceToHost); cublasGetVector(dim_s, sizeof(T), d_S, 1, s, 1); cublasGetMatrix(m, m, sizeof(T), d_U, m, u, m); @@ -415,6 +478,7 @@ inline int svd_factor(cusolverDnHandle_t solverHandle, bool compute_vectors, int cudaFree(d_V); cudaFree(work); cudaFree(rwork); + cudaFree(d_info); return info; } @@ -447,8 +511,12 @@ inline int complex_svd_factor(cusolverDnHandle_t solverHandle, bool compute_vect R* rwork = NULL; cudaMalloc((void**)&rwork, 5 * dim_s * sizeof(R)); + int* d_info = NULL; + cudaMalloc((void**)&d_info, sizeof(int)); + char job = compute_vectors ? 'A' : 'N'; - gesvd(solverHandle, job, job, m, n, d_A, m, d_S, d_U, m, d_V, n, work, lWork, rwork, &info); + gesvd(solverHandle, job, job, m, n, d_A, m, d_S, d_U, m, d_V, n, work, lWork, rwork, d_info); + cudaMemcpy(&info, d_info, 1, cudaMemcpyDeviceToHost); cublasGetVector(dim_s, sizeof(T), d_S, 1, s_local, 1); cublasGetMatrix(m, m, sizeof(T), d_U, m, u, m); @@ -466,6 +534,7 @@ inline int complex_svd_factor(cusolverDnHandle_t solverHandle, bool compute_vect cudaFree(d_V); cudaFree(work); cudaFree(rwork); + cudaFree(d_info); return info; } @@ -643,22 +712,22 @@ inline int complex_svd_factor(cusolverDnHandle_t solverHandle, bool compute_vect #define zgesvdbsize cusolverDnZgesvd_bufferSize -inline int sgetri(cublasHandle_t handle, int n, const float a[], int lda, const int *ipiv, float c[], int ldc, int *info) +inline int sgetri(cublasHandle_t handle, int n, const float a[], int lda, const int ipiv[], float c[], int ldc, int *info) { return cublasSgetriBatched(handle, n, &a, lda, ipiv, &c, ldc, info, 1); } -inline int dgetri(cublasHandle_t handle, int n, const double a[], int lda, const int *ipiv, double c[], int ldc, int *info) +inline int dgetri(cublasHandle_t handle, int n, const double a[], int lda, const int ipiv[], double c[], int ldc, int *info) { return cublasDgetriBatched(handle, n, &a, lda, ipiv, &c, ldc, info, 1); } -inline int cgetri(cublasHandle_t handle, int n, const cuComplex a[], int lda, const int *ipiv, cuComplex c[], int ldc, int *info) +inline int cgetri(cublasHandle_t handle, int n, const cuComplex a[], int lda, const int ipiv[], cuComplex c[], int ldc, int *info) { return cublasCgetriBatched(handle, n, &a, lda, ipiv, &c, ldc, info, 1); } -inline int zgetri(cublasHandle_t handle, int n, const cuDoubleComplex a[], int lda, const int *ipiv, cuDoubleComplex c[], int ldc, int *info) +inline int zgetri(cublasHandle_t handle, int n, const cuDoubleComplex a[], int lda, const int ipiv[], cuDoubleComplex c[], int ldc, int *info) { return cublasZgetriBatched(handle, n, &a, lda, ipiv, &c, ldc, info, 1); } diff --git a/src/NativeProviders/Windows/CUDA/CUDAWrapper.vcxproj b/src/NativeProviders/Windows/CUDA/CUDAWrapper.vcxproj index a9786588..393b8c88 100644 --- a/src/NativeProviders/Windows/CUDA/CUDAWrapper.vcxproj +++ b/src/NativeProviders/Windows/CUDA/CUDAWrapper.vcxproj @@ -24,6 +24,7 @@ + @@ -111,14 +112,20 @@ Level3 Disabled - true $(CUDA_PATH)\include;$(ProjectDir)..\..\Common;$(ProjectDir)..\..\CUDA;%(AdditionalIncludeDirectories) + _WINDOWS;%(PreprocessorDefinitions) + MultiThreadedDebug true cudart.lib;cusolver.lib;cublas.lib;cublas_device.lib;%(AdditionalDependencies) $(CUDA_PATH)\lib\x64;%(AdditionalLibraryDirectories) + + copy "$(CUDA_PATH)\bin\cublas64_70.dll" $(OutputPath) +copy "$(CUDA_PATH)\bin\cusolver64_70.dll" $(OutputPath) +copy "$(CUDA_PATH)\bin\cudart64_70.dll" $(OutputPath) + @@ -143,8 +150,10 @@ MaxSpeed true true - true $(CUDA_PATH)\include;$(ProjectDir)..\..\Common;$(ProjectDir)..\..\CUDA;%(AdditionalIncludeDirectories) + _WINDOWS;%(PreprocessorDefinitions) + /Qvec-report:1 + MultiThreaded true @@ -153,6 +162,11 @@ cudart.lib;cusolver.lib;cublas.lib;%(AdditionalDependencies) $(CUDA_PATH)\lib\x64;%(AdditionalLibraryDirectories) + + copy "$(CUDA_PATH)\bin\cublas64_70.dll" $(OutputPath) +copy "$(CUDA_PATH)\bin\cusolver64_70.dll" $(OutputPath) +copy "$(CUDA_PATH)\bin\cudart64_70.dll" $(OutputPath) + diff --git a/src/NativeProviders/Windows/CUDA/CUDAWrapper.vcxproj.filters b/src/NativeProviders/Windows/CUDA/CUDAWrapper.vcxproj.filters index 9bda620d..3c9be8c5 100644 --- a/src/NativeProviders/Windows/CUDA/CUDAWrapper.vcxproj.filters +++ b/src/NativeProviders/Windows/CUDA/CUDAWrapper.vcxproj.filters @@ -29,5 +29,8 @@ Source Files + + Source Files + \ No newline at end of file diff --git a/src/Numerics/Control.cs b/src/Numerics/Control.cs index fa1028fc..0e9570fd 100644 --- a/src/Numerics/Control.cs +++ b/src/Numerics/Control.cs @@ -78,6 +78,10 @@ namespace MathNet.Numerics case "MKL": LinearAlgebraProvider = new Providers.LinearAlgebra.Mkl.MklLinearAlgebraProvider(); break; + + case "CUDA": + LinearAlgebraProvider = new Providers.LinearAlgebra.Cuda.CudaLinearAlgebraProvider(); + break; #endif default: LinearAlgebraProvider = new ManagedLinearAlgebraProvider(); @@ -127,6 +131,11 @@ namespace MathNet.Numerics { LinearAlgebraProvider = new Providers.LinearAlgebra.Mkl.MklLinearAlgebraProvider(consistency, precision, accuracy); } + + public static void UseNativeCUDA() + { + LinearAlgebraProvider = new Providers.LinearAlgebra.Cuda.CudaLinearAlgebraProvider(); + } #endif /// diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index b783234e..b451cf78 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -157,6 +157,12 @@ + + + + + + diff --git a/src/Numerics/Properties/AssemblyInfo.cs b/src/Numerics/Properties/AssemblyInfo.cs index 2281025c..930f5179 100644 --- a/src/Numerics/Properties/AssemblyInfo.cs +++ b/src/Numerics/Properties/AssemblyInfo.cs @@ -76,6 +76,7 @@ using System.Runtime.InteropServices; #else [assembly: InternalsVisibleTo("MathNet.Numerics.UnitTests")] [assembly: InternalsVisibleTo("MathNet.Numerics.UnitTestsMKL")] +[assembly: InternalsVisibleTo("MathNet.Numerics.UnitTestsCUDA")] [assembly: InternalsVisibleTo("Performance")] #endif diff --git a/src/Numerics/Properties/Resources.Designer.cs b/src/Numerics/Properties/Resources.Designer.cs index 0c7fae4d..443ba150 100644 --- a/src/Numerics/Properties/Resources.Designer.cs +++ b/src/Numerics/Properties/Resources.Designer.cs @@ -1,15 +1,13 @@ //------------------------------------------------------------------------------ // // This code was generated by a tool. -// Runtime Version:4.0.30319.34209 +// Runtime Version:4.0.30319.34014 // // Changes to this file may cause incorrect behavior and will be lost if // the code is regenerated. // //------------------------------------------------------------------------------ -using System.Reflection; - namespace MathNet.Numerics.Properties { using System; @@ -40,18 +38,10 @@ namespace MathNet.Numerics.Properties { [global::System.ComponentModel.EditorBrowsableAttribute(global::System.ComponentModel.EditorBrowsableState.Advanced)] public static global::System.Resources.ResourceManager ResourceManager { get { -#if NET45REFLECTION - if (object.ReferenceEquals(resourceMan, null)) - { - global::System.Resources.ResourceManager temp = new global::System.Resources.ResourceManager("MathNet.Numerics.Properties.Resources", typeof(Resources).GetTypeInfo().Assembly); - resourceMan = temp; - } -#else if (object.ReferenceEquals(resourceMan, null)) { global::System.Resources.ResourceManager temp = new global::System.Resources.ResourceManager("MathNet.Numerics.Properties.Resources", typeof(Resources).Assembly); resourceMan = temp; } -#endif return resourceMan; } } @@ -935,6 +925,15 @@ namespace MathNet.Numerics.Properties { } } + /// + /// Looks up a localized string similar to User work buffers are not supported by this provider.. + /// + public static string UserWorkBufferNotSupported { + get { + return ResourceManager.GetString("UserWorkBufferNotSupported", resourceCulture); + } + } + /// /// Looks up a localized string similar to Vectors can not be empty and must have at least one element.. /// diff --git a/src/Numerics/Properties/Resources.resx b/src/Numerics/Properties/Resources.resx index d46d44c8..e1360503 100644 --- a/src/Numerics/Properties/Resources.resx +++ b/src/Numerics/Properties/Resources.resx @@ -412,4 +412,7 @@ Vectors can not be empty and must have at least one element. + + User work buffers are not supported by this provider. + \ No newline at end of file diff --git a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex.cs b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex.cs new file mode 100644 index 00000000..bad3b091 --- /dev/null +++ b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex.cs @@ -0,0 +1,703 @@ +// +// 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-2013 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// 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; +using System.Numerics; +using System.Security; +using MathNet.Numerics.LinearAlgebra.Factorization; +using MathNet.Numerics.Properties; + +namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda +{ + /// + /// Intel's Math Kernel Library (MKL) linear algebra provider. + /// + public partial class CudaLinearAlgebraProvider + { + /// + /// 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(_blasHandle, x.Length, x, y); + } + + /// + /// Adds a scaled vector to another: result = y + alpha*x. + /// + /// The vector to update. + /// The value to scale by. + /// The vector to add to . + /// The result of the addition. + /// This is similar to the AXPY BLAS routine. + [SecuritySafeCritical] + public override void AddVectorToScaledVector(Complex[] y, Complex alpha, Complex[] x, Complex[] result) + { + if (y == null) + { + throw new ArgumentNullException("y"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (y.Length != x.Length) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength); + } + + if (!ReferenceEquals(y, result)) + { + Array.Copy(y, 0, result, 0, y.Length); + } + + if (alpha == Complex.Zero) + { + return; + } + + SafeNativeMethods.z_axpy(_blasHandle, y.Length, alpha, x, result); + } + + /// + /// Scales an array. Can be used to scale a vector and a matrix. + /// + /// The scalar. + /// The values to scale. + /// This result of the scaling. + /// This is similar to the SCAL BLAS routine. + [SecuritySafeCritical] + public override void ScaleArray(Complex alpha, Complex[] x, Complex[] result) + { + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (!ReferenceEquals(x, result)) + { + Array.Copy(x, 0, result, 0, x.Length); + } + + if (alpha == Complex.One) + { + return; + } + + SafeNativeMethods.z_scale(_blasHandle, x.Length, alpha, result); + } + + /// + /// Multiples two matrices. result = x * y + /// + /// The x matrix. + /// The number of rows in the x matrix. + /// The number of columns in the x matrix. + /// The y matrix. + /// The number of rows in the y matrix. + /// The number of columns in the y matrix. + /// Where to store the result of the multiplication. + /// This is a simplified version of the BLAS GEMM routine with alpha + /// set to Complex.One and beta set to Complex.Zero, and x and y are not transposed. + public override void MatrixMultiply(Complex[] x, int rowsX, int columnsX, Complex[] y, int rowsY, int columnsY, Complex[] result) + { + MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, Complex.One, x, rowsX, columnsX, y, rowsY, columnsY, Complex.Zero, result); + } + + /// + /// Multiplies two matrices and updates another with the result. c = alpha*op(a)*op(b) + beta*c + /// + /// How to transpose the matrix. + /// How to transpose the matrix. + /// The value to scale matrix. + /// The a matrix. + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The b matrix + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The value to scale the matrix. + /// The c matrix. + [SecuritySafeCritical] + public override void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] a, int rowsA, int columnsA, Complex[] b, int rowsB, int columnsB, Complex beta, Complex[] c) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (c == null) + { + throw new ArgumentNullException("c"); + } + + var m = transposeA == Transpose.DontTranspose ? rowsA : columnsA; + var n = transposeB == Transpose.DontTranspose ? columnsB : rowsB; + var k = transposeA == Transpose.DontTranspose ? columnsA : rowsA; + var l = transposeB == Transpose.DontTranspose ? rowsB : columnsB; + + if (c.Length != m*n) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + if (k != l) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + SafeNativeMethods.z_matrix_multiply(_blasHandle, transposeA.ToCUDA(), transposeB.ToCUDA(), m, n, k, alpha, a, b, beta, c); + } + + /// + /// Computes the LUP factorization of A. P*A = L*U. + /// + /// An by matrix. The matrix is overwritten with the + /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always Complex.One + /// for the L factor). The upper triangular factor U is stored on and above the diagonal of . + /// The order of the square matrix . + /// On exit, it contains the pivot indices. The size of the array must be . + /// This is equivalent to the GETRF LAPACK routine. + [SecuritySafeCritical] + public override void LUFactor(Complex[] data, int order, int[] ipiv) + { + if (data == null) + { + throw new ArgumentNullException("data"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (data.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "data"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + Solver(SafeNativeMethods.z_lu_factor(_solverHandle, order, data, ipiv)); + } + + /// + /// Computes the inverse of matrix using LU factorization. + /// + /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . + /// This is equivalent to the GETRF and GETRI LAPACK routines. + [SecuritySafeCritical] + public override void LUInverse(Complex[] a, int order) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + Solver(SafeNativeMethods.z_lu_inverse(_solverHandle, _blasHandle, order, a)); + } + + /// + /// Computes the inverse of a previously factored matrix. + /// + /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . + /// The pivot indices of . + /// This is equivalent to the GETRI LAPACK routine. + [SecuritySafeCritical] + public override void LUInverseFactored(Complex[] a, int order, int[] ipiv) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + BLAS(SafeNativeMethods.z_lu_inverse_factored(_blasHandle, order, a, ipiv)); + } + + /// + /// Computes the inverse of matrix using LU factorization. + /// + /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . + /// Not supported. Should be left null. + /// This is equivalent to the GETRF and GETRI LAPACK routines. + [SecuritySafeCritical] + public override void LUInverse(Complex[] a, int order, Complex[] work) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (work != null) + { + throw new NotSupportedException(Resources.UserWorkBufferNotSupported); + } + + Solver(SafeNativeMethods.z_lu_inverse(_solverHandle, _blasHandle, order, a)); + } + + /// + /// Computes the inverse of a previously factored matrix. + /// + /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . + /// The pivot indices of . + /// Not supported. Should be left null. + /// This is equivalent to the GETRI LAPACK routine. + [SecuritySafeCritical] + public override void LUInverseFactored(Complex[] a, int order, int[] ipiv, Complex[] work) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + if (work != null) + { + throw new NotSupportedException(Resources.UserWorkBufferNotSupported); + } + + BLAS(SafeNativeMethods.z_lu_inverse_factored(_blasHandle, order, a, ipiv)); + } + + /// + /// Solves A*X=B for X using LU factorization. + /// + /// The number of columns of B. + /// The square matrix A. + /// The order of the square matrix . + /// On entry the B matrix; on exit the X matrix. + /// This is equivalent to the GETRF and GETRS LAPACK routines. + [SecuritySafeCritical] + public override void LUSolve(int columnsOfB, Complex[] a, int order, Complex[] b) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (b.Length != columnsOfB*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.z_lu_solve(_solverHandle, order, columnsOfB, a, b)); + } + + /// + /// Solves A*X=B for X using a previously factored A matrix. + /// + /// The number of columns of B. + /// The factored A matrix. + /// The order of the square matrix . + /// The pivot indices of . + /// On entry the B matrix; on exit the X matrix. + /// This is equivalent to the GETRS LAPACK routine. + [SecuritySafeCritical] + public override void LUSolveFactored(int columnsOfB, Complex[] a, int order, int[] ipiv, Complex[] b) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + if (b.Length != columnsOfB*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.z_lu_solve_factored(_solverHandle, order, columnsOfB, a, ipiv, b)); + } + + /// + /// Computes the Cholesky factorization of A. + /// + /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the + /// the Cholesky factorization. + /// The number of rows or columns in the matrix. + /// This is equivalent to the POTRF LAPACK routine. + [SecuritySafeCritical] + public override void CholeskyFactor(Complex[] a, int order) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (order < 1) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + Solver(SafeNativeMethods.z_cholesky_factor(_solverHandle, order, a)); + } + + /// + /// Solves A*X=B for X using Cholesky factorization. + /// + /// The square, positive definite matrix A. + /// The number of rows and columns in A. + /// On entry the B matrix; on exit the X matrix. + /// The number of columns in the B matrix. + /// This is equivalent to the POTRF add POTRS LAPACK routines. + /// + [SecuritySafeCritical] + public override void CholeskySolve(Complex[] a, int orderA, Complex[] b, int columnsB) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (b.Length != orderA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.z_cholesky_solve(_solverHandle, orderA, columnsB, a, b)); + } + + /// + /// Solves A*X=B for X using a previously factored A matrix. + /// + /// The square, positive definite matrix A. + /// The number of rows and columns in A. + /// On entry the B matrix; on exit the X matrix. + /// The number of columns in the B matrix. + /// This is equivalent to the POTRS LAPACK routine. + [SecuritySafeCritical] + public override void CholeskySolveFactored(Complex[] a, int orderA, Complex[] b, int columnsB) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (b.Length != orderA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.z_cholesky_solve_factored(_solverHandle, orderA, columnsB, a, b)); + } + + /// + /// Computes the singular value decomposition of A. + /// + /// Compute the singular U and VT vectors or not. + /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The singular values of A in ascending value. + /// If is true, on exit U contains the left + /// singular vectors. + /// If is true, on exit VT contains the transposed + /// right singular vectors. + /// This is equivalent to the GESVD LAPACK routine. + [SecuritySafeCritical] + public override void SingularValueDecomposition(bool computeVectors, Complex[] a, int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (s == null) + { + throw new ArgumentNullException("s"); + } + + if (u == null) + { + throw new ArgumentNullException("u"); + } + + if (vt == null) + { + throw new ArgumentNullException("vt"); + } + + if (u.Length != rowsA*rowsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); + } + + if (vt.Length != columnsA*columnsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); + } + + if (s.Length != Math.Min(rowsA, columnsA)) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); + } + + SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, null); + } + + /// + /// Solves A*X=B for X using the singular value decomposition of A. + /// + /// On entry, the M by N matrix to decompose. + /// 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. + public override void SvdSolve(Complex[] a, int rowsA, int columnsA, Complex[] b, int columnsB, Complex[] x) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (b.Length != rowsA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (x.Length != columnsA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + var s = new Complex[Math.Min(rowsA, columnsA)]; + var u = new Complex[rowsA*rowsA]; + var vt = new Complex[columnsA*columnsA]; + + var clone = new Complex[a.Length]; + a.Copy(clone); + SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt, null); + SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); + } + + /// + /// Computes the singular value decomposition of A. + /// + /// Compute the singular U and VT vectors or not. + /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The singular values of A in ascending value. + /// If is true, on exit U contains the left + /// singular vectors. + /// If is true, on exit VT contains the transposed + /// right singular vectors. + /// User work buffers are not supported. Should be null. + /// This is equivalent to the GESVD LAPACK routine. + [SecuritySafeCritical] + public override void SingularValueDecomposition(bool computeVectors, Complex[] a, int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt, Complex[] work) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (s == null) + { + throw new ArgumentNullException("s"); + } + + if (u == null) + { + throw new ArgumentNullException("u"); + } + + if (vt == null) + { + throw new ArgumentNullException("vt"); + } + + if (work != null) + { + throw new NotSupportedException(Resources.UserWorkBufferNotSupported); + } + + if (u.Length != rowsA*rowsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); + } + + if (vt.Length != columnsA*columnsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); + } + + if (s.Length != Math.Min(rowsA, columnsA)) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); + } + + Solver(SafeNativeMethods.z_svd_factor(_solverHandle, computeVectors, rowsA, columnsA, a, s, u, vt)); + } + } +} + +#endif diff --git a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex32.cs b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex32.cs new file mode 100644 index 00000000..3e06022c --- /dev/null +++ b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex32.cs @@ -0,0 +1,703 @@ +// +// 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-2013 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// 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; +using System.Numerics; +using System.Security; +using MathNet.Numerics.LinearAlgebra.Factorization; +using MathNet.Numerics.Properties; + +namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda +{ + /// + /// Intel's Math Kernel Library (MKL) linear algebra provider. + /// + public partial class CudaLinearAlgebraProvider + { + /// + /// 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(_blasHandle, x.Length, x, y); + } + + /// + /// Adds a scaled vector to another: result = y + alpha*x. + /// + /// The vector to update. + /// The value to scale by. + /// The vector to add to . + /// The result of the addition. + /// This is similar to the AXPY BLAS routine. + [SecuritySafeCritical] + public override void AddVectorToScaledVector(Complex32[] y, Complex32 alpha, Complex32[] x, Complex32[] result) + { + if (y == null) + { + throw new ArgumentNullException("y"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (y.Length != x.Length) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength); + } + + if (!ReferenceEquals(y, result)) + { + Array.Copy(y, 0, result, 0, y.Length); + } + + if (alpha == Complex32.Zero) + { + return; + } + + SafeNativeMethods.c_axpy(_blasHandle, y.Length, alpha, x, result); + } + + /// + /// Scales an array. Can be used to scale a vector and a matrix. + /// + /// The scalar. + /// The values to scale. + /// This result of the scaling. + /// This is similar to the SCAL BLAS routine. + [SecuritySafeCritical] + public override void ScaleArray(Complex32 alpha, Complex32[] x, Complex32[] result) + { + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (!ReferenceEquals(x, result)) + { + Array.Copy(x, 0, result, 0, x.Length); + } + + if (alpha == Complex32.One) + { + return; + } + + SafeNativeMethods.c_scale(_blasHandle, x.Length, alpha, result); + } + + /// + /// Multiples two matrices. result = x * y + /// + /// The x matrix. + /// The number of rows in the x matrix. + /// The number of columns in the x matrix. + /// The y matrix. + /// The number of rows in the y matrix. + /// The number of columns in the y matrix. + /// Where to store the result of the multiplication. + /// This is a simplified version of the BLAS GEMM routine with alpha + /// set to Complex32.One and beta set to Complex32.Zero, and x and y are not transposed. + public override void MatrixMultiply(Complex32[] x, int rowsX, int columnsX, Complex32[] y, int rowsY, int columnsY, Complex32[] result) + { + MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, Complex32.One, x, rowsX, columnsX, y, rowsY, columnsY, Complex32.Zero, result); + } + + /// + /// Multiplies two matrices and updates another with the result. c = alpha*op(a)*op(b) + beta*c + /// + /// How to transpose the matrix. + /// How to transpose the matrix. + /// The value to scale matrix. + /// The a matrix. + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The b matrix + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The value to scale the matrix. + /// The c matrix. + [SecuritySafeCritical] + public override void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex32 alpha, Complex32[] a, int rowsA, int columnsA, Complex32[] b, int rowsB, int columnsB, Complex32 beta, Complex32[] c) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (c == null) + { + throw new ArgumentNullException("c"); + } + + var m = transposeA == Transpose.DontTranspose ? rowsA : columnsA; + var n = transposeB == Transpose.DontTranspose ? columnsB : rowsB; + var k = transposeA == Transpose.DontTranspose ? columnsA : rowsA; + var l = transposeB == Transpose.DontTranspose ? rowsB : columnsB; + + if (c.Length != m*n) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + if (k != l) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + SafeNativeMethods.c_matrix_multiply(_blasHandle, transposeA.ToCUDA(), transposeB.ToCUDA(), m, n, k, alpha, a, b, beta, c); + } + + /// + /// Computes the LUP factorization of A. P*A = L*U. + /// + /// An by matrix. The matrix is overwritten with the + /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always Complex32.One + /// for the L factor). The upper triangular factor U is stored on and above the diagonal of . + /// The order of the square matrix . + /// On exit, it contains the pivot indices. The size of the array must be . + /// This is equivalent to the GETRF LAPACK routine. + [SecuritySafeCritical] + public override void LUFactor(Complex32[] data, int order, int[] ipiv) + { + if (data == null) + { + throw new ArgumentNullException("data"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (data.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "data"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + Solver(SafeNativeMethods.c_lu_factor(_solverHandle, order, data, ipiv)); + } + + /// + /// Computes the inverse of matrix using LU factorization. + /// + /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . + /// This is equivalent to the GETRF and GETRI LAPACK routines. + [SecuritySafeCritical] + public override void LUInverse(Complex32[] a, int order) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + Solver(SafeNativeMethods.c_lu_inverse(_solverHandle, _blasHandle, order, a)); + } + + /// + /// Computes the inverse of a previously factored matrix. + /// + /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . + /// The pivot indices of . + /// This is equivalent to the GETRI LAPACK routine. + [SecuritySafeCritical] + public override void LUInverseFactored(Complex32[] a, int order, int[] ipiv) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + BLAS(SafeNativeMethods.c_lu_inverse_factored(_blasHandle, order, a, ipiv)); + } + + /// + /// Computes the inverse of matrix using LU factorization. + /// + /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . + /// Not supported. Should be left null. + /// This is equivalent to the GETRF and GETRI LAPACK routines. + [SecuritySafeCritical] + public override void LUInverse(Complex32[] a, int order, Complex32[] work) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (work != null) + { + throw new NotSupportedException(Resources.UserWorkBufferNotSupported); + } + + Solver(SafeNativeMethods.c_lu_inverse(_solverHandle, _blasHandle, order, a)); + } + + /// + /// Computes the inverse of a previously factored matrix. + /// + /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . + /// The pivot indices of . + /// Not supported. Should be left null. + /// This is equivalent to the GETRI LAPACK routine. + [SecuritySafeCritical] + public override void LUInverseFactored(Complex32[] a, int order, int[] ipiv, Complex32[] work) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + if (work != null) + { + throw new NotSupportedException(Resources.UserWorkBufferNotSupported); + } + + BLAS(SafeNativeMethods.c_lu_inverse_factored(_blasHandle, order, a, ipiv)); + } + + /// + /// Solves A*X=B for X using LU factorization. + /// + /// The number of columns of B. + /// The square matrix A. + /// The order of the square matrix . + /// On entry the B matrix; on exit the X matrix. + /// This is equivalent to the GETRF and GETRS LAPACK routines. + [SecuritySafeCritical] + public override void LUSolve(int columnsOfB, Complex32[] a, int order, Complex32[] b) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (b.Length != columnsOfB*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.c_lu_solve(_solverHandle, order, columnsOfB, a, b)); + } + + /// + /// Solves A*X=B for X using a previously factored A matrix. + /// + /// The number of columns of B. + /// The factored A matrix. + /// The order of the square matrix . + /// The pivot indices of . + /// On entry the B matrix; on exit the X matrix. + /// This is equivalent to the GETRS LAPACK routine. + [SecuritySafeCritical] + public override void LUSolveFactored(int columnsOfB, Complex32[] a, int order, int[] ipiv, Complex32[] b) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + if (b.Length != columnsOfB*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.c_lu_solve_factored(_solverHandle, order, columnsOfB, a, ipiv, b)); + } + + /// + /// Computes the Cholesky factorization of A. + /// + /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the + /// the Cholesky factorization. + /// The number of rows or columns in the matrix. + /// This is equivalent to the POTRF LAPACK routine. + [SecuritySafeCritical] + public override void CholeskyFactor(Complex32[] a, int order) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (order < 1) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + Solver(SafeNativeMethods.c_cholesky_factor(_solverHandle, order, a)); + } + + /// + /// Solves A*X=B for X using Cholesky factorization. + /// + /// The square, positive definite matrix A. + /// The number of rows and columns in A. + /// On entry the B matrix; on exit the X matrix. + /// The number of columns in the B matrix. + /// This is equivalent to the POTRF add POTRS LAPACK routines. + /// + [SecuritySafeCritical] + public override void CholeskySolve(Complex32[] a, int orderA, Complex32[] b, int columnsB) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (b.Length != orderA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.c_cholesky_solve(_solverHandle, orderA, columnsB, a, b)); + } + + /// + /// Solves A*X=B for X using a previously factored A matrix. + /// + /// The square, positive definite matrix A. + /// The number of rows and columns in A. + /// On entry the B matrix; on exit the X matrix. + /// The number of columns in the B matrix. + /// This is equivalent to the POTRS LAPACK routine. + [SecuritySafeCritical] + public override void CholeskySolveFactored(Complex32[] a, int orderA, Complex32[] b, int columnsB) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (b.Length != orderA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.c_cholesky_solve_factored(_solverHandle, orderA, columnsB, a, b)); + } + + /// + /// Computes the singular value decomposition of A. + /// + /// Compute the singular U and VT vectors or not. + /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The singular values of A in ascending value. + /// If is true, on exit U contains the left + /// singular vectors. + /// If is true, on exit VT contains the transposed + /// right singular vectors. + /// This is equivalent to the GESVD LAPACK routine. + [SecuritySafeCritical] + public override void SingularValueDecomposition(bool computeVectors, Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (s == null) + { + throw new ArgumentNullException("s"); + } + + if (u == null) + { + throw new ArgumentNullException("u"); + } + + if (vt == null) + { + throw new ArgumentNullException("vt"); + } + + if (u.Length != rowsA*rowsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); + } + + if (vt.Length != columnsA*columnsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); + } + + if (s.Length != Math.Min(rowsA, columnsA)) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); + } + + SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt); + } + + /// + /// Solves A*X=B for X using the singular value decomposition of A. + /// + /// On entry, the M by N matrix to decompose. + /// 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. + public override void SvdSolve(Complex32[] a, int rowsA, int columnsA, Complex32[] b, int columnsB, Complex32[] x) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (b.Length != rowsA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (x.Length != columnsA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + var s = new Complex32[Math.Min(rowsA, columnsA)]; + var u = new Complex32[rowsA*rowsA]; + var vt = new Complex32[columnsA*columnsA]; + + var clone = new Complex32[a.Length]; + a.Copy(clone); + SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt, null); + SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); + } + + /// + /// Computes the singular value decomposition of A. + /// + /// Compute the singular U and VT vectors or not. + /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The singular values of A in ascending value. + /// If is true, on exit U contains the left + /// singular vectors. + /// If is true, on exit VT contains the transposed + /// right singular vectors. + /// Not supported. Should be left null. + /// This is equivalent to the GESVD LAPACK routine. + [SecuritySafeCritical] + public override void SingularValueDecomposition(bool computeVectors, Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] work) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (s == null) + { + throw new ArgumentNullException("s"); + } + + if (u == null) + { + throw new ArgumentNullException("u"); + } + + if (vt == null) + { + throw new ArgumentNullException("vt"); + } + + if (work != null) + { + throw new NotSupportedException(Resources.UserWorkBufferNotSupported); + } + + if (u.Length != rowsA*rowsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); + } + + if (vt.Length != columnsA*columnsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); + } + + if (s.Length != Math.Min(rowsA, columnsA)) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); + } + + Solver(SafeNativeMethods.c_svd_factor(_solverHandle, computeVectors, rowsA, columnsA, a, s, u, vt)); + } + } +} + +#endif diff --git a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Double.cs b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Double.cs new file mode 100644 index 00000000..8213d8a5 --- /dev/null +++ b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Double.cs @@ -0,0 +1,704 @@ +// +// 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-2013 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// 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; +using System.Numerics; +using System.Security; +using MathNet.Numerics.LinearAlgebra.Factorization; +using MathNet.Numerics.Properties; + +namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda +{ + /// + /// Intel's Math Kernel Library (MKL) linear algebra provider. + /// + public partial class CudaLinearAlgebraProvider + { + /// + /// 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(_blasHandle, x.Length, x, y); + } + + /// + /// Adds a scaled vector to another: result = y + alpha*x. + /// + /// The vector to update. + /// The value to scale by. + /// The vector to add to . + /// The result of the addition. + /// This is similar to the AXPY BLAS routine. + [SecuritySafeCritical] + public override void AddVectorToScaledVector(double[] y, double alpha, double[] x, double[] result) + { + if (y == null) + { + throw new ArgumentNullException("y"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (y.Length != x.Length) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength); + } + + if (!ReferenceEquals(y, result)) + { + Array.Copy(y, 0, result, 0, y.Length); + } + + if (alpha == 0.0) + { + return; + } + + SafeNativeMethods.d_axpy(_blasHandle, y.Length, alpha, x, result); + } + + /// + /// Scales an array. Can be used to scale a vector and a matrix. + /// + /// The scalar. + /// The values to scale. + /// This result of the scaling. + /// This is similar to the SCAL BLAS routine. + [SecuritySafeCritical] + public override void ScaleArray(double alpha, double[] x, double[] result) + { + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (!ReferenceEquals(x, result)) + { + Array.Copy(x, 0, result, 0, x.Length); + } + + if (alpha == 1.0) + { + return; + } + + SafeNativeMethods.d_scale(_blasHandle, x.Length, alpha, result); + } + + /// + /// Multiples two matrices. result = x * y + /// + /// The x matrix. + /// The number of rows in the x matrix. + /// The number of columns in the x matrix. + /// The y matrix. + /// The number of rows in the y matrix. + /// The number of columns in the y matrix. + /// Where to store the result of the multiplication. + /// This is a simplified version of the BLAS GEMM routine with alpha + /// set to 1.0 and beta set to 0.0, and x and y are not transposed. + public override void MatrixMultiply(double[] x, int rowsX, int columnsX, double[] y, int rowsY, int columnsY, double[] result) + { + MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, 1.0, x, rowsX, columnsX, y, rowsY, columnsY, 0.0, result); + } + + /// + /// Multiplies two matrices and updates another with the result. c = alpha*op(a)*op(b) + beta*c + /// + /// How to transpose the matrix. + /// How to transpose the matrix. + /// The value to scale matrix. + /// The a matrix. + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The b matrix + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The value to scale the matrix. + /// The c matrix. + [SecuritySafeCritical] + public override void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, int rowsA, int columnsA, double[] b, int rowsB, int columnsB, double beta, double[] c) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (c == null) + { + throw new ArgumentNullException("c"); + } + + var m = transposeA == Transpose.DontTranspose ? rowsA : columnsA; + var n = transposeB == Transpose.DontTranspose ? columnsB : rowsB; + var k = transposeA == Transpose.DontTranspose ? columnsA : rowsA; + var l = transposeB == Transpose.DontTranspose ? rowsB : columnsB; + + if (c.Length != m*n) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + if (k != l) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + SafeNativeMethods.d_matrix_multiply(_blasHandle, transposeA.ToCUDA(), transposeB.ToCUDA(), m, n, k, alpha, a, b, beta, c); + } + + /// + /// Computes the LUP factorization of A. P*A = L*U. + /// + /// An by matrix. The matrix is overwritten with the + /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0 + /// for the L factor). The upper triangular factor U is stored on and above the diagonal of . + /// The order of the square matrix . + /// On exit, it contains the pivot indices. The size of the array must be . + /// This is equivalent to the GETRF LAPACK routine. + [SecuritySafeCritical] + public override void LUFactor(double[] data, int order, int[] ipiv) + { + if (data == null) + { + throw new ArgumentNullException("data"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (data.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "data"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + Solver(SafeNativeMethods.d_lu_factor(_solverHandle, order, data, ipiv)); + } + + /// + /// Computes the inverse of matrix using LU factorization. + /// + /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . + /// This is equivalent to the GETRF and GETRI LAPACK routines. + [SecuritySafeCritical] + public override void LUInverse(double[] a, int order) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + Solver(SafeNativeMethods.d_lu_inverse(_solverHandle, _blasHandle, order, a)); + } + + /// + /// Computes the inverse of a previously factored matrix. + /// + /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . + /// The pivot indices of . + /// This is equivalent to the GETRI LAPACK routine. + [SecuritySafeCritical] + public override void LUInverseFactored(double[] a, int order, int[] ipiv) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + BLAS(SafeNativeMethods.d_lu_inverse_factored(_blasHandle, order, a, ipiv)); + } + + /// + /// Computes the inverse of matrix using LU factorization. + /// + /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . + /// Not supported. Should be left null. + /// This is equivalent to the GETRF and GETRI LAPACK routines. + [SecuritySafeCritical] + public override void LUInverse(double[] a, int order, double[] work) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (work != null) + { + throw new NotSupportedException(Resources.UserWorkBufferNotSupported); + } + + Solver(SafeNativeMethods.d_lu_inverse(_solverHandle, _blasHandle, order, a)); + } + + /// + /// Computes the inverse of a previously factored matrix. + /// + /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . + /// The pivot indices of . + /// Not supported. Should be left null. + /// This is equivalent to the GETRI LAPACK routine. + [SecuritySafeCritical] + public override void LUInverseFactored(double[] a, int order, int[] ipiv, double[] work) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + if (work != null) + { + throw new NotSupportedException(Resources.UserWorkBufferNotSupported); + } + + BLAS(SafeNativeMethods.d_lu_inverse_factored(_blasHandle, order, a, ipiv)); + } + + /// + /// Solves A*X=B for X using LU factorization. + /// + /// The number of columns of B. + /// The square matrix A. + /// The order of the square matrix . + /// On entry the B matrix; on exit the X matrix. + /// This is equivalent to the GETRF and GETRS LAPACK routines. + [SecuritySafeCritical] + public override void LUSolve(int columnsOfB, double[] a, int order, double[] b) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (b.Length != columnsOfB*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.d_lu_solve(_solverHandle, order, columnsOfB, a, b)); + } + + /// + /// Solves A*X=B for X using a previously factored A matrix. + /// + /// The number of columns of B. + /// The factored A matrix. + /// The order of the square matrix . + /// The pivot indices of . + /// On entry the B matrix; on exit the X matrix. + /// This is equivalent to the GETRS LAPACK routine. + [SecuritySafeCritical] + public override void LUSolveFactored(int columnsOfB, double[] a, int order, int[] ipiv, double[] b) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + if (b.Length != columnsOfB*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.d_lu_solve_factored(_solverHandle, order, columnsOfB, a, ipiv, b)); + } + + /// + /// Computes the Cholesky factorization of A. + /// + /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the + /// the Cholesky factorization. + /// The number of rows or columns in the matrix. + /// This is equivalent to the POTRF LAPACK routine. + [SecuritySafeCritical] + public override void CholeskyFactor(double[] a, int order) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (order < 1) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + Solver(SafeNativeMethods.d_cholesky_factor(_solverHandle, order, a)); + } + + /// + /// Solves A*X=B for X using Cholesky factorization. + /// + /// The square, positive definite matrix A. + /// The number of rows and columns in A. + /// On entry the B matrix; on exit the X matrix. + /// The number of columns in the B matrix. + /// This is equivalent to the POTRF add POTRS LAPACK routines. + /// + [SecuritySafeCritical] + public override void CholeskySolve(double[] a, int orderA, double[] b, int columnsB) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (b.Length != orderA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.d_cholesky_solve(_solverHandle, orderA, columnsB, a, b)); + } + + /// + /// Solves A*X=B for X using a previously factored A matrix. + /// + /// The square, positive definite matrix A. + /// The number of rows and columns in A. + /// On entry the B matrix; on exit the X matrix. + /// The number of columns in the B matrix. + /// This is equivalent to the POTRS LAPACK routine. + [SecuritySafeCritical] + public override void CholeskySolveFactored(double[] a, int orderA, double[] b, int columnsB) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (b.Length != orderA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.d_cholesky_solve_factored(_solverHandle, orderA, columnsB, a, b)); + } + + /// + /// Computes the singular value decomposition of A. + /// + /// Compute the singular U and VT vectors or not. + /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The singular values of A in ascending value. + /// If is true, on exit U contains the left + /// singular vectors. + /// If is true, on exit VT contains the transposed + /// right singular vectors. + /// This is equivalent to the GESVD LAPACK routine. + [SecuritySafeCritical] + public override void SingularValueDecomposition(bool computeVectors, double[] a, int rowsA, int columnsA, double[] s, double[] u, double[] vt) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (s == null) + { + throw new ArgumentNullException("s"); + } + + if (u == null) + { + throw new ArgumentNullException("u"); + } + + if (vt == null) + { + throw new ArgumentNullException("vt"); + } + + if (u.Length != rowsA*rowsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); + } + + if (vt.Length != columnsA*columnsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); + } + + if (s.Length != Math.Min(rowsA, columnsA)) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); + } + + SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt); + } + + /// + /// Solves A*X=B for X using the singular value decomposition of A. + /// + /// On entry, the M by N matrix to decompose. + /// 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. + public override void SvdSolve(double[] a, int rowsA, int columnsA, double[] b, int columnsB, double[] x) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (b.Length != rowsA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (x.Length != columnsA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + var s = new double[Math.Min(rowsA, columnsA)]; + var u = new double[rowsA*rowsA]; + var vt = new double[columnsA*columnsA]; + + var clone = new double[a.Length]; + a.Copy(clone); + SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt); + SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); + } + + /// + /// Computes the singular value decomposition of A. + /// + /// Compute the singular U and VT vectors or not. + /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The singular values of A in ascending value. + /// If is true, on exit U contains the left + /// singular vectors. + /// If is true, on exit VT contains the transposed + /// right singular vectors. + /// Not supported. Should be left null. + /// This is equivalent to the GESVD LAPACK routine. + [SecuritySafeCritical] + public override void SingularValueDecomposition(bool computeVectors, double[] a, int rowsA, int columnsA, double[] s, double[] u, double[] vt, double[] work) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (s == null) + { + throw new ArgumentNullException("s"); + } + + if (u == null) + { + throw new ArgumentNullException("u"); + } + + if (vt == null) + { + throw new ArgumentNullException("vt"); + } + + if (work != null) + { + throw new NotSupportedException(Resources.UserWorkBufferNotSupported); + } + + if (u.Length != rowsA*rowsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); + } + + if (vt.Length != columnsA*columnsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); + } + + if (s.Length != Math.Min(rowsA, columnsA)) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); + } + + + Solver (SafeNativeMethods.d_svd_factor(_solverHandle, computeVectors, rowsA, columnsA, a, s, u, vt)); + } + } +} + +#endif diff --git a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Single.cs b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Single.cs new file mode 100644 index 00000000..f964e286 --- /dev/null +++ b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Single.cs @@ -0,0 +1,703 @@ +// +// 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-2013 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// 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; +using System.Numerics; +using System.Security; +using MathNet.Numerics.LinearAlgebra.Factorization; +using MathNet.Numerics.Properties; + +namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda +{ + /// + /// Intel's Math Kernel Library (MKL) linear algebra provider. + /// + public partial class CudaLinearAlgebraProvider + { + /// + /// 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(_blasHandle, x.Length, x, y); + } + + /// + /// Adds a scaled vector to another: result = y + alpha*x. + /// + /// The vector to update. + /// The value to scale by. + /// The vector to add to . + /// The result of the addition. + /// This is similar to the AXPY BLAS routine. + [SecuritySafeCritical] + public override void AddVectorToScaledVector(float[] y, float alpha, float[] x, float[] result) + { + if (y == null) + { + throw new ArgumentNullException("y"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (y.Length != x.Length) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength); + } + + if (!ReferenceEquals(y, result)) + { + Array.Copy(y, 0, result, 0, y.Length); + } + + if (alpha == 0.0f) + { + return; + } + + SafeNativeMethods.s_axpy(_blasHandle, y.Length, alpha, x, result); + } + + /// + /// Scales an array. Can be used to scale a vector and a matrix. + /// + /// The scalar. + /// The values to scale. + /// This result of the scaling. + /// This is similar to the SCAL BLAS routine. + [SecuritySafeCritical] + public override void ScaleArray(float alpha, float[] x, float[] result) + { + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (!ReferenceEquals(x, result)) + { + Array.Copy(x, 0, result, 0, x.Length); + } + + if (alpha == 1.0f) + { + return; + } + + SafeNativeMethods.s_scale(_blasHandle, x.Length, alpha, result); + } + + /// + /// Multiples two matrices. result = x * y + /// + /// The x matrix. + /// The number of rows in the x matrix. + /// The number of columns in the x matrix. + /// The y matrix. + /// The number of rows in the y matrix. + /// The number of columns in the y matrix. + /// Where to store the result of the multiplication. + /// This is a simplified version of the BLAS GEMM routine with alpha + /// set to 1.0f and beta set to 0.0f, and x and y are not transposed. + public override void MatrixMultiply(float[] x, int rowsX, int columnsX, float[] y, int rowsY, int columnsY, float[] result) + { + MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, 1.0f, x, rowsX, columnsX, y, rowsY, columnsY, 0.0f, result); + } + + /// + /// Multiplies two matrices and updates another with the result. c = alpha*op(a)*op(b) + beta*c + /// + /// How to transpose the matrix. + /// How to transpose the matrix. + /// The value to scale matrix. + /// The a matrix. + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The b matrix + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The value to scale the matrix. + /// The c matrix. + [SecuritySafeCritical] + public override void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, float alpha, float[] a, int rowsA, int columnsA, float[] b, int rowsB, int columnsB, float beta, float[] c) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (c == null) + { + throw new ArgumentNullException("c"); + } + + var m = transposeA == Transpose.DontTranspose ? rowsA : columnsA; + var n = transposeB == Transpose.DontTranspose ? columnsB : rowsB; + var k = transposeA == Transpose.DontTranspose ? columnsA : rowsA; + var l = transposeB == Transpose.DontTranspose ? rowsB : columnsB; + + if (c.Length != m*n) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + if (k != l) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + SafeNativeMethods.s_matrix_multiply(_blasHandle, transposeA.ToCUDA(), transposeB.ToCUDA(), m, n, k, alpha, a, b, beta, c); + } + + /// + /// Computes the LUP factorization of A. P*A = L*U. + /// + /// An by matrix. The matrix is overwritten with the + /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0f + /// for the L factor). The upper triangular factor U is stored on and above the diagonal of . + /// The order of the square matrix . + /// On exit, it contains the pivot indices. The size of the array must be . + /// This is equivalent to the GETRF LAPACK routine. + [SecuritySafeCritical] + public override void LUFactor(float[] data, int order, int[] ipiv) + { + if (data == null) + { + throw new ArgumentNullException("data"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (data.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "data"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + Solver(SafeNativeMethods.s_lu_factor(_solverHandle, order, data, ipiv)); + } + + /// + /// Computes the inverse of matrix using LU factorization. + /// + /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . + /// This is equivalent to the GETRF and GETRI LAPACK routines. + [SecuritySafeCritical] + public override void LUInverse(float[] a, int order) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + Solver(SafeNativeMethods.s_lu_inverse(_solverHandle, _blasHandle, order, a)); + } + + /// + /// Computes the inverse of a previously factored matrix. + /// + /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . + /// The pivot indices of . + /// This is equivalent to the GETRI LAPACK routine. + [SecuritySafeCritical] + public override void LUInverseFactored(float[] a, int order, int[] ipiv) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + BLAS(SafeNativeMethods.s_lu_inverse_factored(_blasHandle, order, a, ipiv)); + } + + /// + /// Computes the inverse of matrix using LU factorization. + /// + /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . + /// Not supported. Should be left null. + /// This is equivalent to the GETRF and GETRI LAPACK routines. + [SecuritySafeCritical] + public override void LUInverse(float[] a, int order, float[] work) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (work != null) + { + throw new NotSupportedException(Resources.UserWorkBufferNotSupported); + } + + Solver(SafeNativeMethods.s_lu_inverse(_solverHandle, _blasHandle, order, a)); + } + + /// + /// Computes the inverse of a previously factored matrix. + /// + /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . + /// The pivot indices of . + /// Not supported. This should be left null. + /// This is equivalent to the GETRI LAPACK routine. + [SecuritySafeCritical] + public override void LUInverseFactored(float[] a, int order, int[] ipiv, float[] work) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + if (work != null) + { + throw new NotSupportedException(Resources.UserWorkBufferNotSupported); + } + + BLAS(SafeNativeMethods.s_lu_inverse_factored(_blasHandle, order, a, ipiv)); + } + + /// + /// Solves A*X=B for X using LU factorization. + /// + /// The number of columns of B. + /// The square matrix A. + /// The order of the square matrix . + /// On entry the B matrix; on exit the X matrix. + /// This is equivalent to the GETRF and GETRS LAPACK routines. + [SecuritySafeCritical] + public override void LUSolve(int columnsOfB, float[] a, int order, float[] b) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (b.Length != columnsOfB*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.s_lu_solve(_solverHandle, order, columnsOfB, a, b)); + } + + /// + /// Solves A*X=B for X using a previously factored A matrix. + /// + /// The number of columns of B. + /// The factored A matrix. + /// The order of the square matrix . + /// The pivot indices of . + /// On entry the B matrix; on exit the X matrix. + /// This is equivalent to the GETRS LAPACK routine. + [SecuritySafeCritical] + public override void LUSolveFactored(int columnsOfB, float[] a, int order, int[] ipiv, float[] b) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + if (b.Length != columnsOfB*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.s_lu_solve_factored(_solverHandle, order, columnsOfB, a, ipiv, b)); + } + + /// + /// Computes the Cholesky factorization of A. + /// + /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the + /// the Cholesky factorization. + /// The number of rows or columns in the matrix. + /// This is equivalent to the POTRF LAPACK routine. + [SecuritySafeCritical] + public override void CholeskyFactor(float[] a, int order) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (order < 1) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + Solver(SafeNativeMethods.s_cholesky_factor(_solverHandle, order, a)); + } + + /// + /// Solves A*X=B for X using Cholesky factorization. + /// + /// The square, positive definite matrix A. + /// The number of rows and columns in A. + /// On entry the B matrix; on exit the X matrix. + /// The number of columns in the B matrix. + /// This is equivalent to the POTRF add POTRS LAPACK routines. + /// + [SecuritySafeCritical] + public override void CholeskySolve(float[] a, int orderA, float[] b, int columnsB) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (b.Length != orderA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.s_cholesky_solve(_solverHandle, orderA, columnsB, a, b)); + } + + /// + /// Solves A*X=B for X using a previously factored A matrix. + /// + /// The square, positive definite matrix A. + /// The number of rows and columns in A. + /// On entry the B matrix; on exit the X matrix. + /// The number of columns in the B matrix. + /// This is equivalent to the POTRS LAPACK routine. + [SecuritySafeCritical] + public override void CholeskySolveFactored(float[] a, int orderA, float[] b, int columnsB) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (b.Length != orderA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.s_cholesky_solve_factored(_solverHandle, orderA, columnsB, a, b)); + } + + /// + /// Computes the singular value decomposition of A. + /// + /// Compute the singular U and VT vectors or not. + /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The singular values of A in ascending value. + /// If is true, on exit U contains the left + /// singular vectors. + /// If is true, on exit VT contains the transposed + /// right singular vectors. + /// This is equivalent to the GESVD LAPACK routine. + [SecuritySafeCritical] + public override void SingularValueDecomposition(bool computeVectors, float[] a, int rowsA, int columnsA, float[] s, float[] u, float[] vt) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (s == null) + { + throw new ArgumentNullException("s"); + } + + if (u == null) + { + throw new ArgumentNullException("u"); + } + + if (vt == null) + { + throw new ArgumentNullException("vt"); + } + + if (u.Length != rowsA*rowsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); + } + + if (vt.Length != columnsA*columnsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); + } + + if (s.Length != Math.Min(rowsA, columnsA)) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); + } + + SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, null); + } + + /// + /// Solves A*X=B for X using the singular value decomposition of A. + /// + /// On entry, the M by N matrix to decompose. + /// 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. + public override void SvdSolve(float[] a, int rowsA, int columnsA, float[] b, int columnsB, float[] x) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (b.Length != rowsA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (x.Length != columnsA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + var s = new float[Math.Min(rowsA, columnsA)]; + var u = new float[rowsA*rowsA]; + var vt = new float[columnsA*columnsA]; + + var clone = new float[a.Length]; + a.Copy(clone); + SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt, null); + SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); + } + + /// + /// Computes the singular value decomposition of A. + /// + /// Compute the singular U and VT vectors or not. + /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The singular values of A in ascending value. + /// If is true, on exit U contains the left + /// singular vectors. + /// If is true, on exit VT contains the transposed + /// right singular vectors. + /// Not supported. Should be left null. + /// This is equivalent to the GESVD LAPACK routine. + [SecuritySafeCritical] + public override void SingularValueDecomposition(bool computeVectors, float[] a, int rowsA, int columnsA, float[] s, float[] u, float[] vt, float[] work) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (s == null) + { + throw new ArgumentNullException("s"); + } + + if (u == null) + { + throw new ArgumentNullException("u"); + } + + if (vt == null) + { + throw new ArgumentNullException("vt"); + } + + if (work != null) + { + throw new ArgumentException(Resources.UserWorkBufferNotSupported); + } + + if (u.Length != rowsA*rowsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); + } + + if (vt.Length != columnsA*columnsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); + } + + if (s.Length != Math.Min(rowsA, columnsA)) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); + } + + Solver(SafeNativeMethods.s_svd_factor(_solverHandle, computeVectors, rowsA, columnsA, a, s, u, vt)); + } + } +} + +#endif diff --git a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.cs new file mode 100644 index 00000000..fe45efa8 --- /dev/null +++ b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.cs @@ -0,0 +1,209 @@ +// +// 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. +// + +using System; + +#if NATIVE + +namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda +{ + /// + /// Consistency vs. performance trade-off between runs on different machines. + /// + + + /// + /// Intel's Math Kernel Library (MKL) linear algebra provider. + /// + public partial class CudaLinearAlgebraProvider : ManagedLinearAlgebraProvider, IDisposable + { + private int _nativeRevision; + private bool _nativeIX86; + private bool _nativeX64; + private bool _nativeIA64; + private IntPtr _blasHandle; + private IntPtr _solverHandle; + + /// + /// Sets the desired bit consistency on repeated identical computations on varying CPU architectures, + /// as a trade-off with performance. + /// + /// VML optimal precision and rounding. + /// VML accuracy mode. + [CLSCompliant(false)] + public CudaLinearAlgebraProvider() + { + } + + /// + /// Initialize and verify that the provided is indeed available. + /// If calling this method fails, consider to fall back to alternatives like the managed provider. + /// + public override void InitializeVerify() + { + int a, b, linearAlgebra; + try + { + // Load the native library + NativeProviderLoader.TryLoad(SafeNativeMethods.DllName); + + a = SafeNativeMethods.query_capability(0); + b = SafeNativeMethods.query_capability(1); + + _nativeIX86 = SafeNativeMethods.query_capability(8) > 0; + _nativeX64 = SafeNativeMethods.query_capability(9) > 0; + _nativeIA64 = SafeNativeMethods.query_capability(10) > 0; + + _nativeRevision = SafeNativeMethods.query_capability(64); + linearAlgebra = SafeNativeMethods.query_capability(128); + } + catch (DllNotFoundException e) + { + throw new NotSupportedException("Cuda Native Provider not found.", e); + } + catch (BadImageFormatException e) + { + throw new NotSupportedException("Cuda Native Provider found but failed to load. Please verify that the platform matches (x64 vs x32, Windows vs Linux).", e); + } + catch (EntryPointNotFoundException e) + { + throw new NotSupportedException("Cuda Native Provider does not support capability querying and is therefore not compatible. Consider upgrading to a newer version.", e); + } + + if (a != 0 || b != -1 || linearAlgebra <=0 || _nativeRevision < 1) + { + throw new NotSupportedException("Cuda Native Provider too old or not compatible. Consider upgrading to a newer version."); + } + + BLAS(SafeNativeMethods.createBLASHandle(ref _blasHandle)); + Solver(SafeNativeMethods.createSolverHandle(ref _solverHandle)); + } + + private void BLAS(int status) + { + switch (status) + { + case 0: // CUBLAS_STATUS_SUCCESS + return; + + case 1: // CUBLAS_STATUS_NOT_INITIALIZED + throw new Exception("The CUDA Runtime initialization failed"); + + case 2: // CUSOLVER_STATUS_ALLOC_FAILED + throw new OutOfMemoryException("The resources could not be allocated"); + + case 7: // CUBLAS_STATUS_INVALID_VALUE + throw new ArgumentException("Invalid value"); + + case 8: // CUBLAS_STATUS_ARCH_MISMATCH + throw new NotSupportedException("The device does not support this opeation."); + + case 11: // CUBLAS_STATUS_MAPPING_ERROR + throw new Exception("Mapping error."); + + case 13: // CUBLAS_STATUS_EXECUTION_FAILED + throw new Exception("Execution failed"); + + case 14: // CUBLAS_STATUS_INTERNAL_ERROR + throw new Exception("Internal error"); + + case 15: // CUBLAS_STATUS_NOT_SUPPORTED + throw new NotSupportedException(); + + case 16: // CUBLAS_STATUS_LICENSE_ERROR + throw new Exception("License error"); + + default: + throw new Exception("Unrecognized cuBLAS status code: " + status); + } + } + + private void Solver(int status) + { + switch (status) + { + case 0: // CUSOLVER_STATUS_SUCCESS + return; + + case 1: // CUSOLVER_STATUS_NOT_INITIALIZED + throw new Exception("The library was not initialized"); + + case 2: // CUSOLVER_STATUS_ALLOC_FAILED + throw new OutOfMemoryException("The resources could not be allocated"); + + case 3: // CUSOLVER_STATUS_INVALID_VALUE + throw new ArgumentException("Invalid value"); + + case 4: // CUSOLVER_STATUS_ARCH_MISMATCH + throw new NotSupportedException("The device does not support compute capability 2.0 and above"); + + case 5: // CUSOLVER_STATUS_MAPPING_ERROR + throw new Exception("Mapping error"); + + case 6: // CUSOLVER_STATUS_EXECUTION_FAILED + throw new Exception("Execution failed"); + + case 7: //CUSOLVER_STATUS_INTERNAL_ERROR + throw new Exception("Internal error"); + + case 8: // CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED + throw new NotSupportedException("Matrix type not supported"); + + case 9: // CUSOLVER_STATUS_NOT_SUPPORTED + throw new NotSupportedException(); + + case 10: // CUSOLVER_STATUS_ZERO_PIVOT + throw new Exception("Zero pivot"); + + case 11: //CUSOLVER_STATUS_INVALID_LICENSE + throw new Exception("Invalid license"); + + default: + throw new Exception("Unrecognized cuSolverDn status code: " + status); + + + } + } + + public override string ToString() + { + return string.Format("Nvidia CUDA ({1}; revision {0})", _nativeRevision, _nativeIX86 ? "x86" : _nativeX64 ? "x64" : _nativeIA64 ? "IA64" : "unknown"); + } + + + public void Dispose() + { + BLAS(SafeNativeMethods.destroyBLASHandle(_blasHandle)); + Solver(SafeNativeMethods.destroySolverHandle(_solverHandle)); + } + } +} + +#endif diff --git a/src/Numerics/Providers/LinearAlgebra/Cuda/SafeNativeMethods.cs b/src/Numerics/Providers/LinearAlgebra/Cuda/SafeNativeMethods.cs new file mode 100644 index 00000000..ba58a483 --- /dev/null +++ b/src/Numerics/Providers/LinearAlgebra/Cuda/SafeNativeMethods.cs @@ -0,0 +1,378 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009-2014 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; +using System.Numerics; +using System.Runtime.InteropServices; +using System.Security; + +namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda +{ + /// + /// P/Invoke methods to the native math libraries. + /// + [SuppressUnmanagedCodeSecurity] + [SecurityCritical] + internal static class SafeNativeMethods + { + // ReSharper disable InconsistentNaming + + /// + /// Name of the native DLL. + /// + const string _DllName = "MathNet.Numerics.CUDA.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 int createBLASHandle(ref IntPtr blasHandle); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int destroyBLASHandle(IntPtr blasHandle); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int createSolverHandle(ref IntPtr solverHandle); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int destroySolverHandle(IntPtr solverHandle); + + #region BLAS + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void s_axpy(IntPtr blasHandle, 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(IntPtr blasHandle, 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(IntPtr blasHandle, 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(IntPtr blasHandle, 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(IntPtr blasHandle, int n, float alpha, [Out] float[] x); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void d_scale(IntPtr blasHandle, int n, double alpha, [Out] double[] x); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void c_scale(IntPtr blasHandle, int n, Complex32 alpha, [In, Out] Complex32[] x); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void z_scale(IntPtr blasHandle, int n, Complex alpha, [In, Out] Complex[] x); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern float s_dot_product(IntPtr blasHandle, int n, float[] x, float[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern double d_dot_product(IntPtr blasHandle, int n, double[] x, double[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern Complex32 c_dot_product(IntPtr blasHandle, int n, Complex32[] x, Complex32[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern Complex z_dot_product(IntPtr blasHandle, int n, Complex[] x, Complex[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void s_matrix_multiply(IntPtr blasHandle, int transA, int 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(IntPtr blasHandle, int transA, int 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(IntPtr blasHandle, int transA, int 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(IntPtr blasHandle, int transA, int transB, int m, int n, int k, Complex alpha, Complex[] x, Complex[] y, Complex beta, [In, Out] Complex[] c); + + internal static int ToCUDA(this Transpose transpose) + { + switch (transpose) + { + case Transpose.DontTranspose: + return 0; + + case Transpose.Transpose: + return 1; + + case Transpose.ConjugateTranspose: + return 2; + + default: + throw new ArgumentException("Unsupported transpose: " + transpose); + } + } + + #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 double 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(IntPtr solverHandle, int n, [In, Out] float[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_cholesky_factor(IntPtr solverHandle, int n, [In, Out] double[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_cholesky_factor(IntPtr solverHandle, int n, [In, Out] Complex32[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_cholesky_factor(IntPtr solverHandle, int n, [In, Out] Complex[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_lu_factor(IntPtr solverHandle, 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(IntPtr solverHandle, 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(IntPtr solverHandle, 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(IntPtr solverHandle, 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(IntPtr solverHandle, IntPtr blasHandle, int n, [In, Out] float[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_lu_inverse(IntPtr solverHandle, IntPtr blasHandle, int n, [In, Out] double[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_lu_inverse(IntPtr solverHandle, IntPtr blasHandle, int n, [In, Out] Complex32[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_lu_inverse(IntPtr solverHandle, IntPtr blasHandle, int n, [In, Out] Complex[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_lu_inverse_factored(IntPtr blasHandle, 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_inverse_factored(IntPtr blasHandle, 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_inverse_factored(IntPtr blasHandle, 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_inverse_factored(IntPtr blasHandle, 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_solve_factored(IntPtr solverHandle, 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(IntPtr solverHandle, 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(IntPtr solverHandle, 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(IntPtr solverHandle, 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(IntPtr solverHandle, 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(IntPtr solverHandle, 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(IntPtr solverHandle, 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(IntPtr solverHandle, 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(IntPtr solverHandle, 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(IntPtr solverHandle, 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(IntPtr solverHandle, 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(IntPtr solverHandle, 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(IntPtr solverHandle, 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(IntPtr solverHandle, 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(IntPtr solverHandle, 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(IntPtr solverHandle, 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_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_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_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_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(IntPtr solverHandle, [MarshalAs(UnmanagedType.U1)] bool computeVectors, int m, int n, [In, Out] float[] a, [In, Out] float[] s, [In, Out] float[] u, [In, Out] float[] v); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_svd_factor(IntPtr solverHandle, [MarshalAs(UnmanagedType.U1)] bool computeVectors, int m, int n, [In, Out] double[] a, [In, Out] double[] s, [In, Out] double[] u, [In, Out] double[] v); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_svd_factor(IntPtr solverHandle, [MarshalAs(UnmanagedType.U1)] bool computeVectors, int m, int n, [In, Out] Complex32[] a, [In, Out] Complex32[] s, [In, Out] Complex32[] u, [In, Out] Complex32[] v); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_svd_factor(IntPtr solverHandle, [MarshalAs(UnmanagedType.U1)] bool computeVectors, int m, int n, [In, Out] Complex[] a, [In, Out] Complex[] s, [In, Out] Complex[] u, [In, Out] Complex[] v); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern int s_eigen([MarshalAs(UnmanagedType.U1)] bool isSymmetric, int n, [In] 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([MarshalAs(UnmanagedType.U1)] bool isSymmetric, int n, [In] 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([MarshalAs(UnmanagedType.U1)] bool isSymmetric, int n, [In] 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([MarshalAs(UnmanagedType.U1)] bool isSymmetric, int n, [In] Complex[] a, [In, Out] Complex[] vectors, [In, Out] Complex[] values, [In, Out] Complex[] d); + + #endregion LAPACK + + #region Vector Functions + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void s_vector_add(int n, float[] x, float[] y, [In, Out] float[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void s_vector_subtract(int n, float[] x, float[] y, [In, Out] float[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void s_vector_multiply(int n, float[] x, float[] y, [In, Out] float[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void s_vector_divide(int n, float[] x, float[] y, [In, Out] float[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void d_vector_add(int n, double[] x, double[] y, [In, Out] double[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void d_vector_subtract(int n, double[] x, double[] y, [In, Out] double[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void d_vector_multiply(int n, double[] x, double[] y, [In, Out] double[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void d_vector_divide(int n, double[] x, double[] y, [In, Out] double[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void c_vector_add(int n, Complex32[] x, Complex32[] y, [In, Out] Complex32[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void c_vector_subtract(int n, Complex32[] x, Complex32[] y, [In, Out] Complex32[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void c_vector_multiply(int n, Complex32[] x, Complex32[] y, [In, Out] Complex32[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void c_vector_divide(int n, Complex32[] x, Complex32[] y, [In, Out] Complex32[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void z_vector_add(int n, Complex[] x, Complex[] y, [In, Out] Complex[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void z_vector_subtract(int n, Complex[] x, Complex[] y, [In, Out] Complex[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void z_vector_multiply(int n, Complex[] x, Complex[] y, [In, Out] Complex[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void z_vector_divide(int n, Complex[] x, Complex[] y, [In, Out] Complex[] result); + + #endregion Vector Functions + + // ReSharper restore InconsistentNaming + } +} + +#endif diff --git a/src/UnitTests/UnitTests-CUDA.csproj b/src/UnitTests/UnitTests-CUDA.csproj new file mode 100644 index 00000000..598ce9a2 --- /dev/null +++ b/src/UnitTests/UnitTests-CUDA.csproj @@ -0,0 +1,348 @@ + + + + 10.0 + Debug + AnyCPU + 8.0.30703 + 2.0 + {E79C0395-01DC-4BC9-B86C-ED45790892C5} + Library + Properties + MathNet.Numerics.UnitTests + MathNet.Numerics.UnitTestsCUDA + v4.5 + 512 + ..\..\ + + + TRACE;CUDA + ..\..\out\CUDA\Windows\ + ..\..\obj\CUDA\Windows\x86\ + ..\..\obj\CUDA\Windows\x86\ + true + pdbonly + prompt + MinimumRecommendedRules.ruleset + 1591 + AnyCPU + + + TRACE;DEBUG;CUDA + ..\..\out\CUDA\Windows\ + ..\..\obj\CUDA\Windows\x86\ + ..\..\obj\CUDA\Windows\x86\ + false + full + true + prompt + 4 + 1591 + AnyCPU + + + + + + + + + + + + + + + + + + + data\Codeplex-5667.csv + Always + + + data\Github-Cureos-1.csv + Always + + + data\Matlab\A.mat + Always + + + data\Matlab\collection-nocompress.mat + Always + + + data\Matlab\collection.mat + Always + + + data\Matlab\complex.mat + Always + + + data\Matlab\sparse-large.mat + Always + + + data\Matlab\sparse-small.mat + Always + + + data\Matlab\sparse_complex.mat + Always + + + data\Matlab\v.mat + Always + + + data\NIST\AtmWtAgt.dat + Always + + + data\NIST\Bennett5.dat + Always + + + data\NIST\BoxBOD.dat + Always + + + data\NIST\Chwirut1.dat + Always + + + data\NIST\Chwirut2.dat + Always + + + data\NIST\DanWood.dat + Always + + + data\NIST\Eckerle4.dat + Always + + + data\NIST\ENSO.dat + Always + + + data\NIST\Filip.dat + Always + + + data\NIST\Gauss1.dat + Always + + + data\NIST\Gauss2.dat + Always + + + data\NIST\Gauss3.dat + Always + + + data\NIST\Hahn1.dat + Always + + + data\NIST\Kirby2.dat + Always + + + data\NIST\Lanczos1.dat + Always + + + data\NIST\Lanczos2.dat + Always + + + data\NIST\Lanczos3.dat + Always + + + data\NIST\Lew.dat + Always + + + data\NIST\Longley.dat + Always + + + data\NIST\Lottery.dat + Always + + + data\NIST\Mavro.dat + Always + + + data\NIST\MGH09.dat + Always + + + data\NIST\MGH10.dat + Always + + + data\NIST\MGH17.dat + Always + + + data\NIST\Michelso.dat + Always + + + data\NIST\Misra1a.dat + Always + + + data\NIST\Misra1b.dat + Always + + + data\NIST\Misra1c.dat + Always + + + data\NIST\Misra1d.dat + Always + + + data\NIST\Nelson.dat + Always + + + data\NIST\NoInt1.dat + Always + + + data\NIST\NoInt2.dat + Always + + + data\NIST\Norris.dat + Always + + + data\NIST\NumAcc1.dat + Always + + + data\NIST\NumAcc2.dat + Always + + + data\NIST\NumAcc3.dat + Always + + + data\NIST\NumAcc4.dat + Always + + + data\NIST\Pontius.dat + Always + + + data\NIST\Rat42.dat + Always + + + data\NIST\Rat43.dat + Always + + + data\NIST\Roszman1.dat + Always + + + data\NIST\SiRstvt.dat + Always + + + data\NIST\SmLs01t.dat + Always + + + data\NIST\SmLs02t.dat + Always + + + data\NIST\SmLs03t.dat + Always + + + data\NIST\SmLs04t.dat + Always + + + data\NIST\SmLs05t.dat + Always + + + data\NIST\SmLs06t.dat + Always + + + data\NIST\SmLs07t.dat + Always + + + data\NIST\SmLs08t.dat + Always + + + data\NIST\SmLs09t.dat + Always + + + data\NIST\Thurber.dat + Always + + + data\NIST\Wampler1.dat + Always + + + data\NIST\Wampler2.dat + Always + + + data\NIST\Wampler3.dat + Always + + + data\NIST\Wampler4.dat + Always + + + data\NIST\Wampler5.dat + Always + + + + data\NIST\Meixner.dat + Always + + + + + + {b7cae5f4-a23f-4438-b5be-41226618b695} + Numerics + + + + + + ..\..\packages\NUnit\lib\nunit.framework.dll + True + True + + + \ No newline at end of file diff --git a/src/UnitTests/UseLinearAlgebraProvider.cs b/src/UnitTests/UseLinearAlgebraProvider.cs index fc4a42d6..d45957e8 100644 --- a/src/UnitTests/UseLinearAlgebraProvider.cs +++ b/src/UnitTests/UseLinearAlgebraProvider.cs @@ -40,6 +40,9 @@ namespace MathNet.Numerics.UnitTests { #if !NET35 && NATIVE Control.UseNativeMKL(); +#endif +#if CUDA + Control.UseNativeCUDA(); #endif }