From a3dee9d03ec9c1b916ffb0b063f2a9c6a0fd319d Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Sat, 18 Aug 2018 01:56:23 +0200 Subject: [PATCH] Add PreconditionIC --- include/deal.II/lac/cuda_precondition.h | 396 +++-- source/lac/CMakeLists.txt | 1 + source/lac/cuda_precondition.cu | 1055 +++++------- tests/cuda/precondition_01.cu | 81 + tests/cuda/precondition_01.output | 168 ++ tests/cuda/precondition_02.cu | 2097 +++++++++++++++++++++++ tests/cuda/precondition_02.output | 87 + 7 files changed, 3059 insertions(+), 826 deletions(-) create mode 100644 tests/cuda/precondition_01.cu create mode 100644 tests/cuda/precondition_01.output create mode 100644 tests/cuda/precondition_02.cu create mode 100644 tests/cuda/precondition_02.output diff --git a/include/deal.II/lac/cuda_precondition.h b/include/deal.II/lac/cuda_precondition.h index 724c5f44b5..076b483ae6 100644 --- a/include/deal.II/lac/cuda_precondition.h +++ b/include/deal.II/lac/cuda_precondition.h @@ -13,222 +13,268 @@ // // --------------------------------------------------------------------- +#ifndef dealii_cuda_precondition_h +#define dealii_cuda_precondition_h + #include +#include + #include +#ifdef DEAL_II_WITH_CUDA + DEAL_II_NAMESPACE_OPEN +// forward-definition +namespace LinearAlgebra +{ namespace CUDAWrappers { + template + class Vector; + } +} // namespace LinearAlgebra + +namespace CUDAWrappers +{ + // forward definition + template + class SparseMatrix; + + /** + * This class implements an incomplete Cholesky factorization (IC) + * preconditioner for @em symmetric CUDAWrappers::SparseMatrix matrices. + * + * The implementation closely follows the one documented in the cuSPARSE + * documentation + * (https://docs.nvidia.com/cuda/cusparse/index.html#cusparse-lt-t-gt-csric02). + * + * @note Instantiations for this template are provided for @ and + * @. + * + * @ingroup Preconditioners CUDAWrappers + * @author Daniel Arndt + * @date 2018 + */ + template + class PreconditionIC + { + public: /** - * This class implements an incomplete Cholesky factorization (IC) - * preconditioner for @em symmetric CUDAWrappers::SparseMatrix matrices. - * - * The implementation closely follows the one documented in the cuSPARSE - * documentation - * (https://docs.nvidia.com/cuda/cusparse/index.html#cusparse-lt-t-gt-csric02). - * - * @note Instantiations for this template are provided for @ and - * @. - * - * @ingroup Preconditioners CUDAWrappers - * @author Daniel Arndt - * @date 2018 + * Declare the type for container size. */ - template - class PreconditionIC + using size_type = int; + + /** + * Standardized data struct to pipe additional flags to the + * preconditioner. + */ + struct AdditionalData { - public: /** - * Declare the type for container size. + * Constructor. cuSPARSE allows to compute and use level information. + * According to the documentation it is this might improve performance. + * It is suggested to try both options. */ - using size_type = int; + AdditionalData(bool use_level_analysis = true); /** - * Standardized data struct to pipe additional flags to the - * preconditioner. + * Flag that determines if level informations are used when creating and + * applying the preconditioner. See the documentation for + * cusparseSolvePolicy_t at + * https://docs.nvidia.com/cuda/cusparse/index.html#cusparsesolvepolicy_t + * for more information. */ - struct AdditionalData - { - /** - * Constructor. cuSPARSE allows to compute and use level information. - * According to the documentation it is this might improve performance. - * It is suggested to try both options. - */ - AdditionalData(bool use_level_analysis = true); - - /** - * Flag that determines if level informations are used when creating and - * applying the preconditioner. See the documentation for - * cusparseSolvePolicy_t at - * https://docs.nvidia.com/cuda/cusparse/index.html#cusparsesolvepolicy_t - * for more information. - */ - bool use_level_analysis; - }; + bool use_level_analysis; + }; - /** - * Constructor. - */ - PreconditionIC(const Utilities::CUDA::Handle &handle); + /** + * Constructor. + */ + PreconditionIC(const Utilities::CUDA::Handle &handle); - /** - * The copy constructor is deleted. - */ - PreconditionIC(const PreconditionIC &) = delete; + /** + * The copy constructor is deleted. + */ + PreconditionIC(const PreconditionIC &) = delete; - /** - * The copy assignment operator is deleted. - */ - PreconditionIC & - operator=(const PreconditionIC &) = delete; + /** + * The copy assignment operator is deleted. + */ + PreconditionIC & + operator=(const PreconditionIC &) = delete; - /** - * Destructor. Free all resources that were initialized in this class. - */ - ~PreconditionIC(); + /** + * Destructor. Free all resources that were initialized in this class. + */ + ~PreconditionIC(); - /** - * Initialize this object. In particular, the given matrix is copied to be - * modified in-place. For the underlying sparsity pattern pointers are - * stored. Specifically, this means - * that the current object can only be used reliably as long as @p matrix is valid - * and has not been changed since calling this function. - * - * The @p additional_data determines if level information are used. - */ - void - initialize(const SparseMatrix &matrix, - const AdditionalData &additional_data = AdditionalData()); + /** + * Initialize this object. In particular, the given matrix is copied to be + * modified in-place. For the underlying sparsity pattern pointers are + * stored. Specifically, this means + * that the current object can only be used reliably as long as @p matrix is valid + * and has not been changed since calling this function. + * + * The @p additional_data determines if level information are used. + */ + void + initialize(const SparseMatrix &matrix, + const AdditionalData & additional_data = AdditionalData()); - /** - * Apply the preconditioner. - */ - void - vmult(LinearAlgebra::CUDAWrappers::Vector & dst, - const LinearAlgebra::CUDAWrappers::Vector &src) const; + /** + * Apply the preconditioner. + */ + void + vmult(LinearAlgebra::CUDAWrappers::Vector & dst, + const LinearAlgebra::CUDAWrappers::Vector &src) const; - /** - * Apply the preconditioner. Since the preconditioner is symmetric, this - * is the same as vmult(). - */ - void - Tvmult(LinearAlgebra::CUDAWrappers::Vector & dst, - const LinearAlgebra::CUDAWrappers::Vector &src) const; + /** + * Apply the preconditioner. Since the preconditioner is symmetric, this + * is the same as vmult(). + */ + void + Tvmult(LinearAlgebra::CUDAWrappers::Vector & dst, + const LinearAlgebra::CUDAWrappers::Vector &src) const; - /** - * Return the dimension of the codomain (or range) space. Note that the - * matrix is square and has dimension $m \times m$. - * - * @note This function should only be called if the preconditioner has been - * initialized. - */ - size_type - m() const; + /** + * Return the dimension of the codomain (or range) space. Note that the + * matrix is square and has dimension $m \times m$. + * + * @note This function should only be called if the preconditioner has been + * initialized. + */ + size_type + m() const; - /** - * Return the dimension of the codomain (or range) space. Note that the - * matrix is square and has dimension $m \times m$. - * - * @note This function should only be called if the preconditioner has been - * initialized. - */ - size_type - n() const; + /** + * Return the dimension of the codomain (or range) space. Note that the + * matrix is square and has dimension $m \times m$. + * + * @note This function should only be called if the preconditioner has been + * initialized. + */ + size_type + n() const; - private: - /** - * cuSPARSE handle used to call cuSPARSE functions. - */ - cusparseHandle_t cusparse_handle; + private: + /** + * cuSPARSE handle used to call cuSPARSE functions. + */ + cusparseHandle_t cusparse_handle; - /** - * cuSPARSE description of the sparse matrix $M=LL^T$. - */ - cusparseMatDescr_t descr_M; + /** + * cuSPARSE description of the sparse matrix $M=LL^T$. + */ + cusparseMatDescr_t descr_M; - /** - * cuSPARSE description of the lower triangular matrix $L$. - */ - cusparseMatDescr_t descr_L; + /** + * cuSPARSE description of the lower triangular matrix $L$. + */ + cusparseMatDescr_t descr_L; - /** - * Solve and analysis structure for $M=LL^T$. - */ - csric02Info_t info_M; + /** + * Solve and analysis structure for $M=LL^T$. + */ + csric02Info_t info_M; - /** - * Solve and analysis structure for the lower triangular matrix $L$. - */ - csrsv2Info_t info_L; + /** + * Solve and analysis structure for the lower triangular matrix $L$. + */ + csrsv2Info_t info_L; - /** - * Solve and analysis structure for the upper triangular matrix $L^T$. - */ - csrsv2Info_t info_Lt; + /** + * Solve and analysis structure for the upper triangular matrix $L^T$. + */ + csrsv2Info_t info_Lt; - /** - * Pointer to the values (on the device) of the computed preconditioning - * matrix. - */ - std::unique_ptr P_val_dev; + /** + * Pointer to the values (on the device) of the computed preconditioning + * matrix. + */ + std::unique_ptr P_val_dev; - /** - * Pointer to the row pointer (on the device) of the sparse matrix this - * object was initialized with. - */ - const int *P_row_ptr_dev; + /** + * Pointer to the row pointer (on the device) of the sparse matrix this + * object was initialized with. + */ + const int *P_row_ptr_dev; - /** - * Pointer to the column indices (on the device) of the sparse matrix this - * object was initialized with. - */ - const int *P_column_index_dev; + /** + * Pointer to the column indices (on the device) of the sparse matrix this + * object was initialized with. + */ + const int *P_column_index_dev; - /** - * Pointer to the value (on the device) for a temporary (helper) vector - * used in vmult(). - */ - std::unique_ptr tmp_dev; + /** + * Pointer to the value (on the device) for a temporary (helper) vector + * used in vmult(). + */ + std::unique_ptr tmp_dev; - /** - * - */ - std::unique_ptr buffer_dev; + /** + * + */ + std::unique_ptr buffer_dev; - /** - * Determine if level information should be generated for the lower - * triangular matrix $L$. This value can be modified through an - * AdditionalData object. - */ - cusparseSolvePolicy_t policy_L; + /** + * Determine if level information should be generated for the lower + * triangular matrix $L$. This value can be modified through an + * AdditionalData object. + */ + cusparseSolvePolicy_t policy_L; - /** - * Determine if level information should be generated for the upper - * triangular matrix $L^T$. This value can be modified through an - * AdditionalData object. - */ - cusparseSolvePolicy_t policy_Lt; + /** + * Determine if level information should be generated for the upper + * triangular matrix $L^T$. This value can be modified through an + * AdditionalData object. + */ + cusparseSolvePolicy_t policy_Lt; - /** - * Determine if level information should be generated for $M=LL^T$. This - * value can be modified through an AdditionalData object. - */ - cusparseSolvePolicy_t policy_M; + /** + * Determine if level information should be generated for $M=LL^T$. This + * value can be modified through an AdditionalData object. + */ + cusparseSolvePolicy_t policy_M; - /** - * The number of rows is the same as for the matrix this object has been - * initialized with. - */ - int n_rows; + /** + * The number of rows is the same as for the matrix this object has been + * initialized with. + */ + int n_rows; - /** - * The number of non-zero elements is the same as for the matrix this - * object has been initialized with. - */ - int n_nonzero_elements; - }; + /** + * The number of non-zero elements is the same as for the matrix this + * object has been initialized with. + */ + int n_nonzero_elements; + }; + + /*--------------------------- inline functions ----------------------------*/ + +# ifndef DOXYGEN + template + inline typename PreconditionIC::size_type + PreconditionIC::m() const + { + return n_rows; } + + + template + inline typename PreconditionIC::size_type + PreconditionIC::n() const + { + return n_rows; + } +# endif // DOXYGEN + +} // namespace CUDAWrappers + DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_WITH_CUDA + +#endif // dealii_cuda_precondition_h diff --git a/source/lac/CMakeLists.txt b/source/lac/CMakeLists.txt index 0abdfcba5d..f08be24285 100644 --- a/source/lac/CMakeLists.txt +++ b/source/lac/CMakeLists.txt @@ -147,6 +147,7 @@ IF(DEAL_II_WITH_CUDA) SET(_separate_src ${_separate_src} cuda_kernels.cu + cuda_precondition.cu cuda_solver_direct.cu cuda_sparse_matrix.cu cuda_vector.cu diff --git a/source/lac/cuda_precondition.cu b/source/lac/cuda_precondition.cu index ab67f7bce8..16da6e7782 100644 --- a/source/lac/cuda_precondition.cu +++ b/source/lac/cuda_precondition.cu @@ -13,7 +13,8 @@ // // --------------------------------------------------------------------- -#include +#include +#include DEAL_II_NAMESPACE_OPEN @@ -27,16 +28,16 @@ namespace */ template cusparseStatus_t - cusparseXcsric02(cusparseHandle_t handle, - int m, - int nnz, - const cusparseMatDescr_t descrA, - Number * csrValA_valM, - const int * csrRowPtrA, - const int * csrColIndA, - csric02Info_t info, - cusparseSolvePolicy_t policy, - void * pBuffer) + cusparseXcsric02(cusparseHandle_t /*handle*/, + int /*m*/, + int /*nnz*/, + const cusparseMatDescr_t /*descrA*/, + Number * /*csrValA_valM*/, + const int * /*csrRowPtrA*/, + const int * /*csrColIndA*/, + csric02Info_t /*info*/, + cusparseSolvePolicy_t /*policy*/, + void * /*pBuffer*/) { AssertThrow(false, ExcNotImplemented()); return CUSPARSE_STATUS_INVALID_VALUE; @@ -92,79 +93,82 @@ namespace pBuffer); } - template <> - cusparseStatus_t - cusparseXcsric02(cusparseHandle_t handle, - int m, - int nnz, - const cusparseMatDescr_t descrA, - cuComplex * csrValA_valM, - const int * csrRowPtrA, - const int * csrColIndA, - csric02Info_t info, - cusparseSolvePolicy_t policy, - void * pBuffer) - { - return cusparseCcsric02(handle, - m, - nnz, - descrA, - csrValA_valM, - csrRowPtrA, - csrColIndA, - info, - policy, - pBuffer); - } + /* + template <> + cusparseStatus_t + cusparseXcsric02(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + cuComplex * csrValA_valM, + const int * csrRowPtrA, + const int * csrColIndA, + csric02Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + return cusparseCcsric02(handle, + m, + nnz, + descrA, + csrValA_valM, + csrRowPtrA, + csrColIndA, + info, + policy, + pBuffer); + } + + template <> + cusparseStatus_t + cusparseXcsric02(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + cuDoubleComplex * csrValA_valM, + const int * csrRowPtrA, + const int * csrColIndA, + csric02Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + return cusparseZcsric02(handle, + m, + nnz, + descrA, + csrValA_valM, + csrRowPtrA, + csrColIndA, + info, + policy, + pBuffer); + } + */ - template <> - cusparseStatus_t - cusparseXcsric02(cusparseHandle_t handle, - int m, - int nnz, - const cusparseMatDescr_t descrA, - cuDoubleComplex * csrValA_valM, - const int * csrRowPtrA, - const int * csrColIndA, - csric02Info_t info, - cusparseSolvePolicy_t policy, - void * pBuffer) - { - return cusparseZcsric02(handle, - m, - nnz, - descrA, - csrValA_valM, - csrRowPtrA, - csrColIndA, - info, - policy, - pBuffer); - } /** * Template wrapper for cusparsecsrsv2_solve - *(https://docs.nvidia.com/cuda/cusparse/index.html#cusparse-lt-t-gt-csrsv2_solve). + * (https://docs.nvidia.com/cuda/cusparse/index.html#cusparse-lt-t-gt-csrsv2_solve). * This function performs the solve phase of csrsv2, a new sparse triangular - *linear system op(A)*y = alpha*x. + * linear system op(A)*y = alpha*x. */ template cusparseStatus_t - cusparseXcsrsv2_solve(cusparseHandle_t handle, - cusparseOperation_t transA, - int m, - int nnz, - const Number * alpha, - const cusparseMatDescr_t descra, - const Number * csrValA, - const int * csrRowPtrA, - const int * csrColIndA, - csrsv2Info_t info, - const Number * x, - Number * y, - cusparseSolvePolicy_t policy, - void * pBuffer) + cusparseXcsrsv2_solve(cusparseHandle_t /*handle*/, + cusparseOperation_t /*transA*/, + int /*m*/, + int /*nnz*/, + const Number * /*alpha*/, + const cusparseMatDescr_t /*descra*/, + const Number * /*csrValA*/, + const int * /*csrRowPtrA*/, + const int * /*csrColIndA*/, + csrsv2Info_t /*info*/, + const Number * /*x*/, + Number * /*y*/, + cusparseSolvePolicy_t /*policy*/, + void * /*pBuffer*/) { AssertThrow(false, ExcNotImplemented()); return CUSPARSE_STATUS_INVALID_VALUE; @@ -236,6 +240,7 @@ namespace pBuffer); } + /* template <> cusparseStatus_t cusparseXcsrsv2_solve(cusparseHandle_t handle, @@ -301,6 +306,8 @@ namespace policy, pBuffer); } +*/ + /** @@ -311,17 +318,17 @@ namespace */ template cusparseStatus_t - cusparseXcsrsv2_analysis(cusparseHandle_t handle, - cusparseOperation_t transA, - int m, - int nnz, - const cusparseMatDescr_t descrA, - const Number * csrValA, - const int * csrRowPtrA, - const int * csrColIndA, - csrsv2Info_t info, - cusparseSolvePolicy_t policy, - void * pBuffer) + cusparseXcsrsv2_analysis(cusparseHandle_t /*handle*/, + cusparseOperation_t /*transA*/, + int /*m*/, + int /*nnz*/, + const cusparseMatDescr_t /*descrA*/, + const Number * /*csrValA*/, + const int * /*csrRowPtrA*/, + const int * /*csrColIndA*/, + csrsv2Info_t /*info*/, + cusparseSolvePolicy_t /*policy*/, + void * /*pBuffer*/) { AssertThrow(false, ExcNotImplemented()); return CUSPARSE_STATUS_INVALID_VALUE; @@ -381,6 +388,7 @@ namespace pBuffer); } + /* template <> cusparseStatus_t cusparseXcsrsv2_analysis(cusparseHandle_t handle, @@ -434,7 +442,7 @@ namespace policy, pBuffer); } - +*/ /** @@ -445,16 +453,16 @@ namespace */ template cusparseStatus_t - cusparseXcsric02_analysis(cusparseHandle_t handle, - int m, - int nnz, - const cusparseMatDescr_t descrA, - const Number * csrValA, - const int * csrRowPtrA, - const int * csrColIndA, - csric02Info_t info, - cusparseSolvePolicy_t policy, - void * pBuffer) + cusparseXcsric02_analysis(cusparseHandle_t /*handle*/, + int /*m*/, + int /*nnz*/, + const cusparseMatDescr_t /*descrA*/, + const Number * /*csrValA*/, + const int * /*csrRowPtrA*/, + const int * /*csrColIndA*/, + csric02Info_t /*info*/, + cusparseSolvePolicy_t /*policy*/, + void * /*pBuffer*/) { AssertThrow(false, ExcNotImplemented()); return CUSPARSE_STATUS_INVALID_VALUE; @@ -510,6 +518,7 @@ namespace pBuffer); } + /* template <> cusparseStatus_t cusparseXcsric02_analysis(cusparseHandle_t handle, @@ -559,6 +568,8 @@ namespace policy, pBuffer); } +*/ + /** @@ -569,16 +580,16 @@ namespace */ template cusparseStatus_t - cusparseXcsrsv2_bufferSize(cusparseHandle_t handle, - cusparseOperation_t transA, - int m, - int nnz, - const cusparseMatDescr_t descrA, - Number * csrValA, - const int * csrRowPtrA, - const int * csrColIndA, - csrsv2Info_t info, - int * pBufferSizeInBytes) + cusparseXcsrsv2_bufferSize(cusparseHandle_t /*handle*/, + cusparseOperation_t /*transA*/, + int /*m*/, + int /*nnz*/, + const cusparseMatDescr_t /*descrA*/, + Number * /*csrValA*/, + const int * /*csrRowPtrA*/, + const int * /*csrColIndA*/, + csrsv2Info_t /*info*/, + int * /*pBufferSizeInBytes*/) { AssertThrow(false, ExcNotImplemented()); return CUSPARSE_STATUS_INVALID_VALUE; @@ -634,6 +645,7 @@ namespace pBufferSizeInBytes); } + /* template <> cusparseStatus_t cusparseXcsrsv2_bufferSize(cusparseHandle_t handle, @@ -683,26 +695,27 @@ namespace info, pBufferSizeInBytes); } +*/ /** * Template wrapper for cusparsecsric02_bufferSize * (https://docs.nvidia.com/cuda/cusparse/index.html#cusparse-lt-t-gt-csric02_bufferSize). - *This function returns size of buffer used in computing the - *incomplete-Cholesky factorization with 0 fill-in and no pivoting. + * This function returns size of buffer used in computing the + * incomplete-Cholesky factorization with 0 fill-in and no pivoting. */ template cusparseStatus_t - cusparseXcsric02_bufferSize(cusparseHandle_t handle, - int m, - int nnz, - const cusparseMatDescr_t descrA, - Number * csrValA, - const int * csrRowPtrA, - const int * csrColIndA, - csric02Info_t info, - int * pBufferSizeInBytes) + cusparseXcsric02_bufferSize(cusparseHandle_t /*handle*/, + int /*m*/, + int /*nnz*/, + const cusparseMatDescr_t /*descrA*/, + Number * /*csrValA*/, + const int * /*csrRowPtrA*/, + const int * /*csrColIndA*/, + csric02Info_t /*info*/, + int * /*pBufferSizeInBytes*/) { AssertThrow(false, ExcNotImplemented()); return CUSPARSE_STATUS_INVALID_VALUE; @@ -754,6 +767,7 @@ namespace pBufferSizeInBytes); } + /* template <> cusparseStatus_t cusparseXcsric02_bufferSize(cusparseHandle_t handle, @@ -799,6 +813,9 @@ namespace info, pBufferSizeInBytes); } + */ + + template void @@ -820,578 +837,314 @@ namespace } } // namespace - namespace CUDAWrappers +namespace CUDAWrappers +{ + template + PreconditionIC::AdditionalData::AdditionalData( + bool use_level_analysis_) + : use_level_analysis(use_level_analysis_) + {} + + + + template + PreconditionIC::PreconditionIC(const Utilities::CUDA::Handle &handle) + : cusparse_handle(handle.cusparse_handle) + , P_val_dev(nullptr, delete_device_vector) + , P_row_ptr_dev(nullptr) + , P_column_index_dev(nullptr) + , tmp_dev(nullptr, delete_device_vector) + , buffer_dev(nullptr, delete_device_vector) + , policy_L(CUSPARSE_SOLVE_POLICY_USE_LEVEL) + , policy_Lt(CUSPARSE_SOLVE_POLICY_USE_LEVEL) + , policy_M(CUSPARSE_SOLVE_POLICY_USE_LEVEL) + , n_rows(0) + , n_nonzero_elements(0) { - template - PreconditionIC::AdditionalData::AdditionalData( - bool use_level_analysis_) - : use_level_analysis(use_level_analysis_) - {} - - - - template - PreconditionIC::PreconditionIC( - const Utilities::CUDA::Handle &handle) - : cusparse_handle(handle.cusparse_handle) - , P_val_dev(nullptr, delete_device_vector) - , P_row_ptr_dev(nullptr) - , P_column_index_dev(nullptr) - , tmp_dev(nullptr, delete_device_vector) - , buffer_dev(nullptr, delete_device_vector) - , policy_L(CUSPARSE_SOLVE_POLICY_USE_LEVEL) - , policy_Lt(CUSPARSE_SOLVE_POLICY_USE_LEVEL) - , policy_M(CUSPARSE_SOLVE_POLICY_USE_LEVEL) - , n_rows(0) - , n_nonzero_elements(0) - { - cusparseStatus_t status; - // step 1: create a descriptor which contains - // - matrix M is base-0 - // - matrix L is base-0 - // - matrix L is lower triangular - // - matrix L has non-unit diagonal - status = cusparseCreateMatDescr(&descr_M); - AssertCusparse(status); - status = cusparseSetMatIndexBase(descr_M, CUSPARSE_INDEX_BASE_ZERO); - AssertCusparse(status); - status = cusparseSetMatType(descr_M, CUSPARSE_MATRIX_TYPE_GENERAL); - AssertCusparse(status); - - status = cusparseCreateMatDescr(&descr_L); - AssertCusparse(status); - status = cusparseSetMatIndexBase(descr_L, CUSPARSE_INDEX_BASE_ZERO); - AssertCusparse(status); - status = cusparseSetMatType(descr_L, CUSPARSE_MATRIX_TYPE_GENERAL); - AssertCusparse(status); - status = cusparseSetMatFillMode(descr_L, CUSPARSE_FILL_MODE_LOWER); - AssertCusparse(status); - status = cusparseSetMatDiagType(descr_L, CUSPARSE_DIAG_TYPE_NON_UNIT); - AssertCusparse(status); - - // step 2: create a empty info structure - // we need one info for csric02 and two info's for csrsv2 - status = cusparseCreateCsric02Info(&info_M); - AssertCusparse(status); - status = cusparseCreateCsrsv2Info(&info_L); - AssertCusparse(status); - status = cusparseCreateCsrsv2Info(&info_Lt); - AssertCusparse(status); - } + cusparseStatus_t status; + // step 1: create a descriptor which contains + // - matrix M is base-0 + // - matrix L is base-0 + // - matrix L is lower triangular + // - matrix L has non-unit diagonal + status = cusparseCreateMatDescr(&descr_M); + AssertCusparse(status); + status = cusparseSetMatIndexBase(descr_M, CUSPARSE_INDEX_BASE_ZERO); + AssertCusparse(status); + status = cusparseSetMatType(descr_M, CUSPARSE_MATRIX_TYPE_GENERAL); + AssertCusparse(status); + + status = cusparseCreateMatDescr(&descr_L); + AssertCusparse(status); + status = cusparseSetMatIndexBase(descr_L, CUSPARSE_INDEX_BASE_ZERO); + AssertCusparse(status); + status = cusparseSetMatType(descr_L, CUSPARSE_MATRIX_TYPE_GENERAL); + AssertCusparse(status); + status = cusparseSetMatFillMode(descr_L, CUSPARSE_FILL_MODE_LOWER); + AssertCusparse(status); + status = cusparseSetMatDiagType(descr_L, CUSPARSE_DIAG_TYPE_NON_UNIT); + AssertCusparse(status); + + // step 2: create a empty info structure + // we need one info for csric02 and two info's for csrsv2 + status = cusparseCreateCsric02Info(&info_M); + AssertCusparse(status); + status = cusparseCreateCsrsv2Info(&info_L); + AssertCusparse(status); + status = cusparseCreateCsrsv2Info(&info_Lt); + AssertCusparse(status); + } - template - PreconditionIC::~PreconditionIC() - { - // step 8: free resources - cusparseStatus_t status = cusparseDestroyMatDescr(descr_M); - AssertNothrowCusparse(status); + template + PreconditionIC::~PreconditionIC() + { + // step 8: free resources + cusparseStatus_t status = cusparseDestroyMatDescr(descr_M); + AssertNothrowCusparse(status); - status = cusparseDestroyMatDescr(descr_L); - AssertNothrowCusparse(status); + status = cusparseDestroyMatDescr(descr_L); + AssertNothrowCusparse(status); - status = cusparseDestroyCsric02Info(info_M); - AssertNothrowCusparse(status); + status = cusparseDestroyCsric02Info(info_M); + AssertNothrowCusparse(status); - status = cusparseDestroyCsrsv2Info(info_L); - AssertNothrowCusparse(status); + status = cusparseDestroyCsrsv2Info(info_L); + AssertNothrowCusparse(status); - status = cusparseDestroyCsrsv2Info(info_Lt); - AssertNothrowCusparse(status); - } + status = cusparseDestroyCsrsv2Info(info_Lt); + AssertNothrowCusparse(status); + } - template - void - PreconditionIC::initialize(const SparseMatrix &A, - const AdditionalData &additional_data) - { - if (additional_data.use_level_analysis) - { - policy_L = CUSPARSE_SOLVE_POLICY_USE_LEVEL; - policy_Lt = CUSPARSE_SOLVE_POLICY_USE_LEVEL; - policy_M = CUSPARSE_SOLVE_POLICY_USE_LEVEL; - } - else - { - policy_L = CUSPARSE_SOLVE_POLICY_NO_LEVEL; - policy_Lt = CUSPARSE_SOLVE_POLICY_NO_LEVEL; - policy_M = CUSPARSE_SOLVE_POLICY_NO_LEVEL; - } - - n_rows = A.m(); - n_nonzero_elements = A.n_nonzero_elements(); - AssertDimension(A.m(), A.n()); - - const auto cusparse_matrix = A.get_cusparse_matrix(); - const Number *const A_val_dev = std::get<0>(cusparse_matrix); - - // create a copy of the matrix entries - P_val_dev.reset(allocate_device_vector(n_nonzero_elements)); - cudaError_t cuda_status = cudaMemcpy(P_val_dev.get(), - A_val_dev, - n_nonzero_elements * sizeof(Number), - cudaMemcpyDeviceToDevice); - P_column_index_dev = std::get<1>(cusparse_matrix); - P_row_ptr_dev = std::get<2>(cusparse_matrix); - const cusparseMatDescr_t mat_descr = std::get<3>(cusparse_matrix); - - // initializa an internal buffer we need later on - tmp_dev.reset(allocate_device_vector(n_rows)); - - // step 3: query how much memory used in csric02 and csrsv2, and allocate - // the buffer - int BufferSize_M; - cusparseStatus_t status = cusparseXcsric02_bufferSize(cusparse_handle, - n_rows, - n_nonzero_elements, - descr_M, - P_val_dev.get(), - P_row_ptr_dev, - P_column_index_dev, - info_M, - &BufferSize_M); - AssertCusparse(status); - - int BufferSize_L; - status = cusparseXcsrsv2_bufferSize(cusparse_handle, - CUSPARSE_OPERATION_NON_TRANSPOSE, - n_rows, - n_nonzero_elements, - descr_L, - P_val_dev.get(), - P_row_ptr_dev, - P_column_index_dev, - info_L, - &BufferSize_L); - AssertCusparse(status); - - int BufferSize_Lt; - status = cusparseXcsrsv2_bufferSize(cusparse_handle, - CUSPARSE_OPERATION_TRANSPOSE, - n_rows, - n_nonzero_elements, - descr_L, - P_val_dev.get(), - P_row_ptr_dev, - P_column_index_dev, - info_Lt, - &BufferSize_Lt); - AssertCusparse(status); - - const int BufferSize = - std::max(BufferSize_M, std::max(BufferSize_L, BufferSize_Lt)); - // workaround: since allocate_device_vector needs a type, we pass char - // which is required to have size 1. - buffer_dev.reset(static_cast( - allocate_device_vector(BufferSize / sizeof(char)))); - - // step 4: perform analysis of incomplete Cholesky on M - // perform analysis of triangular solve on L - // perform analysis of triangular solve on L' - // The lower triangular part of M has the same sparsity pattern as L, so - // we can do analysis of csric02 and csrsv2 simultaneously. - - status = cusparseXcsric02_analysis(cusparse_handle, - n_rows, - n_nonzero_elements, - descr_M, - P_val_dev.get(), - P_row_ptr_dev, - P_column_index_dev, - info_M, - policy_M, - buffer_dev.get()); - AssertCusparse(status); - - int structural_zero; - status = - cusparseXcsric02_zeroPivot(cusparse_handle, info_M, &structural_zero); - AssertCusparse(status); - - status = cusparseXcsrsv2_analysis(cusparse_handle, - CUSPARSE_OPERATION_TRANSPOSE, + template + void + PreconditionIC::initialize(const SparseMatrix &A, + const AdditionalData &additional_data) + { + if (additional_data.use_level_analysis) + { + policy_L = CUSPARSE_SOLVE_POLICY_USE_LEVEL; + policy_Lt = CUSPARSE_SOLVE_POLICY_USE_LEVEL; + policy_M = CUSPARSE_SOLVE_POLICY_USE_LEVEL; + } + else + { + policy_L = CUSPARSE_SOLVE_POLICY_NO_LEVEL; + policy_Lt = CUSPARSE_SOLVE_POLICY_NO_LEVEL; + policy_M = CUSPARSE_SOLVE_POLICY_NO_LEVEL; + } + + n_rows = A.m(); + n_nonzero_elements = A.n_nonzero_elements(); + AssertDimension(A.m(), A.n()); + + const auto cusparse_matrix = A.get_cusparse_matrix(); + const Number *const A_val_dev = std::get<0>(cusparse_matrix); + + // create a copy of the matrix entries + P_val_dev.reset(allocate_device_vector(n_nonzero_elements)); + cudaError_t cuda_status = cudaMemcpy(P_val_dev.get(), + A_val_dev, + n_nonzero_elements * sizeof(Number), + cudaMemcpyDeviceToDevice); + P_column_index_dev = std::get<1>(cusparse_matrix); + P_row_ptr_dev = std::get<2>(cusparse_matrix); + const cusparseMatDescr_t mat_descr = std::get<3>(cusparse_matrix); + + // initializa an internal buffer we need later on + tmp_dev.reset(allocate_device_vector(n_rows)); + + // step 3: query how much memory used in csric02 and csrsv2, and allocate + // the buffer + int BufferSize_M; + cusparseStatus_t status = cusparseXcsric02_bufferSize(cusparse_handle, + n_rows, + n_nonzero_elements, + descr_M, + P_val_dev.get(), + P_row_ptr_dev, + P_column_index_dev, + info_M, + &BufferSize_M); + AssertCusparse(status); + + int BufferSize_L; + status = cusparseXcsrsv2_bufferSize(cusparse_handle, + CUSPARSE_OPERATION_NON_TRANSPOSE, n_rows, n_nonzero_elements, descr_L, P_val_dev.get(), P_row_ptr_dev, P_column_index_dev, - info_Lt, - policy_Lt, - buffer_dev.get()); - AssertCusparse(status); + info_L, + &BufferSize_L); + AssertCusparse(status); - status = cusparseXcsrsv2_analysis(cusparse_handle, - CUSPARSE_OPERATION_NON_TRANSPOSE, + int BufferSize_Lt; + status = cusparseXcsrsv2_bufferSize(cusparse_handle, + CUSPARSE_OPERATION_TRANSPOSE, n_rows, n_nonzero_elements, descr_L, P_val_dev.get(), P_row_ptr_dev, P_column_index_dev, - info_L, - policy_L, - buffer_dev.get()); - AssertCusparse(status); - - // step 5: M = L * L' - status = cusparseXcsric02(cusparse_handle, - n_rows, - n_nonzero_elements, - descr_M, - P_val_dev.get(), - P_row_ptr_dev, - P_column_index_dev, - info_M, - policy_M, - buffer_dev.get()); - AssertCusparse(status); - - int numerical_zero; - status = - cusparseXcsric02_zeroPivot(cusparse_handle, info_M, &numerical_zero); - AssertCusparse(status); - } - - - - template - void - PreconditionIC::vmult( - LinearAlgebra::CUDAWrappers::Vector & dst, - const LinearAlgebra::CUDAWrappers::Vector &src) const - { - Assert(P_val_dev != nullptr, ExcNotInitialized()); - Assert(P_row_ptr_dev != nullptr, ExcNotInitialized()); - Assert(P_column_index_dev != nullptr, ExcNotInitialized()); - AssertDimension(dst.size(), static_cast(n_rows)); - AssertDimension(src.size(), static_cast(n_rows)); - Assert(tmp_dev != nullptr, ExcInternalError()); - - const Number *const src_dev = src.get_values(); - Number *const dst_dev = dst.get_values(); - // step 6: solve L*z = alpha*x - const double alpha = 1.; - cusparseStatus_t status = - cusparseXcsrsv2_solve(cusparse_handle, - CUSPARSE_OPERATION_NON_TRANSPOSE, + info_Lt, + &BufferSize_Lt); + AssertCusparse(status); + + const int BufferSize = + std::max(BufferSize_M, std::max(BufferSize_L, BufferSize_Lt)); + // workaround: since allocate_device_vector needs a type, we pass char + // which is required to have size 1. + buffer_dev.reset(static_cast( + allocate_device_vector(BufferSize / sizeof(char)))); + + // step 4: perform analysis of incomplete Cholesky on M + // perform analysis of triangular solve on L + // perform analysis of triangular solve on L' + // The lower triangular part of M has the same sparsity pattern as L, so + // we can do analysis of csric02 and csrsv2 simultaneously. + + status = cusparseXcsric02_analysis(cusparse_handle, + n_rows, + n_nonzero_elements, + descr_M, + P_val_dev.get(), + P_row_ptr_dev, + P_column_index_dev, + info_M, + policy_M, + buffer_dev.get()); + AssertCusparse(status); + + int structural_zero; + status = + cusparseXcsric02_zeroPivot(cusparse_handle, info_M, &structural_zero); + AssertCusparse(status); + + status = cusparseXcsrsv2_analysis(cusparse_handle, + CUSPARSE_OPERATION_TRANSPOSE, + n_rows, + n_nonzero_elements, + descr_L, + P_val_dev.get(), + P_row_ptr_dev, + P_column_index_dev, + info_Lt, + policy_Lt, + buffer_dev.get()); + AssertCusparse(status); + + status = cusparseXcsrsv2_analysis(cusparse_handle, + CUSPARSE_OPERATION_NON_TRANSPOSE, + n_rows, + n_nonzero_elements, + descr_L, + P_val_dev.get(), + P_row_ptr_dev, + P_column_index_dev, + info_L, + policy_L, + buffer_dev.get()); + AssertCusparse(status); + + // step 5: M = L * L' + status = cusparseXcsric02(cusparse_handle, n_rows, n_nonzero_elements, - &alpha, - descr_L, + descr_M, P_val_dev.get(), P_row_ptr_dev, P_column_index_dev, - info_L, - src_dev, - tmp_dev.get(), - policy_L, + info_M, + policy_M, buffer_dev.get()); - AssertCusparse(status); - - // step 7: solve L'*y = alpha*z - status = cusparseXcsrsv2_solve(cusparse_handle, - CUSPARSE_OPERATION_TRANSPOSE, - n_rows, - n_nonzero_elements, - &alpha, - descr_L, - P_val_dev.get(), - P_row_ptr_dev, - P_column_index_dev, - info_Lt, - tmp_dev.get(), - dst_dev, - policy_Lt, - buffer_dev.get()); - AssertCusparse(status); - } - + AssertCusparse(status); - - template - void - PreconditionIC::Tvmult( - LinearAlgebra::CUDAWrappers::Vector & dst, - const LinearAlgebra::CUDAWrappers::Vector &src) const - { - // the constructed preconditioner is symmetric - vmult(dst, src); - } - - - - template - PreconditionIC::size_type - PreconditionIC::m() const - { - return n_rows; - } + int numerical_zero; + status = + cusparseXcsric02_zeroPivot(cusparse_handle, info_M, &numerical_zero); + AssertCusparse(status); + } - template - PreconditionIC::size_type - PreconditionIC::n() const - { - return n_rows; - } + template + void + PreconditionIC::vmult( + LinearAlgebra::CUDAWrappers::Vector & dst, + const LinearAlgebra::CUDAWrappers::Vector &src) const + { + Assert(P_val_dev != nullptr, ExcNotInitialized()); + Assert(P_row_ptr_dev != nullptr, ExcNotInitialized()); + Assert(P_column_index_dev != nullptr, ExcNotInitialized()); + AssertDimension(dst.size(), static_cast(n_rows)); + AssertDimension(src.size(), static_cast(n_rows)); + Assert(tmp_dev != nullptr, ExcInternalError()); + + const Number *const src_dev = src.get_values(); + Number *const dst_dev = dst.get_values(); + // step 6: solve L*z = alpha*x + const Number alpha = internal::NumberType::value(1.); + cusparseStatus_t status = + cusparseXcsrsv2_solve(cusparse_handle, + CUSPARSE_OPERATION_NON_TRANSPOSE, + n_rows, + n_nonzero_elements, + &alpha, + descr_L, + P_val_dev.get(), + P_row_ptr_dev, + P_column_index_dev, + info_L, + src_dev, + tmp_dev.get(), + policy_L, + buffer_dev.get()); + AssertCusparse(status); + + // step 7: solve L'*y = alpha*z + status = cusparseXcsrsv2_solve(cusparse_handle, + CUSPARSE_OPERATION_TRANSPOSE, + n_rows, + n_nonzero_elements, + &alpha, + descr_L, + P_val_dev.get(), + P_row_ptr_dev, + P_column_index_dev, + info_Lt, + tmp_dev.get(), + dst_dev, + policy_Lt, + buffer_dev.get()); + AssertCusparse(status); + } - template - void - apply_preconditioner(const SparseMatrix &A, - const cusparseHandle_t cusparse_handle, - LinearAlgebra::CUDAWrappers::Vector & dst, - const LinearAlgebra::CUDAWrappers::Vector &src) - { - const Number *const src_dev = src.get_values(); - Number * dst_dev = dst.get_values(); - const cusparseHandle_t handle = cusparse_handle; - - const auto cusparse_matrix = A.get_cusparse_matrix(); - Number * A_val_dev = std::get<0>(cusparse_matrix); - const int *const A_row_ptr_dev = std::get<2>(cusparse_matrix); - const int *const A_column_index_dev = std::get<1>(cusparse_matrix); - const cusparseMatDescr_t mat_descr = std::get<3>(cusparse_matrix); - - const unsigned int n_rows = A.m(); - const unsigned int n_nonzero_elements = A.n_nonzero_elements(); - - AssertDimension(dst.size(), src.size()); - AssertDimension(A.m(), src.size()); - AssertDimension(A.n(), src.size()); - - std::unique_ptr tmp_dev( - allocate_device_vector(dst.size()), - delete_device_vector); - - // Suppose that A is a m x m sparse matrix represented by CSR format, - // Assumption: - // - handle is already created by cusparseCreate(), - // - (A_row_ptr_dev, A_column_index_dev, A_val_dev) is CSR of A on device - // memory, - // - src_dev is right hand side vector on device memory, - // - dst_dev is solution vector on device memory. - // - tmp_dev is intermediate result on device memory. - - cusparseMatDescr_t descr_M = mat_descr; - cusparseMatDescr_t descr_L = mat_descr; - csric02Info_t info_M = 0; - csrsv2Info_t info_L = 0; - csrsv2Info_t info_Lt = 0; - int BufferSize_M; - int BufferSize_L; - int BufferSize_Lt; - int BufferSize; - void * buffer_dev = 0; - int structural_zero; - int numerical_zero; - const double alpha = 1.; - const cusparseSolvePolicy_t policy_M = CUSPARSE_SOLVE_POLICY_NO_LEVEL; - const cusparseSolvePolicy_t policy_L = CUSPARSE_SOLVE_POLICY_NO_LEVEL; - const cusparseSolvePolicy_t policy_Lt = CUSPARSE_SOLVE_POLICY_USE_LEVEL; - - cusparseStatus_t status; - // step 1: create a descriptor which contains - // - matrix M is base-0 - // - matrix L is base-0 - // - matrix L is lower triangular - // - matrix L has non-unit diagonal - status = cusparseCreateMatDescr(&descr_M); - AssertCusparse(status); - status = cusparseSetMatIndexBase(descr_M, CUSPARSE_INDEX_BASE_ZERO); - AssertCusparse(status); - status = cusparseSetMatType(descr_M, CUSPARSE_MATRIX_TYPE_GENERAL); - AssertCusparse(status); - - status = cusparseCreateMatDescr(&descr_L); - AssertCusparse(status); - status = cusparseSetMatIndexBase(descr_L, CUSPARSE_INDEX_BASE_ZERO); - AssertCusparse(status); - status = cusparseSetMatType(descr_L, CUSPARSE_MATRIX_TYPE_GENERAL); - AssertCusparse(status); - status = cusparseSetMatFillMode(descr_L, CUSPARSE_FILL_MODE_LOWER); - AssertCusparse(status); - status = cusparseSetMatDiagType(descr_L, CUSPARSE_DIAG_TYPE_NON_UNIT); - AssertCusparse(status); - - // step 2: create a empty info structure - // we need one info for csric02 and two info's for csrsv2 - status = cusparseCreateCsric02Info(&info_M); - AssertCusparse(status); - status = cusparseCreateCsrsv2Info(&info_L); - AssertCusparse(status); - status = cusparseCreateCsrsv2Info(&info_Lt); - AssertCusparse(status); - - // step 3: query how much memory used in csric02 and csrsv2, and allocate - // the buffer - status = cusparseXcsric02_bufferSize(handle, - n_rows, - n_nonzero_elements, - descr_M, - A_val_dev, - A_row_ptr_dev, - A_column_index_dev, - info_M, - &BufferSize_M); - AssertCusparse(status); - status = cusparseXcsrsv2_bufferSize(handle, - CUSPARSE_OPERATION_NON_TRANSPOSE, - n_rows, - n_nonzero_elements, - descr_L, - A_val_dev, - A_row_ptr_dev, - A_column_index_dev, - info_L, - &BufferSize_L); - AssertCusparse(status); - status = cusparseXcsrsv2_bufferSize(handle, - CUSPARSE_OPERATION_TRANSPOSE, - n_rows, - n_nonzero_elements, - descr_L, - A_val_dev, - A_row_ptr_dev, - A_column_index_dev, - info_Lt, - &BufferSize_Lt); - AssertCusparse(status); - - BufferSize = max(BufferSize_M, max(BufferSize_L, BufferSize_Lt)); - - // buffer_dev returned by cudaMalloc is automatically aligned to 128 - // bytes. - cudaError_t status_cuda = cudaMalloc((void **)&buffer_dev, BufferSize); - Assert(cudaSuccess == status_cuda, ExcInternalError()); - - // step 4: perform analysis of incomplete Cholesky on M - // perform analysis of triangular solve on L - // perform analysis of triangular solve on L' - // The lower triangular part of M has the same sparsity pattern as L, so - // we can do analysis of csric02 and csrsv2 simultaneously. - - status = cusparseXcsric02_analysis(handle, - n_rows, - n_nonzero_elements, - descr_M, - A_val_dev, - A_row_ptr_dev, - A_column_index_dev, - info_M, - policy_M, - buffer_dev); - AssertCusparse(status); - status = cusparseXcsric02_zeroPivot(handle, info_M, &structural_zero); - if (CUSPARSE_STATUS_ZERO_PIVOT == status) - { - printf("A(%d,%d) is missing\n", structural_zero, structural_zero); - } - - status = cusparseXcsrsv2_analysis(handle, - CUSPARSE_OPERATION_TRANSPOSE, - n_rows, - n_nonzero_elements, - descr_L, - A_val_dev, - A_row_ptr_dev, - A_column_index_dev, - info_Lt, - policy_Lt, - buffer_dev); - AssertCusparse(status); - - status = cusparseXcsrsv2_analysis(handle, - CUSPARSE_OPERATION_NON_TRANSPOSE, - n_rows, - n_nonzero_elements, - descr_L, - A_val_dev, - A_row_ptr_dev, - A_column_index_dev, - info_L, - policy_L, - buffer_dev); - AssertCusparse(status); - - // step 5: M = L * L' - status = cusparseXcsric02(handle, - n_rows, - n_nonzero_elements, - descr_M, - A_val_dev, - A_row_ptr_dev, - A_column_index_dev, - info_M, - policy_M, - buffer_dev); - AssertCusparse(status); - status = cusparseXcsric02_zeroPivot(handle, info_M, &numerical_zero); - if (CUSPARSE_STATUS_ZERO_PIVOT == status) - { - printf("L(%d,%d) is zero\n", numerical_zero, numerical_zero); - } - - // step 6: solve L*z = x - status = cusparseXcsrsv2_solve(handle, - CUSPARSE_OPERATION_NON_TRANSPOSE, - n_rows, - n_nonzero_elements, - &alpha, - descr_L, - A_val_dev, - A_row_ptr_dev, - A_column_index_dev, - info_L, - src_dev, - tmp_dev.get(), - policy_L, - buffer_dev); - AssertCusparse(status); - - // step 7: solve L'*y = z - status = cusparseXcsrsv2_solve(handle, - CUSPARSE_OPERATION_TRANSPOSE, - n_rows, - n_nonzero_elements, - &alpha, - descr_L, - A_val_dev, - A_row_ptr_dev, - A_column_index_dev, - info_Lt, - tmp_dev.get(), - dst_dev, - policy_Lt, - buffer_dev); - AssertCusparse(status); - - // step 8: free resources - status_cuda = cudaFree(buffer_dev); - AssertCuda(status_cuda); - status = cusparseDestroyMatDescr(descr_M); - AssertCusparse(status); - status = cusparseDestroyMatDescr(descr_L); - AssertCusparse(status); - status = cusparseDestroyCsric02Info(info_M); - AssertCusparse(status); - status = cusparseDestroyCsrsv2Info(info_L); - AssertCusparse(status); - status = cusparseDestroyCsrsv2Info(info_Lt); - AssertCusparse(status); - } + template + void + PreconditionIC::Tvmult( + LinearAlgebra::CUDAWrappers::Vector & dst, + const LinearAlgebra::CUDAWrappers::Vector &src) const + { + // the constructed preconditioner is symmetric + vmult(dst, src); + } - // explicit instantiations - template class PreconditionIC; - template class PreconditionIC; - } // namespace CUDAWrappers + // explicit instantiations + template class PreconditionIC; + template class PreconditionIC; + // template class PreconditionIC; + // template class PreconditionIC; +} // namespace CUDAWrappers DEAL_II_NAMESPACE_CLOSE diff --git a/tests/cuda/precondition_01.cu b/tests/cuda/precondition_01.cu new file mode 100644 index 0000000000..05cbe128d1 --- /dev/null +++ b/tests/cuda/precondition_01.cu @@ -0,0 +1,81 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// Check that dealii::SolverCG works with CUDAWrappers::SparseMatrix + +#include +#include + +#include +#include +#include +#include +#include + +#include "../testmatrix.h" +#include "../tests.h" + +template +void +test(Utilities::CUDA::Handle &cuda_handle) +{ + // Build the sparse matrix on the host + const unsigned int problem_size = 10; + unsigned int size = (problem_size - 1) * (problem_size - 1); + FDMatrix testproblem(problem_size, problem_size); + SparsityPattern structure(size, size, 5); + SparseMatrix A; + testproblem.five_point_structure(structure); + structure.compress(); + A.reinit(structure); + testproblem.five_point(A); + A.print(std::cout); + + // Solve on the device + CUDAWrappers::SparseMatrix A_dev(cuda_handle, A); + LinearAlgebra::CUDAWrappers::Vector sol_dev(size); + LinearAlgebra::CUDAWrappers::Vector rhs_dev(size); + LinearAlgebra::ReadWriteVector rw_vector(size); + for (unsigned int i = 0; i < size; ++i) + rw_vector[i] = static_cast(i); + rhs_dev.import(rw_vector, VectorOperation::insert); + SolverControl control(100, 1.e-10); + SolverCG> cg_dev(control); + + CUDAWrappers::PreconditionIC prec_ic(cuda_handle); + prec_ic.initialize(A_dev); + + cg_dev.solve(A_dev, sol_dev, rhs_dev, prec_ic); + + // Check the result + rw_vector.import(sol_dev, VectorOperation::insert); + for (unsigned int i = 0; i < size; ++i) + deallog << rw_vector[i] << std::endl; +} + +int +main() +{ + initlog(); + deallog.depth_console(0); + + Utilities::CUDA::Handle cuda_handle; + test(cuda_handle); + test(cuda_handle); + + deallog << "OK" << std::endl; + + return 0; +} diff --git a/tests/cuda/precondition_01.output b/tests/cuda/precondition_01.output new file mode 100644 index 0000000000..0f5ff01738 --- /dev/null +++ b/tests/cuda/precondition_01.output @@ -0,0 +1,168 @@ + +DEAL:cg::Starting value 416.989 +DEAL:cg::Convergence step 19 value 1.41136e-11 +DEAL::20.9607 +DEAL::38.8073 +DEAL::52.3525 +DEAL::61.2758 +DEAL::65.4369 +DEAL::64.6135 +DEAL::58.3945 +DEAL::46.1455 +DEAL::27.0190 +DEAL::45.0353 +DEAL::80.9161 +DEAL::107.327 +DEAL::124.314 +DEAL::131.858 +DEAL::129.623 +DEAL::116.819 +DEAL::92.1685 +DEAL::53.9305 +DEAL::69.2645 +DEAL::122.495 +DEAL::160.725 +DEAL::184.794 +DEAL::195.060 +DEAL::191.200 +DEAL::172.090 +DEAL::135.779 +DEAL::79.5345 +DEAL::91.5276 +DEAL::160.074 +DEAL::208.285 +DEAL::238.076 +DEAL::250.389 +DEAL::245.025 +DEAL::220.563 +DEAL::174.324 +DEAL::102.428 +DEAL::109.771 +DEAL::189.991 +DEAL::245.262 +DEAL::278.838 +DEAL::292.394 +DEAL::285.950 +DEAL::257.813 +DEAL::204.524 +DEAL::120.855 +DEAL::121.567 +DEAL::207.855 +DEAL::265.937 +DEAL::300.618 +DEAL::314.399 +DEAL::307.567 +DEAL::278.215 +DEAL::222.104 +DEAL::132.468 +DEAL::123.643 +DEAL::207.924 +DEAL::263.011 +DEAL::295.300 +DEAL::308.015 +DEAL::301.706 +DEAL::274.376 +DEAL::221.208 +DEAL::133.913 +DEAL::111.079 +DEAL::182.188 +DEAL::226.884 +DEAL::252.556 +DEAL::262.656 +DEAL::257.865 +DEAL::236.376 +DEAL::193.440 +DEAL::119.974 +DEAL::75.4858 +DEAL::118.864 +DEAL::144.783 +DEAL::159.382 +DEAL::165.189 +DEAL::162.720 +DEAL::150.825 +DEAL::126.202 +DEAL::81.5441 +DEAL:cg::Starting value 416.989 +DEAL:cg::Convergence step 17 value 5.54040e-11 +DEAL::20.9607 +DEAL::38.8073 +DEAL::52.3525 +DEAL::61.2757 +DEAL::65.4369 +DEAL::64.6135 +DEAL::58.3945 +DEAL::46.1455 +DEAL::27.0190 +DEAL::45.0353 +DEAL::80.9161 +DEAL::107.327 +DEAL::124.314 +DEAL::131.858 +DEAL::129.623 +DEAL::116.819 +DEAL::92.1685 +DEAL::53.9305 +DEAL::69.2645 +DEAL::122.495 +DEAL::160.725 +DEAL::184.794 +DEAL::195.060 +DEAL::191.200 +DEAL::172.090 +DEAL::135.779 +DEAL::79.5345 +DEAL::91.5276 +DEAL::160.074 +DEAL::208.285 +DEAL::238.076 +DEAL::250.389 +DEAL::245.025 +DEAL::220.563 +DEAL::174.324 +DEAL::102.428 +DEAL::109.771 +DEAL::189.991 +DEAL::245.262 +DEAL::278.838 +DEAL::292.394 +DEAL::285.950 +DEAL::257.813 +DEAL::204.524 +DEAL::120.855 +DEAL::121.567 +DEAL::207.855 +DEAL::265.937 +DEAL::300.618 +DEAL::314.399 +DEAL::307.567 +DEAL::278.215 +DEAL::222.104 +DEAL::132.468 +DEAL::123.643 +DEAL::207.924 +DEAL::263.011 +DEAL::295.300 +DEAL::308.015 +DEAL::301.706 +DEAL::274.376 +DEAL::221.208 +DEAL::133.913 +DEAL::111.079 +DEAL::182.188 +DEAL::226.884 +DEAL::252.556 +DEAL::262.656 +DEAL::257.865 +DEAL::236.376 +DEAL::193.440 +DEAL::119.974 +DEAL::75.4858 +DEAL::118.864 +DEAL::144.783 +DEAL::159.382 +DEAL::165.189 +DEAL::162.720 +DEAL::150.825 +DEAL::126.202 +DEAL::81.5441 +DEAL::OK diff --git a/tests/cuda/precondition_02.cu b/tests/cuda/precondition_02.cu new file mode 100644 index 0000000000..c5f68ab960 --- /dev/null +++ b/tests/cuda/precondition_02.cu @@ -0,0 +1,2097 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// Check that dealii::SolverCG works with CUDAWrappers::SparseMatrix + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include + +#include "../testmatrix.h" +#include "../tests.h" + +DEAL_II_NAMESPACE_OPEN + +namespace CUDAWrappers +{ + /** \addtogroup CUDAWrappers + * @{ + */ + + /** + * Template wrapper for cusparsecsrilu02. + * (https://docs.nvidia.com/cuda/cusparse/index.html#cusparse-lt-t-gt-csrilu02). + * function performs the solve phase of the incomplete-LU factorization with + * 0 fill-in and no pivoting. + */ + template + cusparseStatus_t + cusparseXcsrilu02(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + Number * csrValA_valM, + const int * csrRowPtrA, + const int * csrColIndA, + csrilu02Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + AssertThrow(false, ExcNotImplemented()); + return CUSPARSE_STATUS_INVALID_VALUE; + } + + template <> + cusparseStatus_t + cusparseXcsrilu02(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + float * csrValA_valM, + const int * csrRowPtrA, + const int * csrColIndA, + csrilu02Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + return cusparseScsrilu02(handle, + m, + nnz, + descrA, + csrValA_valM, + csrRowPtrA, + csrColIndA, + info, + policy, + pBuffer); + } + + template <> + cusparseStatus_t + cusparseXcsrilu02(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + double * csrValA_valM, + const int * csrRowPtrA, + const int * csrColIndA, + csrilu02Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + return cusparseDcsrilu02(handle, + m, + nnz, + descrA, + csrValA_valM, + csrRowPtrA, + csrColIndA, + info, + policy, + pBuffer); + } + + template <> + cusparseStatus_t + cusparseXcsrilu02(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + cuComplex * csrValA_valM, + const int * csrRowPtrA, + const int * csrColIndA, + csrilu02Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + return cusparseCcsrilu02(handle, + m, + nnz, + descrA, + csrValA_valM, + csrRowPtrA, + csrColIndA, + info, + policy, + pBuffer); + } + + template <> + cusparseStatus_t + cusparseXcsrilu02(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + cuDoubleComplex * csrValA_valM, + const int * csrRowPtrA, + const int * csrColIndA, + csrilu02Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + return cusparseZcsrilu02(handle, + m, + nnz, + descrA, + csrValA_valM, + csrRowPtrA, + csrColIndA, + info, + policy, + pBuffer); + } + + + + /** + * Template wrapper for cusparsecsrilu02_analysis. + * (https://docs.nvidia.com/cuda/cusparse/index.html#cusparse-lt-t-gt-csrilu02_analysis). + * This function performs the analysis phase of the incomplete-LU + * factorization with 0 fill-in and no pivoting. + */ + template + cusparseStatus_t + cusparseXcsrilu02_analysis(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + const Number * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csrilu02Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + AssertThrow(false, ExcNotImplemented()); + return CUSPARSE_STATUS_INVALID_VALUE; + } + + template <> + cusparseStatus_t + cusparseXcsrilu02_analysis(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + const float * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csrilu02Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + return cusparseScsrilu02_analysis(handle, + m, + nnz, + descrA, + csrValA, + csrRowPtrA, + csrColIndA, + info, + policy, + pBuffer); + } + + template <> + cusparseStatus_t + cusparseXcsrilu02_analysis(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + const double * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csrilu02Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + return cusparseDcsrilu02_analysis(handle, + m, + nnz, + descrA, + csrValA, + csrRowPtrA, + csrColIndA, + info, + policy, + pBuffer); + } + + template <> + cusparseStatus_t + cusparseXcsrilu02_analysis(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + const cuComplex * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csrilu02Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + return cusparseCcsrilu02_analysis(handle, + m, + nnz, + descrA, + csrValA, + csrRowPtrA, + csrColIndA, + info, + policy, + pBuffer); + } + + template <> + cusparseStatus_t + cusparseXcsrilu02_analysis(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + const cuDoubleComplex * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csrilu02Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + return cusparseZcsrilu02_analysis(handle, + m, + nnz, + descrA, + csrValA, + csrRowPtrA, + csrColIndA, + info, + policy, + pBuffer); + } + + + + /** + * Template wrapper for cusparsecsrilu02_bufferSize. + * (https://docs.nvidia.com/cuda/cusparse/index.html#cusparse-lt-t-gt-csrilu02_bufferSize). + * This function returns size of the buffer used in computing the + * incomplete-LU factorization with 0 fill-in and no pivoting. + */ + template + cusparseStatus_t + cusparseXcsrilu02_bufferSize(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + Number * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csrilu02Info_t info, + int * pBufferSizeInBytes) + { + AssertThrow(false, ExcNotImplemented()); + return CUSPARSE_STATUS_INVALID_VALUE; + } + + template <> + cusparseStatus_t + cusparseXcsrilu02_bufferSize(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + float * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csrilu02Info_t info, + int *pBufferSizeInBytes) + { + return cusparseScsrilu02_bufferSize(handle, + m, + nnz, + descrA, + csrValA, + csrRowPtrA, + csrColIndA, + info, + pBufferSizeInBytes); + } + + template <> + cusparseStatus_t + cusparseXcsrilu02_bufferSize(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + double * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csrilu02Info_t info, + int *pBufferSizeInBytes) + { + return cusparseDcsrilu02_bufferSize(handle, + m, + nnz, + descrA, + csrValA, + csrRowPtrA, + csrColIndA, + info, + pBufferSizeInBytes); + } + + + template <> + cusparseStatus_t + cusparseXcsrilu02_bufferSize(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + cuComplex * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csrilu02Info_t info, + int *pBufferSizeInBytes) + { + return cusparseCcsrilu02_bufferSize(handle, + m, + nnz, + descrA, + csrValA, + csrRowPtrA, + csrColIndA, + info, + pBufferSizeInBytes); + } + + template <> + cusparseStatus_t + cusparseXcsrilu02_bufferSize(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + cuDoubleComplex *csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csrilu02Info_t info, + int *pBufferSizeInBytes) + { + return cusparseZcsrilu02_bufferSize(handle, + m, + nnz, + descrA, + csrValA, + csrRowPtrA, + csrColIndA, + info, + pBufferSizeInBytes); + } + + + /** + * Template wrapper for cusparsecsric02 + * (https://docs.nvidia.com/cuda/cusparse/index.html#cusparse-lt-t-gt-csric02). + * This function performs the solve phase of the computing the + * incomplete-Cholesky factorization with 0 fill-in and no pivoting. + */ + template + cusparseStatus_t + cusparseXcsric02(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + Number * csrValA_valM, + const int * csrRowPtrA, + const int * csrColIndA, + csric02Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + AssertThrow(false, ExcNotImplemented()); + return CUSPARSE_STATUS_INVALID_VALUE; + } + + template <> + cusparseStatus_t + cusparseXcsric02(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + float * csrValA_valM, + const int * csrRowPtrA, + const int * csrColIndA, + csric02Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + return cusparseScsric02(handle, + m, + nnz, + descrA, + csrValA_valM, + csrRowPtrA, + csrColIndA, + info, + policy, + pBuffer); + } + + template <> + cusparseStatus_t + cusparseXcsric02(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + double * csrValA_valM, + const int * csrRowPtrA, + const int * csrColIndA, + csric02Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + return cusparseDcsric02(handle, + m, + nnz, + descrA, + csrValA_valM, + csrRowPtrA, + csrColIndA, + info, + policy, + pBuffer); + } + + template <> + cusparseStatus_t + cusparseXcsric02(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + cuComplex * csrValA_valM, + const int * csrRowPtrA, + const int * csrColIndA, + csric02Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + return cusparseCcsric02(handle, + m, + nnz, + descrA, + csrValA_valM, + csrRowPtrA, + csrColIndA, + info, + policy, + pBuffer); + } + + template <> + cusparseStatus_t + cusparseXcsric02(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + cuDoubleComplex * csrValA_valM, + const int * csrRowPtrA, + const int * csrColIndA, + csric02Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + return cusparseZcsric02(handle, + m, + nnz, + descrA, + csrValA_valM, + csrRowPtrA, + csrColIndA, + info, + policy, + pBuffer); + } + + + /** + * Template wrapper for cusparsecsrsv2_solve + *(https://docs.nvidia.com/cuda/cusparse/index.html#cusparse-lt-t-gt-csrsv2_solve). + * This function performs the solve phase of csrsv2, a new sparse triangular + *linear system op(A)*y = alpha*x. + */ + template + cusparseStatus_t + cusparseXcsrsv2_solve(cusparseHandle_t handle, + cusparseOperation_t transA, + int m, + int nnz, + const Number * alpha, + const cusparseMatDescr_t descra, + const Number * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csrsv2Info_t info, + const Number * x, + Number * y, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + AssertThrow(false, ExcNotImplemented()); + return CUSPARSE_STATUS_INVALID_VALUE; + } + + template <> + cusparseStatus_t + cusparseXcsrsv2_solve(cusparseHandle_t handle, + cusparseOperation_t transA, + int m, + int nnz, + const float * alpha, + const cusparseMatDescr_t descra, + const float * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csrsv2Info_t info, + const float * x, + float * y, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + return cusparseScsrsv2_solve(handle, + transA, + m, + nnz, + alpha, + descra, + csrValA, + csrRowPtrA, + csrColIndA, + info, + x, + y, + policy, + pBuffer); + } + + template <> + cusparseStatus_t + cusparseXcsrsv2_solve(cusparseHandle_t handle, + cusparseOperation_t transA, + int m, + int nnz, + const double * alpha, + const cusparseMatDescr_t descra, + const double * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csrsv2Info_t info, + const double * x, + double * y, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + return cusparseDcsrsv2_solve(handle, + transA, + m, + nnz, + alpha, + descra, + csrValA, + csrRowPtrA, + csrColIndA, + info, + x, + y, + policy, + pBuffer); + } + + template <> + cusparseStatus_t + cusparseXcsrsv2_solve(cusparseHandle_t handle, + cusparseOperation_t transA, + int m, + int nnz, + const cuComplex * alpha, + const cusparseMatDescr_t descra, + const cuComplex * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csrsv2Info_t info, + const cuComplex * x, + cuComplex * y, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + return cusparseCcsrsv2_solve(handle, + transA, + m, + nnz, + alpha, + descra, + csrValA, + csrRowPtrA, + csrColIndA, + info, + x, + y, + policy, + pBuffer); + } + + template <> + cusparseStatus_t + cusparseXcsrsv2_solve(cusparseHandle_t handle, + cusparseOperation_t transA, + int m, + int nnz, + const cuDoubleComplex * alpha, + const cusparseMatDescr_t descra, + const cuDoubleComplex * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csrsv2Info_t info, + const cuDoubleComplex * x, + cuDoubleComplex * y, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + return cusparseZcsrsv2_solve(handle, + transA, + m, + nnz, + alpha, + descra, + csrValA, + csrRowPtrA, + csrColIndA, + info, + x, + y, + policy, + pBuffer); + } + + + /** + * Template wrapper for cusparsecsrsv2_analysis + * (https://docs.nvidia.com/cuda/cusparse/index.html#cusparse-lt-t-gt-csrsv2_analysis). + * This function performs the analysis phase of csrsv2, a new sparse + * triangular linear system op(A)*y = alpha*x. + */ + template + cusparseStatus_t + cusparseXcsrsv2_analysis(cusparseHandle_t handle, + cusparseOperation_t transA, + int m, + int nnz, + const cusparseMatDescr_t descrA, + const Number * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csrsv2Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + AssertThrow(false, ExcNotImplemented()); + return CUSPARSE_STATUS_INVALID_VALUE; + } + + template <> + cusparseStatus_t + cusparseXcsrsv2_analysis(cusparseHandle_t handle, + cusparseOperation_t transA, + int m, + int nnz, + const cusparseMatDescr_t descrA, + const float * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csrsv2Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + return cusparseScsrsv2_analysis(handle, + transA, + m, + nnz, + descrA, + csrValA, + csrRowPtrA, + csrColIndA, + info, + policy, + pBuffer); + } + + template <> + cusparseStatus_t + cusparseXcsrsv2_analysis(cusparseHandle_t handle, + cusparseOperation_t transA, + int m, + int nnz, + const cusparseMatDescr_t descrA, + const double * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csrsv2Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + return cusparseDcsrsv2_analysis(handle, + transA, + m, + nnz, + descrA, + csrValA, + csrRowPtrA, + csrColIndA, + info, + policy, + pBuffer); + } + + template <> + cusparseStatus_t + cusparseXcsrsv2_analysis(cusparseHandle_t handle, + cusparseOperation_t transA, + int m, + int nnz, + const cusparseMatDescr_t descrA, + const cuComplex * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csrsv2Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + return cusparseCcsrsv2_analysis(handle, + transA, + m, + nnz, + descrA, + csrValA, + csrRowPtrA, + csrColIndA, + info, + policy, + pBuffer); + } + + template <> + cusparseStatus_t + cusparseXcsrsv2_analysis(cusparseHandle_t handle, + cusparseOperation_t transA, + int m, + int nnz, + const cusparseMatDescr_t descrA, + const cuDoubleComplex * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csrsv2Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + return cusparseZcsrsv2_analysis(handle, + transA, + m, + nnz, + descrA, + csrValA, + csrRowPtrA, + csrColIndA, + info, + policy, + pBuffer); + } + + + + /** + * Template wrapper for cusparsecsric02_analysis + * (https://docs.nvidia.com/cuda/cusparse/index.html#cusparse-lt-t-gt-csric02_analysis). + * This function performs the analysis phase of the incomplete-Cholesky + * factorization with 0 fill-in and no pivoting. + */ + template + cusparseStatus_t + cusparseXcsric02_analysis(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + const Number * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csric02Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + AssertThrow(false, ExcNotImplemented()); + return CUSPARSE_STATUS_INVALID_VALUE; + } + + template <> + cusparseStatus_t + cusparseXcsric02_analysis(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + const float * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csric02Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + return cusparseScsric02_analysis(handle, + m, + nnz, + descrA, + csrValA, + csrRowPtrA, + csrColIndA, + info, + policy, + pBuffer); + } + + template <> + cusparseStatus_t + cusparseXcsric02_analysis(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + const double * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csric02Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + return cusparseDcsric02_analysis(handle, + m, + nnz, + descrA, + csrValA, + csrRowPtrA, + csrColIndA, + info, + policy, + pBuffer); + } + + template <> + cusparseStatus_t + cusparseXcsric02_analysis(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + const cuComplex * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csric02Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + return cusparseCcsric02_analysis(handle, + m, + nnz, + descrA, + csrValA, + csrRowPtrA, + csrColIndA, + info, + policy, + pBuffer); + } + + template <> + cusparseStatus_t + cusparseXcsric02_analysis(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + const cuDoubleComplex * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csric02Info_t info, + cusparseSolvePolicy_t policy, + void * pBuffer) + { + return cusparseZcsric02_analysis(handle, + m, + nnz, + descrA, + csrValA, + csrRowPtrA, + csrColIndA, + info, + policy, + pBuffer); + } + + + /** + * Template wrapper for cusparsecsrsv2_bufferSize + * (https://docs.nvidia.com/cuda/cusparse/index.html#cusparse-lt-t-gt-csrsv2_bufferSize). + * This function returns the size of the buffer used in csrsv2, a new sparse + * triangular linear system op(A)*y = alpha*x. + */ + template + cusparseStatus_t + cusparseXcsrsv2_bufferSize(cusparseHandle_t handle, + cusparseOperation_t transA, + int m, + int nnz, + const cusparseMatDescr_t descrA, + Number * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csrsv2Info_t info, + int * pBufferSizeInBytes) + { + AssertThrow(false, ExcNotImplemented()); + return CUSPARSE_STATUS_INVALID_VALUE; + } + + template <> + cusparseStatus_t + cusparseXcsrsv2_bufferSize(cusparseHandle_t handle, + cusparseOperation_t transA, + int m, + int nnz, + const cusparseMatDescr_t descrA, + float * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csrsv2Info_t info, + int * pBufferSizeInBytes) + { + return cusparseScsrsv2_bufferSize(handle, + transA, + m, + nnz, + descrA, + csrValA, + csrRowPtrA, + csrColIndA, + info, + pBufferSizeInBytes); + } + + template <> + cusparseStatus_t + cusparseXcsrsv2_bufferSize(cusparseHandle_t handle, + cusparseOperation_t transA, + int m, + int nnz, + const cusparseMatDescr_t descrA, + double * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csrsv2Info_t info, + int *pBufferSizeInBytes) + { + return cusparseDcsrsv2_bufferSize(handle, + transA, + m, + nnz, + descrA, + csrValA, + csrRowPtrA, + csrColIndA, + info, + pBufferSizeInBytes); + } + + template <> + cusparseStatus_t + cusparseXcsrsv2_bufferSize(cusparseHandle_t handle, + cusparseOperation_t transA, + int m, + int nnz, + const cusparseMatDescr_t descrA, + cuComplex * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csrsv2Info_t info, + int *pBufferSizeInBytes) + { + return cusparseCcsrsv2_bufferSize(handle, + transA, + m, + nnz, + descrA, + csrValA, + csrRowPtrA, + csrColIndA, + info, + pBufferSizeInBytes); + } + + template <> + cusparseStatus_t + cusparseXcsrsv2_bufferSize(cusparseHandle_t handle, + cusparseOperation_t transA, + int m, + int nnz, + const cusparseMatDescr_t descrA, + cuDoubleComplex * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csrsv2Info_t info, + int * pBufferSizeInBytes) + { + return cusparseZcsrsv2_bufferSize(handle, + transA, + m, + nnz, + descrA, + csrValA, + csrRowPtrA, + csrColIndA, + info, + pBufferSizeInBytes); + } + + + + /** + * Template wrapper for cusparsecsric02_bufferSize + * (https://docs.nvidia.com/cuda/cusparse/index.html#cusparse-lt-t-gt-csric02_bufferSize). + *This function returns size of buffer used in computing the + *incomplete-Cholesky factorization with 0 fill-in and no pivoting. + */ + template + cusparseStatus_t + cusparseXcsric02_bufferSize(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + Number * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csric02Info_t info, + int * pBufferSizeInBytes) + { + AssertThrow(false, ExcNotImplemented()); + return CUSPARSE_STATUS_INVALID_VALUE; + } + + template <> + cusparseStatus_t + cusparseXcsric02_bufferSize(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + float * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csric02Info_t info, + int *pBufferSizeInBytes) + { + return cusparseScsric02_bufferSize(handle, + m, + nnz, + descrA, + csrValA, + csrRowPtrA, + csrColIndA, + info, + pBufferSizeInBytes); + } + + template <> + cusparseStatus_t + cusparseXcsric02_bufferSize(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + double * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csric02Info_t info, + int *pBufferSizeInBytes) + { + return cusparseDcsric02_bufferSize(handle, + m, + nnz, + descrA, + csrValA, + csrRowPtrA, + csrColIndA, + info, + pBufferSizeInBytes); + } + + template <> + cusparseStatus_t + cusparseXcsric02_bufferSize(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + cuComplex * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csric02Info_t info, + int *pBufferSizeInBytes) + { + return cusparseCcsric02_bufferSize(handle, + m, + nnz, + descrA, + csrValA, + csrRowPtrA, + csrColIndA, + info, + pBufferSizeInBytes); + } + + template <> + cusparseStatus_t + cusparseXcsric02_bufferSize(cusparseHandle_t handle, + int m, + int nnz, + const cusparseMatDescr_t descrA, + cuDoubleComplex * csrValA, + const int * csrRowPtrA, + const int * csrColIndA, + csric02Info_t info, + int * pBufferSizeInBytes) + { + return cusparseZcsric02_bufferSize(handle, + m, + nnz, + descrA, + csrValA, + csrRowPtrA, + csrColIndA, + info, + pBufferSizeInBytes); + } + /** + * @} + */ +} // namespace CUDAWrappers + +DEAL_II_NAMESPACE_CLOSE + +namespace +{ + template + void + delete_device_vector(Number *device_ptr) noexcept + { + const cudaError_t error_code = cudaFree(device_ptr); + (void)error_code; + AssertNothrow(error_code == cudaSuccess, + dealii::ExcCudaError(cudaGetErrorString(error_code))); + } + template + Number * + allocate_device_vector(const std::size_t size) + { + Number *device_ptr; + Utilities::CUDA::malloc(device_ptr, size); + return device_ptr; + } +} // namespace + +namespace dealii +{ + namespace CUDAWrappers + { + /** + * This class implements an incomplete Cholesky factorization (IC) + * preconditioner for @em symmetric CUDAWrappers::SparseMatrix matrices. + * + * The implementation closely follows the one documented in the cuSPARSE + * documentation + * (https://docs.nvidia.com/cuda/cusparse/index.html#cusparse-lt-t-gt-csrilu02). + * + * @note Instantiations for this template are provided for @ and + * @. + * + * @ingroup Preconditioners CUDAWrappers + * @author Daniel Arndt + * @date 2018 + */ + template + class PreconditionILU + { + public: + /** + * Declare the type for container size. + */ + using size_type = int; + + /** + * Standardized data struct to pipe additional flags to the + * preconditioner. + */ + struct AdditionalData + { + /** + * Constructor. cuSPARSE allows to compute and use level information. + * According to the documentation it is this might improve performance. + * It is suggested to try both options. + */ + AdditionalData(bool use_level_analysis = true); + + /** + * Flag that determines if level informations are used when creating and + * applying the preconditioner. See the documentation for + * cusparseSolvePolicy_t at + * https://docs.nvidia.com/cuda/cusparse/index.html#cusparsesolvepolicy_t + * for more information. + */ + bool use_level_analysis; + }; + + /** + * Constructor. + */ + PreconditionILU(const Utilities::CUDA::Handle &handle); + + /** + * The copy constructor is deleted. + */ + PreconditionILU(const PreconditionILU &) = delete; + + /** + * The copy assignment operator is deleted. + */ + PreconditionILU & + operator=(const PreconditionILU &) = delete; + + /** + * Destructor. Free all resources that were initialized in this class. + */ + ~PreconditionILU(); + + /** + * Initialize this object. In particular, the given matrix is copied to be + * modified in-place. For the underlying sparsity pattern pointers are + * stored. Specifically, this means + * that the current object can only be used reliably as long as @p matrix is valid + * and has not been changed since calling this function. + * + * The @p additional_data determines if level information are used. + */ + void + initialize(const SparseMatrix &matrix, + const AdditionalData &additional_data = AdditionalData()); + + /** + * Apply the preconditioner. + */ + void + vmult(LinearAlgebra::CUDAWrappers::Vector & dst, + const LinearAlgebra::CUDAWrappers::Vector &src) const; + + /** + * Apply the preconditioner. Since the preconditioner is symmetric, this + * is the same as vmult(). + */ + void + Tvmult(LinearAlgebra::CUDAWrappers::Vector & dst, + const LinearAlgebra::CUDAWrappers::Vector &src) const; + + /** + * Return the dimension of the codomain (or range) space. Note that the + * matrix is square and has dimension $m \times m$. + * + * @note This function should only be called if the preconditioner has been + * initialized. + */ + size_type + m() const; + + /** + * Return the dimension of the codomain (or range) space. Note that the + * matrix is square and has dimension $m \times m$. + * + * @note This function should only be called if the preconditioner has been + * initialized. + */ + size_type + n() const; + + private: + /** + * cuSPARSE handle used to call cuSPARSE functions. + */ + cusparseHandle_t cusparse_handle; + + /** + * cuSPARSE description of the sparse matrix $M=LU$. + */ + cusparseMatDescr_t descr_M; + + /** + * cuSPARSE description of the lower triangular matrix $L$. + */ + cusparseMatDescr_t descr_L; + + /** + * cuSPARSE description of the upper triangular matrix $U$. + */ + cusparseMatDescr_t descr_U; + + /** + * Solve and analysis structure for $M=LL^T$. + */ + csrilu02Info_t info_M; + + /** + * Solve and analysis structure for the lower triangular matrix $L$. + */ + csrsv2Info_t info_L; + + /** + * Solve and analysis structure for the upper triangular matrix $U$. + */ + csrsv2Info_t info_U; + + /** + * Pointer to the values (on the device) of the computed preconditioning + * matrix. + */ + std::unique_ptr P_val_dev; + + /** + * Pointer to the row pointer (on the device) of the sparse matrix this + * object was initialized with. + */ + const int *P_row_ptr_dev; + + /** + * Pointer to the column indices (on the device) of the sparse matrix this + * object was initialized with. + */ + const int *P_column_index_dev; + + /** + * Pointer to the value (on the device) for a temporary (helper) vector + * used in vmult(). + */ + std::unique_ptr tmp_dev; + + /** + * + */ + std::unique_ptr buffer_dev; + + /** + * Determine if level information should be generated for the lower + * triangular matrix $L$. This value can be modified through an + * AdditionalData object. + */ + cusparseSolvePolicy_t policy_L; + + /** + * Determine if level information should be generated for the upper + * triangular matrix $L^T$. This value can be modified through an + * AdditionalData object. + */ + cusparseSolvePolicy_t policy_U; + + /** + * Determine if level information should be generated for $M=LL^T$. This + * value can be modified through an AdditionalData object. + */ + cusparseSolvePolicy_t policy_M; + + /** + * The number of rows is the same as for the matrix this object has been + * initialized with. + */ + int n_rows; + + /** + * The number of non-zero elements is the same as for the matrix this + * object has been initialized with. + */ + int n_nonzero_elements; + }; + + template + PreconditionILU::AdditionalData::AdditionalData( + bool use_level_analysis_) + : use_level_analysis(use_level_analysis_) + {} + + + + template + PreconditionILU::PreconditionILU( + const Utilities::CUDA::Handle &handle) + : cusparse_handle(handle.cusparse_handle) + , P_val_dev(nullptr, delete_device_vector) + , P_row_ptr_dev(nullptr) + , P_column_index_dev(nullptr) + , tmp_dev(nullptr, delete_device_vector) + , buffer_dev(nullptr, delete_device_vector) + , policy_L(CUSPARSE_SOLVE_POLICY_USE_LEVEL) + , policy_U(CUSPARSE_SOLVE_POLICY_USE_LEVEL) + , policy_M(CUSPARSE_SOLVE_POLICY_USE_LEVEL) + , n_rows(0) + , n_nonzero_elements(0) + { + cusparseStatus_t status; + // step 1: create a descriptor which contains + // - matrix M is base-0 + // - matrix L is base-0 + // - matrix L is lower triangular + // - matrix L has unit diagonal + // - matrix U is base-0 + // - matrix U is upper triangular + // - matrix U has non-unit diagonal + status = cusparseCreateMatDescr(&descr_M); + AssertCusparse(status); + status = cusparseSetMatIndexBase(descr_M, CUSPARSE_INDEX_BASE_ZERO); + AssertCusparse(status); + status = cusparseSetMatType(descr_M, CUSPARSE_MATRIX_TYPE_GENERAL); + AssertCusparse(status); + + status = cusparseCreateMatDescr(&descr_L); + AssertCusparse(status); + status = cusparseSetMatIndexBase(descr_L, CUSPARSE_INDEX_BASE_ZERO); + AssertCusparse(status); + status = cusparseSetMatType(descr_L, CUSPARSE_MATRIX_TYPE_GENERAL); + AssertCusparse(status); + status = cusparseSetMatFillMode(descr_L, CUSPARSE_FILL_MODE_LOWER); + AssertCusparse(status); + status = cusparseSetMatDiagType(descr_L, CUSPARSE_DIAG_TYPE_UNIT); + AssertCusparse(status); + + status = cusparseCreateMatDescr(&descr_U); + AssertCusparse(status); + status = cusparseSetMatIndexBase(descr_U, CUSPARSE_INDEX_BASE_ZERO); + AssertCusparse(status); + status = cusparseSetMatType(descr_U, CUSPARSE_MATRIX_TYPE_GENERAL); + AssertCusparse(status); + status = cusparseSetMatFillMode(descr_U, CUSPARSE_FILL_MODE_UPPER); + AssertCusparse(status); + status = cusparseSetMatDiagType(descr_U, CUSPARSE_DIAG_TYPE_NON_UNIT); + AssertCusparse(status); + + // step 2: create a empty info structure + // we need one info for csrilu02 and two info's for csrsv2 + status = cusparseCreateCsrilu02Info(&info_M); + AssertCusparse(status); + status = cusparseCreateCsrsv2Info(&info_L); + AssertCusparse(status); + status = cusparseCreateCsrsv2Info(&info_U); + AssertCusparse(status); + } + + template + PreconditionILU::~PreconditionILU() + { + // step 8: free resources + cusparseStatus_t status = cusparseDestroyMatDescr(descr_M); + AssertNothrowCusparse(status); + + status = cusparseDestroyMatDescr(descr_L); + AssertNothrowCusparse(status); + + status = cusparseDestroyMatDescr(descr_U); + AssertNothrowCusparse(status); + + status = cusparseDestroyCsrilu02Info(info_M); + AssertNothrowCusparse(status); + + status = cusparseDestroyCsrsv2Info(info_L); + AssertNothrowCusparse(status); + + status = cusparseDestroyCsrsv2Info(info_U); + AssertNothrowCusparse(status); + } + + + + template + void + PreconditionILU::initialize(const SparseMatrix &A, + const AdditionalData &additional_data) + { + if (additional_data.use_level_analysis) + { + policy_L = CUSPARSE_SOLVE_POLICY_USE_LEVEL; + policy_U = CUSPARSE_SOLVE_POLICY_USE_LEVEL; + policy_M = CUSPARSE_SOLVE_POLICY_USE_LEVEL; + } + else + { + policy_L = CUSPARSE_SOLVE_POLICY_NO_LEVEL; + policy_U = CUSPARSE_SOLVE_POLICY_NO_LEVEL; + policy_M = CUSPARSE_SOLVE_POLICY_NO_LEVEL; + } + + + n_rows = A.m(); + n_nonzero_elements = A.n_nonzero_elements(); + AssertDimension(A.m(), A.n()); + + const auto cusparse_matrix = A.get_cusparse_matrix(); + const Number *const A_val_dev = std::get<0>(cusparse_matrix); + + // create a copy of the matrix entries + P_val_dev.reset(allocate_device_vector(n_nonzero_elements)); + cudaError_t cuda_status = cudaMemcpy(P_val_dev.get(), + A_val_dev, + n_nonzero_elements * sizeof(Number), + cudaMemcpyDeviceToDevice); + P_column_index_dev = std::get<1>(cusparse_matrix); + P_row_ptr_dev = std::get<2>(cusparse_matrix); + const cusparseMatDescr_t mat_descr = std::get<3>(cusparse_matrix); + + // initializa an internal buffer we need later on + tmp_dev.reset(allocate_device_vector(n_rows)); + + // step 3: query how much memory used in csrilu02 and csrsv2, and allocate + // the buffer + int BufferSize_M; + cusparseStatus_t status = cusparseXcsrilu02_bufferSize(cusparse_handle, + n_rows, + n_nonzero_elements, + descr_M, + P_val_dev.get(), + P_row_ptr_dev, + P_column_index_dev, + info_M, + &BufferSize_M); + AssertCusparse(status); + + int BufferSize_L; + status = cusparseXcsrsv2_bufferSize(cusparse_handle, + CUSPARSE_OPERATION_NON_TRANSPOSE, + n_rows, + n_nonzero_elements, + descr_L, + P_val_dev.get(), + P_row_ptr_dev, + P_column_index_dev, + info_L, + &BufferSize_L); + AssertCusparse(status); + + int BufferSize_U; + status = cusparseXcsrsv2_bufferSize(cusparse_handle, + CUSPARSE_OPERATION_NON_TRANSPOSE, + n_rows, + n_nonzero_elements, + descr_U, + P_val_dev.get(), + P_row_ptr_dev, + P_column_index_dev, + info_U, + &BufferSize_U); + AssertCusparse(status); + + const int BufferSize = + std::max(BufferSize_M, std::max(BufferSize_L, BufferSize_U)); + // workaround: since allocate_device_vector needs a type, we pass char + // which is required to have size 1. + buffer_dev.reset(static_cast( + allocate_device_vector(BufferSize / sizeof(char)))); + + // step 4: perform analysis of incomplete Cholesky on M + // perform analysis of triangular solve on L + // perform analysis of triangular solve on U + // The lower(upper) triangular part of M has the same sparsity pattern as + // L(U), we can do analysis of csrilu0 and csrsv2 simultaneously. + + status = cusparseXcsrilu02_analysis(cusparse_handle, + n_rows, + n_nonzero_elements, + descr_M, + P_val_dev.get(), + P_row_ptr_dev, + P_column_index_dev, + info_M, + policy_M, + buffer_dev.get()); + AssertCusparse(status); + + int structural_zero; + status = + cusparseXcsrilu02_zeroPivot(cusparse_handle, info_M, &structural_zero); + AssertCusparse(status); + + status = cusparseXcsrsv2_analysis(cusparse_handle, + CUSPARSE_OPERATION_NON_TRANSPOSE, + n_rows, + n_nonzero_elements, + descr_L, + P_val_dev.get(), + P_row_ptr_dev, + P_column_index_dev, + info_L, + policy_L, + buffer_dev.get()); + AssertCusparse(status); + + status = cusparseXcsrsv2_analysis(cusparse_handle, + CUSPARSE_OPERATION_NON_TRANSPOSE, + n_rows, + n_nonzero_elements, + descr_U, + P_val_dev.get(), + P_row_ptr_dev, + P_column_index_dev, + info_U, + policy_U, + buffer_dev.get()); + + // step 5: M = L * U + status = cusparseXcsrilu02(cusparse_handle, + n_rows, + n_nonzero_elements, + descr_M, + P_val_dev.get(), + P_row_ptr_dev, + P_column_index_dev, + info_M, + policy_M, + buffer_dev.get()); + AssertCusparse(status); + + int numerical_zero; + status = + cusparseXcsrilu02_zeroPivot(cusparse_handle, info_M, &numerical_zero); + AssertCusparse(status); + } + + + + template + void + PreconditionILU::vmult( + LinearAlgebra::CUDAWrappers::Vector & dst, + const LinearAlgebra::CUDAWrappers::Vector &src) const + { + Assert(P_val_dev != nullptr, ExcNotInitialized()); + Assert(P_row_ptr_dev != nullptr, ExcNotInitialized()); + Assert(P_column_index_dev != nullptr, ExcNotInitialized()); + AssertDimension(dst.size(), static_cast(n_rows)); + AssertDimension(src.size(), static_cast(n_rows)); + Assert(tmp_dev != nullptr, ExcInternalError()); + + const Number *const src_dev = src.get_values(); + Number *const dst_dev = dst.get_values(); + + // step 6: solve L*z = alpha*x + const Number alpha = 1.; + cusparseStatus_t status = + cusparseXcsrsv2_solve(cusparse_handle, + CUSPARSE_OPERATION_NON_TRANSPOSE, + n_rows, + n_nonzero_elements, + &alpha, + descr_L, + P_val_dev.get(), + P_row_ptr_dev, + P_column_index_dev, + info_L, + src_dev, + tmp_dev.get(), + policy_L, + buffer_dev.get()); + AssertCusparse(status); + + // step 7: solve U*y = alpha*z + status = cusparseXcsrsv2_solve(cusparse_handle, + CUSPARSE_OPERATION_NON_TRANSPOSE, + n_rows, + n_nonzero_elements, + &alpha, + descr_U, + P_val_dev.get(), + P_row_ptr_dev, + P_column_index_dev, + info_U, + tmp_dev.get(), + dst_dev, + policy_U, + buffer_dev.get()); + AssertCusparse(status); + } + + + + template + void + PreconditionILU::Tvmult( + LinearAlgebra::CUDAWrappers::Vector & dst, + const LinearAlgebra::CUDAWrappers::Vector &src) const + { + // the constructed preconditioner is symmetric + vmult(dst, src); + } + + template + PreconditionILU::size_type + PreconditionILU::m() const + { + return n_rows; + } + + + template + PreconditionILU::size_type + PreconditionILU::n() const + { + return n_rows; + } + + + + template + void + apply_preconditioner(const SparseMatrix &A, + const cusparseHandle_t cusparse_handle, + LinearAlgebra::CUDAWrappers::Vector & dst, + const LinearAlgebra::CUDAWrappers::Vector &src) + { + const Number *const src_dev = src.get_values(); + Number * dst_dev = dst.get_values(); + const cusparseHandle_t handle = cusparse_handle; + + const auto cusparse_matrix = A.get_cusparse_matrix(); + Number * A_val_dev = std::get<0>(cusparse_matrix); + const int *const A_row_ptr_dev = std::get<2>(cusparse_matrix); + const int *const A_column_index_dev = std::get<1>(cusparse_matrix); + const cusparseMatDescr_t mat_descr = std::get<3>(cusparse_matrix); + + const unsigned int n_rows = A.m(); + const unsigned int n_nonzero_elements = A.n_nonzero_elements(); + + AssertDimension(dst.size(), src.size()); + AssertDimension(A.m(), src.size()); + AssertDimension(A.n(), src.size()); + + std::unique_ptr tmp_dev( + allocate_device_vector(dst.size()), + delete_device_vector); + + // Suppose that A is a m x m sparse matrix represented by CSR format, + // Assumption: + // - handle is already created by cusparseCreate(), + // - (A_row_ptr_dev, A_column_index_dev, A_val_dev) is CSR of A on device + // memory, + // - src_dev is right hand side vector on device memory, + // - dst_dev is solution vector on device memory. + // - tmp_dev is intermediate result on device memory. + + cusparseMatDescr_t descr_M = mat_descr; + cusparseMatDescr_t descr_L = mat_descr; + cusparseMatDescr_t descr_U = mat_descr; + csrilu02Info_t info_M = 0; + csrsv2Info_t info_L = 0; + csrsv2Info_t info_U = 0; + int BufferSize_M; + int BufferSize_L; + int BufferSize_U; + int BufferSize; + void * buffer_dev = 0; + int structural_zero; + int numerical_zero; + const double alpha = 1.; + const cusparseSolvePolicy_t policy_M = CUSPARSE_SOLVE_POLICY_NO_LEVEL; + const cusparseSolvePolicy_t policy_L = CUSPARSE_SOLVE_POLICY_NO_LEVEL; + const cusparseSolvePolicy_t policy_U = CUSPARSE_SOLVE_POLICY_USE_LEVEL; + + // step 1: create a descriptor which contains + // - matrix M is base-0 + // - matrix L is base-0 + // - matrix L is lower triangular + // - matrix L has unit diagonal + // - matrix U is base-0 + // - matrix U is upper triangular + // - matrix U has non-unit diagonal + cusparseStatus_t status = cusparseCreateMatDescr(&descr_M); + AssertCusparse(status); + status = cusparseSetMatIndexBase(descr_M, CUSPARSE_INDEX_BASE_ZERO); + AssertCusparse(status); + status = cusparseSetMatType(descr_M, CUSPARSE_MATRIX_TYPE_GENERAL); + AssertCusparse(status); + + status = cusparseCreateMatDescr(&descr_L); + AssertCusparse(status); + status = cusparseSetMatIndexBase(descr_L, CUSPARSE_INDEX_BASE_ZERO); + AssertCusparse(status); + status = cusparseSetMatType(descr_L, CUSPARSE_MATRIX_TYPE_GENERAL); + AssertCusparse(status); + status = cusparseSetMatFillMode(descr_L, CUSPARSE_FILL_MODE_LOWER); + AssertCusparse(status); + status = cusparseSetMatDiagType(descr_L, CUSPARSE_DIAG_TYPE_UNIT); + AssertCusparse(status); + + status = cusparseCreateMatDescr(&descr_U); + AssertCusparse(status); + status = cusparseSetMatIndexBase(descr_U, CUSPARSE_INDEX_BASE_ZERO); + AssertCusparse(status); + status = cusparseSetMatType(descr_U, CUSPARSE_MATRIX_TYPE_GENERAL); + AssertCusparse(status); + status = cusparseSetMatFillMode(descr_U, CUSPARSE_FILL_MODE_UPPER); + AssertCusparse(status); + status = cusparseSetMatDiagType(descr_U, CUSPARSE_DIAG_TYPE_NON_UNIT); + AssertCusparse(status); + + // step 2: create a empty info structure + // we need one info for csrilu02 and two info's for csrsv2 + status = cusparseCreateCsrilu02Info(&info_M); + AssertCusparse(status); + status = cusparseCreateCsrsv2Info(&info_L); + AssertCusparse(status); + status = cusparseCreateCsrsv2Info(&info_U); + AssertCusparse(status); + + // step 3: query how much memory used in csrilu02 and csrsv2, and allocate + // the buffer + status = cusparseXcsrilu02_bufferSize(handle, + n_rows, + n_nonzero_elements, + descr_M, + A_val_dev, + A_row_ptr_dev, + A_column_index_dev, + info_M, + &BufferSize_M); + AssertCusparse(status); + + status = cusparseXcsrsv2_bufferSize(handle, + CUSPARSE_OPERATION_NON_TRANSPOSE, + n_rows, + n_nonzero_elements, + descr_L, + A_val_dev, + A_row_ptr_dev, + A_column_index_dev, + info_L, + &BufferSize_L); + AssertCusparse(status); + + status = cusparseXcsrsv2_bufferSize(handle, + CUSPARSE_OPERATION_NON_TRANSPOSE, + n_rows, + n_nonzero_elements, + descr_U, + A_val_dev, + A_row_ptr_dev, + A_column_index_dev, + info_U, + &BufferSize_U); + AssertCusparse(status); + + BufferSize = max(BufferSize_M, max(BufferSize_L, BufferSize_U)); + + // Buffer returned by cudaMalloc is automatically aligned to 128 bytes. + cudaMalloc((void **)&buffer_dev, BufferSize); + + // step 4: perform analysis of incomplete Cholesky on M + // perform analysis of triangular solve on L + // perform analysis of triangular solve on U + // The lower(upper) triangular part of M has the same sparsity pattern as + // L(U), we can do analysis of csrilu0 and csrsv2 simultaneously. + + status = cusparseXcsrilu02_analysis(handle, + n_rows, + n_nonzero_elements, + descr_M, + A_val_dev, + A_row_ptr_dev, + A_column_index_dev, + info_M, + policy_M, + buffer_dev); + status = cusparseXcsrilu02_zeroPivot(handle, info_M, &structural_zero); + AssertCusparse(status); + if (CUSPARSE_STATUS_ZERO_PIVOT == status) + { + printf("A(%d,%d) is missing\n", structural_zero, structural_zero); + } + + status = cusparseXcsrsv2_analysis(handle, + CUSPARSE_OPERATION_NON_TRANSPOSE, + n_rows, + n_nonzero_elements, + descr_L, + A_val_dev, + A_row_ptr_dev, + A_column_index_dev, + info_L, + policy_L, + buffer_dev); + AssertCusparse(status); + + status = cusparseXcsrsv2_analysis(handle, + CUSPARSE_OPERATION_NON_TRANSPOSE, + n_rows, + n_nonzero_elements, + descr_U, + A_val_dev, + A_row_ptr_dev, + A_column_index_dev, + info_U, + policy_U, + buffer_dev); + AssertCusparse(status); + + // step 5: M = L * U + status = cusparseXcsrilu02(handle, + n_rows, + n_nonzero_elements, + descr_M, + A_val_dev, + A_row_ptr_dev, + A_column_index_dev, + info_M, + policy_M, + buffer_dev); + status = cusparseXcsrilu02_zeroPivot(handle, info_M, &numerical_zero); + AssertCusparse(status); + if (CUSPARSE_STATUS_ZERO_PIVOT == status) + { + printf("U(%d,%d) is zero\n", numerical_zero, numerical_zero); + } + + // step 6: solve L*z = x + status = cusparseXcsrsv2_solve(handle, + CUSPARSE_OPERATION_NON_TRANSPOSE, + n_rows, + n_nonzero_elements, + &alpha, + descr_L, + A_val_dev, + A_row_ptr_dev, + A_column_index_dev, + info_L, + src_dev, + tmp_dev.get(), + policy_L, + buffer_dev); + AssertCusparse(status); + + // step 7: solve U*y = z + status = cusparseXcsrsv2_solve(handle, + CUSPARSE_OPERATION_NON_TRANSPOSE, + n_rows, + n_nonzero_elements, + &alpha, + descr_U, + A_val_dev, + A_row_ptr_dev, + A_column_index_dev, + info_U, + tmp_dev.get(), + dst_dev, + policy_U, + buffer_dev); + AssertCusparse(status); + + // step 8: free resources + cudaFree(buffer_dev); + status = cusparseDestroyMatDescr(descr_M); + AssertCusparse(status); + status = cusparseDestroyMatDescr(descr_L); + AssertCusparse(status); + status = cusparseDestroyMatDescr(descr_U); + AssertCusparse(status); + status = cusparseDestroyCsrilu02Info(info_M); + AssertCusparse(status); + status = cusparseDestroyCsrsv2Info(info_L); + AssertCusparse(status); + status = cusparseDestroyCsrsv2Info(info_U); + AssertCusparse(status); + } + } // namespace CUDAWrappers +} // namespace dealii + +void +test(Utilities::CUDA::Handle &cuda_handle) +{ + // Build the sparse matrix on the host + const unsigned int problem_size = 10; + unsigned int size = (problem_size - 1) * (problem_size - 1); + FDMatrix testproblem(problem_size, problem_size); + SparsityPattern structure(size, size, 5); + SparseMatrix A; + testproblem.five_point_structure(structure); + structure.compress(); + A.reinit(structure); + testproblem.five_point(A); + A.print(std::cout); + + // Solve on the host + PreconditionIdentity prec_no; + SolverControl control(100, 1.e-10); + SolverCG<> cg_host(control); + Vector sol_host(size); + Vector rhs_host(size); + for (unsigned int i = 0; i < size; ++i) + rhs_host[i] = static_cast(i); + cg_host.solve(A, sol_host, rhs_host, prec_no); + + // Solve on the device + CUDAWrappers::SparseMatrix A_dev(cuda_handle, A); + LinearAlgebra::CUDAWrappers::Vector sol_dev(size); + LinearAlgebra::CUDAWrappers::Vector rhs_dev(size); + LinearAlgebra::ReadWriteVector rw_vector(size); + for (unsigned int i = 0; i < size; ++i) + rw_vector[i] = static_cast(i); + rhs_dev.import(rw_vector, VectorOperation::insert); + SolverCG> cg_dev(control); + + A_dev.print(std::cout); + A_dev.print_formatted(std::cout); + CUDAWrappers::PreconditionILU prec_double(cuda_handle); + CUDAWrappers::PreconditionILU prec_float(cuda_handle); + CUDAWrappers::PreconditionILU prec_complex_float(cuda_handle); + CUDAWrappers::PreconditionILU prec_complex_double( + cuda_handle); + + // apply_preconditioner(A_dev, cuda_handle.cusparse_handle, sol_dev, rhs_dev); + // A_dev.print_formatted(std::cout); + prec_double.initialize(A_dev); + // A_dev.print_formatted(std::cout); + // prec_double.vmult(sol_dev, rhs_dev); + // A_dev.print_formatted(std::cout); + cg_dev.solve(A_dev, sol_dev, rhs_dev, prec_double); + + // Check the result + rw_vector.import(sol_dev, VectorOperation::insert); + for (unsigned int i = 0; i < size; ++i) + deallog << rw_vector[i] << " " << sol_host[i] << std::endl; +} + +int +main() +{ + initlog(); + deallog.depth_console(0); + + Utilities::CUDA::Handle cuda_handle; + test(cuda_handle); + + deallog << "OK" << std::endl; + + return 0; +} diff --git a/tests/cuda/precondition_02.output b/tests/cuda/precondition_02.output new file mode 100644 index 0000000000..d6d49d1ae1 --- /dev/null +++ b/tests/cuda/precondition_02.output @@ -0,0 +1,87 @@ + +DEAL:cg::Starting value 416.989 +DEAL:cg::Convergence step 31 value 8.71925e-12 +DEAL:cg::Starting value 416.989 +DEAL:cg::Convergence step 17 value 5.54040e-11 +DEAL::20.9607 20.9607 +DEAL::38.8073 38.8073 +DEAL::52.3525 52.3525 +DEAL::61.2757 61.2757 +DEAL::65.4369 65.4369 +DEAL::64.6135 64.6135 +DEAL::58.3945 58.3945 +DEAL::46.1455 46.1455 +DEAL::27.0190 27.0190 +DEAL::45.0353 45.0353 +DEAL::80.9161 80.9161 +DEAL::107.327 107.327 +DEAL::124.314 124.314 +DEAL::131.858 131.858 +DEAL::129.623 129.623 +DEAL::116.819 116.819 +DEAL::92.1685 92.1685 +DEAL::53.9305 53.9305 +DEAL::69.2645 69.2645 +DEAL::122.495 122.495 +DEAL::160.725 160.725 +DEAL::184.794 184.794 +DEAL::195.060 195.060 +DEAL::191.200 191.200 +DEAL::172.090 172.090 +DEAL::135.779 135.779 +DEAL::79.5345 79.5345 +DEAL::91.5276 91.5276 +DEAL::160.074 160.074 +DEAL::208.285 208.285 +DEAL::238.076 238.076 +DEAL::250.389 250.389 +DEAL::245.025 245.025 +DEAL::220.563 220.563 +DEAL::174.324 174.324 +DEAL::102.428 102.428 +DEAL::109.771 109.771 +DEAL::189.991 189.991 +DEAL::245.262 245.262 +DEAL::278.838 278.838 +DEAL::292.394 292.394 +DEAL::285.950 285.950 +DEAL::257.813 257.813 +DEAL::204.524 204.524 +DEAL::120.855 120.855 +DEAL::121.567 121.567 +DEAL::207.855 207.855 +DEAL::265.937 265.937 +DEAL::300.618 300.618 +DEAL::314.399 314.399 +DEAL::307.567 307.567 +DEAL::278.215 278.215 +DEAL::222.104 222.104 +DEAL::132.468 132.468 +DEAL::123.643 123.643 +DEAL::207.924 207.924 +DEAL::263.011 263.011 +DEAL::295.300 295.300 +DEAL::308.015 308.015 +DEAL::301.706 301.706 +DEAL::274.376 274.376 +DEAL::221.208 221.208 +DEAL::133.913 133.913 +DEAL::111.079 111.079 +DEAL::182.188 182.188 +DEAL::226.884 226.884 +DEAL::252.556 252.556 +DEAL::262.656 262.656 +DEAL::257.865 257.865 +DEAL::236.376 236.376 +DEAL::193.440 193.440 +DEAL::119.974 119.974 +DEAL::75.4858 75.4858 +DEAL::118.864 118.864 +DEAL::144.783 144.783 +DEAL::159.382 159.382 +DEAL::165.189 165.189 +DEAL::162.720 162.720 +DEAL::150.825 150.825 +DEAL::126.202 126.202 +DEAL::81.5441 81.5441 +DEAL::OK -- 2.39.5