Browse Source

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.

pull/306/head
Matthew Johnson 12 years ago
parent
commit
8b7b617481
  1. 44
      MathNet.Numerics.NativeProviders.sln
  2. 6
      src/NativeProviders/CUDA/blas.cpp
  3. 79
      src/NativeProviders/CUDA/capabilities.cpp
  4. 125
      src/NativeProviders/CUDA/lapack.cpp
  5. 18
      src/NativeProviders/Windows/CUDA/CUDAWrapper.vcxproj
  6. 3
      src/NativeProviders/Windows/CUDA/CUDAWrapper.vcxproj.filters
  7. 9
      src/Numerics/Control.cs
  8. 6
      src/Numerics/Numerics.csproj
  9. 1
      src/Numerics/Properties/AssemblyInfo.cs
  10. 21
      src/Numerics/Properties/Resources.Designer.cs
  11. 3
      src/Numerics/Properties/Resources.resx
  12. 703
      src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex.cs
  13. 703
      src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex32.cs
  14. 704
      src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Double.cs
  15. 703
      src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Single.cs
  16. 209
      src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.cs
  17. 378
      src/Numerics/Providers/LinearAlgebra/Cuda/SafeNativeMethods.cs
  18. 348
      src/UnitTests/UnitTests-CUDA.csproj
  19. 3
      src/UnitTests/UseLinearAlgebraProvider.cs

44
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

6
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

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

125
src/NativeProviders/CUDA/lapack.cpp

@ -10,8 +10,6 @@ template<typename T, typename GETRF, typename GETRFBSIZE>
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<typename T, typename GETRF, typename GETRI, typename GETRFBSIZE>
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<typename T, typename GETRI>
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<typename T, typename GETRS>
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<typename T, typename GETRF, typename GETRS, typename GETRFBSIZE>
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);
}

18
src/NativeProviders/Windows/CUDA/CUDAWrapper.vcxproj

@ -24,6 +24,7 @@
<ItemGroup>
<ClCompile Include="..\..\Common\WindowsDLL.cpp" />
<ClCompile Include="..\..\CUDA\blas.cpp" />
<ClCompile Include="..\..\CUDA\capabilities.cpp" />
<ClCompile Include="..\..\CUDA\lapack.cpp" />
</ItemGroup>
<PropertyGroup Label="Globals">
@ -111,14 +112,20 @@
<ClCompile>
<WarningLevel>Level3</WarningLevel>
<Optimization>Disabled</Optimization>
<SDLCheck>true</SDLCheck>
<AdditionalIncludeDirectories>$(CUDA_PATH)\include;$(ProjectDir)..\..\Common;$(ProjectDir)..\..\CUDA;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
<PreprocessorDefinitions>_WINDOWS;%(PreprocessorDefinitions)</PreprocessorDefinitions>
<RuntimeLibrary>MultiThreadedDebug</RuntimeLibrary>
</ClCompile>
<Link>
<GenerateDebugInformation>true</GenerateDebugInformation>
<AdditionalDependencies>cudart.lib;cusolver.lib;cublas.lib;cublas_device.lib;%(AdditionalDependencies)</AdditionalDependencies>
<AdditionalLibraryDirectories>$(CUDA_PATH)\lib\x64;%(AdditionalLibraryDirectories)</AdditionalLibraryDirectories>
</Link>
<PostBuildEvent>
<Command>copy "$(CUDA_PATH)\bin\cublas64_70.dll" $(OutputPath)
copy "$(CUDA_PATH)\bin\cusolver64_70.dll" $(OutputPath)
copy "$(CUDA_PATH)\bin\cudart64_70.dll" $(OutputPath)</Command>
</PostBuildEvent>
</ItemDefinitionGroup>
<ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">
<ClCompile>
@ -143,8 +150,10 @@
<Optimization>MaxSpeed</Optimization>
<FunctionLevelLinking>true</FunctionLevelLinking>
<IntrinsicFunctions>true</IntrinsicFunctions>
<SDLCheck>true</SDLCheck>
<AdditionalIncludeDirectories>$(CUDA_PATH)\include;$(ProjectDir)..\..\Common;$(ProjectDir)..\..\CUDA;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
<PreprocessorDefinitions>_WINDOWS;%(PreprocessorDefinitions)</PreprocessorDefinitions>
<AdditionalOptions>/Qvec-report:1 </AdditionalOptions>
<RuntimeLibrary>MultiThreaded</RuntimeLibrary>
</ClCompile>
<Link>
<GenerateDebugInformation>true</GenerateDebugInformation>
@ -153,6 +162,11 @@
<AdditionalDependencies>cudart.lib;cusolver.lib;cublas.lib;%(AdditionalDependencies)</AdditionalDependencies>
<AdditionalLibraryDirectories>$(CUDA_PATH)\lib\x64;%(AdditionalLibraryDirectories)</AdditionalLibraryDirectories>
</Link>
<PostBuildEvent>
<Command>copy "$(CUDA_PATH)\bin\cublas64_70.dll" $(OutputPath)
copy "$(CUDA_PATH)\bin\cusolver64_70.dll" $(OutputPath)
copy "$(CUDA_PATH)\bin\cudart64_70.dll" $(OutputPath)</Command>
</PostBuildEvent>
</ItemDefinitionGroup>
<Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />
<ImportGroup Label="ExtensionTargets">

