cleaning up docs and sharp edges

src/asgard_discretization.cpp

@@ -155,16 +155,16 @@ void discretization_manager<precision>::start_cold()
 
   degree_ = options.degree.value();
 
-  if (not stop_verbosity())
-    std::cout << "\n -- ASGarD discretization options --\n";
-
   if (high_verbosity()) {
     std::cout << "Branch: " << GIT_BRANCH << '\n';
     std::cout << "Commit Summary: " << GIT_COMMIT_HASH
                     << GIT_COMMIT_SUMMARY << '\n';
-    std::cout << "This executable was built on " << BUILD_TIME << '\n';
+    std::cout << "The library was built on " << BUILD_TIME << '\n';
   }
 
+  if (not stop_verbosity())
+    std::cout << "\n -- ASGarD discretization options --\n";
+
   sgrid = sparse_grid(options);
 
   if (not stop_verbosity()) {
@@ -572,7 +572,7 @@ void discretization_manager<precision>::ode_sv(imex_flag imflag,
       kronops.make(imflag, *pde, matrices, grid);
       precision const tolerance = *options.isolver_tolerance;
       int const restart         = *options.isolver_iterations;
-      int const max_iter        = *options.isolver_outer_iterations;
+      int const max_iter        = *options.isolver_inner_iterations;
       sol.resize(static_cast<int>(x.size()));
       std::copy(x.begin(), x.end(), sol.begin());
       if (solver == solve_opts::gmres)

src/asgard_discretization.hpp

@@ -175,8 +175,7 @@ class discretization_manager
     terms.apply_all(sgrid, conn, alpha, x, beta, y);
   }
   //! applies ADI preconditioner for all terms
-  void terms_apply_adi(std::vector<precision> const &x,
-                       std::vector<precision> &y) const
+  void terms_apply_adi(precision const x[], precision y[]) const
   {
     tools::time_event performance_("terms_apply_adi kronmult");
     terms.apply_all_adi(sgrid, conn, x, y);

src/asgard_io.cpp

@@ -141,7 +141,7 @@ void write_output(PDE<P> const &pde, // std::vector<moment<P>> const &moments,
 
   H5Easy::dump(file, "isolver_tolerance", options.isolver_tolerance.value());
   H5Easy::dump(file, "isolver_iterations", options.isolver_iterations.value());
-  H5Easy::dump(file, "isolver_outer_iterations", options.isolver_outer_iterations.value());
+  H5Easy::dump(file, "isolver_inner_iterations", options.isolver_inner_iterations.value());
 
   // save some basic build info
   H5Easy::dump(file, "GIT_BRANCH", std::string(GIT_BRANCH));
@@ -281,7 +281,7 @@ void h5writer<P>::write(PDEv2<P> const &pde, int degree, sparse_grid const &grid
     H5Easy::dump(file, "solver_method", static_cast<int>(options.solver.value_or(solve_opts::direct)));
     H5Easy::dump(file, "solver_itol", options.isolver_tolerance.value_or(-1));
     H5Easy::dump(file, "solver_iter", options.isolver_iterations.value_or(-1));
-    H5Easy::dump(file, "solver_oiter", options.isolver_outer_iterations.value_or(-1));
+    H5Easy::dump(file, "solver_inner", options.isolver_inner_iterations.value_or(-1));
   }
 
   H5Easy::dump(file, "timer_report", tools::timer.report());
@@ -498,12 +498,21 @@ void h5writer<P>::read(std::string const &filename, bool silent, PDEv2<P> &pde,
   { // solver data section
     if (not pde.options_.solver)
       pde.options_.solver = static_cast<solve_opts>(H5Easy::load<int>(file, "solver_method"));
-    if (not pde.options_.isolver_tolerance)
+    if (not pde.options_.isolver_tolerance) {
       pde.options_.isolver_tolerance = H5Easy::load<double>(file, "solver_itol");
-    if (not pde.options_.isolver_iterations)
+      if (pde.options_.isolver_tolerance.value() < 0)
+        pde.options_.isolver_tolerance = pde.options_.default_isolver_tolerance;
+    }
+    if (not pde.options_.isolver_iterations) {
       pde.options_.isolver_iterations = H5Easy::load<int>(file, "solver_iter");
-    if (not pde.options_.isolver_outer_iterations)
-      pde.options_.isolver_outer_iterations = H5Easy::load<int>(file, "solver_oiter");
+      if (pde.options_.isolver_iterations.value() < 0)
+        pde.options_.isolver_iterations = pde.options_.default_isolver_iterations;
+    }
+    if (not pde.options_.isolver_inner_iterations) {
+      pde.options_.isolver_inner_iterations = H5Easy::load<int>(file, "solver_inner");
+      if (pde.options_.isolver_inner_iterations.value() < 0)
+        pde.options_.isolver_inner_iterations = pde.options_.default_isolver_inner_iterations;
+    }
   }
 
   state = H5Easy::load<std::vector<P>>(file, "state");

src/asgard_program_options.cpp

@@ -113,8 +113,9 @@ Options          Short   Value      Description
                                     adi - very experimental, not very stable (yet)
 -isolve-tol      -ist    double     Iterative solver tolerance, applies to GMRES and BICG.
 -isolve-iter     -isi    int        Iterative solver maximum number of iterations,
-                                    for GMRES this is the number of inner iterations.
--isolve-outer    -iso    int        (GMRES only) The maximum number of outer GMRES iterations.
+                                    for GMRES this is the number of outer iterations.
+-isolve-inner    -isn    int        (GMRES only) The maximum number of inner GMRES iterations,
+                                    this is ignored by BiCGSTAB.
 
 Leaving soon:
 -memory                  int        Memory limit for the GPU, applied to the earlier versions
@@ -222,8 +223,8 @@ void prog_opts::process_inputs(std::vector<std::string_view> const &argv,
       {"-kron-mode", optentry::kron_mode},
       {"-isolve-tol", optentry::isol_tolerance}, {"-ist", optentry::isol_tolerance},
       {"-isolve-iter", optentry::isol_iterations}, {"-isi", optentry::isol_iterations},
-      {"-isolve-outer", optentry::isol_outer_iterations},
-      {"-iso", optentry::isol_outer_iterations},
+      {"-isolve-inner", optentry::isol_inner_iterations},
+      {"-isn", optentry::isol_inner_iterations},
       {"-restart", optentry::restart_file},
   };
 
@@ -558,12 +559,12 @@ void prog_opts::process_inputs(std::vector<std::string_view> const &argv,
       }
     }
     break;
-    case optentry::isol_outer_iterations: {
+    case optentry::isol_inner_iterations: {
       auto selected = move_process_next();
       if (not selected)
         throw std::runtime_error(report_no_value());
       try {
-        isolver_outer_iterations = std::stoi(selected->data());
+        isolver_inner_iterations = std::stoi(selected->data());
       } catch(std::invalid_argument &) {
         throw std::runtime_error(report_wrong_value());
       } catch(std::out_of_range &) {

src/asgard_program_options.hpp

@@ -69,7 +69,7 @@ enum class preconditioner_opts
   none = 0,
   //! diagonal Jacobi preconditioner
   jacobi,
-  //! using alternating direction pseudoinverse
+  //! using alternating direction pseudoinverse (experimental)
   adi
 };
 
@@ -430,7 +430,7 @@ struct prog_opts
   //! max number of iterations (inner iterations for gmres)
   std::optional<int> isolver_iterations;
   //! max number of output gmres iterations
-  std::optional<int> isolver_outer_iterations;
+  std::optional<int> isolver_inner_iterations;
 
   //! local kron method only, mode dense or sparse (faster but memory hungry)
   std::optional<kronmult_mode> kron_mode;
@@ -626,6 +626,12 @@ struct prog_opts
   std::optional<solve_opts> default_solver;
   //! used in place of the preconditioner type, if preconditioner is not specified
   std::optional<preconditioner_opts> default_precon;
+  //! used in place of the tolerance, if tolerance is not specified
+  std::optional<double> default_isolver_tolerance;
+  //! max number of iterations (inner iterations for gmres)
+  std::optional<int> default_isolver_iterations;
+  //! max number of outer gmres iterations
+  std::optional<int> default_isolver_inner_iterations;
 
   //! returns the first available from stop-time, default-stop-time or -1
   double get_stop_time() const { return stop_time.value_or(default_stop_time.value_or(-1)); }
@@ -683,7 +689,7 @@ struct prog_opts
     kron_mode,
     isol_tolerance,
     isol_iterations,
-    isol_outer_iterations,
+    isol_inner_iterations,
     restart_file,
     set_verbosity
   };

src/asgard_program_options_tests.cpp

@@ -245,15 +245,15 @@ TEST_CASE("new program options", "[single options]")
     REQUIRE_THROWS_WITH(prog_opts(vecstrview({"exe", "-isi", "dummy"})),
                         "invalid value for -isi, see exe -help");
   }
-  SECTION("-isolve_outer")
+  SECTION("-isolve_inner")
   {
-    prog_opts prog(vecstrview({"", "-isolve-outer", "200"}));
-    REQUIRE(prog.isolver_outer_iterations);
-    REQUIRE(prog.isolver_outer_iterations.value() == 200);
-    REQUIRE_THROWS_WITH(prog_opts(vecstrview({"exe", "-isolve-outer"})),
-                        "-isolve-outer must be followed by a value, see exe -help");
-    REQUIRE_THROWS_WITH(prog_opts(vecstrview({"exe", "-iso", "dummy"})),
-                        "invalid value for -iso, see exe -help");
+    prog_opts prog(vecstrview({"", "-isolve-inner", "200"}));
+    REQUIRE(prog.isolver_inner_iterations);
+    REQUIRE(prog.isolver_inner_iterations.value() == 200);
+    REQUIRE_THROWS_WITH(prog_opts(vecstrview({"exe", "-isolve-inner"})),
+                        "-isolve-inner must be followed by a value, see exe -help");
+    REQUIRE_THROWS_WITH(prog_opts(vecstrview({"exe", "-isn", "dummy"})),
+                        "invalid value for -isn, see exe -help");
   }
   SECTION("-title")
   {

src/asgard_solver.cpp

@@ -657,15 +657,15 @@ int bicgstab<P>::solve(
   auto dot = [&](std::vector<P> const &a, std::vector<P> const &b)
     -> P {
       P r = 0;
-//ASGARD_OMP_PARFOR_SIMD_EXTRA(reduction(+:r))
+ASGARD_OMP_PARFOR_SIMD_EXTRA(reduction(+:r))
       for (int64_t i = 0; i < n; i++)
         r += a[i] * b[i];
       return r;
     };
   auto dot1 = [&](std::vector<P> const &a)
     -> P {
       P r = 0;
-//ASGARD_OMP_PARFOR_SIMD_EXTRA(reduction(+:r))
+ASGARD_OMP_PARFOR_SIMD_EXTRA(reduction(+:r))
       for (int64_t i = 0; i < n; i++)
         r += a[i] * a[i];
       return r;
@@ -684,7 +684,7 @@ ASGARD_OMP_PARFOR_SIMD
   r = rhs;
 
   int num_appy = 1;
-  apply_lhs(-1, x, 1, r); // r0 = b - A * x0
+  apply_lhs(-1, x.data(), 1, r.data()); // r0 = b - A * x0
 
   P rho = dot1(r);
 
@@ -693,7 +693,7 @@ ASGARD_OMP_PARFOR_SIMD
 
   for (int i = 0; i < max_iter_; i++) {
     ++num_appy;
-    apply_lhs(1, p, 0, v); // v = A * p
+    apply_lhs(1, p.data(), 0, v.data()); // v = A * p
 
     P const alpha = rho / dot(rref, v);
 
@@ -705,7 +705,7 @@ ASGARD_OMP_PARFOR_SIMD
     }
 
     ++num_appy;
-    apply_lhs(1, r, 0, t); // t = A * p
+    apply_lhs(1, r.data(), 0, t.data()); // t = A * p
 
     P const omega = dot(r, t) / dot1(t);
 
@@ -732,8 +732,8 @@ ASGARD_OMP_PARFOR_SIMD
 
 template<typename P>
 int gmres<P>::solve(
-    operatoin_apply_precon_arr<P> apply_precon,
-    operatoin_apply_lhs_arr<P> apply_lhs, std::vector<P> const &rhs,
+    operatoin_apply_precon<P> apply_precon,
+    operatoin_apply_lhs<P> apply_lhs, std::vector<P> const &rhs,
     std::vector<P> &x) const
 {
   int const n = static_cast<int>(rhs.size());
@@ -936,6 +936,13 @@ ASGARD_OMP_PARFOR_SIMD
   }
 }
 
+template<typename P>
+void solver_manager<P>::xpby(std::vector<P> const &x, P beta, P y[]) {
+ASGARD_OMP_PARFOR_SIMD
+  for (size_t i = 0; i < x.size(); i++)
+    y[i] = x[i] + beta * y[i];
+}
+
 template<typename P>
 void solver_manager<P>::print_opts(std::ostream &os) const
 {
@@ -946,18 +953,18 @@ void solver_manager<P>::print_opts(std::ostream &os) const
       os << "  direct\n";
       break;
     case 1:
-      os << "  bicgstab\n";
-      os << "  tolerance:      " << std::get<solvers::bicgstab<P>>(var).tolerance() << '\n';
-      os << "  max iterations: " << std::get<solvers::bicgstab<P>>(var).max_iter() << '\n';
-      has_precon = true;
-      break;
-    case 2:
       os << "  gmres\n";
       os << "  tolerance: " << std::get<solvers::gmres<P>>(var).tolerance() << '\n';
       os << "  max inner: " << std::get<solvers::gmres<P>>(var).max_inner() << '\n';
       os << "  max outer: " << std::get<solvers::gmres<P>>(var).max_outer() << '\n';
       has_precon = true;
       break;
+    case 2:
+      os << "  bicgstab\n";
+      os << "  tolerance:      " << std::get<solvers::bicgstab<P>>(var).tolerance() << '\n';
+      os << "  max iterations: " << std::get<solvers::bicgstab<P>>(var).max_iter() << '\n';
+      has_precon = true;
+      break;
     default:
       break;
   }