Adding cusparse ILU and incomplete Cholesky solvers on GPU.

examples/ex1.cpp

@@ -175,7 +175,8 @@ int main(int argc, char *argv[])
    //    domain integrator.
    BilinearForm a(&fespace);
    if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
-   a.AddDomainIntegrator(new DiffusionIntegrator(one));
+   //a.AddDomainIntegrator(new DiffusionIntegrator(one));
+   a.AddDomainIntegrator(new MassIntegrator(one));
 
    // 10. Assemble the bilinear form and the corresponding linear system,
    //     applying any necessary transformations such as: eliminating boundary
@@ -184,19 +185,40 @@ int main(int argc, char *argv[])
    if (static_cond) { a.EnableStaticCondensation(); }
    a.Assemble();
 
-   OperatorPtr A;
+   OperatorPtr A, As;
    Vector B, X;
-   a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+   Array<int> empty_list;
+   a.FormSystemMatrix(empty_list, As);
+   //a.FormLinearSystem(empty_list, x, b, A, X, B);
+   //a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
 
-   cout << "Size of linear system: " << A->Height() << endl;
+   //cout << "Size of linear system: " << A->Height() << endl;
 
    // 11. Solve the linear system A X = B.
    if (!pa)
    {
 #ifndef MFEM_USE_SUITESPARSE
       // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
-      GSSmoother M((SparseMatrix&)(*A));
-      PCG(*A, M, B, X, 1, 200, 1e-12, 0.0);
+      //GSSmoother M((SparseMatrix&)(*A));
+
+      //SparseMatrix &Asp = *As.As<SparseMatrix>();
+      SparseMatrix &Asp = a.SpMat();
+
+      Asp.Finalize();
+      Asp.SortColumnIndices();
+
+      Vector tmpx(B.Size());
+      Vector tmpy(B.Size());
+      tmpx = 1.0;
+      tmpy = 0.0;
+
+      //As.As<SparseMatrix>()->Mult(tmpx, tmpy);
+      Asp.Mult(tmpx, tmpy);
+
+      //IncompleteCholesky M(*As.As<SparseMatrix>());
+      IncompleteCholesky M(Asp);
+      //ILUcusparse M(*A.As<SparseMatrix>());
+      PCG(*As, M, B, X, 1, 200, 1e-12, 0.0);
 #else
       // If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
       UMFPackSolver umf_solver;

examples/ex1p.cpp

@@ -122,7 +122,7 @@ int main(int argc, char *argv[])
    {
       int ref_levels =
          (int)floor(log(10000./mesh.GetNE())/log(2.)/dim);
-      for (int l = 0; l < ref_levels; l++)
+      for (int l = 0; l < ref_levels-1; l++)
       {
          mesh.UniformRefinement();
       }
@@ -134,7 +134,7 @@ int main(int argc, char *argv[])
    ParMesh pmesh(MPI_COMM_WORLD, mesh);
    mesh.Clear();
    {
-      int par_ref_levels = 2;
+      int par_ref_levels = 1;
       for (int l = 0; l < par_ref_levels; l++)
       {
          pmesh.UniformRefinement();
@@ -216,6 +216,13 @@ int main(int argc, char *argv[])
    Vector B, X;
    a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
 
+   SparseMatrix Asp;
+   A.As<HypreParMatrix>()->GetDiag(Asp);
+   Vector diag;
+
+   StopWatch sw;
+   sw.Start();
+
    // 13. Solve the linear system A X = B.
    //     * With full assembly, use the BoomerAMG preconditioner from hypre.
    //     * With partial assembly, use Jacobi smoothing, for now.
@@ -229,7 +236,14 @@ int main(int argc, char *argv[])
    }
    else
    {
-      prec = new HypreBoomerAMG;
+      //prec = new HypreBoomerAMG;
+      Asp.Finalize();
+      Asp.SortColumnIndices();
+
+      Asp.GetDiag(diag);
+      prec = new OperatorJacobiSmoother(diag, ess_tdof_list);
+      //prec = new IncompleteCholesky(Asp);
+      //prec = new ILUcusparse(Asp);
    }
    CGSolver cg(MPI_COMM_WORLD);
    cg.SetRelTol(1e-12);
@@ -240,6 +254,9 @@ int main(int argc, char *argv[])
    cg.Mult(B, X);
    delete prec;
 
+   sw.Stop();
+   cout << "Step 13 solve time " << sw.RealTime() << endl;
+
    // 14. Recover the parallel grid function corresponding to X. This is the
    //     local finite element solution on each processor.
    a.RecoverFEMSolution(X, b, x);

linalg/solvers.cpp

@@ -2972,4 +2972,36 @@ KLUSolver::~KLUSolver()
 
 #endif // MFEM_USE_SUITESPARSE
 
+IncompleteCholesky::IncompleteCholesky(SparseMatrix &A_) : A(&A_)
+{
+#ifdef MFEM_USE_CUDA
+   A->IncompleteCholeskySetup();
+#endif
+}
+
+void IncompleteCholesky::Mult(const Vector &b, Vector &x) const
+{
+#ifdef MFEM_USE_CUDA
+   A->IncompleteCholeskyMult(b, x);
+#else
+   x = b;
+#endif
+}
+
+ILUcusparse::ILUcusparse(SparseMatrix &A_) : A(&A_)
+{
+#ifdef MFEM_USE_CUDA
+   A->ILUSetup();
+#endif
+}
+
+void ILUcusparse::Mult(const Vector &b, Vector &x) const
+{
+#ifdef MFEM_USE_CUDA
+   A->ILUMult(b, x);
+#else
+   x = b;
+#endif
+}
+
 }

linalg/solvers.hpp

@@ -799,6 +799,28 @@ public:
 
 #endif // MFEM_USE_SUITESPARSE
 
+class IncompleteCholesky : public IterativeSolver
+{
+private:
+   SparseMatrix *A;
+
+public:
+   IncompleteCholesky(SparseMatrix &A_);
+
+   virtual void Mult(const Vector &b, Vector &x) const;
+};
+
+class ILUcusparse : public IterativeSolver
+{
+private:
+   SparseMatrix *A;
+
+public:
+   ILUcusparse(SparseMatrix &A_);
+
+   virtual void Mult(const Vector &b, Vector &x) const;
+};
+
 }
 
 #endif // MFEM_SOLVERS

linalg/sparsemat.cpp

@@ -3899,4 +3899,323 @@ void SparseMatrix::Swap(SparseMatrix &other)
    mfem::Swap(isSorted, other.isSorted);
 }
 
