]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add PreconditionILU
authorDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Sun, 19 Aug 2018 21:47:38 +0000 (23:47 +0200)
committerDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Mon, 20 Aug 2018 18:04:51 +0000 (20:04 +0200)
include/deal.II/lac/cuda_precondition.h
source/lac/cuda_precondition.cu
tests/cuda/precondition_01.cu
tests/cuda/precondition_02.cu
tests/cuda/precondition_02.output

index 076b483ae6df0efc589110420c25ea71e93998a7..73ce2a644060b003d32413a27bb25fd235272373 100644 (file)
@@ -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 <tt>@<float@> and
+   * @<double@></tt>.
+   *
+   * @ingroup Preconditioners CUDAWrappers
+   * @author Daniel Arndt
+   * @date 2018
+   */
+  template <typename Number>
+  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<Number> &) = delete;
+
+    /**
+     * The copy assignment operator is deleted.
+     */
+    PreconditionILU &
+    operator=(const PreconditionILU<Number> &) = 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<Number> &matrix,
+               const AdditionalData &      additional_data = AdditionalData());
+
+    /**
+     * Apply the preconditioner.
+     */
+    void
+    vmult(LinearAlgebra::CUDAWrappers::Vector<Number> &      dst,
+          const LinearAlgebra::CUDAWrappers::Vector<Number> &src) const;
+
+    /**
+     * Apply the preconditioner. Since the preconditioner is symmetric, this
+     * is the same as vmult().
+     */
+    void
+    Tvmult(LinearAlgebra::CUDAWrappers::Vector<Number> &      dst,
+           const LinearAlgebra::CUDAWrappers::Vector<Number> &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<Number[], void (*)(Number *)> 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<Number[], void (*)(Number *)> tmp_dev;
+
+    /**
+     *
+     */
+    std::unique_ptr<void, void (*)(void *)> 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 <typename Number>
+  inline typename PreconditionILU<Number>::size_type
+  PreconditionILU<Number>::m() const
+  {
+    return n_rows;
+  }
+
+
+
+  template <typename Number>
+  inline typename PreconditionILU<Number>::size_type
+  PreconditionILU<Number>::n() const
+  {
+    return n_rows;
+  }
 #  endif // DOXYGEN
 
 } // namespace CUDAWrappers
index 16da6e7782fc5a53042b7b2fae073e68da710998..26ab5d58320ea80bb2ec716bb2c9a38ac2354432 100644 (file)
@@ -20,6 +20,376 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace
 {
+  /**
+   * Template wrapper for cusparse<t>csrilu02.
+   * (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 <typename Number>
+  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<float>(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<double>(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<cuComplex>(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<cuDoubleComplex>(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 cusparse<t>csrilu02_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 <typename Number>
+  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<float>(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<double>(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<cuComplex>(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<cuDoubleComplex>(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 cusparse<t>csrilu02_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 <typename Number>
+  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<float>(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<double>(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<cuComplex>(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<cuDoubleComplex>(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 cusparse<t>csric02
    * (https://docs.nvidia.com/cuda/cusparse/index.html#cusparse-lt-t-gt-csric02).
@@ -1140,11 +1510,334 @@ namespace CUDAWrappers
 
 
 
+  template <typename Number>
+  PreconditionILU<Number>::AdditionalData::AdditionalData(
+    bool use_level_analysis_)
+    : use_level_analysis(use_level_analysis_)
+  {}
+
+
+
+  template <typename Number>
+  PreconditionILU<Number>::PreconditionILU(
+    const Utilities::CUDA::Handle &handle)
+    : cusparse_handle(handle.cusparse_handle)
+    , P_val_dev(nullptr, delete_device_vector<Number>)
+    , P_row_ptr_dev(nullptr)
+    , P_column_index_dev(nullptr)
+    , tmp_dev(nullptr, delete_device_vector<Number>)
+    , buffer_dev(nullptr, delete_device_vector<void>)
+    , 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 <typename Number>
+  PreconditionILU<Number>::~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 <typename Number>
+  void
+  PreconditionILU<Number>::initialize(const SparseMatrix<Number> &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<Number>(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<Number>(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<void *>(
+      allocate_device_vector<char>(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 <typename Number>
+  void
+  PreconditionILU<Number>::vmult(
+    LinearAlgebra::CUDAWrappers::Vector<Number> &      dst,
+    const LinearAlgebra::CUDAWrappers::Vector<Number> &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<unsigned int>(n_rows));
+    AssertDimension(src.size(), static_cast<unsigned int>(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<Number>::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 <typename Number>
+  void
+  PreconditionILU<Number>::Tvmult(
+    LinearAlgebra::CUDAWrappers::Vector<Number> & /*dst*/,
+    const LinearAlgebra::CUDAWrappers::Vector<Number> & /*src*/) const
+  {
+    Assert(false, ExcNotImplemented());
+  }
+
+
+
   // explicit instantiations
   template class PreconditionIC<float>;
   template class PreconditionIC<double>;
   // template class PreconditionIC<cuComplex>;
   // template class PreconditionIC<cuDoubleComplex>;
