Squashing commits on AIR_relax

examples/ex9p.cpp

@@ -58,6 +58,81 @@ double inflow_function(const Vector &x);
 // Mesh bounding box
 Vector bb_min, bb_max;
 
+struct AIR_parameters
+{
+   int blocksize;
+   int distanceR;
+   std::string prerelax;
+   std::string postrelax;
+   int interp_type;
+   int relax_type;
+   int coarsen_type;
+   double strength_tolC;
+   double strength_tolR;
+   double filter_tolR;
+   double filterA_tol;
+};
+
+class AIR_prec : public Solver
+{
+private:
+   const HypreParMatrix *A;
+   HypreParMatrix A_s;
+
+   // Preconditioner/solvers for A
+   HypreBoomerAMG *AIR_solver;
+   const AIR_parameters &AIR;
+   int blocksize;
+
+public:
+
+   AIR_prec(const AIR_parameters &_AIR) :
+      AIR_solver(NULL), AIR(_AIR)
+   {
+      blocksize = AIR.blocksize;
+   }
+
+   void SetOperator(const Operator &op)
+   {
+      A = dynamic_cast<const HypreParMatrix *>(&op);
+      delete AIR_solver;
+
+      // Scale A by block-diagonal inverse
+#if MFEM_HYPRE_VERSION >= 21800
+      BlockInverseScale(A, &A_s, NULL, NULL, blocksize, 0);
+      AIR_solver = new HypreBoomerAMG(A_s);
+      AIR_solver->SetLAIROptions(AIR.distanceR, AIR.prerelax,
+                                 AIR.postrelax, AIR.strength_tolC,
+                                 AIR.strength_tolR, AIR.filter_tolR,
+                                 AIR.interp_type, AIR.relax_type,
+                                 AIR.filterA_tol, AIR.coarsen_type,
+                                 -1, 1);
+#else
+      MFEM_ABORT("Must have MFEM_HYPRE_VERSION >= 21800 to use AIR.\n");
+#endif
+      AIR_solver->SetPrintLevel(0);
+      AIR_solver->SetMaxLevels(50);
+   }
+
+   virtual void Mult(const Vector &x, Vector &y) const
+   {
+      // scale the rhs by block inverse and solve system
+      HypreParVector z_s;
+#if MFEM_HYPRE_VERSION >= 21800
+      BlockInverseScale(A, NULL, &x, &z_s, blocksize, 2);
+#endif
+      AIR_solver->Mult(z_s, y);
+   }
+
+   ~AIR_prec()
+   {
+#if MFEM_HYPRE_VERSION >= 21800
+      BlockInverseScale(NULL, NULL, NULL, NULL, 0, -1);
+#endif
+      delete AIR_solver;
+   }
+};
+
 class DG_Solver : public Solver
 {
 private:
@@ -65,7 +140,7 @@ private:
    SparseMatrix M_diag;
    HypreParMatrix *A;
    GMRESSolver linear_solver;
-   BlockILU prec;
+   Solver *prec;
    double dt;
 public:
    DG_Solver(HypreParMatrix &M_, HypreParMatrix &K_, const FiniteElementSpace &fes)
@@ -73,20 +148,40 @@ public:
         K(K_),
         A(NULL),
         linear_solver(M.GetComm()),
-        prec(fes.GetFE(0)->GetDof(),
-             BlockILU::Reordering::MINIMUM_DISCARDED_FILL),
         dt(-1.0)
    {
+      prec = new BlockILU(fes.GetFE(0)->GetDof(),
+                          BlockILU::Reordering::MINIMUM_DISCARDED_FILL);
       linear_solver.iterative_mode = false;
       linear_solver.SetRelTol(1e-9);
       linear_solver.SetAbsTol(0.0);
       linear_solver.SetMaxIter(100);
       linear_solver.SetPrintLevel(0);
-      linear_solver.SetPreconditioner(prec);
+      linear_solver.SetPreconditioner(*prec);
 
       M.GetDiag(M_diag);
    }
 
+   DG_Solver(HypreParMatrix &M_, HypreParMatrix &K_, const FiniteElementSpace &fes,
+             const AIR_parameters &_AIR)
+      : M(M_),
+        K(K_),
+        A(NULL),
+        linear_solver(M.GetComm()),
+        dt(-1.0)
+   {
+      prec = new AIR_prec(_AIR);
+      linear_solver.iterative_mode = false;
+      linear_solver.SetRelTol(1e-9);
+      linear_solver.SetAbsTol(0.0);
+      linear_solver.SetMaxIter(100);
+      linear_solver.SetPrintLevel(0);
+      linear_solver.SetPreconditioner(*prec);
+
+      M.GetDiag(M_diag);
+   }
+
+
    void SetTimeStep(double dt_)
    {
       if (dt_ != dt)
@@ -115,10 +210,12 @@ public:
 
    ~DG_Solver()
    {
+      delete prec;
       delete A;
    }
 };
 
+
 /** A time-dependent operator for the right-hand side of the ODE. The DG weak
     form of du/dt = -v.grad(u) is M du/dt = K u + b, where M and K are the mass
     and advection matrices, and b describes the flow on the boundary. This can
@@ -136,6 +233,8 @@ private:
    mutable Vector z;
 
 public:
+   FE_Evolution(ParBilinearForm &_M, ParBilinearForm &_K, const Vector &_b,
+                const AIR_parameters &_AIR);
    FE_Evolution(ParBilinearForm &_M, ParBilinearForm &_K, const Vector &_b);
 
    virtual void Mult(const Vector &x, Vector &y) const;
@@ -169,7 +268,10 @@ int main(int argc, char *argv[])
    bool paraview = false;
    bool binary = false;
    int vis_steps = 5;
-
+   int solver_type = 1;
+   AIR_parameters AIR0 = {-1, 1, "", "FA", 100, 10, 10,
+                          0.1, 0.01, 0.0, 1e-4
+                         };
    int precision = 8;
    cout.precision(precision);
 
@@ -199,6 +301,8 @@ int main(int argc, char *argv[])
                   "Final time; start time is 0.");
    args.AddOption(&dt, "-dt", "--time-step",
                   "Time step.");
+   args.AddOption(&solver_type, "-st", "--solver-type",
+                  "Solver for implicit solves. 0 for ILU, 1 for pAIR-AMG.");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
@@ -300,6 +404,9 @@ int main(int argc, char *argv[])
       cout << "Number of unknowns: " << global_vSize << endl;
    }
 
+   // Get DG blocksize for AIR
+   AIR0.blocksize = fes->GetFE(0)->GetDof();
+
    // 8. Set up and assemble the parallel bilinear and linear forms (and the
    //    parallel hypre matrices) corresponding to the DG discretization. The
    //    DGTraceIntegrator involves integrals over mesh interior faces.
@@ -427,11 +534,19 @@ int main(int argc, char *argv[])
    // 10. Define the time-dependent evolution operator describing the ODE
    //     right-hand side, and perform time-integration (looping over the time
    //     iterations, ti, with a time-step dt).
-   FE_Evolution adv(*m, *k, *B);
+   FE_Evolution *adv;
+   if (solver_type == 1)
+   {
+      adv = new FE_Evolution(*m, *k, *B, AIR0);
+   }
+   else
+   {
+      adv = new FE_Evolution(*m, *k, *B);
+   }
 
    double t = 0.0;
-   adv.SetTime(t);
-   ode_solver->Init(adv);
+   adv->SetTime(t);
+   ode_solver->Init(*adv);
 
    bool done = false;
    for (int ti = 0; !done; )
@@ -498,6 +613,7 @@ int main(int argc, char *argv[])
    delete ode_solver;
    delete pd;
    delete dc;
+   delete adv;
 
    MPI_Finalize();
    return 0;
@@ -507,8 +623,7 @@ int main(int argc, char *argv[])
 // Implementation of class FE_Evolution
 FE_Evolution::FE_Evolution(ParBilinearForm &_M, ParBilinearForm &_K,
                            const Vector &_b)
-   : TimeDependentOperator(_M.Height()),
-     b(_b),
+   : TimeDependentOperator(_M.Height()), b(_b),
      M_solver(_M.ParFESpace()->GetComm()),
      z(_M.Height())
 {
@@ -551,6 +666,45 @@ FE_Evolution::FE_Evolution(ParBilinearForm &_M, ParBilinearForm &_K,
    M_solver.SetPrintLevel(0);
 }
 
+// Implementation of class FE_Evolution
+FE_Evolution::FE_Evolution(ParBilinearForm &_M, ParBilinearForm &_K,
+                           const Vector &_b, const AIR_parameters &_AIR)
+   : TimeDependentOperator(_M.Height()), b(_b),
+     M_solver(_M.ParFESpace()->GetComm()),
+     z(_M.Height())
+{
+   bool pa = _M.GetAssemblyLevel()==AssemblyLevel::PARTIAL;
+
+   if (pa)
+   {
+      MFEM_ABORT("AIR solver not available for partial assembly.\n");
+   }
+   else
+   {
+      M.Reset(_M.ParallelAssemble(), true);
+      K.Reset(_K.ParallelAssemble(), true);
+   }
+
+   HypreParMatrix &M_mat = *M.As<HypreParMatrix>();
+   HypreParMatrix &K_mat = *K.As<HypreParMatrix>();
+   HypreSmoother *hypre_prec = new HypreSmoother(M_mat, HypreSmoother::GS);
+   M_prec = hypre_prec;
+
+   dg_solver = new DG_Solver(M_mat, K_mat, *_M.FESpace(), _AIR);
+
+   M_solver.SetPreconditioner(*M_prec);
+   M_solver.SetOperator(*M);
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+// Solve the equation:
+//    u_t = M^{-1}(Ku + b),
+// by solving associated linear system
+//    (M - dt*K) d = K*u + b
 void FE_Evolution::ImplicitSolve(const double dt, const Vector &x, Vector &k)
 {
    K->Mult(x, z);

linalg/hypre.cpp

@@ -173,6 +173,13 @@ HypreParVector::HypreParVector(ParFiniteElementSpace *pfes)
    own_ParVector = 1;
 }
 
+void HypreParVector::WrapHypreParVector(hypre_ParVector *y)
+{
+   x = y;
+   _SetDataAndSize_();
+   own_ParVector = 0;
+}
+
 Vector * HypreParVector::GlobalVector() const
 {
    hypre_Vector *hv = hypre_ParVectorToVectorAll(*this);
@@ -978,6 +985,11 @@ void HypreParMatrix::GetOffd(SparseMatrix &offd, HYPRE_Int* &cmap) const
    cmap = A->col_map_offd;
 }
 
+void HypreParMatrix::GetProcRows(SparseMatrix &colCSRMat)
+{
+   MakeWrapper(hypre_MergeDiagAndOffd(A), colCSRMat);
+}
+
 void HypreParMatrix::GetBlocks(Array2D<HypreParMatrix*> &blocks,
                                bool interleaved_rows,
                                bool interleaved_cols) const
@@ -1018,6 +1030,47 @@ HypreParMatrix * HypreParMatrix::Transpose() const
    return new HypreParMatrix(At);
 }
 
+#if MFEM_HYPRE_VERSION >= 21800
+HypreParMatrix * HypreParMatrix::ExtractSubmatrix(Array<int> &indices,
+                                                  double threshhold) const
+{
+   if (!(A->comm))
+   {
+      BuildComm();
+   }
+
+   hypre_ParCSRMatrix *submat;
+
+   // Get number of rows stored on this processor
+   int local_num_vars = hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(A));
+
+   // Form hypre CF-splitting array designating submatrix as F-points (-1)
+   int *CF_marker = new int[local_num_vars];
+   std::fill_n(CF_marker, local_num_vars, 1);
+   for (int j=0; j<indices.Size(); j++)
+   {
+      if (indices[j] > local_num_vars)
+      {
+         MFEM_WARNING("WARNING : " << indices[j] << " > " << local_num_vars);
+      }
+      CF_marker[indices[j]] = -1;
+   }
+
+   // Construct cpts_global array on hypre matrix structure
+   int *cpts_global;
+   hypre_BoomerAMGCoarseParms(MPI_COMM_WORLD, local_num_vars, 1, NULL,
+                              CF_marker, NULL, &cpts_global);
+
+   // Extract submatrix into *submat
+   hypre_ParCSRMatrixExtractSubmatrixFC(A, CF_marker, cpts_global,
+                                        "FF", &submat, threshhold);
+
+   delete[] CF_marker;
+   free(cpts_global);
+   return new HypreParMatrix(submat);
+}
+#endif
+
 HYPRE_Int HypreParMatrix::Mult(HypreParVector &x, HypreParVector &y,
                                double a, double b)
 {
@@ -1587,6 +1640,41 @@ void HypreParMatrix::Destroy()
    }
 }
 
+/* job =  0, extract block diagonal of A and scale A into C
+ * job =  1, job 0 + scale b into d
+ * job =  2, use A to scale b only
+ */
+#if MFEM_HYPRE_VERSION >= 21800
+int BlockInverseScale(const HypreParMatrix *A, HypreParMatrix *C,
+                      const Vector *b, HypreParVector *d, int block, int job)
+{
+   if (0 == job || 1 == job)
+   {
+      hypre_ParCSRMatrix *C_hypre;
+      hypre_ParcsrBdiagInvScal(*A, block, &C_hypre);
+      hypre_ParCSRMatrixDropSmallEntries(C_hypre, 1e-15, 1);
+      (*C).WrapHypreParCSRMatrix(C_hypre);
+   }
+
+   if (1 == job || 2 == job)
+   {
+      HypreParVector *b_Hypre = new HypreParVector(A->GetComm(),
+                                                   A->GetGlobalNumRows(),
+                                                   b->GetData(), A->GetRowStarts());
+      hypre_ParVector *d_hypre;
+      hypre_ParvecBdiagInvScal(*b_Hypre, block, &d_hypre, *A);
+      delete b_Hypre;
+
+      d->WrapHypreParVector(d_hypre);
+      d->SetOwnership(true);
+
+      return 0;
+   }
+
+   return -1;
+}
+#endif
+
 HypreParMatrix *Add(double alpha, const HypreParMatrix &A,
                     double beta,  const HypreParMatrix &B)
 {
@@ -2480,6 +2568,12 @@ void HypreGMRES::SetTol(double tol)
    HYPRE_GMRESSetTol(gmres_solver, tol);
 }
 
+void HypreGMRES::SetAbsTol(double tol)
+{
+   HYPRE_GMRESSetTol(gmres_solver, 0.0);
+   HYPRE_GMRESSetAbsoluteTol(gmres_solver, tol);
+}
+
 void HypreGMRES::SetMaxIter(int max_iter)
 {
    HYPRE_GMRESSetMaxIter(gmres_solver, max_iter);
@@ -2778,6 +2872,75 @@ HypreBoomerAMG::HypreBoomerAMG(HypreParMatrix &A) : HypreSolver(&A)
    SetDefaultOptions();
 }
 
+void HypreBoomerAMG::Mult(const HypreParVector &b, HypreParVector &x) const
+{
+   int myid;
+   HYPRE_Int time_index = 0;
+   HYPRE_Int num_iterations;
+   double final_res_norm;
+   MPI_Comm comm;
+   HYPRE_Int print_level;
+
+   HYPRE_BoomerAMGGetPrintLevel(amg_precond, &print_level);
+
+   HYPRE_ParCSRMatrixGetComm(*A, &comm);
+
+   if (!setup_called)
+   {
+      if (print_level > 0)
+      {
+         time_index = hypre_InitializeTiming("BoomerAMG Setup");
+         hypre_BeginTiming(time_index);
+      }
+
+      HYPRE_BoomerAMGSetup(amg_precond, *A, b, x);
+      setup_called = 1;
+
+      if (print_level > 0)
+      {
+         hypre_EndTiming(time_index);
+         hypre_PrintTiming("Setup phase times", comm);
+         hypre_FinalizeTiming(time_index);
+         hypre_ClearTiming();
+      }
+   }
+
+   if (print_level > 0)
+   {
+      time_index = hypre_InitializeTiming("BoomerAMG Solve");
+      hypre_BeginTiming(time_index);
+   }
+
+   if (!iterative_mode)
+   {
+      x = 0.0;
+   }
+
+   HYPRE_BoomerAMGSolve(amg_precond, *A, b, x);
+
+   if (print_level > 0)
+   {
+      hypre_EndTiming(time_index);
+      hypre_PrintTiming("Solve phase times", comm);
+      hypre_FinalizeTiming(time_index);
+      hypre_ClearTiming();
+
+      HYPRE_BoomerAMGGetNumIterations(amg_precond, &num_iterations);
+      HYPRE_BoomerAMGGetFinalRelativeResidualNorm(amg_precond,
+                                                  &final_res_norm);
+
+      MPI_Comm_rank(comm, &myid);
+
+      if (myid == 0)
+      {
+         mfem::out << "BoomerAMG Iterations = " << num_iterations << endl
+                   << "Final Relative Residual Norm = " << final_res_norm
+                   << endl;
+      }
+   }
+}
+
+
 void HypreBoomerAMG::SetDefaultOptions()
 {
    // AMG coarsening options:
@@ -3025,6 +3188,243 @@ void HypreBoomerAMG::SetElasticityOptions(ParFiniteElementSpace *fespace)
    error_mode = IGNORE_HYPRE_ERRORS;
 }
 
+#if MFEM_HYPRE_VERSION >= 21800
+
+void HypreBoomerAMG::SetLAIROptions(int distance,
+                                    std::string prerelax,
+                                    std::string postrelax,
+                                    double strength_tolC,
+                                    double strength_tolR,
+                                    double filter_tolR,
+                                    int interp_type,
+                                    int relax_type,
+                                    double filterA_tol,
+                                    int splitting,
+                                    int blksize,
+                                    int Sabs)
+{
+   int ns_down, ns_up, ns_coarse;
+   if (distance > 0)
+   {
+      ns_down = prerelax.length();
+      ns_up = postrelax.length();
+      ns_coarse = 1;
+      std::string F("F");
+      std::string C("C");
+      std::string A("A");
+
+      // Array to store relaxation scheme and pass to Hypre
+      int **grid_relax_points = (int **) malloc(4*sizeof(int *));
+      grid_relax_points[0] = NULL;
+      grid_relax_points[1] = (int *) malloc(sizeof(int)*ns_down);
+      grid_relax_points[2] = (int *) malloc(sizeof(int)*ns_up);
+      grid_relax_points[3] = (int *) malloc(sizeof(int));
+      grid_relax_points[3][0] = 0;
+
+      // set down relax scheme
+      for (unsigned int i = 0; i<ns_down; i++)
+      {
+         if (prerelax.compare(i,1,F) == 0)
+         {
+            grid_relax_points[1][i] = -1;
+         }
+         else if (prerelax.compare(i,1,C) == 0)
+         {
+            grid_relax_points[1][i] = 1;
+         }
+         else if (prerelax.compare(i,1,A) == 0)
+         {
+            grid_relax_points[1][i] = 0;
+         }
+      }
+
+      // set up relax scheme
+      for (unsigned int i = 0; i<ns_up; i++)
+      {
+         if (postrelax.compare(i,1,F) == 0)
+         {
+            grid_relax_points[2][i] = -1;
+         }
+         else if (postrelax.compare(i,1,C) == 0)
+         {
+            grid_relax_points[2][i] = 1;
+         }
+         else if (postrelax.compare(i,1,A) == 0)
+         {
+            grid_relax_points[2][i] = 0;
+         }
+      }
+
+      HYPRE_BoomerAMGSetRestriction(amg_precond, distance);
+
+      HYPRE_BoomerAMGSetGridRelaxPoints(amg_precond, grid_relax_points);
+
+      HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type);
+   }
+
+   if (Sabs)
+   {
+      HYPRE_BoomerAMGSetSabs(amg_precond, Sabs);
+   }
+
+   if (blksize > 0)
+   {
+      HYPRE_BoomerAMGSetNumFunctions(amg_precond, blksize);
+      HYPRE_BoomerAMGSetNodal(amg_precond, 1);
+   }
+
+   HYPRE_BoomerAMGSetCoarsenType(amg_precond, splitting);
+
+   /* does not support aggressive coarsening */
+   HYPRE_BoomerAMGSetAggNumLevels(amg_precond, 0);
+
+   HYPRE_BoomerAMGSetStrongThreshold(amg_precond, strength_tolC);
+
+   if (distance > 0)
+   {
+      HYPRE_BoomerAMGSetStrongThresholdR(amg_precond, strength_tolR);
+      HYPRE_BoomerAMGSetFilterThresholdR(amg_precond, filter_tolR);
+   }
+
+   if (relax_type > -1)
+   {
+      HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type);
+   }
+
+   if (distance > 0)
+   {
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_coarse, 3);
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_down,   1);
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_up,     2);
+
+      HYPRE_BoomerAMGSetADropTol(amg_precond, filterA_tol);
+      /* type = -1: drop based on row inf-norm */
+      HYPRE_BoomerAMGSetADropType(amg_precond, -1);
+   }
+}
+
+void HypreBoomerAMG::SetNAIROptions(int neumann_degree,
+                                    std::string prerelax,
+                                    std::string postrelax,
+                                    double strength_tolC,
+                                    double strength_tolR,
+                                    double filter_tolR,
+                                    int interp_type,
+                                    int relax_type,
+                                    double filterA_tol,
+                                    int splitting,
+                                    int blksize,
+                                    int Sabs)
+{
+   int ns_down, ns_up, ns_coarse;
+   if (neumann_degree > 0)
+   {
+      ns_down = prerelax.length();
+      ns_up = postrelax.length();
+      ns_coarse = 1;
+      std::string F("F");
+      std::string C("C");
+      std::string A("A");
+
+      // Array to store relaxation scheme and pass to Hypre
+      int **grid_relax_points = (int **) malloc(4*sizeof(int *));
+      grid_relax_points[0] = NULL;
+      grid_relax_points[1] = (int *) malloc(sizeof(int)*ns_down);
+      grid_relax_points[2] = (int *) malloc(sizeof(int)*ns_up);
+      grid_relax_points[3] = (int *) malloc(sizeof(int));
+      grid_relax_points[3][0] = 0;
+
+      // set down relax scheme
+      for (unsigned int i = 0; i<ns_down; i++)
+      {
+         if (prerelax.compare(i,1,F) == 0)
+         {
+            grid_relax_points[1][i] = -1;
+         }
+         else if (prerelax.compare(i,1,C) == 0)
+         {
+            grid_relax_points[1][i] = 1;
+         }
+         else if (prerelax.compare(i,1,A) == 0)
+         {
+            grid_relax_points[1][i] = 0;
+         }
+      }
+
+      // set up relax scheme
+      for (unsigned int i = 0; i<ns_up; i++)
+      {
+         if (postrelax.compare(i,1,F) == 0)
+         {
+            grid_relax_points[2][i] = -1;
+         }
+         else if (postrelax.compare(i,1,C) == 0)
+         {
+            grid_relax_points[2][i] = 1;
+         }
+         else if (postrelax.compare(i,1,A) == 0)
+         {
+            grid_relax_points[2][i] = 0;
+         }
+      }
+
+      HYPRE_BoomerAMGSetRestriction(amg_precond, 3+neumann_degree);
+
+      HYPRE_BoomerAMGSetGridRelaxPoints(amg_precond, grid_relax_points);
+
+      HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type);
+   }
+
+   if (Sabs)
+   {
+      HYPRE_BoomerAMGSetSabs(amg_precond, Sabs);
+   }
+
+   if (blksize > 0)
+   {
+      HYPRE_BoomerAMGSetNumFunctions(amg_precond, blksize);
+      HYPRE_BoomerAMGSetNodal(amg_precond, 1);
+   }
+
+   HYPRE_BoomerAMGSetCoarsenType(amg_precond, splitting);
+
+   /* does not support aggressive coarsening */
+   HYPRE_BoomerAMGSetAggNumLevels(amg_precond, 0);
+
+   HYPRE_BoomerAMGSetStrongThreshold(amg_precond, strength_tolC);
+
+   if (neumann_degree > 0)
+   {
+      HYPRE_BoomerAMGSetStrongThresholdR(amg_precond, strength_tolR);
+   }
+
+   if (relax_type > -1)
+   {
+      HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type);
+   }
+
+   if (neumann_degree > 0)
+   {
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_coarse, 3);
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_down,   1);
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_up,     2);
+
+      HYPRE_BoomerAMGSetADropTol(amg_precond, filterA_tol);
+      /* type = -1: drop based on row inf-norm */
+      HYPRE_BoomerAMGSetADropType(amg_precond, -1);
+   }
+}
+
+
+void HypreBoomerAMG::SetRelaxationOrdering(int *ordering)
+{
+   HYPRE_BoomerAMGSetRelaxType(amg_precond, 10);
+   relax_ordering = (HYPRE_Int *) ordering;
+   hypre_ParCSRMatrixProcOrdering(A) = relax_ordering;
+}
+
+#endif
+
 HypreBoomerAMG::~HypreBoomerAMG()
 {
    for (int i = 0; i < rbms.Size(); i++)
@@ -3032,6 +3432,12 @@ HypreBoomerAMG::~HypreBoomerAMG()
       HYPRE_ParVectorDestroy(rbms[i]);
    }
 