+#ifdef MFEM_USE_CUDA
+void SparseMatrix::IncompleteCholeskyMult(const Vector &x, Vector &y) const
+{
+   if (!(Device::Allows(Backend::CUDA_MASK) && useCuSparse))
+   {
+      y = x;
+      return;
+   }
+
+   MFEM_VERIFY(initCholesky, "Setup not done");
+
+   const double alpha = 1.0;
+
+   auto d_x = x.Read();
+   auto d_y = y.ReadWrite();
+   auto d_z = vecZ.ReadWrite();
+
+   const int height = this->height;
+   const int nnz = J.Capacity();
+   int64_t m = height;
+
+   auto d_csrRowPtr = Read(I, height+1);
+   auto d_csrColInd = Read(J, nnz);
+   auto d_csrVal = Read(A, nnz);
+
+   cusparseDnVecSetValues(vecX_descr, const_cast<double *>(d_x));
+   cusparseDnVecSetValues(vecY_descr, d_y);
+   cusparseDnVecSetValues(vecZ_descr, d_z);
+
+   const cusparseOperation_t trans_L  = CUSPARSE_OPERATION_NON_TRANSPOSE;
+   const cusparseOperation_t trans_Lt = CUSPARSE_OPERATION_TRANSPOSE;
+
+   const cusparseSolvePolicy_t policy_L  = CUSPARSE_SOLVE_POLICY_NO_LEVEL;
+   const cusparseSolvePolicy_t policy_Lt = CUSPARSE_SOLVE_POLICY_USE_LEVEL;
+
+   // Solve L*z = x
+   cusparseDcsrsv2_solve(handle, trans_L, m, nnz, &alpha, descr_L,
+                         d_csrVal, d_csrRowPtr, d_csrColInd, info_L,
+                         d_x, d_z, policy_L, pBuffer);
+
+   // Solve L'*y = z
+   cusparseDcsrsv2_solve(handle, trans_Lt, m, nnz, &alpha, descr_L,
+                         d_csrVal, d_csrRowPtr, d_csrColInd, info_Lt,
+                         d_z, d_y, policy_Lt, pBuffer);
+}
+
+void SparseMatrix::IncompleteCholeskySetup()
+{
+   if (!(Device::Allows(Backend::CUDA_MASK) && useCuSparse))
+   {
+      return;
+   }
+
+   MFEM_VERIFY(!initILU && !initCholesky, "");
+
+   const int height = this->height;
+   const int nnz = J.Capacity();
+
+   auto d_csrRowPtr = Read(I, height+1);
+   auto d_csrColInd = Read(J, nnz);
+   auto d_csrVal = Read(A, nnz);
+
+   //MFEM_VERIFY(I[0] == 0, "cusparse thinks this is not zero based");
+
+   csric02Info_t info_M  = 0;
+   int bufferSize_M;
+   int pBufferSize_L;
+   int pBufferSize_Lt;
+   int pBufferSize;
+   int structural_zero;
+   int numerical_zero;
+
+   const cusparseSolvePolicy_t policy_M  = CUSPARSE_SOLVE_POLICY_NO_LEVEL;
+   const cusparseSolvePolicy_t policy_L  = CUSPARSE_SOLVE_POLICY_NO_LEVEL;
+   const cusparseSolvePolicy_t policy_Lt = CUSPARSE_SOLVE_POLICY_USE_LEVEL;
+   const cusparseOperation_t trans_L  = CUSPARSE_OPERATION_NON_TRANSPOSE;
+   const cusparseOperation_t trans_Lt = CUSPARSE_OPERATION_TRANSPOSE;
+
+   // step 1: create a descriptor which contains
+   // - matrix M is base-0
+   // - matrix L is base-0
+   // - matrix L is lower triangular
+   // - matrix L has non-unit diagonal
+   cusparseCreateMatDescr(&descr_M);
+   cusparseSetMatIndexBase(descr_M, CUSPARSE_INDEX_BASE_ZERO);
+   cusparseSetMatType(descr_M, CUSPARSE_MATRIX_TYPE_GENERAL);
+
+   cusparseCreateMatDescr(&descr_L);
+   cusparseSetMatIndexBase(descr_L, CUSPARSE_INDEX_BASE_ZERO);
+   cusparseSetMatType(descr_L, CUSPARSE_MATRIX_TYPE_GENERAL);
+   cusparseSetMatFillMode(descr_L, CUSPARSE_FILL_MODE_LOWER);
+   cusparseSetMatDiagType(descr_L, CUSPARSE_DIAG_TYPE_NON_UNIT);
+
+   // step 2: create a empty info structure
+   // we need one info for csric02 and two info's for csrsv2
+   cusparseCreateCsric02Info(&info_M);
+   cusparseCreateCsrsv2Info(&info_L);
+   cusparseCreateCsrsv2Info(&info_Lt);
+
+   int64_t m = height;
+
+   // step 3: query how much memory used in csric02 and csrsv2, and allocate the buffer
+   cusparseDcsric02_bufferSize(handle, m, nnz,
+                               descr_M, const_cast<double *>(d_csrVal), const_cast<int *>(d_csrRowPtr),
+                               const_cast<int *>(d_csrColInd), info_M, &bufferSize_M);
+   cusparseDcsrsv2_bufferSize(handle, trans_L, m, nnz,
+                              descr_L, const_cast<double *>(d_csrVal), const_cast<int *>(d_csrRowPtr),
+                              const_cast<int *>(d_csrColInd), info_L, &pBufferSize_L);
+   cusparseDcsrsv2_bufferSize(handle, trans_Lt, m, nnz,
+                              descr_L, const_cast<double *>(d_csrVal), const_cast<int *>(d_csrRowPtr),
+                              const_cast<int *>(d_csrColInd), info_Lt,&pBufferSize_Lt);
+
+   pBufferSize = max(bufferSize_M, max(pBufferSize_L, pBufferSize_Lt));
+
+   // pBuffer returned by cudaMalloc is automatically aligned to 128 bytes.
+   cudaMalloc((void**)&pBuffer, pBufferSize);
+
+   // step 4: perform analysis of incomplete Cholesky on M
+   //         perform analysis of triangular solve on L
+   //         perform analysis of triangular solve on L'
+   // The lower triangular part of M has the same sparsity pattern as L, so
+   // we can do analysis of csric02 and csrsv2 simultaneously.
+
+   cusparseDcsric02_analysis(handle, m, nnz, descr_M,
+                             d_csrVal, d_csrRowPtr, d_csrColInd, info_M,
+                             policy_M, pBuffer);
+   status = cusparseXcsric02_zeroPivot(handle, info_M, &structural_zero);
+   if (CUSPARSE_STATUS_ZERO_PIVOT == status)
+   {
+      printf("A(%d,%d) is missing\n", structural_zero, structural_zero);
+   }
+
+   cusparseDcsrsv2_analysis(handle, trans_L, m, nnz, descr_L,
+                            d_csrVal, d_csrRowPtr, d_csrColInd,
+                            info_L, policy_L, pBuffer);
+
+   cusparseDcsrsv2_analysis(handle, trans_Lt, m, nnz, descr_L,
+                            d_csrVal, d_csrRowPtr, d_csrColInd,
+                            info_Lt, policy_Lt, pBuffer);
+
+   // step 5: M = L * L'
+   cusparseDcsric02(handle, m, nnz, descr_M,
+                    const_cast<double *>(d_csrVal), const_cast<int *>(d_csrRowPtr),
+                    const_cast<int *>(d_csrColInd), info_M, policy_M, pBuffer);
+   status = cusparseXcsric02_zeroPivot(handle, info_M, &numerical_zero);
+   if (CUSPARSE_STATUS_ZERO_PIVOT == status)
+   {
+      printf("L(%d,%d) is zero\n", numerical_zero, numerical_zero);
+   }
+
+   vecZ.SetSize(height);
+   vecZ = 0.0;
+   auto d_z = vecZ.ReadWrite();
+   cusparseCreateDnVec(&vecZ_descr, vecZ.Size(), d_z, CUDA_R_64F);
+
+   initCholesky = true;
+}
+
+void SparseMatrix::ILUMult(const Vector &x, Vector &y) const
+{
+   MFEM_VERIFY(initILU, "Setup not done");
+
+   const double alpha = 1.0;
+
+   auto d_x = x.Read();
+   auto d_y = y.ReadWrite();
+   auto d_z = vecZ.ReadWrite();
+
+   const int height = this->height;
+   const int nnz = J.Capacity();
+   int64_t m = height;
+
+   auto d_csrRowPtr = Read(I, height+1);
+   auto d_csrColInd = Read(J, nnz);
+   auto d_csrVal = Read(A, nnz);
+
+   cusparseDnVecSetValues(vecX_descr, const_cast<double *>(d_x));
+   cusparseDnVecSetValues(vecY_descr, d_y);
+   cusparseDnVecSetValues(vecZ_descr, d_z);
+
+   const cusparseSolvePolicy_t policy_L = CUSPARSE_SOLVE_POLICY_NO_LEVEL;
+   const cusparseSolvePolicy_t policy_U = CUSPARSE_SOLVE_POLICY_USE_LEVEL;
+   const cusparseOperation_t trans_L  = CUSPARSE_OPERATION_NON_TRANSPOSE;
+   const cusparseOperation_t trans_U  = CUSPARSE_OPERATION_NON_TRANSPOSE;
+
+   // Solve L*z = x
+   cusparseDcsrsv2_solve(handle, trans_L, m, nnz, &alpha, descr_L,
+                         d_csrVal, d_csrRowPtr, d_csrColInd, info_L,
+                         d_x, d_z, policy_L, pBuffer);
+
+   // Solve U*y = z
+   cusparseDcsrsv2_solve(handle, trans_U, m, nnz, &alpha, descr_U,
+                         d_csrVal, d_csrRowPtr, d_csrColInd, info_U,
+                         d_z, d_y, policy_U, pBuffer);
+
+   // TODO: destructor
+}
+
+void SparseMatrix::ILUSetup()
+{
+   if (!(Device::Allows(Backend::CUDA_MASK) && useCuSparse))
+   {
+      return;
+   }
+
+   MFEM_VERIFY(!initILU && !initCholesky, "");
+
+   const int height = this->height;
+   const int nnz = J.Capacity();
+
+   auto d_csrRowPtr = Read(I, height+1);
+   auto d_csrColInd = Read(J, nnz);
+   auto d_csrVal = Read(A, nnz);
+
+   csrilu02Info_t info_M  = 0;
+   int pBufferSize_M;
+   int pBufferSize_L;
+   int pBufferSize_U;
+   int pBufferSize;
+   int structural_zero;
+   int numerical_zero;
+   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;
+   const cusparseOperation_t trans_L  = CUSPARSE_OPERATION_NON_TRANSPOSE;
+   const cusparseOperation_t trans_U  = CUSPARSE_OPERATION_NON_TRANSPOSE;
+
+   // 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
+   cusparseCreateMatDescr(&descr_M);
+   cusparseSetMatIndexBase(descr_M, CUSPARSE_INDEX_BASE_ZERO);
+   cusparseSetMatType(descr_M, CUSPARSE_MATRIX_TYPE_GENERAL);
+
+   cusparseCreateMatDescr(&descr_L);
+   cusparseSetMatIndexBase(descr_L, CUSPARSE_INDEX_BASE_ZERO);
+   cusparseSetMatType(descr_L, CUSPARSE_MATRIX_TYPE_GENERAL);
+   cusparseSetMatFillMode(descr_L, CUSPARSE_FILL_MODE_LOWER);
+   cusparseSetMatDiagType(descr_L, CUSPARSE_DIAG_TYPE_UNIT);
+
+   cusparseCreateMatDescr(&descr_U);
+   cusparseSetMatIndexBase(descr_U, CUSPARSE_INDEX_BASE_ZERO);
+   cusparseSetMatType(descr_U, CUSPARSE_MATRIX_TYPE_GENERAL);
+   cusparseSetMatFillMode(descr_U, CUSPARSE_FILL_MODE_UPPER);
+   cusparseSetMatDiagType(descr_U, CUSPARSE_DIAG_TYPE_NON_UNIT);
+
+   // step 2: create a empty info structure
+   // we need one info for csrilu02 and two info's for csrsv2
+   cusparseCreateCsrilu02Info(&info_M);
+   cusparseCreateCsrsv2Info(&info_L);
+   cusparseCreateCsrsv2Info(&info_U);
+
+   // step 3: query how much memory used in csrilu02 and csrsv2, and allocate the buffer
+   int64_t m = height;
+
+   cusparseDcsrilu02_bufferSize(handle, m, nnz, descr_M,
+                                const_cast<double *>(d_csrVal), const_cast<int *>(d_csrRowPtr),
+                                const_cast<int *>(d_csrColInd), info_M, &pBufferSize_M);
+   cusparseDcsrsv2_bufferSize(handle, trans_L, m, nnz, descr_L,
+                              const_cast<double *>(d_csrVal), const_cast<int *>(d_csrRowPtr),
+                              const_cast<int *>(d_csrColInd), info_L, &pBufferSize_L);
+   cusparseDcsrsv2_bufferSize(handle, trans_U, m, nnz, descr_U,
+                              const_cast<double *>(d_csrVal), const_cast<int *>(d_csrRowPtr),
+                              const_cast<int *>(d_csrColInd), info_U, &pBufferSize_U);
+
+   pBufferSize = max(pBufferSize_M, max(pBufferSize_L, pBufferSize_U));
+
+   // pBuffer returned by cudaMalloc is automatically aligned to 128 bytes.
+   cudaMalloc((void**)&pBuffer, pBufferSize);
+
+   // 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.
+
+   cusparseDcsrilu02_analysis(handle, m, nnz, descr_M,
+                              d_csrVal, d_csrRowPtr, d_csrColInd, info_M,
+                              policy_M, pBuffer);
+   status = cusparseXcsrilu02_zeroPivot(handle, info_M, &structural_zero);
+   if (CUSPARSE_STATUS_ZERO_PIVOT == status)
+   {
+      printf("A(%d,%d) is missing\n", structural_zero, structural_zero);
+   }
+
+   cusparseDcsrsv2_analysis(handle, trans_L, m, nnz, descr_L,
+                            d_csrVal, d_csrRowPtr, d_csrColInd,
+                            info_L, policy_L, pBuffer);
+
+   cusparseDcsrsv2_analysis(handle, trans_U, m, nnz, descr_U,
+                            d_csrVal, d_csrRowPtr, d_csrColInd,
+                            info_U, policy_U, pBuffer);  // bug?
+
+   // step 5: M = L * U
+   cusparseDcsrilu02(handle, m, nnz, descr_M,
+                     const_cast<double *>(d_csrVal), const_cast<int *>(d_csrRowPtr),
+                     const_cast<int *>(d_csrColInd), info_M, policy_M, pBuffer);
+   status = cusparseXcsrilu02_zeroPivot(handle, info_M, &numerical_zero);
+   if (CUSPARSE_STATUS_ZERO_PIVOT == status)
+   {
+      printf("U(%d,%d) is zero\n", numerical_zero, numerical_zero);
+   }
+
+   vecZ.SetSize(height);
+   vecZ = 0.0;
+   auto d_z = vecZ.ReadWrite();
+   cusparseCreateDnVec(&vecZ_descr, vecZ.Size(), d_z, CUDA_R_64F);
+
+   initILU = true;
+
+   // TODO: destructor
+}
+#endif // MFEM_USE_CUDA
+
 }