+  template class PreconditionILU<float>;
+  template class PreconditionILU<double>;
+  // template class PreconditionILU<cuComplex>;
+  // template class PreconditionILU<cuDoubleComplex>;
 } // namespace CUDAWrappers
 
 DEAL_II_NAMESPACE_CLOSE
index 05cbe128d13743a0abef47ee29d3ce32e5b91270..b238ad379a74e5a6ffbf9c737318ea502f5a5fa2 100644 (file)
@@ -14,6 +14,7 @@
 // ---------------------------------------------------------------------
 
 // Check that dealii::SolverCG works with CUDAWrappers::SparseMatrix
+// and PreconditionIC
 
 #include <deal.II/base/cuda.h>
 #include <deal.II/base/exceptions.h>
index c5f68ab9600cfe5437a7571b1268ecd230e21b58..990d4a92b93b747d04c27664a06ad2684f9c4d7b 100644 (file)
 // ---------------------------------------------------------------------
 
 // Check that dealii::SolverCG works with CUDAWrappers::SparseMatrix
+// and PreconditionILU
 
 #include <deal.II/base/cuda.h>
 #include <deal.II/base/exceptions.h>
 
+#include <deal.II/lac/cuda_precondition.h>
 #include <deal.II/lac/cuda_sparse_matrix.h>
-#include <deal.II/lac/precondition.h>
 #include <deal.II/lac/read_write_vector.h>
 #include <deal.II/lac/solver_cg.h>
 #include <deal.II/lac/solver_control.h>
-#include <deal.II/lac/vector.h>
-
-#include <memory>
 
 #include "../testmatrix.h"
 #include "../tests.h"
 