+   // If relaxation ordering provided to hypre, set internal
+   // hypre pointer to NULL. Memory is managed in MFEM.
+   if (relax_ordering) {
+      hypre_ParCSRMatrixProcOrdering(A) = NULL;
+   }
+
    HYPRE_BoomerAMGDestroy(amg_precond);
 }
 
@@ -4088,4 +4494,4 @@ HypreAME::StealEigenvectors()
 
 }
 
-#endif
+#endif

linalg/hypre.hpp

@@ -84,6 +84,9 @@ private:
    inline void _SetDataAndSize_();
 
 public:
+
+   HypreParVector() {}
+
    /** @brief Creates vector with given global size and parallel partitioning of
        the rows/columns given by @a col. */
    /** @anchor hypre_partitioning_descr
@@ -116,9 +119,9 @@ public:
    /// MPI communicator
    MPI_Comm GetComm() { return x->comm; }
 
-   /// Returns the parallel row/column partitioning
-   /** See @ref hypre_partitioning_descr "here" for a description of the
-       partitioning array. */
+   void WrapHypreParVector(hypre_ParVector *y);
+
+   /// Returns the row partitioning
    inline HYPRE_Int *Partitioning() { return x->partitioning; }
 
    /// Returns the global number of rows
@@ -234,15 +237,20 @@ public:
    /// An empty matrix to be used as a reference to an existing matrix
    HypreParMatrix();
 
