Some utilities (#710)

* some utilities and minor fixes

src/asgard.cpp

@@ -8,6 +8,12 @@ void simulate(parser const &cli_input, std::unique_ptr<PDE<precision>> &pde)
 {
   options const opts(cli_input);
 
+  if (cli_input.show_libinfo())
+  {
+    print_info();
+    return;
+  }
+
   node_out() << "Branch: " << GIT_BRANCH << '\n';
   node_out() << "Commit Summary: " << GIT_COMMIT_HASH
                      << GIT_COMMIT_SUMMARY << '\n';
@@ -395,4 +401,49 @@ template void simulate(parser const &cli_input, std::unique_ptr<PDE<double>> &pd
 template void simulate(parser const &cli_input, std::unique_ptr<PDE<float>> &pde);
 #endif
 
+void print_info(std::ostream &os)
+{
+  os << "\nASGarD v" << ASGARD_VERSION << "  git-hash: " << GIT_COMMIT_HASH << "\n";
+  os << "git-branch (" << GIT_BRANCH << ")\n";
+#ifdef KRON_MODE_GLOBAL
+#ifdef KRON_MODE_GLOBAL_BLOCK
+  os << "Kronmult method          Block-Global\n";
+#else
+  os << "Kronmult method          Global\n";
+#endif
+#else
+  os << "Kronmult method          Local\n";
+#endif
+#ifdef ASGARD_USE_CUDA
+  os << "GPU Acceleration         CUDA\n";
+#else
+  os << "GPU Acceleration         Disabled\n";
+#endif
+#ifdef ASGARD_USE_OPENMP
+  os << "OpenMP multithreading    Enablded\n";
+#else
+  os << "OpenMP multithreading    Disabled\n";
+#endif
+#ifdef ASGARD_USE_MPI
+  os << "MPI distributed grid     Enabled\n";
+#else
+  os << "MPI distributed grid     Disabled\n";
+#endif
+#ifdef ASGARD_IO_HIGHFIVE
+  os << "HDF5 - HighFive I/O      Enabled\n";
+#else
+  os << "HDF5 - HighFive I/O      Disabled\n";
+#endif
+#ifdef ASGARD_ENABLE_DOUBLE
+#ifdef ASGARD_ENABLE_FLOAT
+  os << "Available precisions     double/float\n";
+#else
+  os << "Available precision      double\n";
+#endif
+#else
+  os << "Available precision      float\n";
+#endif
+  os << '\n';
+}
+
 } // namespace asgard

src/asgard.hpp

@@ -1,11 +1,6 @@
 #pragma once
 
-#include "batch.hpp"
-
-#include "asgard_tools.hpp"
 #include "coefficients.hpp"
-#include "distribution.hpp"
-#include "elements.hpp"
 
 #ifdef ASGARD_IO_HIGHFIVE
 #include "io.hpp"
@@ -19,12 +14,7 @@
 #include "matlab_plot.hpp"
 #endif
 
-#include "asgard_vector.hpp"
-#include "pde.hpp"
-#include "program_options.hpp"
 #include "time_advance.hpp"
-#include "transformations.hpp"
-#include <numeric>
 
 namespace asgard
 {
@@ -40,4 +30,6 @@ void simulate(parser const &cli_input)
   simulate(cli_input, pde);
 }
 
+void print_info(std::ostream &os = std::cout);
+
 }

src/asgard_interpolation.hpp

@@ -19,7 +19,7 @@ class interpolation
 {
 public:
   //! constructs and empty interpolation class
-  interpolation() : num_dimenisons_(0), block_size(0), workspace_(nullptr)
+  interpolation()
   {}
   //! convert to true if the class has been initialized
   operator bool() const { return (num_dimenisons_ > 0); }
@@ -158,14 +158,14 @@ class interpolation
   }
 
 private:
-  int num_dimenisons_;
+  int num_dimenisons_         = 0;
   static constexpr int pterms = 2; // remove static later
-  int64_t block_size;
+  int64_t block_size          = 0;
   kronmult::permutes perms;
 
   wavelet_interp1d<1, precision> wav1d;
 
-  mutable kronmult::block_global_workspace<precision> *workspace_;
+  mutable kronmult::block_global_workspace<precision> *workspace_ = nullptr;
 };
 
 #endif

src/asgard_interpolation1d.hpp

