Draft implementation of som_core::compute_final_solution_cc()

This method computes the final solution using the consistent
constraints approach as described in Section II.B of
O. Goulko et al. Phys. Rev. B 95, 014102 (2017).

The new method is also exposed in Python, although I could not mark it
as `release_GIL_and_enable_signal = True` because of cpp2py issue
TRIQS/cpp2py#44.

c++/som/kernels/mesh_traits.hpp

@@ -21,6 +21,7 @@
 #pragma once
 
 #include <string>
+#include <variant>
 
 #include <triqs/gfs.hpp>
 
@@ -49,4 +50,8 @@ template <> struct mesh_traits<triqs::gfs::refreq> {
   inline static std::string name() { return "real frequency"; }
 };
 
+using mesh_variant_t = std::variant<triqs::gfs::gf_mesh<triqs::gfs::imtime>,
+                                    triqs::gfs::gf_mesh<triqs::gfs::imfreq>,
+                                    triqs::gfs::gf_mesh<triqs::gfs::legendre>>;
+
 } // namespace som

c++/som/som_core/accumulate.cpp

@@ -26,7 +26,8 @@
 #include <utility>
 #include <vector>
 
-#include <boost/preprocessor/seq/for_each_product.hpp>
+#include <boost/preprocessor/seq/elem.hpp>
+#include <boost/preprocessor/seq/enum.hpp>
 
 #include <triqs/utility/signal_handler.hpp>
 
@@ -226,16 +227,12 @@ void som_core::accumulate(accumulate_parameters_t const& p) {
   triqs::signal_handler::start();
   accumulate_status = 0;
   try {
-#define RUN_IMPL_CASE(r, okmk)                                                 \
-  case(int(BOOST_PP_SEQ_ELEM(0, okmk)) +                                       \
-       n_observable_kinds * mesh_traits<BOOST_PP_SEQ_ELEM(1, okmk)>::index):   \
-    accumulate_impl<kernel<BOOST_PP_SEQ_ENUM(okmk)>>();                        \
-    break;
-    switch(int(kind) + n_observable_kinds * mesh.index()) {
-      BOOST_PP_SEQ_FOR_EACH_PRODUCT(RUN_IMPL_CASE,
-                                    (ALL_OBSERVABLES)(ALL_INPUT_MESHES))
-    }
-#undef RUN_IMPL_CASE
+#define IMPL_CASE(r, okmk)                                                     \
+    case(kernel_id(BOOST_PP_SEQ_ELEM(0, okmk), BOOST_PP_SEQ_ELEM(1, okmk){})): \
+      accumulate_impl<kernel<BOOST_PP_SEQ_ENUM(okmk)>>(); break;
+
+    SELECT_KERNEL(IMPL_CASE, som_core::accumulate())
+#undef IMPL_CASE
   } catch(stopped& e) {
     accumulate_status = e.code;
     triqs::signal_handler::received(true);

c++/som/som_core/adjust_f.cpp

@@ -168,18 +168,12 @@ int som_core::adjust_f(adjust_f_parameters_t const& p) {
               << " particular solutions ..." << std::endl;
   }
 
-#define RUN_IMPL_CASE(r, okmk)                                                 \
-  case(int(BOOST_PP_SEQ_ELEM(0, okmk)) +                                       \
-       n_observable_kinds * mesh_traits<BOOST_PP_SEQ_ELEM(1, okmk)>::index):   \
+#define IMPL_CASE(r, okmk)                                                     \
+  case(kernel_id(BOOST_PP_SEQ_ELEM(0, okmk), BOOST_PP_SEQ_ELEM(1, okmk){})):   \
     return adjust_f_impl<kernel<BOOST_PP_SEQ_ENUM(okmk)>>(p);
-  switch(int(kind) + n_observable_kinds * mesh.index()) {
-    BOOST_PP_SEQ_FOR_EACH_PRODUCT(RUN_IMPL_CASE,
-                                  (ALL_OBSERVABLES)(ALL_INPUT_MESHES))
-    default:
-      fatal_error("Unknown observable kind " + to_string(kind) +
-                  " in adjust_f()");
-  }
-#undef RUN_IMPL_CASE
+
+  SELECT_KERNEL(IMPL_CASE, som_core::adjust_f())
+#undef IMPL_CASE
 }
 
 } // namespace som

c++/som/som_core/common.hxx

@@ -44,7 +44,7 @@ inline void warning(std::string const& message) {
   std::cout << "WARNING: " << message << std::endl;
 }
 
-inline std::ostream & mpi_cout(mpi::communicator const& comm) {
+inline std::ostream& mpi_cout(mpi::communicator const& comm) {
   return (std::cout << "[Rank " << comm.rank() << "] ");
 }
 
@@ -74,4 +74,24 @@ void check_gf_stat(gf_view<MeshType> g, statistic_enum expected_stat) {
   return check_gf_stat(make_const_view(g), expected_stat);
 }
 
+template <typename MeshType>
+inline constexpr int kernel_id(observable_kind kind, MeshType) {
+  return int(kind) + n_observable_kinds * mesh_traits<MeshType>::index;
+}
+
+inline int kernel_id(observable_kind kind, mesh_variant_t const& mesh) {
+  return int(kind) + n_observable_kinds * mesh.index();
+}
+
+#define FOR_EACH_KERNEL(F)                                                     \
+  BOOST_PP_SEQ_FOR_EACH_PRODUCT(F, (ALL_OBSERVABLES)(ALL_INPUT_MESHES))
+
+#define SELECT_KERNEL(F, F_NAME)                                               \
+  switch(kernel_id(kind, mesh)) {                                              \
+    FOR_EACH_KERNEL(F)                                                         \
+    default:                                                                   \
+      fatal_error("Unknown observable kind " + std::to_string(kind) +          \
+                  " in " #F_NAME);                                             \
+  }
+
 } // namespace som