+   void WrapHypreParCSRMatrix(hypre_ParCSRMatrix *a, bool owner = true)
+   {
+      A = a;
+      if (!owner) { ParCSROwner = 0; }
+      height = GetNumRows();
+      width = GetNumCols();
+   }
+
    /// Converts hypre's format to HypreParMatrix
    /** If @a owner is false, ownership of @a a is not transferred */
    explicit HypreParMatrix(hypre_ParCSRMatrix *a, bool owner = true)
    {
       Init();
-      A = a;
-      if (!owner) { ParCSROwner = 0; }
-      height = GetNumRows();
-      width = GetNumCols();
+      WrapHypreParCSRMatrix(a, owner);
    }
 
    /// Creates block-diagonal square parallel matrix.
@@ -334,6 +342,7 @@ public:
 
    /// MPI communicator
    MPI_Comm GetComm() const { return A->comm; }
+   void BuildComm() const { hypre_MatvecCommPkgCreate(A); }
 
    /// Typecasting to hypre's hypre_ParCSRMatrix*
    operator hypre_ParCSRMatrix*() const { return A; }
@@ -393,6 +402,8 @@ public:
    void GetDiag(SparseMatrix &diag) const;
    /// Get the local off-diagonal block. NOTE: 'offd' will not own any data.
    void GetOffd(SparseMatrix &offd, HYPRE_Int* &cmap) const;
