From: Daniel Arndt Date: Sun, 19 Aug 2018 21:47:38 +0000 (+0200) Subject: Add PreconditionILU X-Git-Tag: v9.1.0-rc1~769^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8d9e9dba39c68f84957860fdb8a2b57f54273ed4;p=dealii.git Add PreconditionILU --- diff --git a/include/deal.II/lac/cuda_precondition.h b/include/deal.II/lac/cuda_precondition.h index 076b483ae6..73ce2a6440 100644 --- a/include/deal.II/lac/cuda_precondition.h +++ b/include/deal.II/lac/cuda_precondition.h @@ -251,6 +251,220 @@ namespace CUDAWrappers int n_nonzero_elements; }; + /** + * This class implements an incomplete LU factorization 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. + * to the documentation 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=LU$. + */ + 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 $U$. This value can be modified through an + * AdditionalData object. + */ + cusparseSolvePolicy_t policy_U; + + /** + * Determine if level information should be generated for $M=LU$. 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; + }; + /*--------------------------- inline functions ----------------------------*/ # ifndef DOXYGEN @@ -269,6 +483,24 @@ namespace CUDAWrappers { return n_rows; } + + + + template + inline typename PreconditionILU::size_type + PreconditionILU::m() const + { + return n_rows; + } + + + + template + inline typename PreconditionILU::size_type + PreconditionILU::n() const + { + return n_rows; + } # endif // DOXYGEN } // namespace CUDAWrappers diff --git a/source/lac/cuda_precondition.cu b/source/lac/cuda_precondition.cu index 16da6e7782..26ab5d5832 100644 --- a/source/lac/cuda_precondition.cu +++ b/source/lac/cuda_precondition.cu @@ -20,6 +20,376 @@ DEAL_II_NAMESPACE_OPEN namespace { + /** + * 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). @@ -1140,11 +1510,334 @@ namespace CUDAWrappers + 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 = 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 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 + { + Assert(false, ExcNotImplemented()); + } + + + // explicit instantiations template class PreconditionIC; template class PreconditionIC; // template class PreconditionIC; // template class PreconditionIC; + template class PreconditionILU; + template class PreconditionILU; + // template class PreconditionILU; + // template class PreconditionILU; } // namespace CUDAWrappers DEAL_II_NAMESPACE_CLOSE diff --git a/tests/cuda/precondition_01.cu b/tests/cuda/precondition_01.cu index 05cbe128d1..b238ad379a 100644 --- a/tests/cuda/precondition_01.cu +++ b/tests/cuda/precondition_01.cu @@ -14,6 +14,7 @@ // --------------------------------------------------------------------- // Check that dealii::SolverCG works with CUDAWrappers::SparseMatrix +// and PreconditionIC #include #include diff --git a/tests/cuda/precondition_02.cu b/tests/cuda/precondition_02.cu index c5f68ab960..990d4a92b9 100644 --- a/tests/cuda/precondition_02.cu +++ b/tests/cuda/precondition_02.cu @@ -14,2017 +14,21 @@ // --------------------------------------------------------------------- // Check that dealii::SolverCG works with CUDAWrappers::SparseMatrix +// and PreconditionILU #include #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 - +template void test(Utilities::CUDA::Handle &cuda_handle) { @@ -2033,53 +37,33 @@ test(Utilities::CUDA::Handle &cuda_handle) unsigned int size = (problem_size - 1) * (problem_size - 1); FDMatrix testproblem(problem_size, problem_size); SparsityPattern structure(size, size, 5); - SparseMatrix A; + 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); + 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); + rw_vector[i] = static_cast(i); rhs_dev.import(rw_vector, VectorOperation::insert); - SolverCG> cg_dev(control); + SolverControl control(100, 1.e-10); + 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); + CUDAWrappers::PreconditionILU prec_ilu(cuda_handle); + prec_ilu.initialize(A_dev); - // 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); + cg_dev.solve(A_dev, sol_dev, rhs_dev, prec_ilu); // 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; + deallog << rw_vector[i] << std::endl; } int @@ -2089,7 +73,8 @@ main() deallog.depth_console(0); Utilities::CUDA::Handle cuda_handle; - test(cuda_handle); + test(cuda_handle); + test(cuda_handle); deallog << "OK" << std::endl; diff --git a/tests/cuda/precondition_02.output b/tests/cuda/precondition_02.output index d6d49d1ae1..913218546b 100644 --- a/tests/cuda/precondition_02.output +++ b/tests/cuda/precondition_02.output @@ -1,87 +1,168 @@ DEAL:cg::Starting value 416.989 -DEAL:cg::Convergence step 31 value 8.71925e-12 +DEAL:cg::Convergence step 19 value 1.40142e-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 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::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