increase maximum number of quad pts in 3d

for EA capabilities

examples/ginkgo/CMakeLists.txt

@@ -12,6 +12,8 @@
 set(GINKGO_EXAMPLES_SRCS)
 list(APPEND GINKGO_EXAMPLES_SRCS
   ex1.cpp
+  ex1lor.cpp
+  ex1lorcoeff.cpp
 )
 
 # Include the source directory where mfem.hpp and mfem-performance.hpp are.
@@ -33,6 +35,8 @@ add_mfem_examples(GINKGO_EXAMPLES_SRCS ${PFX} "" test_ginkgo)
 
 # Command line options for the tests.
 set(EX1_COMMON_OPTS ex1 -m ../data/star.mesh --use_gko_solver)
+set(EX1LOR_COMMON_OPTS ex1lor -m ../../data/beam-hex.mesh -pc)
+
 set(EX1_TEST_OPTS ${EX9_COMMON_OPTS})
 
 # Add the tests: one test per source file.

examples/ginkgo/ex1lor.cpp

@@ -0,0 +1,721 @@
+//                          MFEM Example 1, modified
+//
+// This code has been modified from `ex1.cpp` provided in the examples of
+// MFEM.  The sections not marked as related to Ginkgo are largely unchanged
+// from the version provided by MFEM.
+//
+// This code also contains portions from `miniapps/performance/ex1.cpp`
+// pertaining to the LOR preconditioner.  The preconditioner types
+// used in this example are applied to the SparseMatrix created from
+// the LOR mesh, so they are all LOR preconditioners with different
+// subtypes.
+//
+// The default mesh option is "beam-hex.mesh", provided by MFEM.
+// Important non-default options:
+//  -m [file] : Mesh file.
+//  -d "cuda" : Use the MFEM cuda backend and Ginkgo CudaExecutor.
+//  -pc-type "gko:ilu" : Use the Ginkgo ILU preconditioner (default is Block
+//                       Jacobi)
+//  -pc-type "none" : No LOR preconditioner
+//
+//  Options only for the Block Jacobi preconditioner (default:)
+//  -pc-so "none" : Don't let Ginkgo automatically pick options for precision
+//                  reduction in the storage of the Block Jacobi preconditioner
+//  -pc-acc [value] : Accuracy parameter.
+//
+// MFEM's provided information about `ex1.cpp`:
+// Description:  This example code demonstrates the use of MFEM to define a
+//               simple finite element discretization of the Laplace problem
+//               -Delta u = 1 with homogeneous Dirichlet boundary conditions.
+//               Specifically, we discretize using a FE space of the specified
+//               order, or if order < 1 using an isoparametric/isogeometric
+//               space (i.e. quadratic for quadratic curvilinear mesh, NURBS for
+//               NURBS mesh, etc.)
+//
+//               The example highlights the use of mesh refinement, finite
+//               element grid functions, as well as linear and bilinear forms
+//               corresponding to the left-hand side and right-hand side of the
+//               discrete linear system. We also cover the explicit elimination
+//               of essential boundary conditions, static condensation, and the
+//               optional connection to the GLVis tool for visualization.
+
+#include "mfem.hpp"
+
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+void PermuteSparseMatrix(SparseMatrix &A, Array<int> &pinv, double shift=0.0)
+{
+   int n = pinv.Size();
+
+   SparseMatrix PA(A);
+
+   Array<int> p(n);
+   for (int i=0; i<n; ++i)
+   {
+      p[pinv[i]] = i;
+   }
+
+   // Set LU = A(P,P) using the permutation generated above
+   const int *IA = A.GetI();
+   const int *JA = A.GetJ();
+   const double *VA = A.GetData();
+
+   int *I = PA.GetI();
+   int *J = PA.GetJ();
+   double *V = PA.GetData();
+
+   I[0] = 0;
+   for (int i=0; i<n; ++i)
+   {
+      int pi = p[i];
+      int nnz_pi = IA[pi+1] - IA[pi];
+      I[i+1] = I[i] + nnz_pi;
+      for (int jj=0; jj<nnz_pi; ++jj)
+      {
+         int pj = JA[IA[pi] + jj];
+         int j = pinv[pj];
+
+         J[I[i] + jj] = j;
+         V[I[i] + jj] = VA[IA[pi] + jj];
+         if (i == j) { V[I[i] + jj] += shift; }
+      }
+   }
+
+   PA.SortColumnIndices();
+   PA.Swap(A);
+}
+
+// helper functions for cg/pcg solve with timer and iter count return
+int cg_solve(const Operator &A, const Vector &b, Vector &x,
+             int print_iter, int max_num_iter,
+             double RTOLERANCE, double ATOLERANCE, double &it_time)
+{
+
+   CGSolver cg;
+   cg.SetPrintLevel(print_iter);
+   cg.SetMaxIter(max_num_iter);
+   cg.SetRelTol(sqrt(RTOLERANCE));
+   cg.SetAbsTol(sqrt(ATOLERANCE));
+   cg.SetOperator(A);
+
+   tic_toc.Clear();
+   tic_toc.Start();
+
+   cg.Mult(b, x);
+
+   tic_toc.Stop();
+   it_time = tic_toc.RealTime();
+
+   return cg.GetNumIterations();
+}
+
+int pcg_solve(const Operator &A, Solver &B, const Vector &b, Vector &x,
+              int print_iter, int max_num_iter,
+              double RTOLERANCE, double ATOLERANCE, double &it_time)
+{
+
+   CGSolver pcg;
+   pcg.SetPrintLevel(print_iter);
+   pcg.SetMaxIter(max_num_iter);
+   pcg.SetRelTol(sqrt(RTOLERANCE));
+   pcg.SetAbsTol(sqrt(ATOLERANCE));
+   pcg.SetOperator(A);
+   pcg.SetPreconditioner(B);
+
+   tic_toc.Clear();
+   tic_toc.Start();
+
+   pcg.Mult(b, x);
+
+   tic_toc.Stop();
+   it_time = tic_toc.RealTime();
+
+   return pcg.GetNumIterations();
+}
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../../data/beam-hex.mesh";
+   int ref_levels = 3;
+   int order = 2;
+   const char *basis_type = "G"; // Gauss-Lobatto
+   bool static_cond = false;
+   bool pa = true;
+   const char *device_config = "cpu";
+   bool visualization = true;
+   const char *pc_type = "gko:bj";
+   const char *pc_storage_opt = "auto";
+   double pc_acc = 1.e-1;
+   int pc_max_bs = 32;
+   int permute = 0;
+   bool output_sol = false;
+   bool output_pc = false;
+   int isai_sparsity_power = 1;
+   int par_ilu_its = 0;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&ref_levels, "-l", "--refinement-levels",
+                  "Number of uniform refinement levels for mesh.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&pc_type, "-pc-type", "--preconditioner-type",
+                  "Type of preconditioner used on LOR matrix.");
+   args.AddOption(&pc_storage_opt, "-pc-so",
+                  "--preconditioner-storage-optimization",
+                  "Type of precision storage optimization to use for Ginkgo BlockJacobi.");
+   args.AddOption(&pc_acc, "-pc-acc", "--preconditioner-accuracy",
+                  "Accuracy parameter for Ginkgo BlockJacobi.");
+   args.AddOption(&pc_max_bs, "-pc-mbs", "--preconditioner-max-block-size",
+                  "Maximum block size for Ginkgo BlockJacobi.");
+   args.AddOption(&permute, "-per", "--permutation",
+                  "Specify preconditioner permutation.");
+   args.AddOption(&output_sol, "-out", "--output-solution-and-mesh", "-no-out",
+                  "--no-solution-and-mesh-output",
+                  "Output mesh and solution for inspection.");
+   args.AddOption(&output_pc, "-out-pc", "--output-lor-matrix-and-mesh",
+                  "-no-out-pc",
+                  "--no-lor-matrix-and-mesh-output",
+                  "Output LOR mesh and sparse matrix for inspection.");
+   args.AddOption(&isai_sparsity_power, "-isai-sp", "--isai-sparsity-power",
+                  "Power to use for sparsity pattern of ISAI in Ginkgo ILU-ISAI.");
+   args.AddOption(&par_ilu_its, "-pilu-its", "--par-ilu-iterations",
+                  "Number of iterations for the Ginkgo ParILU algorithm.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   enum PCType { NONE, GKO_BLOCK_JACOBI, GKO_ILU, GKO_ILU_ISAI, GKO_CUILU, GKO_CUILU_ISAI, MFEM_GS, MFEM_UMFPACK };
+   PCType pc_choice;
+   bool pc = true;
+   const char *trisolve_type = "exact"; //only used for ILU
+   if (!strcmp(pc_type, "gko:bj")) { pc_choice = GKO_BLOCK_JACOBI; }
+   else if (!strcmp(pc_type, "gko:ilu")) { pc_choice = GKO_ILU; }
+   else if (!strcmp(pc_type, "gko:ilu-isai"))
+   {
+      pc_choice = GKO_ILU_ISAI;
+      trisolve_type = "isai";
+   }
+   else if (!strcmp(pc_type, "gko:cuilu")) { pc_choice = GKO_CUILU; }
+   else if (!strcmp(pc_type, "gko:cuilu-isai"))
+   {
+      pc_choice = GKO_CUILU_ISAI;
+      trisolve_type = "isai";
+   }
+   else if (!strcmp(pc_type, "mfem:gs")) { pc_choice = MFEM_GS; }
+   else if (!strcmp(pc_type, "mfem:umf"))
+   {
+#ifdef MFEM_USE_SUITESPARSE
+      pc_choice = MFEM_UMFPACK;
+#else
+      mfem_error("Preconditioner requires SuiteSparse");
+#endif
+   }
+   else if (!strcmp(pc_type, "none"))
+   {
+      pc_choice = NONE;
+      pc = false;
+   }
+   else
+   {
+      mfem_error("Invalid Preconditioner specified");
+      return 3;
+   }
+
+   // ---------------------------------------------------------------
+   // -------------------- Start Ginkgo set-up ----------------------
+
+   // Create Ginkgo executor.
+
+   // This will point to the selected executor default executor
+   std::shared_ptr<gko::Executor> executor;
+
+   // We will always need an OpenMP executor.
+   auto omp_executor = gko::OmpExecutor::create();
+
+   // If the user has requested to use CUDA, then build a
+   // CudaExecutor and set `executor` to it; otherwise,
+   // use the OmpExecutor
+   if (!strcmp(device_config, "cuda"))
+   {
+      auto cuda_executor =
+         gko::CudaExecutor::create(0, gko::OmpExecutor::create());
+      executor = cuda_executor;
+   }
+   else
+   {
+      executor = omp_executor;
+   }
+
+   // --------------------- End Ginkgo set-up -----------------------
+   // ---------------------------------------------------------------
+
+   // See class BasisType in fem/fe_coll.hpp for available basis types
+   int basis = BasisType::GetType(basis_type[0]);
+   cout << "Using " << BasisType::Name(basis) << " basis ..." << endl;
+
+   // 2. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   device.Print();
+
+   // 3. Read the mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
+   //    the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 4. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement.
+   {
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh->UniformRefinement();
+      }
+   }
+   cout << "Total elements in refined mesh: " <<  mesh->GetNE() << std::endl;
+
+   // 5. Define a finite element space on the mesh. Here we use continuous
+   //    Lagrange finite elements of the specified order. If order < 1, we
+   //    instead use an isoparametric/isogeometric space.
+   FiniteElementCollection *fec;
+   if (order > 0)
+   {
+      fec = new H1_FECollection(order, dim, basis);
+   }
+   else if (mesh->GetNodes())
+   {
+      fec = mesh->GetNodes()->OwnFEC();
+      cout << "Using isoparametric FEs: " << fec->Name() << endl;
+   }
+   else
+   {
+      fec = new H1_FECollection(order = 1, dim, basis);
+   }
+   FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
+   cout << "Number of finite element unknowns: " << fespace->GetTrueVSize()
+        << endl;
+
+   // Create the LOR mesh and finite element space. In the settings of this
+   // example, we can transfer between HO and LOR with the identity operator.
+   Mesh *mesh_lor = NULL;
+   FiniteElementCollection *fec_lor = NULL;
+   FiniteElementSpace *fespace_lor = NULL;
+   Array<int> *inv_reordering = NULL;
+   if (pc)
+   {
+      int basis_lor = basis;
+      if (basis == BasisType::Positive) { basis_lor=BasisType::ClosedUniform; }
+      mesh_lor = new Mesh(mesh, order, basis_lor);
+      fec_lor = new H1_FECollection(1, dim);
+      fespace_lor = new FiniteElementSpace(mesh_lor, fec_lor);
+
+      if (permute == 1)
+      {
+
+         tic_toc.Clear();
+         tic_toc.Start();
+
+         const Table &pre_reorder_dofs = fespace_lor->GetElementToDofTable();
+         const Table pre_reorder_dofs_copy(pre_reorder_dofs);
+         fespace_lor->ReorderElementToDofTable();
+         const Table &post_reorder_dofs = fespace_lor->GetElementToDofTable();
+
+         inv_reordering = new Array<int>(fespace_lor->GetTrueVSize());
+         for (int i = 0; i < pre_reorder_dofs.Size(); i++)
+         {
+
+            Array<int> old_row;
+            Array<int> new_row;
+            pre_reorder_dofs_copy.GetRow(i, old_row);
+            post_reorder_dofs.GetRow(i, new_row);
+            for (int j = 0; j < pre_reorder_dofs_copy.RowSize(i); j++)
+            {
+               int new_dof = new_row[j];
+               int old_dof = old_row[j];
+               (*inv_reordering)[old_dof] = new_dof;
+            }
+         }
+
+         tic_toc.Stop();
+         cout << "Real time spent reordering: " <<
+              tic_toc.RealTime() << "\n";
+      }
+   }
+
+   // 6. Determine the list of true (i.e. conforming) essential boundary dofs.
+   //    In this example, the boundary conditions are defined by marking all
+   //    the boundary attributes from the mesh as essential (Dirichlet) and
+   //    converting them to a list of true dofs.
+   Array<int> ess_tdof_list;
+   if (mesh->bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(mesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   Array<int> ess_pc_tdof_list(ess_tdof_list.Size());
+
+   if (permute == 1)
+   {
+      for (int i = 0; i < ess_tdof_list.Size(); i++)
+      {
+         ess_pc_tdof_list.operator[](i) = inv_reordering->operator[](
+                                             ess_tdof_list.operator[](i));
+      }
+   }
+   else
+   {
+
+      Array<int> ess_bdr(mesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_pc_tdof_list);
+   }
+
+   // Array<int> ess_bdr(mesh->bdr_attributes.Max());
+   // ess_bdr = 1;
+   // fespace->GetEssentialTrueDofs(ess_bdr, ess_pc_tdof_list);
+
+   // 7. Set up the linear form b(.) which corresponds to the right-hand side
+   // of
+   //    the FEM linear system, which in this case is (1,phi_i) where phi_i are
+   //    the basis functions in the finite element fespace.
+   LinearForm *b = new LinearForm(fespace);
+
+   ConstantCoefficient one(1.0);
+   b->AddDomainIntegrator(new DomainLFIntegrator(one));
+   b->Assemble();
+
+   // 8. Define the solution vector x as a finite element grid function
+   //    corresponding to fespace. Initialize x with initial guess of zero,
+   //    which satisfies the boundary conditions.
+   GridFunction x(fespace);
+   x = 0.0;
+
+   // 9. Set up the bilinear form a(.,.) on the finite element space
+   //    corresponding to the Laplacian operator -Delta, by adding the
+   //    Diffusion domain integrator.
+   BilinearForm *a = new BilinearForm(fespace);
+   BilinearForm *a_pc = NULL;
+   if (pc) { a_pc = new BilinearForm(fespace_lor); }
+   if (pa)
+   {
+      a->SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   }
+   a->AddDomainIntegrator(new DiffusionIntegrator(one));
+
+   // 10. Assemble the bilinear form and the corresponding linear system,
+   //     applying any necessary transformations such as: eliminating boundary
+   //     conditions, applying conforming constraints for non-conforming AMR,
+   //     static condensation, etc.
+   if (static_cond)
+   {
+      a->EnableStaticCondensation();
+   }
+   a->Assemble();
+
+   OperatorPtr A;
+   Vector B, X;
+
+   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
+
+   // 11. Solve the linear system A X = B.
+   double it_time = 0.;
+   int total_its = 0;
+
+   SparseMatrix A_pc;
+   if (pc)
+   {
+
+      tic_toc.Clear();
+      tic_toc.Start();
+
+      a_pc->AddDomainIntegrator(new DiffusionIntegrator(one));
+      a_pc->UsePrecomputedSparsity();
+      a_pc->Assemble();
+
+      tic_toc.Stop();
+      cout << "Real time assembling A_pc SparseMatrix: " <<
+           tic_toc.RealTime() << "\n";
+
+      tic_toc.Clear();
+      tic_toc.Start();
+
+      a_pc->FormSystemMatrix(ess_pc_tdof_list, A_pc);
+
+      tic_toc.Stop();
+      cout << "Real time adjusting A_pc for essential BC: " <<
+           tic_toc.RealTime() << "\n";
+
+
+      if (permute == 2)
+      {
+         tic_toc.Clear();
+         tic_toc.Start();
+
+         Array<int> perm(fespace_lor->GetTrueVSize());
+         SparseMatrix A_pc_tmp(A_pc);
+         MinimumDiscardedFillOrdering(A_pc_tmp, perm);
+         inv_reordering = new Array<int>(fespace_lor->GetTrueVSize());
+         for (int i=0; i<perm.Size(); ++i)
+         {
+            (*inv_reordering)[perm[i]] = i;
+         }
+         PermuteSparseMatrix(A_pc, *inv_reordering);
+
+         tic_toc.Stop();
+         cout << "Real time spent reordering: " <<
+              tic_toc.RealTime() << "\n";
+      }
+
+
+      if (pc_choice == GKO_BLOCK_JACOBI)
+      {
+
+         // Create Ginkgo Jacobi preconditioner
+         if (permute)
+         {
+            tic_toc.Clear();
+            tic_toc.Start();
+            GinkgoWrappers::GinkgoJacobiPreconditioner M(executor, A_pc, *inv_reordering,
+                                                         pc_storage_opt,
+                                                         pc_acc, pc_max_bs);
+            tic_toc.Stop();
+            cout << "Real time creating Ginkgo BlockJacobi preconditioner: " <<
+                 tic_toc.RealTime() << "\n";
+
+            // Use preconditioned CG
+            total_its = pcg_solve(*A, M, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+            cout << "Real time in PCG: " << it_time << "\n";
+         }
+         else
+         {
+            tic_toc.Clear();
+            tic_toc.Start();
+            GinkgoWrappers::GinkgoJacobiPreconditioner M(executor, A_pc, pc_storage_opt,
+                                                         pc_acc, pc_max_bs);
+            tic_toc.Stop();
+            cout << "Real time creating Ginkgo BlockJacobi preconditioner: " <<
+                 tic_toc.RealTime() << "\n";
+
+            // Use preconditioned CG
+            total_its = pcg_solve(*A, M, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+            cout << "Real time in PCG: " << it_time << "\n";
+         }
+
+      }
+      else if (pc_choice == GKO_ILU || pc_choice == GKO_ILU_ISAI)
+      {
+
+         // Create Ginkgo ILU preconditioner
+
+         if (permute)
+         {
+
+            tic_toc.Clear();
+            tic_toc.Start();
+
+            GinkgoWrappers::GinkgoIluPreconditioner M(executor, A_pc, *inv_reordering,
+                                                      trisolve_type, isai_sparsity_power, par_ilu_its);
+
+            tic_toc.Stop();
+            cout << "Real time creating Ginkgo Ilu preconditioner: " <<
+                 tic_toc.RealTime() << "\n";
+
+            // Use preconditioned CG
+            total_its = pcg_solve(*A, M, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+            cout << "Real time in PCG: " << it_time << "\n";
+
+         }
+         else
+         {
+
+            tic_toc.Clear();
+            tic_toc.Start();
+
+            GinkgoWrappers::GinkgoIluPreconditioner M(executor, A_pc, trisolve_type,
+                                                      isai_sparsity_power, par_ilu_its);
+
+            tic_toc.Stop();
+            cout << "Real time creating Ginkgo Ilu preconditioner: " <<
+                 tic_toc.RealTime() << "\n";
+
+            // Use preconditioned CG
+            total_its = pcg_solve(*A, M, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+            cout << "Real time in PCG: " << it_time << "\n";
+
+         }
+      }
+      else if (pc_choice == GKO_CUILU || pc_choice == GKO_CUILU_ISAI)
+      {
+
+         // Create Ginkgo CuILU preconditioner (uses cuSPARSE for factorization)
+
+         if (permute)
+         {
+
+            tic_toc.Clear();
+            tic_toc.Start();
+
+            GinkgoWrappers::GinkgoCuIluPreconditioner M(executor, A_pc, *inv_reordering,
+                                                      trisolve_type, isai_sparsity_power);
+
+            tic_toc.Stop();
+            cout << "Real time creating Ginkgo CuIlu preconditioner: " <<
+                 tic_toc.RealTime() << "\n";
+
+            // Use preconditioned CG
+            total_its = pcg_solve(*A, M, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+            cout << "Real time in PCG: " << it_time << "\n";
+
+         }
+         else
+         {
+
+            tic_toc.Clear();
+            tic_toc.Start();
+
+            GinkgoWrappers::GinkgoCuIluPreconditioner M(executor, A_pc, trisolve_type,
+                                                      isai_sparsity_power);
+
+            tic_toc.Stop();
+            cout << "Real time creating Ginkgo CuIlu preconditioner: " <<
+                 tic_toc.RealTime() << "\n";
+
+            // Use preconditioned CG
+            total_its = pcg_solve(*A, M, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+            cout << "Real time in PCG: " << it_time << "\n";
+
+         }
+      }
+      else if (pc_choice == MFEM_GS)
+      {
+
+         // Create MFEM preconditioner
+         tic_toc.Clear();
+         tic_toc.Start();
+
+         GSSmoother M(A_pc);
+
+         tic_toc.Stop();
+         cout << "Real time creating MFEM GS preconditioner: " <<
+              tic_toc.RealTime() << "\n";
+
+         // Use preconditioned CG
+         total_its = pcg_solve(*A, M, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+         cout << "Real time in PCG: " << it_time << "\n";
+
+      }
+      else if (pc_choice == MFEM_UMFPACK)
+      {
+
+#ifdef MFEM_USE_SUITESPARSE
+         // Create MFEM preconditioner
+         tic_toc.Clear();
+         tic_toc.Start();
+
+         UMFPackSolver M;
+         M.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
+         M.SetOperator(A_pc);
+
+         tic_toc.Stop();
+         cout << "Real time creating MFEM UMFPACK preconditioner: " <<
+              tic_toc.RealTime() << "\n";
+
+         // Use preconditioned CG
+         total_its = pcg_solve(*A, M, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+         cout << "Real time in PCG: " << it_time << "\n";
+#endif
+      }
+   }
+   else
+   {
+
+      total_its = cg_solve(*A, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+      cout << "Real time in CG: " << it_time << "\n";
+   }
+
+   cout << "Total iterations: " << total_its << "\n";
+   cout << "Avg time per iteration: " << it_time/double(total_its) << "\n";
+
+   // 12. Recover the solution as a finite element grid function.
+   a->RecoverFEMSolution(X, *b, x);
+
+   // 13. Save the refined mesh and the solution. This output can be viewed
+   // later
+   //     using GLVis: "glvis -m refined.mesh -g sol.gf".
+
+   if (output_sol)
+   {
+      ofstream mesh_ofs("refined.mesh");
+      mesh_ofs.precision(8);
+      mesh->Print(mesh_ofs);
+      ofstream sol_ofs("sol.gf");
+      sol_ofs.precision(8);
+      x.Save(sol_ofs);
+   }
+
+   if (pc && output_pc)
+   {
+      ofstream mesh_lor_ofs("lor-refined.mesh");
+      mesh_lor_ofs.precision(8);
+      mesh_lor->Print(mesh_lor_ofs);
+
+      ofstream apc_lor_ofs("lor-mat.dat");
+      mesh_lor_ofs.precision(8);
+      A_pc.PrintCSR(apc_lor_ofs);
+   }
+
+   // 14. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int visport = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *mesh << x << flush;
+   }
+
+   // 15. Free the used memory.
+   delete a;
+   delete b;
+   delete fespace;
+   delete fespace_lor;
+   delete fec_lor;
+   delete mesh_lor;
+   if (order > 0)
+   {
+      delete fec;
+   }
+   delete mesh;
+}
+

examples/ginkgo/ex1lorcoeff.cpp

@@ -0,0 +1,750 @@
+//                          MFEM Example 1, modified
+//
+// This code has been modified from `ex1.cpp` provided in the examples of
+// MFEM.  The sections not marked as related to Ginkgo are largely unchanged
+// from the version provided by MFEM.
+//
+// This code also contains portions from `miniapps/performance/ex1.cpp`
+// pertaining to the LOR preconditioner.  The preconditioner types
+// used in this example are applied to the SparseMatrix created from
+// the LOR mesh, so they are all LOR preconditioners with different
+// subtypes.
+//
+// The default mesh option is "star.mesh", provided by MFEM.
+// Important non-default options:
+//  -m [file] : Mesh file.
+//  -d "cuda" : Use the MFEM cuda backend and Ginkgo CudaExecutor.
+//  -pc-type "gko:ilu" : Use the Ginkgo ILU preconditioner (default is Block
+//                       Jacobi)
+//  -pc-type "none" : No LOR preconditioner
+//
+//  Options only for the Block Jacobi preconditioner (default:)
+//  -pc-so "none" : Don't let Ginkgo automatically pick options for precision
+//                  reduction in the storage of the Block Jacobi preconditioner
+//  -pc-acc [value] : Accuracy parameter.
+//
+// MFEM's provided information about `ex1.cpp`:
+// Description:  This example code demonstrates the use of MFEM to define a
+//               simple finite element discretization of the Laplace problem
+//               -Delta u = 1 with homogeneous Dirichlet boundary conditions.
+//               Specifically, we discretize using a FE space of the specified
+//               order, or if order < 1 using an isoparametric/isogeometric
+//               space (i.e. quadratic for quadratic curvilinear mesh, NURBS for
+//               NURBS mesh, etc.)
+//
+//               The example highlights the use of mesh refinement, finite
+//               element grid functions, as well as linear and bilinear forms
+//               corresponding to the left-hand side and right-hand side of the
+//               discrete linear system. We also cover the explicit elimination
+//               of essential boundary conditions, static condensation, and the
+//               optional connection to the GLVis tool for visualization.
+
+#include "mfem.hpp"
+#include "../../general/forall.hpp"
+
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+double coeff_func(const Vector & x);
+
+// helper functions for cg/pcg solve with timer and iter count return
+int cg_solve(const Operator &A, const Vector &b, Vector &x,
+             int print_iter, int max_num_iter,
+             double RTOLERANCE, double ATOLERANCE, double &it_time)
+{
+
+   CGSolver cg;
+   cg.SetPrintLevel(print_iter);
+   cg.SetMaxIter(max_num_iter);
+   cg.SetRelTol(sqrt(RTOLERANCE));
+   cg.SetAbsTol(sqrt(ATOLERANCE));
+   cg.SetOperator(A);
+
+   tic_toc.Clear();
+   tic_toc.Start();
+
+   cg.Mult(b, x);
+
+   tic_toc.Stop();
+   it_time = tic_toc.RealTime();
+
+   return cg.GetNumIterations();
+}
+
+int pcg_solve(const Operator &A, Solver &B, const Vector &b, Vector &x,
+              int print_iter, int max_num_iter,
+              double RTOLERANCE, double ATOLERANCE, double &it_time)
+{
+
+   CGSolver pcg;
+   pcg.SetPrintLevel(print_iter);
+   pcg.SetMaxIter(max_num_iter);
+   pcg.SetRelTol(sqrt(RTOLERANCE));
+   pcg.SetAbsTol(sqrt(ATOLERANCE));
+   pcg.SetOperator(A);
+   pcg.SetPreconditioner(B);
+
+   tic_toc.Clear();
+   tic_toc.Start();
+
+   pcg.Mult(b, x);
+
+   tic_toc.Stop();
+   it_time = tic_toc.RealTime();
+
+   return pcg.GetNumIterations();
+}
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../../data/star.mesh";
+   int ref_levels = 3;
+   const char *coeff_name = "var";
+   int order = 2;
+   const char *basis_type = "G"; // Gauss-Lobatto
+   bool static_cond = false;
+   bool pa = true;
+   const char *device_config = "cpu";
+   bool visualization = false;
+   const char *pc_type = "gko:bj";
+   const char *pc_storage_opt = "auto";
+   double pc_acc = 1.e-1;
+   int pc_max_bs = 32;
+   bool permute = false;
+   bool output_sol = false;
+   bool output_pc = false;
+   int isai_sparsity_power = 1;
+   int par_ilu_its = 0;
+   double sigma_val = 1.0;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&ref_levels, "-l", "--refinement-levels",
+                  "Number of uniform refinement levels for mesh.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&coeff_name, "-c", "--coeff",
+                  "Type of coefficient for Laplace operator.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&pc_type, "-pc-type", "--preconditioner-type",
+                  "Type of preconditioner used on LOR matrix.");
+   args.AddOption(&pc_storage_opt, "-pc-so",
+                  "--preconditioner-storage-optimization",
+                  "Type of precision storage optimization to use for Ginkgo BlockJacobi.");
+   args.AddOption(&pc_acc, "-pc-acc", "--preconditioner-accuracy",
+                  "Accuracy parameter for Ginkgo BlockJacobi.");
+   args.AddOption(&pc_max_bs, "-pc-mbs", "--preconditioner-max-block-size",
+                  "Maximum block size for Ginkgo BlockJacobi.");
+   args.AddOption(&permute, "-per", "--permutation", "-no-per",
+                  "--no-permutation", "Enable preconditioner permutation.");
+   args.AddOption(&output_sol, "-out", "--output-solution-and-mesh", "-no-out",
+                  "--no-solution-and-mesh-output",
+                  "Output mesh and solution for inspection.");
+   args.AddOption(&output_pc, "-out-pc", "--output-lor-matrix-and-mesh",
+                  "-no-out-pc",
+                  "--no-lor-matrix-and-mesh-output",
+                  "Output LOR mesh and sparse matrix for inspection.");
+   args.AddOption(&isai_sparsity_power, "-isai-sp", "--isai-sparsity-power",
+                  "Power to use for sparsity pattern of ISAI in Ginkgo ILU-ISAI.");
+   args.AddOption(&par_ilu_its, "-pilu-its", "--par-ilu-iterations",
+                  "Number of iterations for the Ginkgo ParILU algorithm.");
+   args.AddOption(&sigma_val, "-sv", "--sigma-value",
+                  "Non-unity value for piecewise discontinuous coefficient.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   enum CoeffType { CONSTANT, VARIABLE, PW_CONSTANT };
+   CoeffType coeff_type;
+   if (!strcmp(coeff_name, "const")) { coeff_type = CONSTANT; }
+   else if (!strcmp(coeff_name, "var")) { coeff_type = VARIABLE; }
+   else if (!strcmp(coeff_name, "pw")) { coeff_type = PW_CONSTANT; }
+   else
+   {
+      mfem_error("Invalid coefficient type specified");
+   }
+
+   enum PCType { NONE, GKO_BLOCK_JACOBI, GKO_ILU, GKO_ILU_ISAI, GKO_CUILU, GKO_CUILU_ISAI, MFEM_GS, MFEM_UMFPACK };
+   PCType pc_choice;
+   bool pc = true;
+   const char *trisolve_type = "exact"; //only used for ILU
+   if (!strcmp(pc_type, "gko:bj")) { pc_choice = GKO_BLOCK_JACOBI; }
+   else if (!strcmp(pc_type, "gko:ilu")) { pc_choice = GKO_ILU; }
+   else if (!strcmp(pc_type, "gko:ilu-isai"))
+   {
+      pc_choice = GKO_ILU_ISAI;
+      trisolve_type = "isai";
+   }
+   else if (!strcmp(pc_type, "gko:cuilu")) { pc_choice = GKO_CUILU; }
+   else if (!strcmp(pc_type, "gko:cuilu-isai"))
+   {
+      pc_choice = GKO_CUILU_ISAI;
+      trisolve_type = "isai";
+   }
+   else if (!strcmp(pc_type, "mfem:gs")) { pc_choice = MFEM_GS; }
+   else if (!strcmp(pc_type, "mfem:umf"))
+   {
+#ifdef MFEM_USE_SUITESPARSE
+      pc_choice = MFEM_UMFPACK;
+#else
+      mfem_error("Preconditioner requires SuiteSparse");
+#endif
+   }
+   else if (!strcmp(pc_type, "none"))
+   {
+      pc_choice = NONE;
+      pc = false;
+   }
+   else
+   {
+      mfem_error("Invalid Preconditioner specified");
+      return 3;
+   }
+
+   // ---------------------------------------------------------------
+   // -------------------- Start Ginkgo set-up ----------------------
+
+   // Create Ginkgo executor.
+
+   // This will point to the selected executor default executor
+   std::shared_ptr<gko::Executor> executor;
+
+   // We will always need an OpenMP executor.
+   auto omp_executor = gko::OmpExecutor::create();
+
+   // If the user has requested to use CUDA, then build a
+   // CudaExecutor and set `executor` to it; otherwise,
+   // use the OmpExecutor
+   if (!strcmp(device_config, "cuda"))
+   {
+      auto cuda_executor =
+         gko::CudaExecutor::create(0, gko::OmpExecutor::create());
+      executor = cuda_executor;
+   }
+   else
+   {
+      executor = omp_executor;
+   }
+
+   // --------------------- End Ginkgo set-up -----------------------
+   // ---------------------------------------------------------------
+
+   // See class BasisType in fem/fe_coll.hpp for available basis types
+   int basis = BasisType::GetType(basis_type[0]);
+   cout << "Using " << BasisType::Name(basis) << " basis ..." << endl;
+
+   // 2. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   device.Print();
+
+   // 3. Read the mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
+   //    the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 4. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement.
+   {
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh->UniformRefinement();
+      }
+   }
+   cout << "Total elements in refined mesh: " <<  mesh->GetNE() << endl;
+
+   // 5. Define a finite element space on the mesh. Here we use continuous
+   //    Lagrange finite elements of the specified order. If order < 1, we
+   //    instead use an isoparametric/isogeometric space.
+   FiniteElementCollection *fec;
+   if (order > 0)
+   {
+      fec = new H1_FECollection(order, dim, basis);
+   }
+   else if (mesh->GetNodes())
+   {
+      fec = mesh->GetNodes()->OwnFEC();
+      cout << "Using isoparametric FEs: " << fec->Name() << endl;
+   }
+   else
+   {
+      fec = new H1_FECollection(order = 1, dim, basis);
+   }
+   FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
+   cout << "Number of finite element unknowns: " << fespace->GetTrueVSize()
+        << endl;
+
+   // Create the LOR mesh and finite element space. In the settings of this
+   // example, we can transfer between HO and LOR with the identity operator.
+   Mesh *mesh_lor = NULL;
+   FiniteElementCollection *fec_lor = NULL;
+   FiniteElementSpace *fespace_lor = NULL;
+   Array<int> *inv_reordering = NULL;
+   if (pc)
+   {
+      int basis_lor = basis;
+      if (basis == BasisType::Positive) { basis_lor=BasisType::ClosedUniform; }
+      mesh_lor = new Mesh(mesh, order, basis_lor);
+      fec_lor = new H1_FECollection(1, dim);
+      fespace_lor = new FiniteElementSpace(mesh_lor, fec_lor);
+
+      if (permute)
+      {
+
+         tic_toc.Clear();
+         tic_toc.Start();
+
+         const Table &pre_reorder_dofs = fespace_lor->GetElementToDofTable();
+         const Table pre_reorder_dofs_copy(pre_reorder_dofs);
+         fespace_lor->ReorderElementToDofTable();
+         const Table &post_reorder_dofs = fespace_lor->GetElementToDofTable();
+
+         inv_reordering = new Array<int>(fespace_lor->GetTrueVSize());
+         for (int i = 0; i < pre_reorder_dofs.Size(); i++)
+         {
+
+            Array<int> old_row;
+            Array<int> new_row;
+            pre_reorder_dofs_copy.GetRow(i, old_row);
+            post_reorder_dofs.GetRow(i, new_row);
+            for (int j = 0; j < pre_reorder_dofs_copy.RowSize(i); j++)
+            {
+               int new_dof = new_row[j];
+               int old_dof = old_row[j];
+               (*inv_reordering)[old_dof] = new_dof;
+            }
+         }
+
+         tic_toc.Stop();
+         cout << "Real time spent reordering: " <<
+              tic_toc.RealTime() << "\n";
+      }
+   }
+
+   // 6. Determine the list of true (i.e. conforming) essential boundary dofs.
+   //    In this example, the boundary conditions are defined by marking all
+   //    the boundary attributes from the mesh as essential (Dirichlet) and
+   //    converting them to a list of true dofs.
+   Array<int> ess_tdof_list;
+   if (mesh->bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(mesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   Array<int> ess_pc_tdof_list(ess_tdof_list.Size());
+   if (permute)
+   {
+      for (int i = 0; i < ess_tdof_list.Size(); i++)
+      {
+         ess_pc_tdof_list.operator[](i) = inv_reordering->operator[](
+                                             ess_tdof_list.operator[](i));
+      }
+   }
+   else
+   {
+
+      Array<int> ess_bdr(mesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_pc_tdof_list);
+   }
+
+   // 7. Set up the linear form b(.) which corresponds to the right-hand side
+   // of
+   //    the FEM linear system, which in this case is (1,phi_i) where phi_i are
+   //    the basis functions in the finite element fespace.
+   LinearForm *b = new LinearForm(fespace);
+
+   ConstantCoefficient one(1.0);
+   b->AddDomainIntegrator(new DomainLFIntegrator(one));
+   b->Assemble();
+
+   // 8. Define the solution vector x as a finite element grid function
+   //    corresponding to fespace. Initialize x with initial guess of zero,
+   //    which satisfies the boundary conditions.
+   GridFunction x(fespace);
+   x = 0.0;
+
+   // 9. Set up the bilinear form a(.,.) on the finite element space
+   //    corresponding to the Laplacian operator -Delta, by adding the
+   //    Diffusion domain integrator.
+   BilinearForm *a = new BilinearForm(fespace);
+   BilinearForm *a_pc = NULL;
+   if (pc) { a_pc = new BilinearForm(fespace_lor); }
+   if (pa)
+   {
+      a->SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   }
+
+   Coefficient *coeff = NULL;
+
+   if (coeff_type == CONSTANT)
+   {
+      coeff = new ConstantCoefficient(1.0);
+   }
+   else if (coeff_type == VARIABLE)
+   {
+      coeff = new FunctionCoefficient(&coeff_func);
+   }
+   else if (coeff_type == PW_CONSTANT)
+   {
+      int num_subregions = mesh->attributes.Max();
+      cout << "Number of subregions in mesh: " << num_subregions << "\n";
+      Vector sigma(num_subregions);
+      sigma = 1;
+      if (num_subregions < 2)
+      {
+         cout << "Warning: PW Constant Coefficient not used, mesh only has one element attribute!\n";
+      }
+      else
+      {
+         sigma(num_subregions-1) = sigma_val;
+      }
+      coeff = new PWConstCoefficient(sigma);
+   }
+   a->AddDomainIntegrator(new DiffusionIntegrator(*coeff));
+
+   // 10. Assemble the bilinear form and the corresponding linear system,
+   //     applying any necessary transformations such as: eliminating boundary
+   //     conditions, applying conforming constraints for non-conforming AMR,
+   //     static condensation, etc.
+   if (static_cond)
+   {
+      a->EnableStaticCondensation();
+   }
+   a->Assemble();
+
+   OperatorPtr A;
+   Vector B, X;
+
+   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
+
+   // 11. Solve the linear system A X = B.
+   double it_time = 0.;
+   int total_its = 0;
+
+   bool on_dev = false;
+   if (!strcmp(device_config, "cuda"))
+   {
+      on_dev = true;
+   }
+   SparseMatrix *A_pc;
+   if (pc)
+   {
+      tic_toc.Clear();
+      tic_toc.Start();
+
+      a_pc->AddDomainIntegrator(new DiffusionIntegrator(*coeff));
+      a_pc->SetAssemblyLevel(AssemblyLevel::FULL);
+      a_pc->Assemble();
+      A_pc = &(a_pc->SpMat());
+
+      tic_toc.Stop();
+      cout << "Real time assembling A_pc SparseMatrix: " <<
+           tic_toc.RealTime() << "\n";
+
+      tic_toc.Clear();
+      tic_toc.Start();
+
+      // Manually set essential BC rows/columns, so we can use device if applicable
+      // Get device or host pointers
+      auto d_I = A_pc->ReadWriteI(on_dev);
+      auto d_J = A_pc->ReadWriteJ(on_dev);
+      auto d_A = A_pc->ReadWriteData(on_dev);
+      auto d_ess_list = ess_pc_tdof_list.Read(on_dev);
+      MFEM_FORALL_SWITCH(on_dev, i, ess_pc_tdof_list.Size(),
+      {
+         int rc = d_ess_list[i];
+         if (rc < 0 ) { rc = -1-rc; }
+         for (int j = d_I[rc]; j < d_I[rc+1]; j++)
+         {
+            const int col = d_J[j];
+            if (col != rc)
+            {
+               d_A[j] = 0.0;
+               for (int k = d_I[col]; 1; k++)
+               {
+                  if (k == d_I[col+1])
+                  {
+                     break;
+                  }
+                  else if (d_J[k] == rc)
+                  {
+                     d_A[k] = 0.0;
+                     break;
+                  }
+               }
+            }
+         }
+      });
+
+      tic_toc.Stop();
+      cout << "Real time adjusting A_pc for essential BC: " <<
+           tic_toc.RealTime() << "\n";
+
+      if (pc_choice == GKO_BLOCK_JACOBI)
+      {
+
+         // Create Ginkgo Jacobi preconditioner
+
+
+         if (permute)
+         {
+            tic_toc.Clear();
+            tic_toc.Start();
+            GinkgoWrappers::GinkgoJacobiPreconditioner M(executor, *A_pc, *inv_reordering,
+                                                         pc_storage_opt,
+                                                         pc_acc, pc_max_bs);
+            tic_toc.Stop();
+            cout << "Real time creating Ginkgo BlockJacobi preconditioner: " <<
+                 tic_toc.RealTime() << endl;
+
+            // Use preconditioned CG
+            total_its = pcg_solve(*A, M, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+            cout << "Real time in PCG: " << it_time << endl;
+         }
+         else
+         {
+            tic_toc.Clear();
+            tic_toc.Start();
+            GinkgoWrappers::GinkgoJacobiPreconditioner M(executor, *A_pc, pc_storage_opt,
+                                                         pc_acc, pc_max_bs);
+            tic_toc.Stop();
+            cout << "Real time creating Ginkgo BlockJacobi preconditioner: " <<
+                 tic_toc.RealTime() << endl;
+
+            // Use preconditioned CG
+            total_its = pcg_solve(*A, M, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+            cout << "Real time in PCG: " << it_time << endl;
+         }
+
+      }
+      else if (pc_choice == GKO_ILU || pc_choice == GKO_ILU_ISAI)
+      {
+
+         // Create Ginkgo ILU preconditioner
+
+         if (permute)
+         {
+
+            tic_toc.Clear();
+            tic_toc.Start();
+
+            GinkgoWrappers::GinkgoIluPreconditioner M(executor, *A_pc, *inv_reordering,
+                                                      trisolve_type, isai_sparsity_power, par_ilu_its);
+
+            tic_toc.Stop();
+            cout << "Real time creating Ginkgo Ilu preconditioner: " <<
+                 tic_toc.RealTime() << endl;
+
+            // Use preconditioned CG
+            total_its = pcg_solve(*A, M, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+            cout << "Real time in PCG: " << it_time << endl;
+
+         }
+         else
+         {
+
+            tic_toc.Clear();
+            tic_toc.Start();
+
+            GinkgoWrappers::GinkgoIluPreconditioner M(executor, *A_pc, trisolve_type,
+                                                      isai_sparsity_power, par_ilu_its);
+
+            tic_toc.Stop();
+            cout << "Real time creating Ginkgo Ilu preconditioner: " <<
+                 tic_toc.RealTime() << endl;
+
+            // Use preconditioned CG
+            total_its = pcg_solve(*A, M, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+            cout << "Real time in PCG: " << it_time << endl;
+
+         }
+      }
+      else if (pc_choice == GKO_CUILU || pc_choice == GKO_CUILU_ISAI)
+      {
+
+         // Create Ginkgo ILU preconditioner using ILU from cuSPARSE
+
+         if (permute)
+         {
+
+            tic_toc.Clear();
+            tic_toc.Start();
+
+            GinkgoWrappers::GinkgoCuIluPreconditioner M(executor, *A_pc, *inv_reordering,
+                                                        trisolve_type, isai_sparsity_power);
+
+            tic_toc.Stop();
+            cout << "Real time creating Ginkgo CuIlu preconditioner: " <<
+                 tic_toc.RealTime() << endl;
+
+            // Use preconditioned CG
+            total_its = pcg_solve(*A, M, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+            cout << "Real time in PCG: " << it_time << endl;
+
+         }
+         else
+         {
+
+            tic_toc.Clear();
+            tic_toc.Start();
+
+            GinkgoWrappers::GinkgoCuIluPreconditioner M(executor, *A_pc, trisolve_type,
+                                                      isai_sparsity_power);
+
+            tic_toc.Stop();
+            cout << "Real time creating Ginkgo CuIlu preconditioner: " <<
+                 tic_toc.RealTime() << endl;
+
+            // Use preconditioned CG
+            total_its = pcg_solve(*A, M, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+            cout << "Real time in PCG: " << it_time << endl;
+
+         }
+      }
+      else if (pc_choice == MFEM_GS)
+      {
+
+         // Create MFEM preconditioner
+         tic_toc.Clear();
+         tic_toc.Start();
+
+         GSSmoother M(*A_pc);
+
+         tic_toc.Stop();
+         cout << "Real time creating MFEM GS preconditioner: " <<
+              tic_toc.RealTime() << endl;
+
+         // Use preconditioned CG
+         total_its = pcg_solve(*A, M, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+         cout << "Real time in PCG: " << it_time << endl;
+
+      }
+      else if (pc_choice == MFEM_UMFPACK)
+      {
+
+#ifdef MFEM_USE_SUITESPARSE
+         // Create MFEM preconditioner
+         tic_toc.Clear();
+         tic_toc.Start();
+
+         UMFPackSolver M;
+         M.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
+         M.SetOperator(*A_pc);
+
+         tic_toc.Stop();
+         cout << "Real time creating MFEM UMFPACK preconditioner: " <<
+              tic_toc.RealTime() << endl;
+
+         // Use preconditioned CG
+         total_its = pcg_solve(*A, M, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+         cout << "Real time in PCG: " << it_time << endl;
+#endif
+      }
+   }
+   else
+   {
+
+      total_its = cg_solve(*A, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+      cout << "Real time in CG: " << it_time << endl;
+   }
+
+   cout << "Total iterations: " << total_its << endl;
+   cout << "Avg time per iteration: " << it_time/double(
+           total_its) << endl;
+
+   // 12. Recover the solution as a finite element grid function.
+   a->RecoverFEMSolution(X, *b, x);
+
+   // 13. Save the refined mesh and the solution. This output can be viewed
+   // later
+   //     using GLVis: "glvis -m refined.mesh -g sol.gf".
+
+   if (output_sol)
+   {
+      ofstream mesh_ofs("refined.mesh");
+      mesh_ofs.precision(8);
+      mesh->Print(mesh_ofs);
+      ofstream sol_ofs("sol.gf");
+      sol_ofs.precision(8);
+      x.Save(sol_ofs);
+   }
+
+   if (pc && output_pc)
+   {
+      ofstream mesh_lor_ofs("lor-refined.mesh");
+      mesh_lor_ofs.precision(8);
+      mesh_lor->Print(mesh_lor_ofs);
+
+      ofstream apc_lor_ofs("lor-mat.dat");
+      mesh_lor_ofs.precision(8);
+      A_pc->PrintCSR(apc_lor_ofs);
+   }
+
+   // 14. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int visport = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *mesh << x << flush;
+   }
+
+   // 15. Free the used memory.
+   delete a;
+   delete b;
+   delete fespace;
+   delete fespace_lor;
+   delete fec_lor;
+   delete mesh_lor;
+   delete coeff;
+   if (order > 0)
+   {
+      delete fec;
+   }
+   delete mesh;
+}
+
+double coeff_func(const Vector & x)
+{
+   double xi(x(0));
+   double yi(x(1));
+   double zi(0.0);
+
+   if (x.Size() == 3)
+   {
+      zi = x(2);
+   }
+
+   return 0.2*(sin(2*xi)*sin(2*yi)*cos(2*zi) + 2.0);
+}

fem/bilinearform.hpp

@@ -299,8 +299,17 @@ public:
    /// Returns a reference to the sparse matrix:  \f$ M \f$
    SparseMatrix &SpMat()
    {
-      MFEM_VERIFY(mat, "mat is NULL and can't be dereferenced");
-      return *mat;
+      FABilinearFormExtension *fa_ext = dynamic_cast<FABilinearFormExtension*>
+                                        (this->ext);
+      if (fa_ext)
+      {
+         return fa_ext->SpMat();
+      }
+      else
+      {
+         MFEM_VERIFY(mat, "mat is NULL and can't be dereferenced");
+         return *mat;
+      }
    }
 
    /**  @brief Nullifies the internal matrix \f$ M \f$ and returns a pointer

fem/bilinearform_ext.hpp

@@ -128,6 +128,10 @@ public:
    void Assemble();
    void Mult(const Vector &x, Vector &y) const;
    void MultTranspose(const Vector &x, Vector &y) const;
+   SparseMatrix &SpMat()
+   {
+      return mat;
+   }
 };
 
 /// Data and methods for matrix-free bilinear forms NOT YET IMPLEMENTED.

fem/quadinterpolator.hpp

@@ -49,7 +49,7 @@ protected:
    static const int MAX_ND2D = 100;
    static const int MAX_VDIM2D = 3;
 
-   static const int MAX_NQ3D = 1000;
+   static const int MAX_NQ3D = 1728;
    static const int MAX_ND3D = 1000;
    static const int MAX_VDIM3D = 3;
 

linalg/ginkgo.cpp

@@ -276,8 +276,8 @@ CGSolver::CGSolver(
                                ATOLERANCE)
 {
    using cg = gko::solver::Cg<double>;
-   this->solver_gen =
-      cg::build().with_criteria(this->combined_factory).on(this->executor);
+   //   this->solver_gen =
+   //      cg::build().with_criteria(this->combined_factory).on(this->executor);
 }
 
 CGSolver::CGSolver(
@@ -292,10 +292,10 @@ CGSolver::CGSolver(
                                ATOLERANCE)
 {
    using cg         = gko::solver::Cg<double>;
-   this->solver_gen = cg::build()
-                      .with_criteria(this->combined_factory)
-                      .with_preconditioner(preconditioner)
-                      .on(this->executor);
+   //   this->solver_gen = cg::build()
+   //                      .with_criteria(this->combined_factory)
+   //                      .with_preconditioner(preconditioner)
+   //                      .on(this->executor);
 }
 
 
@@ -311,9 +311,9 @@ BICGSTABSolver::BICGSTABSolver(
                                ATOLERANCE)
 {
    using bicgstab   = gko::solver::Bicgstab<double>;
-   this->solver_gen = bicgstab::build()
-                      .with_criteria(this->combined_factory)
-                      .on(this->executor);
+   //   this->solver_gen = bicgstab::build()
+   //                      .with_criteria(this->combined_factory)
+   //                      .on(this->executor);
 }
 
 BICGSTABSolver::BICGSTABSolver(
@@ -328,10 +328,10 @@ BICGSTABSolver::BICGSTABSolver(
                                ATOLERANCE)
 {
    using bicgstab   = gko::solver::Bicgstab<double>;
-   this->solver_gen = bicgstab::build()
-                      .with_criteria(this->combined_factory)
-                      .with_preconditioner(preconditioner)
-                      .on(this->executor);
+   //   this->solver_gen = bicgstab::build()
+   //                      .with_criteria(this->combined_factory)
+   //                      .with_preconditioner(preconditioner)
+   //                      .on(this->executor);
 }
 
 
@@ -347,8 +347,8 @@ CGSSolver::CGSSolver(
                                ATOLERANCE)
 {
    using cgs = gko::solver::Cgs<double>;
-   this->solver_gen =
-      cgs::build().with_criteria(this->combined_factory).on(this->executor);
+   //   this->solver_gen =
+   //      cgs::build().with_criteria(this->combined_factory).on(this->executor);
 }
 
 CGSSolver::CGSSolver(
@@ -363,10 +363,10 @@ CGSSolver::CGSSolver(
                                ATOLERANCE)
 {
    using cgs        = gko::solver::Cgs<double>;
-   this->solver_gen = cgs::build()
-                      .with_criteria(this->combined_factory)
-                      .with_preconditioner(preconditioner)
-                      .on(this->executor);
+   //   this->solver_gen = cgs::build()
+   //                      .with_criteria(this->combined_factory)
+   //                      .with_preconditioner(preconditioner)
+   //                      .on(this->executor);
 }
 
 
@@ -382,8 +382,8 @@ FCGSolver::FCGSolver(
                                ATOLERANCE)
 {
    using fcg = gko::solver::Fcg<double>;
-   this->solver_gen =
-      fcg::build().with_criteria(this->combined_factory).on(this->executor);
+   //   this->solver_gen =
+   //      fcg::build().with_criteria(this->combined_factory).on(this->executor);
 }
 
 FCGSolver::FCGSolver(
@@ -398,10 +398,10 @@ FCGSolver::FCGSolver(
                                ATOLERANCE)
 {
    using fcg        = gko::solver::Fcg<double>;
-   this->solver_gen = fcg::build()
-                      .with_criteria(this->combined_factory)
-                      .with_preconditioner(preconditioner)
-                      .on(this->executor);
+   //   this->solver_gen = fcg::build()
+   //                      .with_criteria(this->combined_factory)
+   //                      .with_preconditioner(preconditioner)
+   //                      .on(this->executor);
 }
 
 
@@ -417,10 +417,10 @@ GMRESSolver::GMRESSolver(
                                ATOLERANCE)
 {
    using gmres      = gko::solver::Gmres<double>;
-   this->solver_gen = gmres::build()
-                      .with_krylov_dim(m)
-                      .with_criteria(this->combined_factory)
-                      .on(this->executor);
+   //   this->solver_gen = gmres::build()
+   //                      .with_krylov_dim(m)
+   //                      .with_criteria(this->combined_factory)
+   //                      .on(this->executor);
 }
 
 GMRESSolver::GMRESSolver(
@@ -435,11 +435,11 @@ GMRESSolver::GMRESSolver(
                                ATOLERANCE)
 {
    using gmres      = gko::solver::Gmres<double>;
-   this->solver_gen = gmres::build()
-                      .with_krylov_dim(m)
-                      .with_criteria(this->combined_factory)
-                      .with_preconditioner(preconditioner)
-                      .on(this->executor);
+   //   this->solver_gen = gmres::build()
+   //                      .with_krylov_dim(m)
+   //                      .with_criteria(this->combined_factory)
+   //                      .with_preconditioner(preconditioner)
+   //                      .on(this->executor);
 }
 
 
@@ -455,8 +455,8 @@ IRSolver::IRSolver(
                                ATOLERANCE)
 {
    using ir = gko::solver::Ir<double>;
-   this->solver_gen =
-      ir::build().with_criteria(this->combined_factory).on(this->executor);
+   //   this->solver_gen =
+   //      ir::build().with_criteria(this->combined_factory).on(this->executor);
 }
 
 IRSolver::IRSolver(
@@ -471,10 +471,99 @@ IRSolver::IRSolver(
                                ATOLERANCE)
 {
    using ir         = gko::solver::Ir<double>;
-   this->solver_gen = ir::build()
-                      .with_criteria(this->combined_factory)
-                      .with_solver(inner_solver)
-                      .on(this->executor);
+   //   this->solver_gen = ir::build()
+   //                      .with_criteria(this->combined_factory)
+   //                      .with_solver(inner_solver)
+   //                      .on(this->executor);
+}
+
+
+/* ---------------------- GinkgoPreconditioner ------------------------ */
+void GinkgoPreconditionerBase::Mult(const Vector &x, Vector &y) const
+{
+
+   if (!iterative_mode)
+   {
+      y = 0.0;
+   }
+
+   //Create Ginkgo wrapped-vectors
+   bool on_device = false;
+   if (exec_ != exec_->get_master())
+   {
+      on_device = true;
+   }
+   using vec = gko::matrix::Dense<double>;
+   auto gko_x = vec::create(exec_, gko::dim<2> {x.Size(), 1},
+                            gko::Array<double>::view(exec_,
+                                                     x.Size(), const_cast<double *>(
+                                                        x.Read(on_device))), 1);
+   auto gko_y = vec::create(exec_, gko::dim<2> {y.Size(), 1},
+                            gko::Array<double>::view(exec_,
+                                                     y.Size(), y.ReadWrite(on_device)), 1);
+
+   if (permute_)
+   {
+
+      auto gko_x_perm = vec::create(exec_, gko::dim<2> {x.Size(), 1});
+      auto gko_y_perm = vec::create(exec_, gko::dim<2> {y.Size(), 1});
+
+      vec_permute_->apply(gko::lend(gko_x), gko::lend(gko_x_perm));
+      vec_permute_->apply(gko::lend(gko_y), gko::lend(gko_y_perm));
+
+      gko_precond_.get()->apply(gko::lend(gko_x_perm), gko::lend(gko_y_perm));
+
+      vec_inv_permute_->apply(gko::lend(gko_y_perm), gko::lend(gko_y));
+
+      // Reset MFEM Vector to use Ginkgo output ?
+      //   y = gko_y->get_values();
+
+   }
+   else
+   {
+
+      gko_precond_.get()->apply(gko::lend(gko_x), gko::lend(gko_y));
+
+   }
+}
+
+void GinkgoPreconditionerBase::SetOperator(const Operator &op)
+{
+
+   // Only accept SparseMatrix for this type (see SparseSmoother::SetOperator)
+   SparseMatrix *op_mat = const_cast<SparseMatrix*>(
+                             dynamic_cast<const SparseMatrix*>(&op));
+   if (op_mat == NULL)
+   {
+      mfem_error("GinkgoPreconditioner::SetOperator : not a SparseMatrix!");
+   }
+   height = op_mat->Height();
+   width = op_mat->Width();
+
+   // Release current preconditioner
+   gko_precond_.release();
+
+   bool on_device = false;
+   if (exec_ != exec_->get_master())
+   {
+      on_device = true;
+   }
+
+   using mtx = gko::matrix::Csr<double, int>;
+   auto gko_sparse = mtx::create(
+                        exec_, gko::dim<2>(op_mat->Height(), op_mat->Width()),
+                        gko::Array<double>::view(exec_,
+                                                 op_mat->NumNonZeroElems(),
+                                                 op_mat->ReadWriteData(on_device)),
+                        gko::Array<int>::view(exec_,
+                                              op_mat->NumNonZeroElems(),
+                                              op_mat->ReadWriteJ(on_device)),
+                        gko::Array<int>::view(exec_,
+                                              op_mat->Height() + 1,
+                                              op_mat->ReadWriteI(on_device)));
+   gko_precond_ = gko_precond_factory_.get()->generate(
+                     gko::give(gko_sparse));
+
 }
 
 

linalg/ginkgo.hpp

@@ -640,6 +640,465 @@ public:
 };
 
 
+class GinkgoPreconditionerBase : public Solver
+{
+protected:
+   GinkgoPreconditionerBase(std::shared_ptr<const gko::Executor> exec,
+                            SparseMatrix &a,
+                            bool iter_mode=false)
+      : Solver(a.Height(), a.Width(), iter_mode)
+   {
+      exec_ = std::move(exec);
+      permute_ = false;
+   }
+
+   GinkgoPreconditionerBase(std::shared_ptr<const gko::Executor> exec,
+                            SparseMatrix &a, Array<int> &inv_permutation_indices,
+                            bool iter_mode=false)
+      : Solver(a.Height(), a.Width(), iter_mode)
+   {
+      exec_ = std::move(exec);
+
+      bool on_device = false;
+      if (exec_->get_master() != exec_)
+      {
+         on_device = true;
+      }
+      permute_ = true;
+      auto gko_inv_perm_ind = gko::Array<int>::view(
+                                 exec_,
+                                 inv_permutation_indices.Size(),
+                                 inv_permutation_indices.ReadWrite(
+                                    on_device));
+      // Note the "forward" permutation uses the inverse flag because
+      //  the indices are for the inverse permutation as defined by Ginkgo
+      vec_permute_ = gko::matrix::Permutation<int>::create(
+                        exec_, gko::dim<2> {inv_permutation_indices.Size()},
+                        gko_inv_perm_ind,
+                        gko::matrix::row_permute | gko::matrix::inverse_permute);
+      vec_inv_permute_ = gko::matrix::Permutation<int>::create(
+                            exec_, gko::dim<2> {inv_permutation_indices.Size()},
+                            gko_inv_perm_ind,
+                            gko::matrix::row_permute);
+   }
+
+public:
+   std::shared_ptr<const gko::Executor> get_exec() {return this->exec_; }
+   const gko::LinOpFactory* get_gko_precond_factory()
+   {
+      return this->gko_precond_factory_.get();
+   }
+   gko::LinOp* get_gko_precond() {return this->gko_precond_.get(); }
+
+   virtual void Mult(const Vector &x, Vector &y) const;
+   virtual void SetOperator(const Operator &op);
+
+protected:
+   std::shared_ptr<const gko::Executor> exec_;
+   std::unique_ptr<const gko::LinOpFactory> gko_precond_factory_;
+   std::unique_ptr<gko::LinOp> gko_precond_;
+   bool permute_;
+   std::unique_ptr<gko::matrix::Permutation<int>> vec_permute_;
+   std::unique_ptr<gko::matrix::Permutation<int>> vec_inv_permute_;
+
+};
+
+
+class GinkgoJacobiPreconditioner : public GinkgoPreconditionerBase
+{
+public:
+   GinkgoJacobiPreconditioner(std::shared_ptr<const gko::Executor> exec,
+                              SparseMatrix &a,
+                              const char *storage_opt="none",
+                              const double accuracy=1.e-1,
+                              const int max_block_size=32,
+                              bool iter_mode=false)
+      : GinkgoPreconditionerBase(exec, a, iter_mode)
+   {
+
+      bool on_device = false;
+      if (exec->get_master() != exec)
+      {
+         on_device = true;
+      }
+
+      using mtx = gko::matrix::Csr<double, int>;
+      const int nnz =  a.GetMemoryData().Capacity();
+      auto gko_sparse = mtx::create(
+                           exec, gko::dim<2>(a.Height(), a.Width()),
+                           gko::Array<double>::view(exec,
+                                                    nnz,
+                                                    a.ReadWriteData(on_device)),
+                           gko::Array<int>::view(exec,
+                                                 nnz,
+                                                 a.ReadWriteJ(on_device)),
+                           gko::Array<int>::view(exec, a.Height() + 1,
+                                                 a.ReadWriteI(on_device)));
+
+      gko_sparse->sort_by_column_index();
+
+      if (storage_opt == "auto")
+      {
+         gko_precond_factory_ = gko::preconditioner::Jacobi<double, int>::build()
+                                .with_storage_optimization(
+                                   gko::precision_reduction::autodetect())
+                                .with_accuracy(accuracy)
+                                .with_max_block_size(max_block_size)
+                                .on(exec);
+      }
+      else
+      {
+         gko_precond_factory_ = gko::preconditioner::Jacobi<double, int>::build()
+                                .with_storage_optimization(
+                                   gko::precision_reduction(0, 0))
+                                .with_accuracy(accuracy)
+                                .with_max_block_size(max_block_size)
+                                .on(exec);
+      }
+
+      gko_precond_ = gko_precond_factory_.get()->generate(
+                        gko::give(gko_sparse));
+   }
+
+   GinkgoJacobiPreconditioner(std::shared_ptr<const gko::Executor> exec,
+                              SparseMatrix &a,
+                              Array<int> &inv_permutation_indices,
+                              const char *storage_opt="none",
+                              const double accuracy=1.e-1,
+                              const int max_block_size=32,
+                              bool iter_mode=false)
+      : GinkgoPreconditionerBase(exec, a, inv_permutation_indices, iter_mode)
+   {
+
+      bool on_device = false;
+      if (exec->get_master() != exec)
+      {
+         on_device = true;
+      }
+
+      using mtx = gko::matrix::Csr<double, int>;
+      const int nnz =  a.GetMemoryData().Capacity();
+      auto gko_sparse = mtx::create(
+                           exec, gko::dim<2>(a.Height(), a.Width()),
+                           gko::Array<double>::view(exec,
+                                                    nnz,
+                                                    a.ReadWriteData(on_device)),
+                           gko::Array<int>::view(exec,
+                                                 nnz,
+                                                 a.ReadWriteJ(on_device)),
+                           gko::Array<int>::view(exec, a.Height() + 1,
+                                                 a.ReadWriteI(on_device)));
+
+
+      gko_sparse->sort_by_column_index();
+
+      if (storage_opt == "auto")
+      {
+         gko_precond_factory_ = gko::preconditioner::Jacobi<double, int>::build()
+                                .with_storage_optimization(
+                                   gko::precision_reduction::autodetect())
+                                .with_accuracy(accuracy)
+                                .with_max_block_size(max_block_size)
+                                .on(exec);
+      }
+      else
+      {
+         gko_precond_factory_ = gko::preconditioner::Jacobi<double, int>::build()
+                                .with_storage_optimization(
+                                   gko::precision_reduction(0, 0))
+                                .with_accuracy(accuracy)
+                                .with_max_block_size(max_block_size)
+                                .on(exec);
+      }
+
+      gko_precond_ = gko_precond_factory_.get()->generate(
+                        gko::give(gko_sparse));
+   }
+};
+
+class GinkgoIluPreconditioner : public GinkgoPreconditionerBase
+{
+public:
+   GinkgoIluPreconditioner(std::shared_ptr<const gko::Executor> exec,
+                           SparseMatrix &a, const char *trisolve_type = "exact",
+                           int sparsity_power=1,
+                           int par_ilu_its=0,
+                           bool iter_mode=false)
+      : GinkgoPreconditionerBase(exec, a, iter_mode)
+   {
+
+      bool on_device = false;
+      if (exec->get_master() != exec)
+      {
+         on_device = true;
+      }
+
+      using mtx = gko::matrix::Csr<double, int>;
+      const int nnz =  a.GetMemoryData().Capacity();
+      auto gko_sparse = mtx::create(
+                           exec, gko::dim<2>(a.Height(), a.Width()),
+                           gko::Array<double>::view(exec,
+                                                    nnz,
+                                                    a.ReadWriteData(on_device)),
+                           gko::Array<int>::view(exec,
+                                                 nnz,
+                                                 a.ReadWriteJ(on_device)),
+                           gko::Array<int>::view(exec, a.Height() + 1,
+                                                 a.ReadWriteI(on_device)));
+
+      using ilu_fact_type = gko::factorization::ParIlu<double, int>;
+      std::shared_ptr<ilu_fact_type::Factory> fact_factory = std::move(
+                                                                ilu_fact_type::build()
+                                                                .with_iterations(par_ilu_its)
+                                                                .on(exec));
+
+
+      if (trisolve_type == "isai")
+      {
+
+         using l_solver_type = gko::preconditioner::LowerIsai<>;
+         using u_solver_type = gko::preconditioner::UpperIsai<>;
+
+         std::shared_ptr<l_solver_type::Factory> l_solver_factory = std::move(
+                                                                       l_solver_type::build()
+                                                                       .with_sparsity_power(sparsity_power)
+                                                                       .on(exec));
+         std::shared_ptr<u_solver_type::Factory> u_solver_factory = std::move(
+                                                                       u_solver_type::build()
+                                                                       .with_sparsity_power(sparsity_power)
+                                                                       .on(exec));
+
+         gko_precond_factory_ = gko::preconditioner::Ilu<l_solver_type,
+         u_solver_type>::build()
+         .with_factorization_factory(fact_factory)
+         .with_l_solver_factory(l_solver_factory)
+         .with_u_solver_factory(u_solver_factory)
+         .on(exec);
+
+      }
+      else
+      {
+
+         gko_precond_factory_ = gko::preconditioner::Ilu<>::build()
+                                .with_factorization_factory(fact_factory)
+                                .on(exec);
+      }
+
+      gko_precond_ = gko_precond_factory_.get()->generate(
+                        gko::give(gko_sparse));
+   }
+
+   GinkgoIluPreconditioner(std::shared_ptr<const gko::Executor> exec,
+                           SparseMatrix &a, Array<int> &inv_permutation_indices,
+                           const char *trisolve_type = "exact",
+                           int sparsity_power=1,
+                           int par_ilu_its=0,
+                           bool iter_mode=false)
+      : GinkgoPreconditionerBase(exec, a, inv_permutation_indices, iter_mode)
+   {
+
+      bool on_device = false;
+      if (exec->get_master() != exec)
+      {
+         on_device = true;
+      }
+
+      using mtx = gko::matrix::Csr<double, int>;
+      const int nnz =  a.GetMemoryData().Capacity();
+      auto gko_sparse = mtx::create(
+                           exec, gko::dim<2>(a.Height(), a.Width()),
+                           gko::Array<double>::view(exec,
+                                                    nnz,
+                                                    a.ReadWriteData(on_device)),
+                           gko::Array<int>::view(exec,
+                                                 nnz,
+                                                 a.ReadWriteJ(on_device)),
+                           gko::Array<int>::view(exec, a.Height() + 1,
+                                                 a.ReadWriteI(on_device)));
+
+      using ilu_fact_type = gko::factorization::ParIlu<double, int>;
+      std::shared_ptr<ilu_fact_type::Factory> fact_factory = std::move(
+                                                                ilu_fact_type::build()
+                                                                .with_iterations(par_ilu_its)
+                                                                .on(exec));
+
+      if (trisolve_type == "isai")
+      {
+
+         using l_solver_type = gko::preconditioner::LowerIsai<>;
+         using u_solver_type = gko::preconditioner::UpperIsai<>;
+
+         std::shared_ptr<l_solver_type::Factory> l_solver_factory = std::move(
+                                                                       l_solver_type::build()
+                                                                       .with_sparsity_power(sparsity_power)
+                                                                       .on(exec));
+         std::shared_ptr<u_solver_type::Factory> u_solver_factory = std::move(
+                                                                       u_solver_type::build()
+                                                                       .with_sparsity_power(sparsity_power)
+                                                                       .on(exec));
+
+         gko_precond_factory_ = gko::preconditioner::Ilu<l_solver_type,
+         u_solver_type>::build()
+         .with_factorization_factory(fact_factory)
+         .with_l_solver_factory(l_solver_factory)
+         .with_u_solver_factory(u_solver_factory)
+         .on(exec);
+      }
+      else
+      {
+
+         gko_precond_factory_ = gko::preconditioner::Ilu<>::build()
+                                .with_factorization_factory(fact_factory)
+                                .on(exec);
+      }
+
+      gko_precond_ = gko_precond_factory_.get()->generate(
+                        gko::give(gko_sparse));
+   }
+};
+
+class GinkgoCuIluPreconditioner : public GinkgoPreconditionerBase
+{
+public:
+   GinkgoCuIluPreconditioner(std::shared_ptr<const gko::Executor> exec,
+                             SparseMatrix &a, const char *trisolve_type = "exact",
+                             int sparsity_power=1,
+                             bool iter_mode=false)
+      : GinkgoPreconditionerBase(exec, a, iter_mode)
+   {
+
+      bool on_device = false;
+      if (exec->get_master() != exec)
+      {
+         on_device = true;
+      }
+
+      using mtx = gko::matrix::Csr<double, int>;
+      const int nnz =  a.GetMemoryData().Capacity();
+      auto gko_sparse = mtx::create(
+                           exec, gko::dim<2>(a.Height(), a.Width()),
+                           gko::Array<double>::view(exec,
+                                                    nnz,
+                                                    a.ReadWriteData(on_device)),
+                           gko::Array<int>::view(exec,
+                                                 nnz,
+                                                 a.ReadWriteJ(on_device)),
+                           gko::Array<int>::view(exec, a.Height() + 1,
+                                                 a.ReadWriteI(on_device)));
+
+      // TEST
+      gko_sparse->sort_by_column_index();
+
+      using ilu_fact_type = gko::factorization::Ilu<double, int>;
+      std::shared_ptr<ilu_fact_type::Factory> fact_factory = std::move(
+                                                                ilu_fact_type::build()
+                                                                .on(exec));
+
+
+      if (trisolve_type == "isai")
+      {
+
+         using l_solver_type = gko::preconditioner::LowerIsai<>;
+         using u_solver_type = gko::preconditioner::UpperIsai<>;
+
+         std::shared_ptr<l_solver_type::Factory> l_solver_factory = std::move(
+                                                                       l_solver_type::build()
+                                                                       .with_sparsity_power(sparsity_power)
+                                                                       .on(exec));
+         std::shared_ptr<u_solver_type::Factory> u_solver_factory = std::move(
+                                                                       u_solver_type::build()
+                                                                       .with_sparsity_power(sparsity_power)
+                                                                       .on(exec));
+
+         gko_precond_factory_ = gko::preconditioner::Ilu<l_solver_type,
+         u_solver_type>::build()
+         .with_factorization_factory(fact_factory)
+         .with_l_solver_factory(l_solver_factory)
+         .with_u_solver_factory(u_solver_factory)
+         .on(exec);
+
+      }
+      else
+      {
+
+         gko_precond_factory_ = gko::preconditioner::Ilu<>::build()
+                                .with_factorization_factory(fact_factory)
+                                .on(exec);
+      }
+
+      gko_precond_ = gko_precond_factory_.get()->generate(
+                        gko::give(gko_sparse));
+   }
+
+   GinkgoCuIluPreconditioner(std::shared_ptr<const gko::Executor> exec,
+                             SparseMatrix &a, Array<int> &inv_permutation_indices,
+                             const char *trisolve_type = "exact",
+                             int sparsity_power=1,
+                             bool iter_mode=false)
+      : GinkgoPreconditionerBase(exec, a, inv_permutation_indices, iter_mode)
+   {
+
+      bool on_device = false;
+      if (exec->get_master() != exec)
+      {
+         on_device = true;
+      }
+
+      using mtx = gko::matrix::Csr<double, int>;
+      const int nnz =  a.GetMemoryData().Capacity();
+      auto gko_sparse = mtx::create(
+                           exec, gko::dim<2>(a.Height(), a.Width()),
+                           gko::Array<double>::view(exec,
+                                                    nnz,
+                                                    a.ReadWriteData(on_device)),
+                           gko::Array<int>::view(exec,
+                                                 nnz,
+                                                 a.ReadWriteJ(on_device)),
+                           gko::Array<int>::view(exec, a.Height() + 1,
+                                                 a.ReadWriteI(on_device)));
+
+      // TEST
+      gko_sparse->sort_by_column_index();
+
+      using ilu_fact_type = gko::factorization::Ilu<double, int>;
+      std::shared_ptr<ilu_fact_type::Factory> fact_factory = std::move(
+                                                                ilu_fact_type::build()
+                                                                .on(exec));
+      if (trisolve_type == "isai")
+      {
+
+         using l_solver_type = gko::preconditioner::LowerIsai<>;
+         using u_solver_type = gko::preconditioner::UpperIsai<>;
+
+         std::shared_ptr<l_solver_type::Factory> l_solver_factory = std::move(
+                                                                       l_solver_type::build()
+                                                                       .with_sparsity_power(sparsity_power)
+                                                                       .on(exec));
+         std::shared_ptr<u_solver_type::Factory> u_solver_factory = std::move(
+                                                                       u_solver_type::build()
+                                                                       .with_sparsity_power(sparsity_power)
+                                                                       .on(exec));
+
+         gko_precond_factory_ = gko::preconditioner::Ilu<l_solver_type,
+         u_solver_type>::build()
+         .with_factorization_factory(fact_factory)
+         .with_l_solver_factory(l_solver_factory)
+         .with_u_solver_factory(u_solver_factory)
+         .on(exec);
+      }
+      else
+      {
+
+         gko_precond_factory_ = gko::preconditioner::Ilu<>::build()
+                                .with_factorization_factory(fact_factory)
+                                .on(exec);
+      }
+
+      gko_precond_ = gko_precond_factory_.get()->generate(
+                        gko::give(gko_sparse));
+   }
+};
+
 } // namespace GinkgoWrappers
 
 }

linalg/solvers.hpp

@@ -28,6 +28,8 @@
 namespace mfem
 {
 
+void MinimumDiscardedFillOrdering(SparseMatrix &C, Array<int> &p);
+
 class BilinearForm;
 
 /// Abstract base class for an iterative solver monitor