+   /// Get on-processor rows as CSR matrix.
+   void GetProcRows(SparseMatrix &colCSRMat);
 
    /** Split the matrix into M x N equally sized blocks of parallel matrices.
        The size of 'blocks' must already be set to M x N. */
@@ -403,6 +414,13 @@ public:
    /// Returns the transpose of *this
    HypreParMatrix * Transpose() const;
 
+   /** Returns principle submatrix given by array of indices of connections
+       with relative size > \@ threshold in *this. */
+#if MFEM_HYPRE_VERSION >= 21800
+   HypreParMatrix * ExtractSubmatrix(Array<int> &indices,
+                                     double threshhold=0.0) const;
+#endif
+
    /// Returns the number of rows in the diagonal block of the ParCSRMatrix
    int GetNumRows() const
    {
@@ -549,6 +567,11 @@ public:
    Type GetType() const { return Hypre_ParCSR; }
 };
 
+#if MFEM_HYPRE_VERSION >= 21800
+int BlockInverseScale(const HypreParMatrix *A, HypreParMatrix *C,
+                      const Vector *b, HypreParVector *d, int block, int job);
+#endif
+
 /** @brief Return a new matrix `C = alpha*A + beta*B`, assuming that both `A`
     and `B` use the same row and column partitions and the same `col_map_offd`
     arrays. */