-DEAL_II_NAMESPACE_OPEN
-
-namespace CUDAWrappers
-{
-  /** \addtogroup CUDAWrappers
-   *  @{
-   */
-
-  /**
-   * Template wrapper for cusparse<t>csrilu02.
-   * (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 <typename Number>
-  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<float>(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<double>(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<cuComplex>(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<cuDoubleComplex>(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 cusparse<t>csrilu02_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 <typename Number>
-  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<float>(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<double>(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<cuComplex>(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<cuDoubleComplex>(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 cusparse<t>csrilu02_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 <typename Number>
-  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<float>(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<double>(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<cuComplex>(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<cuDoubleComplex>(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 cusparse<t>csric02
-   * (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 <typename Number>
-  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<float>(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<double>(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<cuComplex>(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<cuDoubleComplex>(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 cusparse<t>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.
-   */
-  template <typename Number>
-  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<float>(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<double>(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<cuComplex>(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<cuDoubleComplex>(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 cusparse<t>csrsv2_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 <typename Number>
-  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<float>(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<double>(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<cuComplex>(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<cuDoubleComplex>(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 cusparse<t>csric02_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 <typename Number>
-  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<float>(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<double>(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<cuComplex>(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<cuDoubleComplex>(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 cusparse<t>csrsv2_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 <typename Number>
-  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<float>(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<double>(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<cuComplex>(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<cuDoubleComplex>(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 cusparse<t>csric02_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 <typename Number>
-  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<float>(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<double>(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<cuComplex>(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<cuDoubleComplex>(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 <typename Number>
-  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 <typename Number>
-  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 <tt>@<float@> and
-     * @<double@></tt>.
-     *
-     * @ingroup Preconditioners CUDAWrappers
-     * @author Daniel Arndt
-     * @date 2018
-     */
-    template <typename Number>
-    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<Number> &) = delete;
-
-      /**
-       * The copy assignment operator is deleted.
-       */
-      PreconditionILU &
-      operator=(const PreconditionILU<Number> &) = 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<Number> &matrix,
-                 const AdditionalData &additional_data = AdditionalData());
-
-      /**
-       * Apply the preconditioner.
-       */
-      void
-      vmult(LinearAlgebra::CUDAWrappers::Vector<Number> &      dst,
-            const LinearAlgebra::CUDAWrappers::Vector<Number> &src) const;
-
-      /**
-       * Apply the preconditioner. Since the preconditioner is symmetric, this
-       * is the same as vmult().
-       */
-      void
-      Tvmult(LinearAlgebra::CUDAWrappers::Vector<Number> &      dst,
-             const LinearAlgebra::CUDAWrappers::Vector<Number> &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<Number[], void (*)(Number *)> 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<Number[], void (*)(Number *)> tmp_dev;
-
-      /**
-       *
-       */
-      std::unique_ptr<void, void (*)(void *)> 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 <typename Number>
-    PreconditionILU<Number>::AdditionalData::AdditionalData(
-      bool use_level_analysis_)
-      : use_level_analysis(use_level_analysis_)
-    {}
-
-
-
-    template <typename Number>
-    PreconditionILU<Number>::PreconditionILU(
-      const Utilities::CUDA::Handle &handle)
-      : cusparse_handle(handle.cusparse_handle)
-      , P_val_dev(nullptr, delete_device_vector<Number>)
-      , P_row_ptr_dev(nullptr)
-      , P_column_index_dev(nullptr)
-      , tmp_dev(nullptr, delete_device_vector<Number>)
-      , buffer_dev(nullptr, delete_device_vector<void>)
-      , 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 <typename Number>
-    PreconditionILU<Number>::~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 <typename Number>
-    void
-    PreconditionILU<Number>::initialize(const SparseMatrix<Number> &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<Number>(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<Number>(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<void *>(
-        allocate_device_vector<char>(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 <typename Number>
-    void
-    PreconditionILU<Number>::vmult(
-      LinearAlgebra::CUDAWrappers::Vector<Number> &      dst,
-      const LinearAlgebra::CUDAWrappers::Vector<Number> &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<unsigned int>(n_rows));
-      AssertDimension(src.size(), static_cast<unsigned int>(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 <typename Number>
-    void
-    PreconditionILU<Number>::Tvmult(
-      LinearAlgebra::CUDAWrappers::Vector<Number> &      dst,
-      const LinearAlgebra::CUDAWrappers::Vector<Number> &src) const
-    {
-      // the constructed preconditioner is symmetric
-      vmult(dst, src);
-    }
-
-    template <typename Number>
-    PreconditionILU<Number>::size_type
-    PreconditionILU<Number>::m() const
-    {
-      return n_rows;
-    }
-
-
-    template <typename Number>
-    PreconditionILU<Number>::size_type
-    PreconditionILU<Number>::n() const
-    {
-      return n_rows;
-    }
-
-
-
-    template <typename Number>
-    void
-    apply_preconditioner(const SparseMatrix<Number> &A,
-                         const cusparseHandle_t      cusparse_handle,
-                         LinearAlgebra::CUDAWrappers::Vector<Number> &      dst,
-                         const LinearAlgebra::CUDAWrappers::Vector<Number> &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<Number[], void (*)(Number *)> tmp_dev(
-        allocate_device_vector<Number>(dst.size()),
-        delete_device_vector<Number>);
-
-      // 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 <typename Number>
 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<double> A;
+  SparseMatrix<Number> 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<double>       sol_host(size);
-  Vector<double>       rhs_host(size);
-  for (unsigned int i = 0; i < size; ++i)
-    rhs_host[i] = static_cast<double>(i);
-  cg_host.solve(A, sol_host, rhs_host, prec_no);
-
   // Solve on the device
-  CUDAWrappers::SparseMatrix<double>          A_dev(cuda_handle, A);
-  LinearAlgebra::CUDAWrappers::Vector<double> sol_dev(size);
-  LinearAlgebra::CUDAWrappers::Vector<double> rhs_dev(size);
-  LinearAlgebra::ReadWriteVector<double>      rw_vector(size);
+  CUDAWrappers::SparseMatrix<Number>          A_dev(cuda_handle, A);
+  LinearAlgebra::CUDAWrappers::Vector<Number> sol_dev(size);
+  LinearAlgebra::CUDAWrappers::Vector<Number> rhs_dev(size);
+  LinearAlgebra::ReadWriteVector<Number>      rw_vector(size);
   for (unsigned int i = 0; i < size; ++i)
-    rw_vector[i] = static_cast<double>(i);
+    rw_vector[i] = static_cast<Number>(i);
   rhs_dev.import(rw_vector, VectorOperation::insert);
-  SolverCG<LinearAlgebra::CUDAWrappers::Vector<double>> cg_dev(control);
+  SolverControl                                         control(100, 1.e-10);
+  SolverCG<LinearAlgebra::CUDAWrappers::Vector<Number>> cg_dev(control);
 
-  A_dev.print(std::cout);
-  A_dev.print_formatted(std::cout);
-  CUDAWrappers::PreconditionILU<double>    prec_double(cuda_handle);
-  CUDAWrappers::PreconditionILU<float>     prec_float(cuda_handle);
-  CUDAWrappers::PreconditionILU<cuComplex> prec_complex_float(cuda_handle);
-  CUDAWrappers::PreconditionILU<cuDoubleComplex> prec_complex_double(
-    cuda_handle);
+  CUDAWrappers::PreconditionILU<Number> 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<float>(cuda_handle);
+  test<double>(cuda_handle);
 
   deallog << "OK" << std::endl;
 
index d6d49d1ae1a4e1e8ffbdc655e1ffcb32aaa9b5c4..913218546b8cade10c97428372926abd0ee38389 100644 (file)
 
 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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.