3
src/NativeProviders/Windows/CUDA/CUDAWrapper.vcxproj.filters

@ -29,5 +29,8 @@
<ClCompile Include="..\..\CUDA\blas.cpp">
<Filter>Source Files</Filter>
</ClCompile>
<ClCompile Include="..\..\CUDA\capabilities.cpp">
<Filter>Source Files</Filter>
</ClCompile>
</ItemGroup>
</Project>

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

6
src/Numerics/Numerics.csproj

@ -157,6 +157,12 @@
<Compile Include="Providers\LinearAlgebra\Acml\AcmlLinearAlgebraProvider.Double.cs" />
<Compile Include="Providers\LinearAlgebra\Acml\AcmlLinearAlgebraProvider.Single.cs" />
<Compile Include="Providers\LinearAlgebra\Acml\SafeNativeMethods.cs" />
<Compile Include="Providers\LinearAlgebra\Cuda\CudaLinearAlgebraProvider.Complex.cs" />
<Compile Include="Providers\LinearAlgebra\Cuda\CudaLinearAlgebraProvider.Complex32.cs" />
<Compile Include="Providers\LinearAlgebra\Cuda\CudaLinearAlgebraProvider.cs" />
<Compile Include="Providers\LinearAlgebra\Cuda\CudaLinearAlgebraProvider.Double.cs" />
<Compile Include="Providers\LinearAlgebra\Cuda\CudaLinearAlgebraProvider.Single.cs" />
<Compile Include="Providers\LinearAlgebra\Cuda\SafeNativeMethods.cs" />
<Compile Include="Providers\LinearAlgebra\GotoBlas\GotoBlasLinearAlgebraProvider.cs" />
<Compile Include="Providers\LinearAlgebra\GotoBlas\GotoBlasLinearAlgebraProvider.Complex.cs" />
<Compile Include="Providers\LinearAlgebra\GotoBlas\GotoBlasLinearAlgebraProvider.Complex32.cs" />

1
src/Numerics/Properties/AssemblyInfo.cs

@ -76,6 +76,7 @@ using System.Runtime.InteropServices;
#else
[assembly: InternalsVisibleTo("MathNet.Numerics.UnitTests")]
[assembly: InternalsVisibleTo("MathNet.Numerics.UnitTestsMKL")]
[assembly: InternalsVisibleTo("MathNet.Numerics.UnitTestsCUDA")]
[assembly: InternalsVisibleTo("Performance")]
#endif

21
src/Numerics/Properties/Resources.Designer.cs

@ -1,15 +1,13 @@
//------------------------------------------------------------------------------
// <auto-generated>
// 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.
// </auto-generated>
//------------------------------------------------------------------------------
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 {
}
}
/// <summary>
/// Looks up a localized string similar to User work buffers are not supported by this provider..
/// </summary>
public static string UserWorkBufferNotSupported {
get {
return ResourceManager.GetString("UserWorkBufferNotSupported", resourceCulture);
}
}
/// <summary>
/// Looks up a localized string similar to Vectors can not be empty and must have at least one element..
/// </summary>

3
src/Numerics/Properties/Resources.resx

@ -412,4 +412,7 @@
<data name="VectorCanNotBeEmpty" xml:space="preserve">
<value>Vectors can not be empty and must have at least one element.</value>
</data>
<data name="UserWorkBufferNotSupported" xml:space="preserve">
<value>User work buffers are not supported by this provider.</value>
</data>
</root>

703
src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex.cs