@@ -28,7 +28,7 @@ class wavelet_interp1d
   static constexpr int matsize = pterms * pterms;
 
   //! null constructor
-  wavelet_interp1d() : conn(nullptr)
+  wavelet_interp1d()
   {}
   //! cache interpolation points and transformation matrices
   wavelet_interp1d(connect_1d const *conn_in) : conn(conn_in)
@@ -94,7 +94,7 @@ class wavelet_interp1d
                          precision mat[]);
 
 private:
-  connect_1d const *conn;
+  connect_1d const *conn = nullptr;
   std::vector<precision> nodes_;
   std::vector<precision> proj2node_;
   std::vector<precision> node2hier_;

src/asgard_kronmult_matrix.hpp

@@ -652,13 +652,13 @@ class local_kronmult_matrix
     flops_ *= int64_t{tensor_size_} * kron_size_;
   }
 
-  int num_dimensions_;
-  int kron_size_; // i.e., n - size of the matrices
-  int num_rows_;
-  int num_cols_;
-  int num_terms_;
-  int64_t tensor_size_;
-  int64_t flops_;
+  int num_dimensions_  = 0;
+  int kron_size_       = 0; // i.e., n - size of the matrices
+  int num_rows_        = 0;
+  int num_cols_        = 0;
+  int num_terms_       = 0;
+  int64_t tensor_size_ = 0;
+  int64_t flops_       = 0;
 
 #ifdef ASGARD_USE_CUDA
   // indicates that the input vectors for single-call-mode will be on the GPU
@@ -704,7 +704,9 @@ class local_kronmult_matrix
   std::vector<fk::vector<precision, mem_type::owner, data_mode>> terms_;
   fk::vector<precision *, mem_type::owner, data_mode> term_pntr_;
   fk::vector<int, mem_type::owner, data_mode> elem_;
-  int row_offset_, col_offset_, num_1d_blocks_;
+  int row_offset_    = 0;
+  int col_offset_    = 0;
+  int num_1d_blocks_ = 0;
 
   // preconditioner
   std::vector<precision> pre_con_;
@@ -1007,9 +1009,9 @@ class global_kron_matrix
 private:
   // description of the multi-indexes and the sparsity pattern
   // global case data
-  int num_dimensions_;
-  int64_t num_active_;
-  int64_t num_padded_;
+  int num_dimensions_ = 0;
+  int64_t num_active_ = 0;
+  int64_t num_padded_ = 0;
   std::vector<kronmult::permutes> perms_;
   std::array<int64_t, num_imex_variants> flops_;
   // data for the 1D tensors
@@ -1020,7 +1022,7 @@ class global_kron_matrix
   std::vector<std::vector<precision>> gvals_;
   // collections of terms
   std::array<std::vector<int>, 3> term_groups;
-  int porder_;
+  int porder_ = 0;
   // preconditioner
   std::vector<precision> pre_con_;
 #ifdef ASGARD_USE_CUDA
@@ -1154,9 +1156,7 @@ template<typename precision>
 class block_global_kron_matrix
 {
 public:
-  block_global_kron_matrix()
-      : num_active_(0), num_padded_(0), num_dimensions_(0), blockn_(0), block_size_(0),
-        conn_volumes_(nullptr), conn_full_(nullptr), workspace_(nullptr) {}
+  block_global_kron_matrix() {}
 
   block_global_kron_matrix(int64_t num_active, int64_t num_padded,
                            int num_dimensions, int blockn, int64_t block_size,
@@ -1261,20 +1261,22 @@ class block_global_kron_matrix
       block_global_kron_matrix<precision> &mat);
 
 private:
-  int64_t num_active_, num_padded_;
-  int num_dimensions_, blockn_;
-  int64_t block_size_;
+  int64_t num_active_ = 0;
+  int64_t num_padded_ = 0;
+  int num_dimensions_ = 0;
+  int blockn_         = 0;
+  int64_t block_size_ = 0;
   vector2d<int> ilist_;
   dimension_sort dsort_;
   std::vector<kronmult::permutes> perms_;
   std::vector<int> flux_dir_;
-  connect_1d const *conn_volumes_;
-  connect_1d const *conn_full_;
+  connect_1d const *conn_volumes_ = nullptr;
+  connect_1d const *conn_full_    = nullptr;
 
   std::vector<std::vector<precision>> gvals_;
   std::array<std::vector<int>, 3> term_groups_;
 
-  mutable kronmult::block_global_workspace<precision> *workspace_;
+  mutable kronmult::block_global_workspace<precision> *workspace_ = nullptr;
 
   mutable std::array<int64_t, num_imex_variants> flops_;
 

src/asgard_reconstruct.hpp

@@ -1,7 +1,6 @@
 #pragma once
 
 #include "asgard_interpolation.hpp"
-#include "elements.hpp"
 
 namespace asgard
 {

src/asgard_tools.cpp

@@ -110,6 +110,4 @@ void simple_timer::get_timing_stats(
   }
 }
 
-simple_timer timer;
-
 } // namespace asgard::tools

