diff --git a/src/NativeProviders/CUDA/blas.c b/src/NativeProviders/CUDA/blas.c deleted file mode 100644 index 18ff81aa..00000000 --- a/src/NativeProviders/CUDA/blas.c +++ /dev/null @@ -1,93 +0,0 @@ -#include "cublas_v2.h" -#include "wrapper_common.h" - -#if GCC -extern "C" { -#endif - DLLEXPORT void s_axpy(const cublasHandle_t handle, const int n, const float alpha, const float x[], float y[]){ - cublasSaxpy(handle, n, &alpha, x, 1, y, 1); - } - - DLLEXPORT void d_axpy(const cublasHandle_t handle, const int n, const double alpha, const double x[], double y[]){ - cublasDaxpy(handle, n, &alpha, x, 1, y, 1); - } - - DLLEXPORT void c_axpy(const cublasHandle_t handle, const int n, const cuComplex alpha, const cuComplex x[], cuComplex y[]){ - cublasCaxpy(handle, n, &alpha, x, 1, y, 1); - } - - DLLEXPORT void z_axpy(const cublasHandle_t handle, const int n, const cuDoubleComplex alpha, const cuDoubleComplex x[], cuDoubleComplex y[]){ - cublasZaxpy(handle, n, &alpha, x, 1, y, 1); - } - - DLLEXPORT void s_scale(const cublasHandle_t handle, const int n, const float alpha, float x[]){ - cublasSscal(handle, n, &alpha, x, 1); - } - - DLLEXPORT void d_scale(const cublasHandle_t handle, const int n, const double alpha, double x[]){ - cublasDscal(handle, n, &alpha, x, 1); - } - - DLLEXPORT void c_scale(const cublasHandle_t handle, const int n, const cuComplex alpha, cuComplex x[]){ - cublasCscal(handle, n, &alpha, x, 1); - } - - DLLEXPORT void z_scale(const cublasHandle_t handle, const int n, const cuDoubleComplex alpha, cuDoubleComplex x[]){ - cublasZscal(handle, n, &alpha, x, 1); - } - - DLLEXPORT float s_dot_product(const cublasHandle_t handle, const int n, const float x[], const float y[]){ - float ret; - cublasSdot(handle, n, x, 1, y, 1, &ret); - return ret; - } - - DLLEXPORT double d_dot_product(const cublasHandle_t handle, const int n, const double x[], const double y[]){ - double ret; - cublasDdot(handle, n, x, 1, y, 1, &ret); - return ret; - } - - DLLEXPORT cuComplex c_dot_product(const cublasHandle_t handle, const int n, const cuComplex x[], const cuComplex y[]){ - cuComplex ret; - cublasCdotu(handle, n, x, 1, y, 1, &ret); - return ret; - } - - DLLEXPORT cuDoubleComplex z_dot_product(const cublasHandle_t handle, const int n, const cuDoubleComplex x[], const cuDoubleComplex y[]){ - cuDoubleComplex ret; - cublasZdotu(handle, n, x, 1, y, 1, &ret); - return ret; - } - - DLLEXPORT void s_matrix_multiply(const cublasHandle_t handle, cublasOperation_t transA, cublasOperation_t transB, const int m, const int n, const int k, const float alpha, const float x[], const float y[], const float beta, float c[]){ - int lda = transA == CUBLAS_OP_N ? m : k; - int ldb = transB == CUBLAS_OP_N ? k : n; - - cublasSgemm(handle, transA, transB, m, n, k, &alpha, x, lda, y, ldb, &beta, c, m); - } - - DLLEXPORT void d_matrix_multiply(const cublasHandle_t handle, cublasOperation_t transA, cublasOperation_t transB, const int m, const int n, const int k, const double alpha, const double x[], const double y[], const double beta, double c[]){ - int lda = transA == CUBLAS_OP_N ? m : k; - int ldb = transB == CUBLAS_OP_N ? k : n; - - cublasDgemm(handle, transA, transB, m, n, k, &alpha, x, lda, y, ldb, &beta, c, m); - } - - DLLEXPORT void c_matrix_multiply(const cublasHandle_t handle, cublasOperation_t transA, cublasOperation_t transB, const int m, const int n, const int k, const cuComplex alpha, const cuComplex x[], const cuComplex y[], const cuComplex beta, cuComplex c[]){ - int lda = transA == CUBLAS_OP_N ? m : k; - int ldb = transB == CUBLAS_OP_N ? k : n; - - cublasCgemm(handle, transA, transB, m, n, k, &alpha, x, lda, y, ldb, &beta, c, m); - } - - DLLEXPORT void z_matrix_multiply(const cublasHandle_t handle, cublasOperation_t transA, cublasOperation_t transB, const int m, const int n, const int k, const cuDoubleComplex alpha, const cuDoubleComplex x[], const cuDoubleComplex y[], const cuDoubleComplex beta, cuDoubleComplex c[]){ - int lda = transA == CUBLAS_OP_N ? m : k; - int ldb = transB == CUBLAS_OP_N ? k : n; - - cublasZgemm(handle, transA, transB, m, n, k, &alpha, x, lda, y, ldb, &beta, c, m); - } - -#if GCC -} -#endif diff --git a/src/NativeProviders/CUDA/blas.cpp b/src/NativeProviders/CUDA/blas.cpp new file mode 100644 index 00000000..47233b72 --- /dev/null +++ b/src/NativeProviders/CUDA/blas.cpp @@ -0,0 +1,167 @@ +#include "cublas_v2.h" +#include "cuda_runtime.h" +#include "wrapper_common.h" + +template +void cuda_axpy(const cublasHandle_t blasHandle, const int n, const T *alpha, const T x[], int incX, T y[], int incY, AXPY axpy) +{ + T *d_X = NULL; + T *d_Y = NULL; + cudaMalloc((void**)&d_X, n*sizeof(T)); + cudaMalloc((void**)&d_Y, n*sizeof(T)); + + cublasSetVector(n, sizeof(T), x, incX, d_X, incX); + cublasSetVector(n, sizeof(T), y, incY, d_Y, incY); + + axpy(blasHandle, n, alpha, d_X, incX, d_Y, incX); + + cublasGetVector(n, sizeof(T), d_Y, incY, y, incY); + + cudaFree(d_X); + cudaFree(d_Y); +} + +template +void cuda_scal(const cublasHandle_t blasHandle, const int n, const T *alpha, T x[], int incX, SCAL scal) +{ + T *d_X = NULL; + cudaMalloc((void**)&d_X, n*sizeof(T)); + + cublasSetVector(n, sizeof(T), x, incX, d_X, incX); + + scal(blasHandle, n, alpha, d_X, incX); + + cublasGetVector(n, sizeof(T), d_X, incX, x, incX); + + cudaFree(d_X); +} + +template +void cuda_dot(const cublasHandle_t blasHandle, const int n, const T x[], int incX, const T y[], int incY, T* result, DOT dot) +{ + T *d_X = NULL; + T *d_Y = NULL; + cudaMalloc((void**)&d_X, n*sizeof(T)); + cudaMalloc((void**)&d_Y, n*sizeof(T)); + + cublasSetVector(n, sizeof(T), x, incX, d_X, incX); + cublasSetVector(n, sizeof(T), y, incY, d_Y, incY); + + dot(blasHandle, n, d_X, incX, d_Y, incY, result); + + cudaFree(d_X); + cudaFree(d_Y); +} + +template +void cuda_gemm(const cublasHandle_t handle, const cublasOperation_t transa, const cublasOperation_t transb, int m, int n, int k, const T *alpha, const T A[], int lda, const T B[], int ldb, const T *beta, T C[], int ldc, GEMM gemm) +{ + T *d_A = NULL; + T *d_B = NULL; + T *d_C = NULL; + cudaMalloc((void**)&d_A, m*k*sizeof(T)); + cudaMalloc((void**)&d_B, k*n*sizeof(T)); + cudaMalloc((void**)&d_C, m*n*sizeof(T)); + + cublasSetMatrix(m, k, sizeof(T), A, m, d_A, m); + cublasSetMatrix(k, n, sizeof(T), B, k, d_B, k); + + gemm(handle, transa, transb, m, n, k, alpha, d_A, lda, d_B, ldb, beta, d_C, ldc); + + cublasGetMatrix(m, n, sizeof(T), d_C, m, C, m); + + cudaFree(d_A); + cudaFree(d_B); + 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); + } + + DLLEXPORT void d_axpy(const cublasHandle_t blasHandle, const int n, const double alpha, const double x[], double y[]){ + cuda_axpy(blasHandle, n, &alpha, x, 1, y, 1, cublasDaxpy); + } + + DLLEXPORT void c_axpy(const cublasHandle_t blasHandle, const int n, const cuComplex alpha, const cuComplex x[], cuComplex y[]){ + cuda_axpy(blasHandle, n, &alpha, x, 1, y, 1, cublasCaxpy); + } + + DLLEXPORT void z_axpy(const cublasHandle_t blasHandle, const int n, const cuDoubleComplex alpha, const cuDoubleComplex x[], cuDoubleComplex y[]){ + cuda_axpy(blasHandle, n, &alpha, x, 1, y, 1, cublasZaxpy); + } + + DLLEXPORT void s_scale(const cublasHandle_t blasHandle, const int n, const float alpha, float x[]){ + cuda_scal(blasHandle, n, &alpha, x, 1, cublasSscal); + } + + DLLEXPORT void d_scale(const cublasHandle_t blasHandle, const int n, const double alpha, double x[]){ + cuda_scal(blasHandle, n, &alpha, x, 1, cublasDscal); + } + + DLLEXPORT void c_scale(const cublasHandle_t blasHandle, const int n, const cuComplex alpha, cuComplex x[]){ + cuda_scal(blasHandle, n, &alpha, x, 1, cublasCscal); + } + + DLLEXPORT void z_scale(const cublasHandle_t blasHandle, const int n, const cuDoubleComplex alpha, cuDoubleComplex x[]){ + cuda_scal(blasHandle, n, &alpha, x, 1, cublasZscal); + } + + DLLEXPORT float s_dot_product(const cublasHandle_t blasHandle, const int n, const float x[], const float y[]){ + float ret; + cuda_dot(blasHandle, n, x, 1, y, 1, &ret, cublasSdot); + return ret; + } + + DLLEXPORT double d_dot_product(const cublasHandle_t blasHandle, const int n, const double x[], const double y[]){ + double ret; + cuda_dot(blasHandle, n, x, 1, y, 1, &ret, cublasDdot); + return ret; + } + + DLLEXPORT cuComplex c_dot_product(const cublasHandle_t blasHandle, const int n, const cuComplex x[], const cuComplex y[]){ + cuComplex ret; + cuda_dot(blasHandle, n, x, 1, y, 1, &ret, cublasCdotu); + return ret; + } + + DLLEXPORT cuDoubleComplex z_dot_product(const cublasHandle_t blasHandle, const int n, const cuDoubleComplex x[], const cuDoubleComplex y[]){ + cuDoubleComplex ret; + cuda_dot(blasHandle, n, x, 1, y, 1, &ret, cublasZdotu); + return ret; + } + + DLLEXPORT void s_matrix_multiply(const cublasHandle_t blasHandle, cublasOperation_t transA, cublasOperation_t transB, const int m, const int n, const int k, const float alpha, const float x[], const float y[], const float beta, float c[]){ + int lda = transA == CUBLAS_OP_N ? m : k; + int ldb = transB == CUBLAS_OP_N ? k : n; + + cuda_gemm(blasHandle, transA, transB, m, n, k, &alpha, x, lda, y, ldb, &beta, c, m, cublasSgemm); + } + + DLLEXPORT void d_matrix_multiply(const cublasHandle_t blasHandle, cublasOperation_t transA, cublasOperation_t transB, const int m, const int n, const int k, const double alpha, const double x[], const double y[], const double beta, double c[]){ + int lda = transA == CUBLAS_OP_N ? m : k; + int ldb = transB == CUBLAS_OP_N ? k : n; + + cuda_gemm(blasHandle, transA, transB, m, n, k, &alpha, x, lda, y, ldb, &beta, c, m, cublasDgemm); + } + + DLLEXPORT void c_matrix_multiply(const cublasHandle_t blasHandle, cublasOperation_t transA, cublasOperation_t transB, const int m, const int n, const int k, const cuComplex alpha, const cuComplex x[], const cuComplex y[], const cuComplex beta, cuComplex c[]){ + int lda = transA == CUBLAS_OP_N ? m : k; + int ldb = transB == CUBLAS_OP_N ? k : n; + + cuda_gemm(blasHandle, transA, transB, m, n, k, &alpha, x, lda, y, ldb, &beta, c, m, cublasCgemm); + } + + DLLEXPORT void z_matrix_multiply(const cublasHandle_t blasHandle, cublasOperation_t transA, cublasOperation_t transB, const int m, const int n, const int k, const cuDoubleComplex alpha, const cuDoubleComplex x[], const cuDoubleComplex y[], const cuDoubleComplex beta, cuDoubleComplex c[]){ + int lda = transA == CUBLAS_OP_N ? m : k; + int ldb = transB == CUBLAS_OP_N ? k : n; + + cuda_gemm(blasHandle, transA, transB, m, n, k, &alpha, x, lda, y, ldb, &beta, c, m, cublasZgemm); + } + +#if GCC +} +#endif diff --git a/src/NativeProviders/CUDA/lapack.cpp b/src/NativeProviders/CUDA/lapack.cpp index 189a0b4c..7798df12 100644 --- a/src/NativeProviders/CUDA/lapack.cpp +++ b/src/NativeProviders/CUDA/lapack.cpp @@ -1,521 +1,976 @@ +#include + #include "lapack_common.h" #include "wrapper_common.h" -#include "cublas.h" +#include "cublas_v2.h" #include "cusolverDn.h" -#include +#include "cuda_runtime.h" -template -inline int lu_factor(int m, T a[], int ipiv[], - int(*getrf) (CBLAS_ORDER, const int, const int, K*, const int, int*)) +template +inline int lu_factor(cusolverDnHandle_t solverHandle, int m, T a[], int ipiv[], GETRF getrf, GETRFBSIZE getrfbsize) { - int info = getrf(CblasColMajor, m, m, a, m, ipiv); + int info = 0; + T* work = NULL; + int lwork = 0; + + T* d_A = NULL; + cudaMalloc((void**)&d_A, m*m*sizeof(T)); + cublasSetMatrix(m, m, sizeof(T), a, m, d_A, m); + + int* d_I = NULL; + cudaMalloc((void**)&d_I, m*sizeof(int)); + + getrfbsize(solverHandle, m, m, a, m, &lwork); + cudaMalloc((void**)lwork, sizeof(T)*lwork); + + getrf(solverHandle, m, m, d_A, m, work, d_I, &info); + + cublasGetMatrix(m, m, sizeof(T), d_A, m, a, m); + cublasGetVector(m, sizeof(T), d_I, 1, ipiv, 1); + shift_ipiv_down(m, ipiv); + + cudaFree(d_A); + cudaFree(d_I); + cudaFree(work); + return info; }; -template -inline int lu_inverse(int n, T a[], - int(*getrf) (CBLAS_ORDER, const int, const int, K*, const int, int*), - int(*getri) (CBLAS_ORDER, const int, K*, const int, const int*)) +template +inline int lu_inverse(cusolverDnHandle_t solverHandle, cublasHandle_t blasHandle, int n, T a[], GETRF getrf, GETRI getri, GETRFBSIZE getrfbsize) { - int* ipiv = new int[n]; - int info = getrf(CblasColMajor, n, n, a, n, ipiv); + int info = 0; + T* work = NULL; + int lwork = 0; + + int* d_I = NULL; + cudaMalloc((void**)&d_I, n*sizeof(T)); + + T* d_A = NULL; + cudaMalloc((void**)&d_A, n*n*sizeof(T)); + cublasSetMatrix(n, n, sizeof(T), a, n, d_A, n); - if (info != 0){ - delete[] ipiv; + getrfbsize(solverHandle, n, n, d_A, n, &lwork); + cudaMalloc((void**)lwork, sizeof(T)*lwork); + + getrf(solverHandle, n, n, d_A, n, work, d_I, &info); + + cudaFree(work); + + if (info != 0) + { + cudaFree(d_A); + cudaFree(d_I); return info; } - info = getri(CblasColMajor, n, a, n, ipiv); - delete[] ipiv; + T* d_C = NULL; + cudaMalloc((void**)&d_C, n*n*sizeof(T)); + + getri(blasHandle, n, d_A, n, d_I, d_C, n, &info); + + cublasGetMatrix(n, n, sizeof(T), d_C, n, a, n); + + cudaFree(d_A); + cudaFree(d_I); + cudaFree(d_C); + return info; }; -template -inline int lu_inverse_factored(int n, T a[], int ipiv[], - int(*getri) (CBLAS_ORDER, const int, K*, const int, const int*)) +template +inline int lu_inverse_factored(cublasHandle_t blasHandle, int n, T a[], int ipiv[], GETRI getri) { shift_ipiv_up(n, ipiv); - int info = getri(CblasColMajor, n, a, n, ipiv); + int info = 0; + + T* d_A = NULL; + cudaMalloc((void**)&d_A, n*n*sizeof(T)); + cublasSetMatrix(n, n, sizeof(T), a, n, d_A, n); + + T* d_C = NULL; + cudaMalloc((void**)&d_C, n*n*sizeof(T)); + + int* d_I = NULL; + 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); + + cublasGetMatrix(n, n, sizeof(T), d_C, n, a, n); + cublasGetVector(n, sizeof(int), d_I, 1, ipiv, 1); + shift_ipiv_down(n, ipiv); + + cudaFree(d_A); + cudaFree(d_I); + cudaFree(d_C); + return info; } -template -inline int lu_solve_factored(int n, int nrhs, T a[], int ipiv[], T b[], - int(*getrs) (CBLAS_ORDER, CBLAS_TRANSPOSE, const int, const int, const K*, const int, const int*, K*, const int)) +template +inline int lu_solve_factored(cusolverDnHandle_t solverHandle, int n, int nrhs, T a[], int ipiv[], T b[], GETRS getrs) { shift_ipiv_up(n, ipiv); - int info = getrs(CblasColMajor, CblasNoTrans, n, nrhs, a, n, ipiv, b, n); + int info = 0; + + T* d_A = NULL; + cudaMalloc((void**)&d_A, n*n*sizeof(T)); + cublasSetMatrix(n, n, sizeof(T), a, n, d_A, n); + + T* d_B = NULL; + cudaMalloc((void**)&d_B, n*nrhs*sizeof(T)); + cublasSetMatrix(n, nrhs, sizeof(T), b, n, d_B, n); + + int* d_I = NULL; + 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); + + cublasGetMatrix(n, nrhs, sizeof(T), d_B, n, b, n); + shift_ipiv_down(n, ipiv); + + cudaFree(d_A); + cudaFree(d_B); + cudaFree(d_I); + return info; } -template -inline int lu_solve(int n, int nrhs, T a[], T b[], - int(*getrf) (CBLAS_ORDER, const int, const int, K*, const int, int*), - int(*getrs) (CBLAS_ORDER, CBLAS_TRANSPOSE, const int, const int, const K*, const int, const int*, K*, const int)) +template +inline int lu_solve(cusolverDnHandle_t solverHandle, int n, int nrhs, T a[], T b[], GETRF getrf, GETRS getrs, GETRFBSIZE getrfbsize) { - T* clone = Clone(n, n, a); - int* ipiv = new int[n]; - int info = getrf(CblasColMajor, n, n, clone, n, ipiv); + int info = 0; + T* work = NULL; + int lwork = 0; + + int* d_I = NULL; + cudaMalloc((void**)&d_I, n*sizeof(T)); + + T* d_A = NULL; + cudaMalloc((void**)&d_A, n*n*sizeof(T)); + cublasSetMatrix(n, n, sizeof(T), a, n, d_A, n); + + getrfbsize(solverHandle, n, n, a, n, &lwork); + cudaMalloc((void**)lwork, sizeof(T)*lwork); + + getrf(solverHandle, n, n, d_A, n, work, d_I, &info); - if (info != 0){ - delete[] ipiv; - delete[] clone; + if (info != 0) + { + cudaFree(d_I); + cudaFree(d_A); return info; } - info = getrs(CblasColMajor, CblasNoTrans, n, nrhs, clone, n, ipiv, b, n); - delete[] ipiv; - delete[] clone; + T* d_B = NULL; + 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); + + cublasGetMatrix(n, nrhs, sizeof(T), d_B, n, b, n); + + cudaFree(d_A); + cudaFree(d_B); + cudaFree(d_I); + return info; } -template -inline int cholesky_factor(int n, T* a, int(*potrf) (CBLAS_ORDER, CBLAS_UPLO, const int, K*, const int)) + +template +inline int cholesky_factor(cusolverDnHandle_t solverHandle, int n, T* a, POTRF potrf, POTRFBSIZE potrfbsize) { - int info = potrf(CblasColMajor, CblasLower, n, a, n); + int info = 0; + + T* d_A = NULL; + cudaMalloc((void**)&d_A, n*n*sizeof(T)); + cublasSetMatrix(n, n, sizeof(T), a, n, d_A, n); + + T* work = NULL; + int lWork = 0; + 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); + + cublasGetMatrix(n, n, sizeof(T), d_A, n, a, n); + T zero = T(); + for (int i = 0; i < n; ++i) { int index = i * n; + for (int j = 0; j < n && i > j; ++j) { a[index + j] = zero; } } - return info; -} - -template -inline int cholesky_solve(int n, int nrhs, T a[], T b[], - int(*potrf) (CBLAS_ORDER, CBLAS_UPLO, const int, K*, const int), - int(*potrs) (CBLAS_ORDER, CBLAS_UPLO, const int, const int, const K*, const int, K*, const int)) -{ - T* clone = Clone(n, n, a); - int info = potrf(CblasColMajor, CblasLower, n, clone, n); - if (info != 0){ - delete[] clone; - return info; - } + cudaFree(d_A); + cudaFree(work); - info = potrs(CblasColMajor, CblasLower, n, nrhs, clone, n, b, n); - delete[] clone; return info; } -template -inline int cholesky_solve_factored(int n, int nrhs, T a[], T b[], - int(*potrs) (CBLAS_ORDER, CBLAS_UPLO, const int, const int, const K*, const int, K*, const int)) +template +inline int cholesky_solve(cusolverDnHandle_t solverHandle, int n, int nrhs, T a[], T b[], POTRF potrf, POTRS potrs, POTRFBSIZE potrfbsize) { - return potrs(CblasColMajor, CblasLower, n, nrhs, a, n, b, n); -} + int info; -template -inline int qr_factor(int m, int n, T r[], T tau[], T q[], T work[], int len, - int(*geqrf) (const int, const int, K*, const int, T*), - int(*orgqr) (const int, const int, const int, K*, const int, const K*)) -{ - int info = geqrf(m, n, r, m, tau); + T* d_A = NULL; + cudaMalloc((void**)&d_A, n*n*sizeof(T)); + cublasSetMatrix(n, n, sizeof(T), a, n, d_A, n); - for (int i = 0; i < m; ++i) - { - for (int j = 0; j < m && j < n; ++j) - { - if (i > j) - { - q[j * m + i] = r[j * m + i]; - } - } - } + T* work = NULL; + int lWork = 0; + potrfbsize(solverHandle, CUBLAS_FILL_MODE_LOWER, n, d_A, n, &lWork); + cudaMalloc((void**)&work, sizeof(T)*lWork); - //compute the q elements explicitly - if (m <= n) - { - info = orgqr(m, m, m, q, m, tau); - } - else + potrf(solverHandle, CUBLAS_FILL_MODE_LOWER, n, d_A, n, work, lWork, &info); + + cudaFree(work); + + if (info != 0) { - info = orgqr(m, m, n, q, m, tau); + cudaFree(d_A); + return info; } - return info; -} + T* d_B = NULL; + cudaMalloc((void**)d_B, n*nrhs*sizeof(T)); + cublasSetMatrix(n, nrhs, sizeof(T), b, n, d_B, n); -template -inline int qr_thin_factor(int m, int n, T q[], T tau[], T r[], T work[], int len, - void(*geqrf) (const int*, const int*, T*, const int*, T*, T*, const int*, int*), - void(*orgqr) (const int*, const int*, const int*, T*, const int*, const T*, T*, const int*, int*)) -{ - int info = 0; - geqrf(&m, &n, q, &m, tau, work, &len, &info); + potrs(solverHandle, CUBLAS_FILL_MODE_LOWER, n, nrhs, d_A, n, d_B, n, &info); - for (int i = 0; i < n; ++i) - { - for (int j = 0; j < n; ++j) - { - if (i <= j) { - r[j * n + i] = q[j * m + i]; - } - } - } + cublasGetMatrix(n, nrhs, sizeof(T), d_B, n, b, n); - orgqr(&m, &n, &n, q, &m, tau, work, &len, &info); + cudaFree(d_A); + cudaFree(d_B); return info; } -template -inline int qr_solve(int m, int n, int bn, T a[], T b[], T x[], T work[], int len, - void(*gels) (const char*, const int*, const int*, const int*, T*, - const int*, T* b, const int*, T*, const int*, int*)) +template +inline int cholesky_solve_factored(cusolverDnHandle_t solverHandle, int n, int nrhs, T a[], T b[], POTRS potrs) { - T* clone_a = new T[m*n]; - std::memcpy(clone_a, a, m*n*sizeof(T)); + int info; - T* clone_b = new T[m*bn]; - std::memcpy(clone_b, b, m*bn*sizeof(T)); + T* d_A = NULL; + cudaMalloc((void**)&d_A, n*n*sizeof(T)); + cublasSetMatrix(n, n, sizeof(T), a, n, d_A, n); - char N = 'N'; - int info = 0; - gels(&N, &m, &n, &bn, clone_a, &m, clone_b, &m, work, &len, &info); - copyBtoX(n, n, bn, clone_b, x); + T* d_B = NULL; + cudaMalloc((void**)d_B, n*nrhs*sizeof(T)); + cublasSetMatrix(n, nrhs, sizeof(T), b, n, d_B, n); - delete[] clone_a; - delete[] clone_b; - return info; -} + potrs(solverHandle, CUBLAS_FILL_MODE_LOWER, n, nrhs, d_A, n, d_B, n, &info); -template -inline int qr_solve_factored(int m, int n, int bn, T r[], T b[], T tau[], T x[], T work[], int len, - void(*ormqr) (const char*, const char*, const int*, const int*, const int*, - const T*, const int*, const T*, T*, const int*, T*, const int*, int* info), - void(*trsm) (const CBLAS_ORDER, const CBLAS_SIDE, const CBLAS_UPLO, const CBLAS_TRANSPOSE, const CBLAS_DIAG, - const int, const int, const T, const T*, const int, T*, const int)) -{ - T* clone_b = new T[m*bn]; - std::memcpy(clone_b, b, m*bn*sizeof(T)); + cublasGetMatrix(n, nrhs, sizeof(T), d_B, n, b, n); - char side = 'L'; - char tran = 'T'; - int info = 0; - ormqr(&side, &tran, &m, &bn, &n, r, &m, tau, clone_b, &m, work, &len, &info); - trsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, 1.0, r, m, clone_b, m); - copyBtoX(n, n, bn, clone_b, x); + cudaFree(d_A); + cudaFree(d_B); - delete[] clone_b; return info; } -template -inline int complex_qr_solve_factored(int m, int n, int bn, T r[], T b[], T tau[], T x[], T work[], int len, - void(*unmqr) (const char*, const char*, const int*, const int*, const int*, - const T*, const int*, const T*, T*, const int*, T*, const int*, int* info), - void(*trsm) (const CBLAS_ORDER, const CBLAS_SIDE, const CBLAS_UPLO, const CBLAS_TRANSPOSE, const CBLAS_DIAG, - const int, const int, const void*, const void*, const int, void*, const int ldb)) +//template +//inline int qr_factor(int m, int n, T r[], T tau[], T q[], T work[], int len, GEQRF geqrf, ORGQR orgqr) +//{ +// int info = 0; +// geqrf(&m, &n, r, &m, tau, work, &len, &info); +// +// for (int i = 0; i < m; ++i) +// { +// for (int j = 0; j < m && j < n; ++j) +// { +// if (i > j) +// { +// q[j * m + i] = r[j * m + i]; +// } +// } +// } +// +// //compute the q elements explicitly +// if (m <= n) +// { +// orgqr(&m, &m, &m, q, &m, tau, work, &len, &info); +// } +// else +// { +// orgqr(&m, &m, &n, q, &m, tau, work, &len, &info); +// } +// +// return info; +//} +// +//template +//inline int qr_thin_factor(int m, int n, T q[], T tau[], T r[], T work[], int len, GEQRF geqrf, ORGQR orgqr) +//{ +// int info = 0; +// geqrf(&m, &n, q, &m, tau, work, &len, &info); +// +// for (int i = 0; i < n; ++i) +// { +// for (int j = 0; j < n; ++j) +// { +// if (i <= j) +// { +// r[j * n + i] = q[j * m + i]; +// } +// } +// } +// +// orgqr(&m, &n, &n, q, &m, tau, work, &len, &info); +// return info; +//} +// +//template +//inline int qr_solve(int m, int n, int bn, T a[], T b[], T x[], T work[], int len, GELS gels) +//{ +// T* clone_a = Clone(m, n, a); +// T* clone_b = Clone(m, bn, b); +// char N = 'N'; +// int info = 0; +// gels(&N, &m, &n, &bn, clone_a, &m, clone_b, &m, work, &len, &info); +// copyBtoX(m, n, bn, clone_b, x); +// delete[] clone_a; +// delete[] clone_b; +// return info; +//} + +//template +//inline int qr_solve_factored(cusolverDnHandle_t solverHandle, cublasHandle_t blasHandle, int m, int n, int bn, T r[], T b[], T tau[], T x[], T work[], int len, ORMQR ormqr, TRSM trsm) +//{ +// T* clone_b = Clone(m, bn, b); +// char side = 'L'; +// char tran = 'T'; +// int info = 0; +// ormqr(solverHandle, &side, &tran, &m, &bn, &n, r, &m, tau, clone_b, &m, work, &len, &info); +// trsm(blasHandle, CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, 1.0, r, m, clone_b, m); +// +// copyBtoX(m, n, bn, clone_b, x); +// delete[] clone_b; +// return info; +//} + +//template +//inline int complex_qr_solve_factored(int m, int n, int bn, T r[], T b[], T tau[], T x[], T work[], int len, UNMQR unmqr, TRSM trsm) +//{ +// T* clone_b = Clone(m, bn, b); +// char side = 'L'; +// char tran = 'C'; +// int info = 0; +// unmqr(&side, &tran, &m, &bn, &n, r, &m, tau, clone_b, &m, work, &len, &info); +// T one = 1.0f; +// trsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, &one, r, m, clone_b, m); +// copyBtoX(m, n, bn, clone_b, x); +// delete[] clone_b; +// return info; +//} + +template +inline int svd_factor(cusolverDnHandle_t solverHandle, bool compute_vectors, int m, int n, T a[], T s[], T u[], T v[], GESVD gesvd, GESVDBSIZE gesvdbsize) { - T* clone_b = new T[m*bn]; - std::memcpy(clone_b, b, m*bn*sizeof(T)); - - char side = 'L'; - char tran = 'C'; int info = 0; - unmqr(&side, &tran, &m, &bn, &n, r, &m, tau, clone_b, &m, work, &len, &info); + int dim_s = std::min(m, n); - T one = { 1.0f, 0.0f }; - trsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, &one, r, m, clone_b, m); - copyBtoX(n, n, bn, clone_b, x); + T* d_A = NULL; + cudaMalloc((void**)&d_A, m*n*sizeof(T)); + cublasSetMatrix(m, n, sizeof(T), a, m, d_A, m); - delete[] clone_b; - return info; -} + T* d_S = NULL; + cudaMalloc((void**)&d_S, dim_s*sizeof(T)); + + T* d_U = NULL; + cudaMalloc((void**)&d_U, m*m*sizeof(T)); + + T* d_V = NULL; + cudaMalloc((void**)&d_V, n*m*sizeof(T)); + + T* work = NULL; + int lWork = 0; + gesvdbsize(solverHandle, m, n, &lWork); + cudaMalloc((void**)&work, lWork*sizeof(T)); + + T* rwork = NULL; + cudaMalloc((void**)&rwork, 5 * dim_s * sizeof(T)); -template -inline int svd_factor(bool compute_vectors, int m, int n, T a[], T s[], T u[], T v[], T work[], int len, - void(*gesvd) (const char*, const char*, const int*, const int*, T*, const int*, - T*, T*, const int*, T*, const int*, T*, const int*, int*)) -{ - int info = 0; char job = compute_vectors ? 'A' : 'N'; - gesvd(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, &info); + gesvd(solverHandle, job, job, m, n, d_A, m, d_S, d_U, m, d_V, n, work, lWork, rwork, &info); + + cublasGetVector(dim_s, sizeof(T), d_S, 1, s, 1); + cublasGetMatrix(m, m, sizeof(T), d_U, m, u, m); + cublasGetMatrix(n, n, sizeof(T), d_V, n, v, n); + + cudaFree(d_A); + cudaFree(d_S); + cudaFree(d_U); + cudaFree(d_V); + cudaFree(work); + cudaFree(rwork); + return info; } - -template -inline int complex_svd_factor(bool compute_vectors, int m, int n, T a[], T s[], T u[], T v[], T work[], int len, - void(*gesvd) (const char*, const char*, const int*, const int*, T*, const int*, - R*, T*, const int*, T*, const int*, T*, const int*, R*, int*)) +template +inline int complex_svd_factor(cusolverDnHandle_t solverHandle, bool compute_vectors, int m, int n, T a[], T s[], T u[], T v[], GESVD gesvd, GESVDBSIZE gesvdbsize) { int info = 0; int dim_s = std::min(m, n); - R* rwork = new R[5 * dim_s]; + + T* d_A = NULL; + cudaMalloc((void**)&d_A, m*n*sizeof(T)); + cublasSetMatrix(m, n, sizeof(T), a, m, d_A, m); + R* s_local = new R[dim_s]; + R* d_S = NULL; + cudaMalloc((void**)&d_S, dim_s*sizeof(R)); + + T* d_U = NULL; + cudaMalloc((void**)&d_U, m*m*sizeof(T)); + + T* d_V = NULL; + cudaMalloc((void**)&d_V, n*m*sizeof(T)); + + T* work = NULL; + int lWork = 0; + gesvdbsize(solverHandle, m, n, &lWork); + cudaMalloc((void**)&work, lWork*sizeof(T)); + + R* rwork = NULL; + cudaMalloc((void**)&rwork, 5 * dim_s * sizeof(R)); + char job = compute_vectors ? 'A' : 'N'; - gesvd(&job, &job, &m, &n, a, &m, s_local, u, &m, v, &n, work, &len, rwork, &info); + gesvd(solverHandle, job, job, m, n, d_A, m, d_S, d_U, m, d_V, n, work, lWork, rwork, &info); + + cublasGetVector(dim_s, sizeof(T), d_S, 1, s_local, 1); + cublasGetMatrix(m, m, sizeof(T), d_U, m, u, m); + cublasGetMatrix(n, n, sizeof(T), d_V, n, v, n); - for (int index = 0; index < dim_s; ++index){ - T value = { s_local[index], 0.0f }; - s[index] = value; + for (int index = 0; index < dim_s; ++index) + { + s[index].x = s_local[index]; } - delete[] rwork; delete[] s_local; + cudaFree(d_A); + cudaFree(d_S); + cudaFree(d_U); + cudaFree(d_V); + cudaFree(work); + cudaFree(rwork); + return info; } +//template +//inline int eigen_factor(int n, T a[], T vectors[], R values[], T d[], GEES gees, TREVC trevc) +//{ +// T* clone_a = Clone(n, n, a); +// T* wr = new T[n]; +// T* wi = new T[n]; +// +// int sdim; +// int info = gees(LAPACK_COL_MAJOR, 'V', 'N', nullptr, n, clone_a, n, &sdim, wr, wi, vectors, n); +// if (info != 0) +// { +// delete[] clone_a; +// delete[] wr; +// delete[] wi; +// return info; +// } +// +// int m; +// info = trevc(LAPACK_COL_MAJOR, 'R', 'B', nullptr, n, clone_a, n, nullptr, n, vectors, n, n, &m); +// if (info != 0) +// { +// delete[] clone_a; +// delete[] wr; +// delete[] wi; +// return info; +// } +// +// for (int index = 0; index < n; ++index) +// { +// values[index] = R(wr[index], wi[index]); +// } +// +// for (int i = 0; i < n; ++i) +// { +// int in = i * n; +// d[in + i] = wr[i]; +// +// if (wi[i] > 0) +// { +// d[in + n + i] = wi[i]; +// } +// else if (wi[i] < 0) +// { +// d[in - n + i] = wi[i]; +// } +// } +// +// delete[] clone_a; +// delete[] wr; +// delete[] wi; +// return info; +//} +// +//template +//inline int eigen_complex_factor(int n, T a[], T vectors[], cuDoubleComplex values[], T d[], GEES gees, TREVC trevc) +//{ +// T* clone_a = Clone(n, n, a); +// T* w = new T[n]; +// +// int sdim; +// int info = gees(LAPACK_COL_MAJOR, 'V', 'N', nullptr, n, clone_a, n, &sdim, w, vectors, n); +// if (info != 0) +// { +// delete[] clone_a; +// delete[] w; +// return info; +// } +// +// int m; +// info = trevc(LAPACK_COL_MAJOR, 'R', 'B', nullptr, n, clone_a, n, nullptr, n, vectors, n, n, &m); +// if (info != 0) +// { +// delete[] clone_a; +// delete[] w; +// return info; +// } +// +// for (int i = 0; i < n; ++i) +// { +// values[i] = w[i]; +// d[i * n + i] = w[i]; +// } +// +// delete[] clone_a; +// delete[] w; +// return info; +//} +// +//template +//inline int sym_eigen_factor(int n, T a[], T vectors[], cuDoubleComplex values[], T d[], SYEV syev) +//{ +// T* clone_a = Clone(n, n, a); +// R* w = new R[n]; +// +// int info = syev(LAPACK_COL_MAJOR, 'V', 'U', n, clone_a, n, w); +// if (info != 0) +// { +// delete[] clone_a; +// delete[] w; +// return info; +// } +// +// memcpy(vectors, clone_a, n*n*sizeof(T)); +// +// for (int index = 0; index < n; ++index) +// { +// values[index] = cuDoubleComplex(w[index]); +// } +// +// for (int j = 0; j < n; ++j) +// { +// int jn = j*n; +// +// for (int i = 0; i < n; ++i) +// { +// if (i == j) +// { +// d[jn + i] = w[i]; +// } +// } +// } +// +// delete[] clone_a; +// delete[] w; +// return info; +//} + +#define sgetrf cusolverDnSgetrf +#define dgetrf cusolverDnDgetrf +#define cgetrf cusolverDnCgetrf +#define zgetrf cusolverDnZgetrf +#define sgetrfbsize cusolverDnSgetrf_bufferSize +#define dgetrfbsize cusolverDnDgetrf_bufferSize +#define cgetrfbsize cusolverDnCgetrf_bufferSize +#define zgetrfbsize cusolverDnZgetrf_bufferSize + +#define sgetrs cusolverDnSgetrs +#define dgetrs cusolverDnDgetrs +#define cgetrs cusolverDnCgetrs +#define zgetrs cusolverDnZgetrs + +#define spotrf cusolverDnSpotrf +#define dpotrf cusolverDnDpotrf +#define cpotrf cusolverDnCpotrf +#define zpotrf cusolverDnZpotrf +#define spotrfbsize cusolverDnSpotrf_bufferSize +#define dpotrfbsize cusolverDnDpotrf_bufferSize +#define cpotrfbsize cusolverDnCpotrf_bufferSize +#define zpotrfbsize cusolverDnZpotrf_bufferSize + +#define spotrs cusolverDnSpotrs +#define dpotrs cusolverDnDpotrs +#define cpotrs cusolverDnCpotrs +#define zpotrs cusolverDnZpotrs + +#define sgeqrf cusolverDnSgeqrf +#define dgeqrf cusolverDnDgeqrf +#define cgeqrf cusolverDnCgeqrf +#define zgeqrf cusolverDnZgeqrf + +#define sormqr cusolverDnSormqr +#define dormqr cusolverDnDormqr + +#define sgesvd cusolverDnSgesvd +#define dgesvd cusolverDnDgesvd +#define cgesvd cusolverDnCgesvd +#define zgesvd cusolverDnZgesvd +#define sgesvdbsize cusolverDnSgesvd_bufferSize +#define dgesvdbsize cusolverDnDgesvd_bufferSize +#define cgesvdbsize cusolverDnCgesvd_bufferSize +#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) +{ + 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) +{ + 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) +{ + 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) +{ + return cublasZgetriBatched(handle, n, &a, lda, ipiv, &c, ldc, info, 1); +} + extern "C" { - DLLEXPORT int s_lu_factor(int m, float a[], int ipiv[]) { - return lu_factor(m, a, ipiv, cusolverDnSgetrf); + + DLLEXPORT int s_lu_factor(cusolverDnHandle_t solverHandle, int m, float a[], int ipiv[]) + { + return lu_factor(solverHandle, m, a, ipiv, sgetrf, sgetrfbsize); } - DLLEXPORT int d_lu_factor(int m, double a[], int ipiv[]) { - return lu_factor(m, a, ipiv, cusolverDnDgetrf); + DLLEXPORT int d_lu_factor(cusolverDnHandle_t solverHandle, int m, double a[], int ipiv[]) + { + return lu_factor(solverHandle, m, a, ipiv, dgetrf, dgetrfbsize); } - DLLEXPORT int c_lu_factor(int m, cuComplex a[], int ipiv[]) { - return lu_factor(m, a, ipiv, cusolverDnCgetrf); + DLLEXPORT int c_lu_factor(cusolverDnHandle_t solverHandle, int m, cuComplex a[], int ipiv[]) + { + return lu_factor(solverHandle, m, a, ipiv, cgetrf, cgetrfbsize); } - DLLEXPORT int z_lu_factor(int m, cuDoubleComplex a[], int ipiv[]) { - return lu_factor(m, a, ipiv, cusolverDnZgetrf); + DLLEXPORT int z_lu_factor(cusolverDnHandle_t solverHandle, int m, cuDoubleComplex a[], int ipiv[]) + { + return lu_factor(solverHandle, m, a, ipiv, zgetrf, zgetrfbsize); } - DLLEXPORT int s_lu_inverse(int n, float a[]) + DLLEXPORT int s_lu_inverse(cusolverDnHandle_t solverHandle, cublasHandle_t blasHandle, int n, float a[]) { - return lu_inverse(n, a, cusolverDnSgetrf, cusolverDnSgetri); + return lu_inverse(solverHandle, blasHandle, n, a, sgetrf, sgetri, sgetrfbsize); } - DLLEXPORT int d_lu_inverse(int n, double a[]) + DLLEXPORT int d_lu_inverse(cusolverDnHandle_t solverHandle, cublasHandle_t blasHandle, int n, double a[]) { - return lu_inverse(n, a, cusolverDnDgetrf, cusolverDnDgetri); + return lu_inverse(solverHandle, blasHandle, n, a, dgetrf, dgetri, dgetrfbsize); } - DLLEXPORT int c_lu_inverse(int n, cuComplex a[]) + DLLEXPORT int c_lu_inverse(cusolverDnHandle_t solverHandle, cublasHandle_t blasHandle, int n, cuComplex a[]) { - return lu_inverse(n, a, cusolverDnCgetrf, cusolverDnCgetri); + return lu_inverse(solverHandle, blasHandle, n, a, cgetrf, cgetri, cgetrfbsize); } - DLLEXPORT int z_lu_inverse(int n, cuDoubleComplex a[]) + DLLEXPORT int z_lu_inverse(cusolverDnHandle_t solverHandle, cublasHandle_t blasHandle, int n, cuDoubleComplex a[]) { - return lu_inverse(n, a, cusolverDnZgetrf, cusolverDnZgetri); + return lu_inverse(solverHandle, blasHandle, n, a, zgetrf, zgetri, zgetrfbsize); } - DLLEXPORT int s_lu_inverse_factored(int n, float a[], int ipiv[], float work[], int lwork) + DLLEXPORT int s_lu_inverse_factored(cublasHandle_t blasHandle, int n, float a[], int ipiv[]) { - return lu_inverse_factored(n, a, ipiv, cusolverDnSgetri); + return lu_inverse_factored(blasHandle, n, a, ipiv, sgetri); } - DLLEXPORT int d_lu_inverse_factored(int n, double a[], int ipiv[], double work[], int lwork) + DLLEXPORT int d_lu_inverse_factored(cublasHandle_t blasHandle, int n, double a[], int ipiv[]) { - return lu_inverse_factored(n, a, ipiv, cusolverDnDgetri); + return lu_inverse_factored(blasHandle, n, a, ipiv, dgetri); } - DLLEXPORT int c_lu_inverse_factored(int n, cuComplex a[], int ipiv[], cuComplex work[], int lwork) + DLLEXPORT int c_lu_inverse_factored(cublasHandle_t blasHandle, int n, cuComplex a[], int ipiv[]) { - return lu_inverse_factored(n, a, ipiv, cusolverDnCgetri); + return lu_inverse_factored(blasHandle, n, a, ipiv, cgetri); } - DLLEXPORT int z_lu_inverse_factored(int n, cuDoubleComplex a[], int ipiv[], cuDoubleComplex work[], int lwork) + DLLEXPORT int z_lu_inverse_factored(cublasHandle_t blasHandle, int n, cuDoubleComplex a[], int ipiv[]) { - return lu_inverse_factored(n, a, ipiv, cusolverDnZgetri); + return lu_inverse_factored(blasHandle, n, a, ipiv, zgetri); } - DLLEXPORT int s_lu_solve_factored(int n, int nrhs, float a[], int ipiv[], float b[]) + DLLEXPORT int s_lu_solve_factored(cusolverDnHandle_t solverHandle, int n, int nrhs, float a[], int ipiv[], float b[]) { - return lu_solve_factored(n, nrhs, a, ipiv, b, cusolverDnSgetrs); + return lu_solve_factored(solverHandle, n, nrhs, a, ipiv, b, sgetrs); } - DLLEXPORT int d_lu_solve_factored(int n, int nrhs, double a[], int ipiv[], double b[]) + DLLEXPORT int d_lu_solve_factored(cusolverDnHandle_t solverHandle, int n, int nrhs, double a[], int ipiv[], double b[]) { - return lu_solve_factored(n, nrhs, a, ipiv, b, cusolverDnDgetrs); + return lu_solve_factored(solverHandle, n, nrhs, a, ipiv, b, dgetrs); } - DLLEXPORT int c_lu_solve_factored(int n, int nrhs, cuComplex a[], int ipiv[], cuComplex b[]) + DLLEXPORT int c_lu_solve_factored(cusolverDnHandle_t solverHandle, int n, int nrhs, cuComplex a[], int ipiv[], cuComplex b[]) { - return lu_solve_factored(n, nrhs, a, ipiv, b, cusolverDnCgetrs); + return lu_solve_factored(solverHandle, n, nrhs, a, ipiv, b, cgetrs); } - DLLEXPORT int z_lu_solve_factored(int n, int nrhs, cuDoubleComplex a[], int ipiv[], cuDoubleComplex b[]) + DLLEXPORT int z_lu_solve_factored(cusolverDnHandle_t solverHandle, int n, int nrhs, cuDoubleComplex a[], int ipiv[], cuDoubleComplex b[]) { - return lu_solve_factored(n, nrhs, a, ipiv, b, cusolverDnZgetrs); + return lu_solve_factored(solverHandle, n, nrhs, a, ipiv, b, zgetrs); } - DLLEXPORT int s_lu_solve(int n, int nrhs, float a[], float b[]) + DLLEXPORT int s_lu_solve(cusolverDnHandle_t solverHandle, int n, int nrhs, float a[], float b[]) { - return lu_solve(n, nrhs, a, b, cusolverDnSgetrf, cusolverDnSgetrs); + return lu_solve(solverHandle, n, nrhs, a, b, sgetrf, sgetrs, sgetrfbsize); } - DLLEXPORT int d_lu_solve(int n, int nrhs, double a[], double b[]) + DLLEXPORT int d_lu_solve(cusolverDnHandle_t solverHandle, int n, int nrhs, double a[], double b[]) { - return lu_solve(n, nrhs, a, b, cusolverDnDgetrf, cusolverDnDgetrs); + return lu_solve(solverHandle, n, nrhs, a, b, dgetrf, dgetrs, dgetrfbsize); } - DLLEXPORT int c_lu_solve(int n, int nrhs, cuComplex a[], cuComplex b[]) + DLLEXPORT int c_lu_solve(cusolverDnHandle_t solverHandle, int n, int nrhs, cuComplex a[], cuComplex b[]) { - return lu_solve(n, nrhs, a, b, cusolverDnCgetrf, cusolverDnCgetrs); + return lu_solve(solverHandle, n, nrhs, a, b, cgetrf, cgetrs, cgetrfbsize); } - DLLEXPORT int z_lu_solve(int n, int nrhs, cuDoubleComplex a[], cuDoubleComplex b[]) + DLLEXPORT int z_lu_solve(cusolverDnHandle_t solverHandle, int n, int nrhs, cuDoubleComplex a[], cuDoubleComplex b[]) { - return lu_solve(n, nrhs, a, b, cusolverDnZgetrf, cusolverDnZgetrs); + return lu_solve(solverHandle, n, nrhs, a, b, zgetrf, zgetrs, zgetrfbsize); } - DLLEXPORT int s_cholesky_factor(int n, float a[]){ - return cholesky_factor(n, a, cusolverDnSpotrf); + DLLEXPORT int s_cholesky_factor(cusolverDnHandle_t solverHandle, int n, float a[]) + { + return cholesky_factor(solverHandle, n, a, spotrf, spotrfbsize); } - DLLEXPORT int d_cholesky_factor(int n, double* a){ - return cholesky_factor(n, a, cusolverDnDpotrf); + DLLEXPORT int d_cholesky_factor(cusolverDnHandle_t solverHandle, int n, double* a) + { + return cholesky_factor(solverHandle, n, a, dpotrf, dpotrfbsize); } - DLLEXPORT int c_cholesky_factor(int n, cuComplex a[]){ - return cholesky_factor(n, a, cusolverDnCpotrf); + DLLEXPORT int c_cholesky_factor(cusolverDnHandle_t solverHandle, int n, cuComplex a[]) + { + return cholesky_factor(solverHandle, n, a, cpotrf, cpotrfbsize); } - DLLEXPORT int z_cholesky_factor(int n, cuDoubleComplex a[]){ - return cholesky_factor(n, a, cusolverDnZpotrf); + DLLEXPORT int z_cholesky_factor(cusolverDnHandle_t solverHandle, int n, cuDoubleComplex a[]) + { + return cholesky_factor(solverHandle, n, a, zpotrf, zpotrfbsize); } - DLLEXPORT int s_cholesky_solve(int n, int nrhs, float a[], float b[]) + DLLEXPORT int s_cholesky_solve(cusolverDnHandle_t solverHandle, int n, int nrhs, float a[], float b[]) { - return cholesky_solve(n, nrhs, a, b, cusolverDnSpotrf, cusolverDnSpotrs); + return cholesky_solve(solverHandle, n, nrhs, a, b, spotrf, spotrs, spotrfbsize); } - DLLEXPORT int d_cholesky_solve(int n, int nrhs, double a[], double b[]) + DLLEXPORT int d_cholesky_solve(cusolverDnHandle_t solverHandle, int n, int nrhs, double a[], double b[]) { - return cholesky_solve(n, nrhs, a, b, cusolverDnDpotrf, cusolverDnDpotrs); + return cholesky_solve(solverHandle, n, nrhs, a, b, dpotrf, dpotrs, dpotrfbsize); } - DLLEXPORT int c_cholesky_solve(int n, int nrhs, cuComplex a[], cuComplex b[]) + DLLEXPORT int c_cholesky_solve(cusolverDnHandle_t solverHandle, int n, int nrhs, cuComplex a[], cuComplex b[]) { - return cholesky_solve(n, nrhs, a, b, cusolverDnCpotrf, cusolverDnCpotrs); + return cholesky_solve(solverHandle, n, nrhs, a, b, cpotrf, cpotrs, cpotrfbsize); } - DLLEXPORT int z_cholesky_solve(int n, int nrhs, cuDoubleComplex a[], cuDoubleComplex b[]) + DLLEXPORT int z_cholesky_solve(cusolverDnHandle_t solverHandle, int n, int nrhs, cuDoubleComplex a[], cuDoubleComplex b[]) { - return cholesky_solve(n, nrhs, a, b, cusolverDnZpotrf, cusolverDnZpotrs); + return cholesky_solve(solverHandle, n, nrhs, a, b, zpotrf, zpotrs, zpotrfbsize); } - DLLEXPORT int s_cholesky_solve_factored(int n, int nrhs, float a[], float b[]) + DLLEXPORT int s_cholesky_solve_factored(cusolverDnHandle_t solverHandle, int n, int nrhs, float a[], float b[]) { - return cholesky_solve_factored(n, nrhs, a, b, cusolverDnSpotrs); + return cholesky_solve_factored(solverHandle, n, nrhs, a, b, spotrs); } - DLLEXPORT int d_cholesky_solve_factored(int n, int nrhs, double a[], double b[]) + DLLEXPORT int d_cholesky_solve_factored(cusolverDnHandle_t solverHandle, int n, int nrhs, double a[], double b[]) { - return cholesky_solve_factored(n, nrhs, a, b, cusolverDnDpotrs); + return cholesky_solve_factored(solverHandle, n, nrhs, a, b, dpotrs); } - DLLEXPORT int c_cholesky_solve_factored(int n, int nrhs, cuComplex a[], cuComplex b[]) + DLLEXPORT int c_cholesky_solve_factored(cusolverDnHandle_t solverHandle, int n, int nrhs, cuComplex a[], cuComplex b[]) { - return cholesky_solve_factored(n, nrhs, a, b, cusolverDnCpotrs); + return cholesky_solve_factored(solverHandle, n, nrhs, a, b, cpotrs); } - DLLEXPORT int z_cholesky_solve_factored(int n, int nrhs, cuDoubleComplex a[], cuDoubleComplex b[]) + DLLEXPORT int z_cholesky_solve_factored(cusolverDnHandle_t solverHandle, int n, int nrhs, cuDoubleComplex a[], cuDoubleComplex b[]) { - return cholesky_solve_factored(n, nrhs, a, b, cusolverDnZpotrs); + return cholesky_solve_factored(solverHandle, n, nrhs, a, b, zpotrs); } + // MJ: I am fairly certain that it would be straightforward to implement ?orgqr and ?gels but I'm focusing on getting the low-hanging fruit working first /*DLLEXPORT int s_qr_factor(int m, int n, float r[], float tau[], float q[], float work[], int len) { - return qr_factor(m, n, r, tau, q, work, len, cusolverDnSgeqrf, cusolverDnSorgqr); + return qr_factor(m, n, r, tau, q, work, len, sgeqrf, sorgqr); } DLLEXPORT int s_qr_thin_factor(int m, int n, float q[], float tau[], float r[], float work[], int len) { - return qr_thin_factor(m, n, q, tau, r, work, len, cusolverDnSgeqrf, cusolverDnSorgqr); + return qr_thin_factor(m, n, q, tau, r, work, len, sgeqrf, sorgqr); } DLLEXPORT int d_qr_factor(int m, int n, double r[], double tau[], double q[], double work[], int len) { - return qr_factor(m, n, r, tau, q, work, len, cusolverDnDgeqrf, cusolverDnDorgqr); + return qr_factor(m, n, r, tau, q, work, len, dgeqrf, dorgqr); } DLLEXPORT int d_qr_thin_factor(int m, int n, double q[], double tau[], double r[], double work[], int len) { - return qr_thin_factor(m, n, q, tau, r, work, len, cusolverDnDgeqrf, cusolverDnDorgqr); + return qr_thin_factor(m, n, q, tau, r, work, len, dgeqrf, dorgqr); } DLLEXPORT int c_qr_factor(int m, int n, cuComplex r[], cuComplex tau[], cuComplex q[], cuComplex work[], int len) { - return qr_factor(m, n, r, tau, q, work, len, cusolverDnCgeqrf, cusolverDnCungqr); + return qr_factor(m, n, r, tau, q, work, len, cgeqrf, cungqr); } DLLEXPORT int c_qr_thin_factor(int m, int n, cuComplex q[], cuComplex tau[], cuComplex r[], cuComplex work[], int len) { - return qr_thin_factor(m, n, q, tau, r, work, len, cusolverDnCgeqrf, cusolverDnCungqr); + return qr_thin_factor(m, n, q, tau, r, work, len, cgeqrf, cungqr); } - DLLEXPORT int z_qr_factor(int m, int n, cuDoubleComplex r[], cuDoubleComplex tau[], cuDoubleComplex q[]) + DLLEXPORT int z_qr_factor(int m, int n, cuDoubleComplex r[], cuDoubleComplex tau[], cuDoubleComplex q[], cuDoubleComplex work[], int len) { - return qr_factor(m, n, r, tau, q, work, len, cusolverDnZgeqrf, cusolverDnZungqr); + return qr_factor(m, n, r, tau, q, work, len, zgeqrf, zungqr); } - DLLEXPORT int z_qr_thin_factor(int m, int n, cuDoubleComplex q[], cuDoubleComplex tau[], cuDoubleComplex r[]) + DLLEXPORT int z_qr_thin_factor(int m, int n, cuDoubleComplex q[], cuDoubleComplex tau[], cuDoubleComplex r[], cuDoubleComplex work[], int len) { - return qr_thin_factor(m, n, q, tau, r, work, len, cusolverDnZgeqrf, cusolverDnZungqr); + return qr_thin_factor(m, n, q, tau, r, work, len, zgeqrf, zungqr); } DLLEXPORT int s_qr_solve(int m, int n, int bn, float a[], float b[], float x[], float work[], int len) { - return qr_solve(m, n, bn, a, b, x, work, len, sgels); + return qr_solve(m, n, bn, a, b, x, work, len, sgels); } DLLEXPORT int d_qr_solve(int m, int n, int bn, double a[], double b[], double x[], double work[], int len) { - return qr_solve(m, n, bn, a, b, x, work, len, dgels); + return qr_solve(m, n, bn, a, b, x, work, len, dgels); } DLLEXPORT int c_qr_solve(int m, int n, int bn, cuComplex a[], cuComplex b[], cuComplex x[], cuComplex work[], int len) { - return qr_solve(m, n, bn, a, b, x, work, len, cgels); + return qr_solve(m, n, bn, a, b, x, work, len, cgels); } DLLEXPORT int z_qr_solve(int m, int n, int bn, cuDoubleComplex a[], cuDoubleComplex b[], cuDoubleComplex x[], cuDoubleComplex work[], int len) { - return qr_solve(m, n, bn, a, b, x, work, len, zgels); - } + return qr_solve(m, n, bn, a, b, x, work, len, zgels); + }*/ + + //DLLEXPORT int s_qr_solve_factored(cusolverDnHandle_t solverHandle, cublasHandle_t blasHandle, int m, int n, int bn, float r[], float b[], float tau[], float x[], float work[], int len) + //{ + // return qr_solve_factored(solverHandle, blasHandle, m, n, bn, r, b, tau, x, work, len, sormqr, cublasStrsm); + //} + + //DLLEXPORT int d_qr_solve_factored(cusolverDnHandle_t solverHandle, cublasHandle_t blasHandle, int m, int n, int bn, double r[], double b[], double tau[], double x[], double work[], int len) + //{ + // return qr_solve_factored(solverHandle, blasHandle, m, n, bn, r, b, tau, x, work, len, dormqr, cublasDtrsm); + //} - DLLEXPORT int s_qr_solve_factored(int m, int n, int bn, float r[], float b[], float tau[], float x[], float work[], int len) + //DLLEXPORT int c_qr_solve_factored(int m, int n, int bn, cuComplex r[], cuComplex b[], cuComplex tau[], cuComplex x[], cuComplex work[], int len) + //{ + // return complex_qr_solve_factored(m, n, bn, r, b, tau, x, work, len, cunmqr, cublasCtrsm); + //} + + //DLLEXPORT int z_qr_solve_factored(int m, int n, int bn, cuDoubleComplex r[], cuDoubleComplex b[], cuDoubleComplex tau[], cuDoubleComplex x[], cuDoubleComplex work[], int len) + //{ + // return complex_qr_solve_factored(m, n, bn, r, b, tau, x, work, len, zunmqr, cublasZtrsm); + //} + + DLLEXPORT int s_svd_factor(cusolverDnHandle_t solverHandle, bool compute_vectors, int m, int n, float a[], float s[], float u[], float v[]) { - return qr_solve_factored(m, n, bn, r, b, tau, x, work, len, sormqr, cblas_strsm); + return svd_factor(solverHandle, compute_vectors, m, n, a, s, u, v, sgesvd, sgesvdbsize); } - DLLEXPORT int d_qr_solve_factored(int m, int n, int bn, double r[], double b[], double tau[], double x[], double work[], int len) + DLLEXPORT int d_svd_factor(cusolverDnHandle_t solverHandle, bool compute_vectors, int m, int n, double a[], double s[], double u[], double v[]) { - return qr_solve_factored(m, n, bn, r, b, tau, x, work, len, dormqr, cblas_dtrsm); + return svd_factor(solverHandle, compute_vectors, m, n, a, s, u, v,dgesvd, dgesvdbsize); } - DLLEXPORT int c_qr_solve_factored(int m, int n, int bn, cuComplex r[], cuComplex b[], cuComplex tau[], cuComplex x[], cuComplex work[], int len) + DLLEXPORT int c_svd_factor(cusolverDnHandle_t solverHandle, bool compute_vectors, int m, int n, cuComplex a[], cuComplex s[], cuComplex u[], cuComplex v[]) { - return complex_qr_solve_factored(m, n, bn, r, b, tau, x, work, len, cunmqr, cblas_ctrsm); + return complex_svd_factor(solverHandle, compute_vectors, m, n, a, s, u, v, cgesvd, cgesvdbsize); } - DLLEXPORT int z_qr_solve_factored(int m, int n, int bn, cuDoubleComplex r[], cuDoubleComplex b[], cuDoubleComplex tau[], cuDoubleComplex x[], cuDoubleComplex work[], int len) + DLLEXPORT int z_svd_factor(cusolverDnHandle_t solverHandle, bool compute_vectors, int m, int n, cuDoubleComplex a[], cuDoubleComplex s[], cuDoubleComplex u[], cuDoubleComplex v[]) { - return complex_qr_solve_factored(m, n, bn, r, b, tau, x, work, len, zunmqr, cblas_ztrsm); + return complex_svd_factor(solverHandle, compute_vectors, m, n, a, s, u, v, zgesvd, zgesvdbsize); } - DLLEXPORT int s_svd_factor(bool compute_vectors, int m, int n, float a[], float s[], float u[], float v[], float work[], int len) + /*DLLEXPORT int s_eigen(bool isSymmetric, int n, float a[], float vectors[], cuDoubleComplex values[], float d[]) { - return svd_factor(compute_vectors, m, n, a, s, u, v, work, len, sgesvd); + if (isSymmetric) + { + return sym_eigen_factor(n, a, vectors, values, d, LAPACKE_ssyev); + } + else + { + return eigen_factor(n, a, vectors, values, d, LAPACKE_sgees, LAPACKE_strevc); + } } - DLLEXPORT int d_svd_factor(bool compute_vectors, int m, int n, double a[], double s[], double u[], double v[], double work[], int len) + DLLEXPORT int d_eigen(bool isSymmetric, int n, double a[], double vectors[], cuDoubleComplex values[], double d[]) { - return svd_factor(compute_vectors, m, n, a, s, u, v, work, len, dgesvd); + if (isSymmetric) + { + return sym_eigen_factor(n, a, vectors, values, d, LAPACKE_dsyev); + } + else + { + return eigen_factor(n, a, vectors, values, d, LAPACKE_dgees, LAPACKE_dtrevc); + } } - DLLEXPORT int c_svd_factor(bool compute_vectors, int m, int n, cuComplex a[], cuComplex s[], cuComplex u[], cuComplex v[], cuComplex work[], int len) + DLLEXPORT int c_eigen(bool isSymmetric, int n, cuComplex a[], cuComplex vectors[], cuDoubleComplex values[], cuComplex d[]) { - return complex_svd_factor(compute_vectors, m, n, a, s, u, v, work, len, cgesvd); + if (isSymmetric) + { + return sym_eigen_factor(n, a, vectors, values, d, LAPACKE_cheev); + } + else + { + return eigen_complex_factor(n, a, vectors, values, d, LAPACKE_cgees, LAPACKE_ctrevc); + } } - DLLEXPORT int z_svd_factor(bool compute_vectors, int m, int n, cuDoubleComplex a[], cuDoubleComplex s[], cuDoubleComplex u[], cuDoubleComplex v[], cuDoubleComplex work[], int len) + DLLEXPORT int z_eigen(bool isSymmetric, int n, cuDoubleComplex a[], cuDoubleComplex vectors[], cuDoubleComplex values[], cuDoubleComplex d[]) { - return complex_svd_factor(compute_vectors, m, n, a, s, u, v, work, len, zgesvd); + if (isSymmetric) + { + return sym_eigen_factor(n, a, vectors, values, d, LAPACKE_zheev); + } + else + { + return eigen_complex_factor(n, a, vectors, values, d, LAPACKE_zgees, LAPACKE_ztrevc); + } }*/ -} \ No newline at end of file +} diff --git a/src/NativeProviders/Windows/CUDA/CUDAWrapper.vcxproj b/src/NativeProviders/Windows/CUDA/CUDAWrapper.vcxproj index ec8944fd..a9786588 100644 --- a/src/NativeProviders/Windows/CUDA/CUDAWrapper.vcxproj +++ b/src/NativeProviders/Windows/CUDA/CUDAWrapper.vcxproj @@ -23,10 +23,8 @@ - - + - {5A52B796-7F41-4C90-8DE2-F3F391C4482C} @@ -118,7 +116,7 @@ true - cusolver.lib;cublas.lib;cublas_device.lib;%(AdditionalDependencies) + cudart.lib;cusolver.lib;cublas.lib;cublas_device.lib;%(AdditionalDependencies) $(CUDA_PATH)\lib\x64;%(AdditionalLibraryDirectories) @@ -152,7 +150,7 @@ true true true - cusolver.lib;cublas.lib;%(AdditionalDependencies) + cudart.lib;cusolver.lib;cublas.lib;%(AdditionalDependencies) $(CUDA_PATH)\lib\x64;%(AdditionalLibraryDirectories) diff --git a/src/NativeProviders/Windows/CUDA/CUDAWrapper.vcxproj.filters b/src/NativeProviders/Windows/CUDA/CUDAWrapper.vcxproj.filters index 0e998fbe..9bda620d 100644 --- a/src/NativeProviders/Windows/CUDA/CUDAWrapper.vcxproj.filters +++ b/src/NativeProviders/Windows/CUDA/CUDAWrapper.vcxproj.filters @@ -23,16 +23,10 @@ Source Files - - Source Files - - - Source Files - Source Files - + Source Files