@@ -633,11 +656,12 @@ public:
        4    = truncated l1-scaled block Gauss-Seidel/SSOR
        5    = lumped Jacobi
        6    = Gauss-Seidel
+       10   = On-processor forward solve for matrix w/ triangular structure
        16   = Chebyshev
        1001 = Taubin polynomial smoother
        1002 = FIR polynomial smoother. */
    enum Type { Jacobi = 0, l1Jacobi = 1, l1GS = 2, l1GStr = 4, lumpedJacobi = 5,
-               GS = 6, Chebyshev = 16, Taubin = 1001, FIR = 1002
+               GS = 6, OPFS = 10, Chebyshev = 16, Taubin = 1001, FIR = 1002
              };
 
    HypreSmoother();
@@ -739,6 +763,26 @@ public:
    virtual ~HypreSolver();
 };
 
+
+#if MFEM_HYPRE_VERSION >= 21800
+/// Abstract class for hypre's solvers and preconditioners
+class HypreTriSolve : public HypreSolver
+{
+public:
+   HypreTriSolve() : HypreSolver() { }
+   explicit HypreTriSolve(HypreParMatrix &A) : HypreSolver(&A) { }
+   virtual operator HYPRE_Solver() const { return NULL; }
+
+   virtual HYPRE_PtrToParSolverFcn SetupFcn() const
+   { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSROnProcTriSetup; }
+   virtual HYPRE_PtrToParSolverFcn SolveFcn() const
+   { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSROnProcTriSolve; }
+
+   HypreParMatrix* GetData() { return A; }
+   virtual ~HypreTriSolve() { }
+};
+#endif
+
 /// PCG solver in hypre
 class HyprePCG : public HypreSolver
 {
@@ -813,6 +857,7 @@ public:
    virtual void SetOperator(const Operator &op);
 
    void SetTol(double tol);
+   void SetAbsTol(double tol);
    void SetMaxIter(int max_iter);
    void SetKDim(int dim);
    void SetLogging(int logging);
@@ -952,6 +997,7 @@ class HypreBoomerAMG : public HypreSolver
 {
 private:
    HYPRE_Solver amg_precond;
+   HYPRE_Int *relax_ordering;
 
    /// Rigid body modes
    Array<HYPRE_ParVector> rbms;
@@ -989,9 +1035,85 @@ public:
        As with SetSystemsOptions(), this solver assumes Ordering::byVDIM. */
    void SetElasticityOptions(ParFiniteElementSpace *fespace);
 
+   void SetRelaxationOrdering(int *ordering);
+
+#if MFEM_HYPRE_VERSION >= 21800
+   /* distance parameter takes on values {1,2,15} for lAIR, meaning R is built using
+   distance 1 neighbors, distance two neighbors, or distance two on processor and
+   distance 1 off processor (i.e., distance 1.5 --> 15).        */
+   void SetLAIROptions(int distance=15,  std::string prerelax="",
+                       std::string postrelax="FFC", double strength_tol=0.1,
+                       double strength_tolR=0.01, double filter_tolR=0.0,
+                       int interp_type=100, int relax_type=3, double filterA_tol=0.0,
+                       int splitting=6, int blksize=0, int Sabs=0);
+
+   void SetNAIROptions(int neumann_degree=2,  std::string prerelax="A",
+                       std::string postrelax="F", double strength_tol=0.1,
+                       double strength_tolR=0.01, double filter_tolR=0.0,
+                       int interp_type=100, int relax_type=10, double filterA_tol=0.0,
+                       int splitting=6, int blksize=0, int Sabs=0);
+
+   void SetStrengthThreshR(double strengthR)
+   { HYPRE_BoomerAMGSetStrongThresholdR(amg_precond, strengthR); }
+
+   void SetFilterThreshR(double filterR)
+   { HYPRE_BoomerAMGSetFilterThresholdR(amg_precond, filterR); }
+
+   void SetRestriction(int restrict_type)
+   { HYPRE_BoomerAMGSetRestriction(amg_precond, restrict_type); }
+
+   void SetTriangular()
+   { HYPRE_BoomerAMGSetIsTriangular(amg_precond, 1); }
+
+   void SetGMRESSwitchR(int gmres_switch)
+   { HYPRE_BoomerAMGSetGMRESSwitchR(amg_precond, gmres_switch); }
+
+   void SetRelaxCycle(int prerelax, int postrelax)
+   {
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, prerelax,  1);
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, postrelax, 2);
+   }
+#endif
+
    void SetPrintLevel(int print_level)
    { HYPRE_BoomerAMGSetPrintLevel(amg_precond, print_level); }
 