src/asgard_tools.hpp

@@ -143,7 +143,7 @@ class simple_timer
       id_to_start_;
 };
 
-extern simple_timer timer;
+inline simple_timer timer;
 
 /*!
  * Allows for RAII style of timing for blocks of code.
@@ -170,3 +170,105 @@ struct time_event
 };
 
 } // namespace asgard::tools
+
+namespace asgard
+{
+/*!
+ * \brief Iterator/generator for a sequence of integers
+ *
+ * This is needed for the indexof template
+ *
+ * Technically satisfies the requirements for legacy iterator
+ * but do not use directly, will be used internally in indexof
+ */
+template<typename idx_type = int64_t>
+struct index_iterator
+{
+  using iterator_category = std::random_access_iterator_tag;
+
+  using value_type      = idx_type;
+  using difference_type = idx_type;
+  using reference       = idx_type &;
+  using pointer         = idx_type *;
+
+  idx_type &operator*() { return value_; }
+  idx_type const &operator*() const { return value_; }
+  bool operator!=(index_iterator const &other) const { return value_ != other.value_; }
+  index_iterator &operator++()
+  {
+    ++value_;
+    return *this;
+  }
+  index_iterator &operator++(int) { return index_iterator{value_++}; }
+  index_iterator &operator--()
+  {
+    --value_;
+    return *this;
+  }
+  index_iterator &operator--(int) { return index_iterator{value_--}; }
+
+  idx_type value_;
+};
+
+/*!
+ * \brief Allows for range for-loops but using indexes
+ *
+ * There is a repeated pattern in coding when cross-referencing entries
+ * between different vectors:
+ * \code
+ *   for (size_t i = 0; i < u.size(); i++)
+ *     u[i] = std::sqrt(x[i]);
+ * \endcode
+ * The operation can be done with a std::transform but it leads to a messy
+ * lambda capture and potential shadow. The index can be used to cross
+ * reference more complex structures where iterators would be messy and
+ * non-trivial, e.g., rows/columns of a matrix, sparse grid indexes, or
+ * entries in a vector2d. The index also helps keep a more expressive
+ * mathematical notation.
+ *
+ * On the other hand, the pattern is tedious to write over and over.
+ *
+ * This template provides an alternative and allows for syntax like:
+ * \code
+ *   for (auto i : indexof(u)) // i is int64_t
+ *     u[i] = std::sqrt(x[i]);
+ *
+ *   for (auto i : indexof<int>(u)) // i is int
+ *     u[i] = std::sqrt(x[i]);
+ *
+ *   for (auto i : indexof<size_t>(1, num_dimensions)) // i is size_t
+ *     u[i] = std::sqrt(x[i]);
+ * \endcode
+ *
+ * At -O3 Godbolt compiler profile yields the same code as for the constructs
+ * for-indexof and the regular for-loop.
+ */
+template<typename idx_type = int64_t>
+struct indexof
+{
+  template<typename vector_type>
+  indexof(vector_type const &f)
+      : beg_(0), end_(static_cast<idx_type>(f.size()))
+  {}
+  indexof(int num)
+      : beg_(0), end_(static_cast<idx_type>(num))
+  {}
+  indexof(int64_t num)
+      : beg_(0), end_(static_cast<idx_type>(num))
+  {}
+  indexof(size_t num)
+      : beg_(0), end_(static_cast<idx_type>(num))
+  {}
+  template<typename cidx_type>
+  indexof(cidx_type b, cidx_type e)
+      : beg_(b), end_(e)
+  {}
+
+  index_iterator<idx_type> begin() const { return index_iterator<idx_type>{beg_}; }
+  index_iterator<idx_type> end() const { return index_iterator<idx_type>{end_}; }
+
+  idx_type beg_;
+  idx_type end_;
+};
+
+} // namespace asgard