@ -0,0 +1,703 @@
// <copyright file="MklLinearAlgebraProvider.Complex.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-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.
// </copyright>
#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
{
/// <summary>
/// Intel's Math Kernel Library (MKL) linear algebra provider.
/// </summary>
public partial class CudaLinearAlgebraProvider
{
/// <summary>
/// Computes the dot product of x and y.
/// </summary>
/// <param name="x">The vector x.</param>
/// <param name="y">The vector y.</param>
/// <returns>The dot product of x and y.</returns>
/// <remarks>This is equivalent to the DOT BLAS routine.</remarks>
[SecuritySafeCritical]
public override Complex DotProduct(Complex[] x, Complex[] y)
{
if (y == null)
{
throw new ArgumentNullException("y");
}
if (x == null)
{
throw new ArgumentNullException("x");
}
if (x.Length != y.Length)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength);
}
return SafeNativeMethods.z_dot_product(_blasHandle, x.Length, x, y);
}
/// <summary>
/// Adds a scaled vector to another: <c>result = y + alpha*x</c>.
/// </summary>
/// <param name="y">The vector to update.</param>
/// <param name="alpha">The value to scale <paramref name="x"/> by.</param>
/// <param name="x">The vector to add to <paramref name="y"/>.</param>
/// <param name="result">The result of the addition.</param>
/// <remarks>This is similar to the AXPY BLAS routine.</remarks>
[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);
}
/// <summary>
/// Scales an array. Can be used to scale a vector and a matrix.
/// </summary>
/// <param name="alpha">The scalar.</param>
/// <param name="x">The values to scale.</param>
/// <param name="result">This result of the scaling.</param>
/// <remarks>This is similar to the SCAL BLAS routine.</remarks>
[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);
}
/// <summary>
/// Multiples two matrices. <c>result = x * y</c>
/// </summary>
/// <param name="x">The x matrix.</param>
/// <param name="rowsX">The number of rows in the x matrix.</param>
/// <param name="columnsX">The number of columns in the x matrix.</param>
/// <param name="y">The y matrix.</param>
/// <param name="rowsY">The number of rows in the y matrix.</param>
/// <param name="columnsY">The number of columns in the y matrix.</param>
/// <param name="result">Where to store the result of the multiplication.</param>
/// <remarks>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.</remarks>
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);
}
/// <summary>
/// Multiplies two matrices and updates another with the result. <c>c = alpha*op(a)*op(b) + beta*c</c>
/// </summary>
/// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
/// <param name="transposeB">How to transpose the <paramref name="b"/> matrix.</param>
/// <param name="alpha">The value to scale <paramref name="a"/> matrix.</param>
/// <param name="a">The a matrix.</param>
/// <param name="rowsA">The number of rows in the <paramref name="a"/> matrix.</param>
/// <param name="columnsA">The number of columns in the <paramref name="a"/> matrix.</param>
/// <param name="b">The b matrix</param>
/// <param name="rowsB">The number of rows in the <paramref name="b"/> matrix.</param>
/// <param name="columnsB">The number of columns in the <paramref name="b"/> matrix.</param>
/// <param name="beta">The value to scale the <paramref name="c"/> matrix.</param>
/// <param name="c">The c matrix.</param>
[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);
}
/// <summary>
/// Computes the LUP factorization of A. P*A = L*U.
/// </summary>
/// <param name="data">An <paramref name="order"/> by <paramref name="order"/> matrix. The matrix is overwritten with the
/// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of <paramref name="data"/> (the diagonal is always Complex.One
/// for the L factor). The upper triangular factor U is stored on and above the diagonal of <paramref name="data"/>.</param>
/// <param name="order">The order of the square matrix <paramref name="data"/>.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size of the array must be <paramref name="order"/>.</param>
/// <remarks>This is equivalent to the GETRF LAPACK routine.</remarks>
[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));
}
/// <summary>
/// Computes the inverse of matrix using LU factorization.
/// </summary>
/// <param name="a">The N by N matrix to invert. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <remarks>This is equivalent to the GETRF and GETRI LAPACK routines.</remarks>
[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));
}
/// <summary>
/// Computes the inverse of a previously factored matrix.
/// </summary>
/// <param name="a">The LU factored N by N matrix. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <remarks>This is equivalent to the GETRI LAPACK routine.</remarks>
[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));
}
/// <summary>
/// Computes the inverse of matrix using LU factorization.
/// </summary>
/// <param name="a">The N by N matrix to invert. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="work">Not supported. Should be left null.</param>
/// <remarks>This is equivalent to the GETRF and GETRI LAPACK routines.</remarks>
[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));
}
/// <summary>
/// Computes the inverse of a previously factored matrix.
/// </summary>
/// <param name="a">The LU factored N by N matrix. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="work">Not supported. Should be left null.</param>
/// <remarks>This is equivalent to the GETRI LAPACK routine.</remarks>
[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));
}
/// <summary>
/// Solves A*X=B for X using LU factorization.
/// </summary>
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The square matrix A.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <remarks>This is equivalent to the GETRF and GETRS LAPACK routines.</remarks>
[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));
}
/// <summary>
/// Solves A*X=B for X using a previously factored A matrix.
/// </summary>
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The factored A matrix.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <remarks>This is equivalent to the GETRS LAPACK routine.</remarks>
[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));
}
/// <summary>
/// Computes the Cholesky factorization of A.
/// </summary>
/// <param name="a">On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the
/// the Cholesky factorization.</param>
/// <param name="order">The number of rows or columns in the matrix.</param>
/// <remarks>This is equivalent to the POTRF LAPACK routine.</remarks>
[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));
}
/// <summary>
/// Solves A*X=B for X using Cholesky factorization.
/// </summary>
/// <param name="a">The square, positive definite matrix A.</param>
/// <param name="orderA">The number of rows and columns in A.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns in the B matrix.</param>
/// <remarks>This is equivalent to the POTRF add POTRS LAPACK routines.
/// </remarks>
[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));
}
/// <summary>
/// Solves A*X=B for X using a previously factored A matrix.
/// </summary>
/// <param name="a">The square, positive definite matrix A.</param>
/// <param name="orderA">The number of rows and columns in A.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns in the B matrix.</param>
/// <remarks>This is equivalent to the POTRS LAPACK routine.</remarks>
[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));
}
/// <summary>
/// Computes the singular value decomposition of A.
/// </summary>
/// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">If <paramref name="computeVectors"/> is <c>true</c>, on exit U contains the left
/// singular vectors.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
[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);
}
/// <summary>
/// Solves A*X=B for X using the singular value decomposition of A.
/// </summary>
/// <param name="a">On entry, the M by N matrix to decompose.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
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);
}
/// <summary>
/// Computes the singular value decomposition of A.
/// </summary>
/// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">If <paramref name="computeVectors"/> is <c>true</c>, on exit U contains the left
/// singular vectors.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <param name="work">User work buffers are not supported. Should be null.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
[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

703
src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex32.cs

@ -0,0 +1,703 @@
// <copyright file="MklLinearAlgebraProvider.Complex32.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-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.
// </copyright>
#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
{
/// <summary>
/// Intel's Math Kernel Library (MKL) linear algebra provider.
/// </summary>
public partial class CudaLinearAlgebraProvider
{
/// <summary>
/// Computes the dot product of x and y.
/// </summary>
/// <param name="x">The vector x.</param>
/// <param name="y">The vector y.</param>
/// <returns>The dot product of x and y.</returns>
/// <remarks>This is equivalent to the DOT BLAS routine.</remarks>
[SecuritySafeCritical]
public override Complex32 DotProduct(Complex32[] x, Complex32[] y)
{
if (y == null)
{
throw new ArgumentNullException("y");
}
if (x == null)
{
throw new ArgumentNullException("x");
}
if (x.Length != y.Length)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength);
}
return SafeNativeMethods.c_dot_product(_blasHandle, x.Length, x, y);
}
/// <summary>
/// Adds a scaled vector to another: <c>result = y + alpha*x</c>.
/// </summary>
/// <param name="y">The vector to update.</param>
/// <param name="alpha">The value to scale <paramref name="x"/> by.</param>
/// <param name="x">The vector to add to <paramref name="y"/>.</param>
/// <param name="result">The result of the addition.</param>
/// <remarks>This is similar to the AXPY BLAS routine.</remarks>
[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);
}
/// <summary>
/// Scales an array. Can be used to scale a vector and a matrix.
/// </summary>
/// <param name="alpha">The scalar.</param>
/// <param name="x">The values to scale.</param>
/// <param name="result">This result of the scaling.</param>
/// <remarks>This is similar to the SCAL BLAS routine.</remarks>
[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);
}
/// <summary>
/// Multiples two matrices. <c>result = x * y</c>
/// </summary>
/// <param name="x">The x matrix.</param>
/// <param name="rowsX">The number of rows in the x matrix.</param>
/// <param name="columnsX">The number of columns in the x matrix.</param>
/// <param name="y">The y matrix.</param>
/// <param name="rowsY">The number of rows in the y matrix.</param>
/// <param name="columnsY">The number of columns in the y matrix.</param>
/// <param name="result">Where to store the result of the multiplication.</param>
/// <remarks>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.</remarks>
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);
}
/// <summary>
/// Multiplies two matrices and updates another with the result. <c>c = alpha*op(a)*op(b) + beta*c</c>
/// </summary>
/// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
/// <param name="transposeB">How to transpose the <paramref name="b"/> matrix.</param>
/// <param name="alpha">The value to scale <paramref name="a"/> matrix.</param>
/// <param name="a">The a matrix.</param>
/// <param name="rowsA">The number of rows in the <paramref name="a"/> matrix.</param>
/// <param name="columnsA">The number of columns in the <paramref name="a"/> matrix.</param>
/// <param name="b">The b matrix</param>
/// <param name="rowsB">The number of rows in the <paramref name="b"/> matrix.</param>
/// <param name="columnsB">The number of columns in the <paramref name="b"/> matrix.</param>
/// <param name="beta">The value to scale the <paramref name="c"/> matrix.</param>
/// <param name="c">The c matrix.</param>
[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);
}
/// <summary>
/// Computes the LUP factorization of A. P*A = L*U.
/// </summary>
/// <param name="data">An <paramref name="order"/> by <paramref name="order"/> matrix. The matrix is overwritten with the
/// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of <paramref name="data"/> (the diagonal is always Complex32.One
/// for the L factor). The upper triangular factor U is stored on and above the diagonal of <paramref name="data"/>.</param>
/// <param name="order">The order of the square matrix <paramref name="data"/>.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size of the array must be <paramref name="order"/>.</param>
/// <remarks>This is equivalent to the GETRF LAPACK routine.</remarks>
[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));
}
/// <summary>
/// Computes the inverse of matrix using LU factorization.
/// </summary>
/// <param name="a">The N by N matrix to invert. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <remarks>This is equivalent to the GETRF and GETRI LAPACK routines.</remarks>
[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));
}
/// <summary>
/// Computes the inverse of a previously factored matrix.
/// </summary>
/// <param name="a">The LU factored N by N matrix. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <remarks>This is equivalent to the GETRI LAPACK routine.</remarks>
[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));
}
/// <summary>
/// Computes the inverse of matrix using LU factorization.
/// </summary>
/// <param name="a">The N by N matrix to invert. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="work">Not supported. Should be left null.</param>
/// <remarks>This is equivalent to the GETRF and GETRI LAPACK routines.</remarks>
[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));
}
/// <summary>
/// Computes the inverse of a previously factored matrix.
/// </summary>
/// <param name="a">The LU factored N by N matrix. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="work">Not supported. Should be left null.</param>
/// <remarks>This is equivalent to the GETRI LAPACK routine.</remarks>
[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));
}
/// <summary>
/// Solves A*X=B for X using LU factorization.
/// </summary>
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The square matrix A.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <remarks>This is equivalent to the GETRF and GETRS LAPACK routines.</remarks>
[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));
}
/// <summary>
/// Solves A*X=B for X using a previously factored A matrix.
/// </summary>
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The factored A matrix.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <remarks>This is equivalent to the GETRS LAPACK routine.</remarks>
[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));
}
/// <summary>
/// Computes the Cholesky factorization of A.
/// </summary>
/// <param name="a">On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the
/// the Cholesky factorization.</param>
/// <param name="order">The number of rows or columns in the matrix.</param>
/// <remarks>This is equivalent to the POTRF LAPACK routine.</remarks>
[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));
}
/// <summary>
/// Solves A*X=B for X using Cholesky factorization.
/// </summary>
/// <param name="a">The square, positive definite matrix A.</param>
/// <param name="orderA">The number of rows and columns in A.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns in the B matrix.</param>
/// <remarks>This is equivalent to the POTRF add POTRS LAPACK routines.
/// </remarks>
[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));
}
/// <summary>
/// Solves A*X=B for X using a previously factored A matrix.
/// </summary>
/// <param name="a">The square, positive definite matrix A.</param>
/// <param name="orderA">The number of rows and columns in A.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns in the B matrix.</param>
/// <remarks>This is equivalent to the POTRS LAPACK routine.</remarks>
[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));
}
/// <summary>
/// Computes the singular value decomposition of A.
/// </summary>
/// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">If <paramref name="computeVectors"/> is <c>true</c>, on exit U contains the left
/// singular vectors.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
[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);
}
/// <summary>
/// Solves A*X=B for X using the singular value decomposition of A.
/// </summary>
/// <param name="a">On entry, the M by N matrix to decompose.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
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);
}
/// <summary>
/// Computes the singular value decomposition of A.
/// </summary>
/// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">If <paramref name="computeVectors"/> is <c>true</c>, on exit U contains the left
/// singular vectors.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <param name="work">Not supported. Should be left null.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
[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

704
src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Double.cs

@ -0,0 +1,704 @@
// <copyright file="MklLinearAlgebraProvider.Double.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-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.
// </copyright>
#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
{
/// <summary>
/// Intel's Math Kernel Library (MKL) linear algebra provider.
/// </summary>
public partial class CudaLinearAlgebraProvider
{
/// <summary>
/// Computes the dot product of x and y.
/// </summary>
/// <param name="x">The vector x.</param>
/// <param name="y">The vector y.</param>
/// <returns>The dot product of x and y.</returns>
/// <remarks>This is equivalent to the DOT BLAS routine.</remarks>
[SecuritySafeCritical]
public override double DotProduct(double[] x, double[] y)
{
if (y == null)
{
throw new ArgumentNullException("y");
}
if (x == null)
{
throw new ArgumentNullException("x");
}
if (x.Length != y.Length)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength);
}
return SafeNativeMethods.d_dot_product(_blasHandle, x.Length, x, y);
}
/// <summary>
/// Adds a scaled vector to another: <c>result = y + alpha*x</c>.
/// </summary>
/// <param name="y">The vector to update.</param>
/// <param name="alpha">The value to scale <paramref name="x"/> by.</param>
/// <param name="x">The vector to add to <paramref name="y"/>.</param>
/// <param name="result">The result of the addition.</param>
/// <remarks>This is similar to the AXPY BLAS routine.</remarks>
[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);
}
/// <summary>
/// Scales an array. Can be used to scale a vector and a matrix.
/// </summary>
/// <param name="alpha">The scalar.</param>
/// <param name="x">The values to scale.</param>
/// <param name="result">This result of the scaling.</param>
/// <remarks>This is similar to the SCAL BLAS routine.</remarks>
[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);
}
/// <summary>
/// Multiples two matrices. <c>result = x * y</c>
/// </summary>
/// <param name="x">The x matrix.</param>
/// <param name="rowsX">The number of rows in the x matrix.</param>
/// <param name="columnsX">The number of columns in the x matrix.</param>
/// <param name="y">The y matrix.</param>
/// <param name="rowsY">The number of rows in the y matrix.</param>
/// <param name="columnsY">The number of columns in the y matrix.</param>
/// <param name="result">Where to store the result of the multiplication.</param>
/// <remarks>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.</remarks>
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);
}
/// <summary>
/// Multiplies two matrices and updates another with the result. <c>c = alpha*op(a)*op(b) + beta*c</c>
/// </summary>
/// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
/// <param name="transposeB">How to transpose the <paramref name="b"/> matrix.</param>
/// <param name="alpha">The value to scale <paramref name="a"/> matrix.</param>
/// <param name="a">The a matrix.</param>
/// <param name="rowsA">The number of rows in the <paramref name="a"/> matrix.</param>
/// <param name="columnsA">The number of columns in the <paramref name="a"/> matrix.</param>
/// <param name="b">The b matrix</param>
/// <param name="rowsB">The number of rows in the <paramref name="b"/> matrix.</param>
/// <param name="columnsB">The number of columns in the <paramref name="b"/> matrix.</param>
/// <param name="beta">The value to scale the <paramref name="c"/> matrix.</param>
/// <param name="c">The c matrix.</param>
[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);
}
/// <summary>
/// Computes the LUP factorization of A. P*A = L*U.
/// </summary>
/// <param name="data">An <paramref name="order"/> by <paramref name="order"/> matrix. The matrix is overwritten with the
/// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of <paramref name="data"/> (the diagonal is always 1.0
/// for the L factor). The upper triangular factor U is stored on and above the diagonal of <paramref name="data"/>.</param>
/// <param name="order">The order of the square matrix <paramref name="data"/>.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size of the array must be <paramref name="order"/>.</param>
/// <remarks>This is equivalent to the GETRF LAPACK routine.</remarks>
[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));
}
/// <summary>
/// Computes the inverse of matrix using LU factorization.
/// </summary>
/// <param name="a">The N by N matrix to invert. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <remarks>This is equivalent to the GETRF and GETRI LAPACK routines.</remarks>
[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));
}
/// <summary>
/// Computes the inverse of a previously factored matrix.
/// </summary>
/// <param name="a">The LU factored N by N matrix. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <remarks>This is equivalent to the GETRI LAPACK routine.</remarks>
[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));
}
/// <summary>
/// Computes the inverse of matrix using LU factorization.
/// </summary>
/// <param name="a">The N by N matrix to invert. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="work">Not supported. Should be left null.</param>
/// <remarks>This is equivalent to the GETRF and GETRI LAPACK routines.</remarks>
[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));
}
/// <summary>
/// Computes the inverse of a previously factored matrix.
/// </summary>
/// <param name="a">The LU factored N by N matrix. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="work">Not supported. Should be left null.</param>
/// <remarks>This is equivalent to the GETRI LAPACK routine.</remarks>
[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));
}
/// <summary>
/// Solves A*X=B for X using LU factorization.
/// </summary>
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The square matrix A.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <remarks>This is equivalent to the GETRF and GETRS LAPACK routines.</remarks>
[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));
}
/// <summary>
/// Solves A*X=B for X using a previously factored A matrix.
/// </summary>
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The factored A matrix.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <remarks>This is equivalent to the GETRS LAPACK routine.</remarks>
[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));
}
/// <summary>
/// Computes the Cholesky factorization of A.
/// </summary>
/// <param name="a">On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the
/// the Cholesky factorization.</param>
/// <param name="order">The number of rows or columns in the matrix.</param>
/// <remarks>This is equivalent to the POTRF LAPACK routine.</remarks>
[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));
}
/// <summary>
/// Solves A*X=B for X using Cholesky factorization.
/// </summary>
/// <param name="a">The square, positive definite matrix A.</param>
/// <param name="orderA">The number of rows and columns in A.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns in the B matrix.</param>
/// <remarks>This is equivalent to the POTRF add POTRS LAPACK routines.
/// </remarks>
[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));
}
/// <summary>
/// Solves A*X=B for X using a previously factored A matrix.
/// </summary>
/// <param name="a">The square, positive definite matrix A.</param>
/// <param name="orderA">The number of rows and columns in A.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns in the B matrix.</param>
/// <remarks>This is equivalent to the POTRS LAPACK routine.</remarks>
[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));
}
/// <summary>
/// Computes the singular value decomposition of A.
/// </summary>
/// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">If <paramref name="computeVectors"/> is <c>true</c>, on exit U contains the left
/// singular vectors.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
[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);
}
/// <summary>
/// Solves A*X=B for X using the singular value decomposition of A.
/// </summary>
/// <param name="a">On entry, the M by N matrix to decompose.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
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);
}
/// <summary>
/// Computes the singular value decomposition of A.
/// </summary>
/// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">If <paramref name="computeVectors"/> is <c>true</c>, on exit U contains the left
/// singular vectors.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <param name="work">Not supported. Should be left null.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
[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

703
src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Single.cs

@ -0,0 +1,703 @@
// <copyright file="MklLinearAlgebraProvider.Single.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-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.
// </copyright>
#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
{
/// <summary>
/// Intel's Math Kernel Library (MKL) linear algebra provider.
/// </summary>
public partial class CudaLinearAlgebraProvider
{
/// <summary>
/// Computes the dot product of x and y.
/// </summary>
/// <param name="x">The vector x.</param>
/// <param name="y">The vector y.</param>
/// <returns>The dot product of x and y.</returns>
/// <remarks>This is equivalent to the DOT BLAS routine.</remarks>
[SecuritySafeCritical]
public override float DotProduct(float[] x, float[] y)
{
if (y == null)
{
throw new ArgumentNullException("y");
}
if (x == null)
{
throw new ArgumentNullException("x");
}
if (x.Length != y.Length)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength);
}
return SafeNativeMethods.s_dot_product(_blasHandle, x.Length, x, y);
}
/// <summary>
/// Adds a scaled vector to another: <c>result = y + alpha*x</c>.
/// </summary>
/// <param name="y">The vector to update.</param>
/// <param name="alpha">The value to scale <paramref name="x"/> by.</param>
/// <param name="x">The vector to add to <paramref name="y"/>.</param>
/// <param name="result">The result of the addition.</param>
/// <remarks>This is similar to the AXPY BLAS routine.</remarks>
[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);
}
/// <summary>
/// Scales an array. Can be used to scale a vector and a matrix.
/// </summary>
/// <param name="alpha">The scalar.</param>
/// <param name="x">The values to scale.</param>
/// <param name="result">This result of the scaling.</param>
/// <remarks>This is similar to the SCAL BLAS routine.</remarks>
[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);
}
/// <summary>
/// Multiples two matrices. <c>result = x * y</c>
/// </summary>
/// <param name="x">The x matrix.</param>
/// <param name="rowsX">The number of rows in the x matrix.</param>
/// <param name="columnsX">The number of columns in the x matrix.</param>
/// <param name="y">The y matrix.</param>
/// <param name="rowsY">The number of rows in the y matrix.</param>
/// <param name="columnsY">The number of columns in the y matrix.</param>
/// <param name="result">Where to store the result of the multiplication.</param>
/// <remarks>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.</remarks>
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);
}
/// <summary>
/// Multiplies two matrices and updates another with the result. <c>c = alpha*op(a)*op(b) + beta*c</c>
/// </summary>
/// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
/// <param name="transposeB">How to transpose the <paramref name="b"/> matrix.</param>
/// <param name="alpha">The value to scale <paramref name="a"/> matrix.</param>
/// <param name="a">The a matrix.</param>
/// <param name="rowsA">The number of rows in the <paramref name="a"/> matrix.</param>
/// <param name="columnsA">The number of columns in the <paramref name="a"/> matrix.</param>
/// <param name="b">The b matrix</param>
/// <param name="rowsB">The number of rows in the <paramref name="b"/> matrix.</param>
/// <param name="columnsB">The number of columns in the <paramref name="b"/> matrix.</param>
/// <param name="beta">The value to scale the <paramref name="c"/> matrix.</param>
/// <param name="c">The c matrix.</param>
[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);
}
/// <summary>
/// Computes the LUP factorization of A. P*A = L*U.
/// </summary>
/// <param name="data">An <paramref name="order"/> by <paramref name="order"/> matrix. The matrix is overwritten with the
/// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of <paramref name="data"/> (the diagonal is always 1.0f
/// for the L factor). The upper triangular factor U is stored on and above the diagonal of <paramref name="data"/>.</param>
/// <param name="order">The order of the square matrix <paramref name="data"/>.</param>
/// <param name="ipiv">On exit, it contains the pivot indices. The size of the array must be <paramref name="order"/>.</param>
/// <remarks>This is equivalent to the GETRF LAPACK routine.</remarks>
[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));
}
/// <summary>
/// Computes the inverse of matrix using LU factorization.
/// </summary>
/// <param name="a">The N by N matrix to invert. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <remarks>This is equivalent to the GETRF and GETRI LAPACK routines.</remarks>
[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));
}
/// <summary>
/// Computes the inverse of a previously factored matrix.
/// </summary>
/// <param name="a">The LU factored N by N matrix. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <remarks>This is equivalent to the GETRI LAPACK routine.</remarks>
[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));
}
/// <summary>
/// Computes the inverse of matrix using LU factorization.
/// </summary>
/// <param name="a">The N by N matrix to invert. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="work">Not supported. Should be left null.</param>
/// <remarks>This is equivalent to the GETRF and GETRI LAPACK routines.</remarks>
[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));
}
/// <summary>
/// Computes the inverse of a previously factored matrix.
/// </summary>
/// <param name="a">The LU factored N by N matrix. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="work">Not supported. This should be left null.</param>
/// <remarks>This is equivalent to the GETRI LAPACK routine.</remarks>
[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));
}
/// <summary>
/// Solves A*X=B for X using LU factorization.
/// </summary>
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The square matrix A.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <remarks>This is equivalent to the GETRF and GETRS LAPACK routines.</remarks>
[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));
}
/// <summary>
/// Solves A*X=B for X using a previously factored A matrix.
/// </summary>
/// <param name="columnsOfB">The number of columns of B.</param>
/// <param name="a">The factored A matrix.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <remarks>This is equivalent to the GETRS LAPACK routine.</remarks>
[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));
}
/// <summary>
/// Computes the Cholesky factorization of A.
/// </summary>
/// <param name="a">On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the
/// the Cholesky factorization.</param>
/// <param name="order">The number of rows or columns in the matrix.</param>
/// <remarks>This is equivalent to the POTRF LAPACK routine.</remarks>
[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));
}
/// <summary>
/// Solves A*X=B for X using Cholesky factorization.
/// </summary>
/// <param name="a">The square, positive definite matrix A.</param>
/// <param name="orderA">The number of rows and columns in A.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns in the B matrix.</param>
/// <remarks>This is equivalent to the POTRF add POTRS LAPACK routines.
/// </remarks>
[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));
}
/// <summary>
/// Solves A*X=B for X using a previously factored A matrix.
/// </summary>
/// <param name="a">The square, positive definite matrix A.</param>
/// <param name="orderA">The number of rows and columns in A.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns in the B matrix.</param>
/// <remarks>This is equivalent to the POTRS LAPACK routine.</remarks>
[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));
}
/// <summary>
/// Computes the singular value decomposition of A.
/// </summary>
/// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">If <paramref name="computeVectors"/> is <c>true</c>, on exit U contains the left
/// singular vectors.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
[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);
}
/// <summary>
/// Solves A*X=B for X using the singular value decomposition of A.
/// </summary>
/// <param name="a">On entry, the M by N matrix to decompose.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
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);
}
/// <summary>
/// Computes the singular value decomposition of A.
/// </summary>
/// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">If <paramref name="computeVectors"/> is <c>true</c>, on exit U contains the left
/// singular vectors.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <param name="work">Not supported. Should be left null.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
[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

209
src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.cs

@ -0,0 +1,209 @@
// <copyright file="MklLinearAlgebraProvider.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2015 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
using System;
#if NATIVE
namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
{
/// <summary>
/// Consistency vs. performance trade-off between runs on different machines.
/// </summary>
/// <summary>
/// Intel's Math Kernel Library (MKL) linear algebra provider.
/// </summary>
public partial class CudaLinearAlgebraProvider : ManagedLinearAlgebraProvider, IDisposable
{
private int _nativeRevision;
private bool _nativeIX86;
private bool _nativeX64;
private bool _nativeIA64;
private IntPtr _blasHandle;
private IntPtr _solverHandle;
/// <param name="consistency">
/// Sets the desired bit consistency on repeated identical computations on varying CPU architectures,
/// as a trade-off with performance.
/// </param>
/// <param name="precision">VML optimal precision and rounding.</param>
/// <param name="accuracy">VML accuracy mode.</param>
[CLSCompliant(false)]
public CudaLinearAlgebraProvider()
{
}
/// <summary>
/// Initialize and verify that the provided is indeed available.
/// If calling this method fails, consider to fall back to alternatives like the managed provider.
/// </summary>
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

378
src/Numerics/Providers/LinearAlgebra/Cuda/SafeNativeMethods.cs

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

348
src/UnitTests/UnitTests-CUDA.csproj

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

3
src/UnitTests/UseLinearAlgebraProvider.cs

@ -40,6 +40,9 @@ namespace MathNet.Numerics.UnitTests
{
#if !NET35 && NATIVE
Control.UseNativeMKL();
#endif
#if CUDA
Control.UseNativeCUDA();
#endif
}

Loading…
Cancel
Save