+   void SetMaxIter(int max_iter)
+   { HYPRE_BoomerAMGSetMaxIter(amg_precond, max_iter); }
+
+   void SetMaxLevels(int max_levels)
+   { HYPRE_BoomerAMGSetMaxLevels(amg_precond, max_levels); }
+
+   void SetTol(double tol)
+   { HYPRE_BoomerAMGSetTol(amg_precond, tol); }
+
+   void SetStrengthThresh(double strength)
+   { HYPRE_BoomerAMGSetStrongThreshold(amg_precond, strength); }
+
+   void SetInterpolation(int interp_type)
+   { HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type); }
+
+   void SetCoarsening(int coarsen_type)
+   { HYPRE_BoomerAMGSetCoarsenType(amg_precond, coarsen_type); }
+
+   void SetRelaxType(int relax_type)
+   { HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type); }
+
+   void SetCycleType(int cycle_type)
+   { HYPRE_BoomerAMGSetCycleType(amg_precond, cycle_type); }
+
+   void GetNumIterations(int &num_it)
+   { HYPRE_BoomerAMGGetNumIterations(amg_precond, &num_it); }
+
+   void SetNodal(int blocksize)
+   {
+      HYPRE_BoomerAMGSetNumFunctions(amg_precond, blocksize);
+      HYPRE_BoomerAMGSetNodal(amg_precond, 1);
+   }
+
+   void SetAggressiveCoarsening(int num_levels)
+   { HYPRE_BoomerAMGSetAggNumLevels(amg_precond, num_levels); }
+
    /// The typecast to HYPRE_Solver returns the internal amg_precond
    virtual operator HYPRE_Solver() const { return amg_precond; }
 
@@ -1000,6 +1122,9 @@ public:
    virtual HYPRE_PtrToParSolverFcn SolveFcn() const
    { return (HYPRE_PtrToParSolverFcn) HYPRE_BoomerAMGSolve; }
 
+   virtual void Mult (const HypreParVector &b, HypreParVector &x) const;
+   using HypreSolver::Mult;
+
    virtual ~HypreBoomerAMG();
 };