linalg/sparsemat.hpp

@@ -103,6 +103,21 @@ protected:
    mutable cusparseSpMatDescr_t matA_descr;
    mutable cusparseDnVecDescr_t vecX_descr;
    mutable cusparseDnVecDescr_t vecY_descr;
+   mutable cusparseDnVecDescr_t vecZ_descr;
+   mutable Vector vecZ;
+
+   cusparseMatDescr_t descr_M = 0;
+   cusparseMatDescr_t descr_L = 0;
+   cusparseMatDescr_t descr_U = 0;
+
+   csrsv2Info_t  info_L  = 0;
+   csrsv2Info_t  info_Lt = 0;
+   csrsv2Info_t  info_U  = 0;
+
+   void *pBuffer = 0;
+
+   bool initILU = false;
+   bool initCholesky = false;
 #endif
 
 public:
@@ -610,6 +625,14 @@ public:
 
    void Swap(SparseMatrix &other);
 
+#ifdef MFEM_USE_CUDA
+   void IncompleteCholeskySetup();
+   void IncompleteCholeskyMult(const Vector &x, Vector &y) const;
+
+   void ILUSetup();
+   void ILUMult(const Vector &x, Vector &y) const;
+#endif
+
    /// Destroys sparse matrix.
    virtual ~SparseMatrix()
    {