Complex triangular

CMakeLists.txt

@@ -122,6 +122,7 @@ set(DCA_LIBS
   ctaux
   ctint
   ss_ct_hyb
+  g0_interpolation
   function_transform
   gaussian_quadrature
   nfft
@@ -140,21 +141,24 @@ set(DCA_LIBS
   quantum_domains
   time_and_frequency_domains
   signals
-  symmetrization
   coarsegraining
   ${DCA_CONCURRENCY_LIB}
   ${DCA_THREADING_LIBS}
   lapack
   cuda_utils
 )
 
+if(TARGET stdc++fs)    # std::filesystem
+  list(APPEND DCA_LIBS stdc++fs)
+endif()
+
 if (DCA_HAVE_CUDA)
   list(APPEND DCA_LIBS
     blas_kernels
-    dnfft_kernels
     lapack_kernels
     mc_kernels
     special_transform_kernels
+    cuda_utils
     ${DCA_CUDA_LIBS})
 endif()
 

applications/cluster_solver_check/cluster_solver_check.cpp

@@ -99,7 +99,7 @@ int main(int argc, char** argv) {
     ed_solver.finalize(dca_loop_data);
 
     const auto Sigma_ed(dca_data_imag.Sigma);
-    const int tested_frequencies = 10;
+    const int tested_frequencies = 5;
     const auto G_ed(dca::math::util::cutFrequency(dca_data_imag.G_k_w, tested_frequencies));
 
     if (concurrency.id() == concurrency.first()) {

applications/dca/main_dca.cpp

@@ -22,6 +22,7 @@
 #include "dca/util/git_version.hpp"
 #include "dca/util/modules.hpp"
 #include "dca/util/signal_handler.hpp"
+#include "dca/application/dca_loop_dispatch.hpp"
 
 int main(int argc, char** argv) {
   if (argc < 2) {

build-aux/summit.cmake

@@ -42,3 +42,6 @@ set(MAGMA_DIR $ENV{OLCF_MAGMA_ROOT} CACHE PATH
 # FFTW paths.
 set(FFTW_INCLUDE_DIR $ENV{OLCF_FFTW_ROOT}/include CACHE PATH "Path to fftw3.h.")
 set(FFTW_LIBRARY $ENV{OLCF_FFTW_ROOT}/lib/libfftw3.so CACHE FILEPATH "The FFTW3(-compatible) library.")
+
+# Flags for summit.
+set(CMAKE_CXX_FLAGS -march=power9)

cmake/dca_config.cmake

@@ -73,9 +73,9 @@ endif()
 
 # Lattice type
 set(DCA_LATTICE "square" CACHE STRING "Lattice type, options are: bilayer | square | triangular |
-    hund | twoband_Cu | threeband | FeAs.")
+    hund | twoband_Cu | threeband | FeAs | Rashba.")
 set_property(CACHE DCA_LATTICE PROPERTY STRINGS bilayer square triangular hund twoband_Cu threeband
-             FeAs)
+             FeAs Rashba)
 
 if (DCA_LATTICE STREQUAL "bilayer")
   set(DCA_LATTICE_TYPE dca::phys::models::bilayer_lattice<PointGroup>)
@@ -99,22 +99,24 @@ elseif (DCA_LATTICE STREQUAL "threeband")
   set(DCA_LATTICE_INCLUDE
     "dca/phys/models/analytic_hamiltonians/threeband_hubbard.hpp")
 
-elseif (DCA_LATTICE STREQUAL "twoband_chain")
-  set(DCA_LATTICE_TYPE dca::phys::models::twoband_chain<dca::phys::domains::no_symmetry<1>>)
-  set(DCA_LATTICE_INCLUDE
-      "dca/phys/models/analytic_hamiltonians/hund_lattice.hpp")
 elseif (DCA_LATTICE STREQUAL "FeAs")
   set(DCA_LATTICE_TYPE dca::phys::models::FeAsLattice<PointGroup>)
   set(DCA_LATTICE_INCLUDE
       "dca/phys/models/analytic_hamiltonians/fe_as_lattice.hpp")
+
 elseif (DCA_LATTICE STREQUAL "twoband_Cu")
   set(DCA_LATTICE_TYPE dca::phys::models::TwoBandCu<PointGroup>)
   set(DCA_LATTICE_INCLUDE
       "dca/phys/models/analytic_hamiltonians/twoband_Cu.hpp")
 
+elseif (DCA_LATTICE STREQUAL "Rashba")
+  set(DCA_LATTICE_TYPE dca::phys::models::RashbaHubbard<PointGroup>)
+  set(DCA_LATTICE_INCLUDE
+          "dca/phys/models/analytic_hamiltonians/rashba_hubbard.hpp")
+
 else()
   message(FATAL_ERROR "Please set DCA_LATTICE to a valid option: bilayer | square | triangular |
-          hund | twoband_Cu | threeband | FeAs.")
+          hund | twoband_Cu | threeband | FeAs | Rashba.")
 endif()
 
 # Model type
@@ -307,22 +309,14 @@ else()
 endif()
 
 option(DCA_WITH_SINGLE_PRECISION_MC "Perform Monte Carlo and measurements in single precision." OFF)
-option(DCA_WITH_SINGLE_PRECISION_TP_MEASUREMENTS "Measure two particle function in single precision." OFF)
 
 if (DCA_WITH_SINGLE_PRECISION_MC)
   set(DCA_WITH_SINGLE_PRECISION_TP_MEASUREMENTS ON CACHE BOOL "Measure two particle function in single precision." FORCE)
-  set(MC_SCALAR float)
-else()
-  set(MC_SCALAR double)
-endif()
-
-if (DCA_WITH_SINGLE_PRECISION_TP_MEASUREMENTS)
-  set(TP_ACCUMULATION_SCALAR float)
+  set(MC_SINGLE_PRECISION true)
 else()
-  set(TP_ACCUMULATION_SCALAR double)
+  set(MC_SINGLE_PRECISION false)
 endif()
 
-
 option(DCA_WITH_MANAGED_MEMORY "Use managed memory allocator." OFF)
 mark_as_advanced(DCA_WITH_MANAGED_MEMORY)
 if (DCA_WITH_MANAGED_MEMORY)
@@ -335,15 +329,6 @@ configure_file("${PROJECT_SOURCE_DIR}/include/dca/config/mc_options.hpp.in"
         "${CMAKE_BINARY_DIR}/include/dca/config/mc_options.hpp" @ONLY)
 
 
-################################################################################
-# Symmetrization
-option(DCA_SYMMETRIZE "Apply cluster, time and frequency symmetries to single particle functions."
-       ON)
-
-if(DCA_SYMMETRIZE)
-  add_compile_definitions(DCA_WITH_SYMMETRIZATION)
-endif()
-
 ################################################################################
 # Generate applications' config files.
 configure_file("${PROJECT_SOURCE_DIR}/include/dca/config/analysis.hpp.in"

cmake/dca_cuda.cmake

@@ -20,6 +20,7 @@ if (CUDA_FOUND)
   list(APPEND DCA_CUDA_LIBS ${CUDA_LIBRARIES} ${CUDA_cusparse_LIBRARY} ${CUDA_cublas_LIBRARY})
   CUDA_INCLUDE_DIRECTORIES(${CUDA_INCLUDE_DIRS})
   set(CUDA_SEPARABLE_COMPILATION ON)
+  list(APPEND CUDA_NVCC_FLAGS "--expt-relaxed-constexpr")
 
   set(CVD_LAUNCHER "" CACHE INTERNAL "launch script for setting the Cuda visible devices.")
   # Use the following script for systems with multiple gpus visible from a rank.

cmake/dca_external_libs.cmake

@@ -27,6 +27,7 @@ if (NOT HDF5_FOUND)
   find_package(HDF5 REQUIRED COMPONENTS C CXX)
 endif()
 
+find_package(HDF5 REQUIRED COMPONENTS C CXX)
 list(APPEND DCA_EXTERNAL_LIBS ${HDF5_LIBRARIES})
 list(APPEND DCA_EXTERNAL_INCLUDE_DIRS ${HDF5_INCLUDE_DIRS})
 

include/dca/config/cmake_options.hpp

@@ -40,12 +40,10 @@ struct CMakeOptions {
   static const std::string dca_with_autotuning;
   static const std::string dca_with_gnuplot;
   static const std::string dca_with_single_precision_mc;
-  static const std::string dca_with_single_precision_tp_measurements;
   static const std::string dca_with_memory_savings;
   static const std::string dca_with_managed_memory;
 
   static const std::string dca_with_qmc_bit;
-
   static void print();
 };
 

include/dca/config/mc_options.hpp.in

@@ -12,6 +12,8 @@
 #ifndef DCA_CONFIG_MC_OPTIONS_HPP
 #define DCA_CONFIG_MC_OPTIONS_HPP
 
+#include <complex>
+#include <type_traits>
 #ifdef DCA_HAVE_CUDA
 #include "dca/linalg/util/allocators/device_allocator.hpp"
 #include "dca/linalg/util/allocators/managed_allocator.hpp"
@@ -22,9 +24,9 @@ namespace config {
 // dca::config::
 
 struct McOptions {
-  using MCScalar = @MC_SCALAR@;
-
-  using TPAccumulationScalar = @TP_ACCUMULATION_SCALAR@;
+  static constexpr bool single_precision = @MC_SINGLE_PRECISION@;
+  using Real = std::conditional_t<single_precision, float, double>;
+  using Complex = std::complex<Real>;
 
   static constexpr bool memory_savings = @MEMORY_SAVINGS@;
 

include/dca/function/domains/dmn.hpp

@@ -22,7 +22,7 @@
 #include <string>
 #include <vector>
 
-#include "dca/util/type_utils.hpp"  // for dca::util::type_name
+#include "dca/util/print_type.hpp"  // for dca::util::type_name
 
 namespace dca {
 namespace func {

include/dca/function/function.hpp

@@ -32,6 +32,7 @@
 #include "dca/function/set_to_zero.hpp"
 #include "dca/util/pack_operations.hpp"
 #include "dca/util/integer_division.hpp"
+#include "dca/util/print_type.hpp"
 #include "dca/util/type_utils.hpp"
 
 namespace dca {
@@ -64,6 +65,8 @@ class function {
   function(const function<scalartype, domain>& other, const std::string& name) : function(other) {
     name_ = name;
   }
+  template <typename Scalar2>
+  function(const function<Scalar2, domain>& other);
 
   // Move constructor
   // Constructs the function with elements and name of other using move semantics.
@@ -82,9 +85,9 @@ class function {
   // Precondition: The other function has been resetted, if the domain had been initialized after
   //               the other function's construction.
   // Postcondition: The function's name is unchanged.
-  function<scalartype, domain>& operator=(const function<scalartype, domain>& other);
+  function& operator=(const function& other);
   template <typename Scalar2>
-  function<scalartype, domain>& operator=(const function<Scalar2, domain>& other);
+  function& operator=(const function<Scalar2, domain>& other);
 
   // Move assignment operator
   // Replaces the function's elements with those of other using move semantics.
@@ -105,7 +108,6 @@ class function {
   const std::string& get_name() const {
     return name_;
   }
-  // TODO: Remove this method and use constructor parameter instead.
   void set_name(const std::string& name) {
     name_ = name;
   }
@@ -120,6 +122,8 @@ class function {
   void resize(std::size_t nb_elements_new) {
     fnc_values_.resize(nb_elements_new);
   }
+
+  // TODO: use a more standard method name.
   // Returns the size of the leaf domain with the given index.
   // Does not return function values!
   int operator[](const int index) const {
@@ -140,10 +144,10 @@ class function {
   auto end() {
     return fnc_values_.end();
   }
-  const auto begin() const {
+  auto begin() const {
     return fnc_values_.begin();
   }
-  const auto end() const {
+  auto end() const {
     return fnc_values_.end();
   }
 
@@ -190,7 +194,7 @@ class function {
   // Enable only if all arguments are integral to prevent subind_to_linind(int*, int) to resolve to
   // subind_to_linind(int...) rather than subind_to_linind(const int* const, int).
   template <typename... Ts>
-  std::enable_if_t<util::if_all<std::is_integral<Ts>::value...>::value, int> subind_2_linind(
+  std::enable_if_t<util::ifAll(std::is_integral<Ts>::value...), int> subind_2_linind(
       const Ts... subindices) const {
     // We need to cast all subindices to the same type for dmn_variadic.
     return dmn(static_cast<int>(subindices)...);
@@ -347,6 +351,24 @@ function<scalartype, domain>::function(const function<scalartype, domain>& other
     throw std::logic_error("Copy construction from a not yet resetted function.");
 }
 
+template <typename scalartype, class domain>
+template <typename Scalar2>
+function<scalartype, domain>::function(const function<Scalar2, domain>& other)
+    : name_(other.get_name()),
+      function_type(__PRETTY_FUNCTION__),
+      dmn(),
+      Nb_sbdms(dmn.get_leaf_domain_sizes().size()),
+      size_sbdm(dmn.get_leaf_domain_sizes()),
+      step_sbdm(dmn.get_leaf_domain_steps()),
+      fnc_values_(dmn.get_size()) {
+  if (size() != other.size()) {
+    // The other function has not been resetted after the domain was initialized.
+    throw std::logic_error("Copy construction from a not yet resetted function.");
+  }
+
+  std::copy(other.begin(), other.end(), begin());
+}
+
 template <typename scalartype, class domain>
 function<scalartype, domain>::function(function<scalartype, domain>&& other)
     : name_(std::move(other.name_)),
@@ -362,8 +384,7 @@ function<scalartype, domain>::function(function<scalartype, domain>&& other)
 }
 
 template <typename scalartype, class domain>
-function<scalartype, domain>& function<scalartype, domain>::operator=(
-    const function<scalartype, domain>& other) {
+function<scalartype, domain>& function<scalartype, domain>::operator=(const function& other) {
   if (this != &other) {
     if (dmn.get_size() != other.dmn.get_size()) {
       // Domain had not been initialized when the functions were created.
@@ -385,10 +406,17 @@ template <typename Scalar, class domain>
 template <typename Scalar2>
 function<Scalar, domain>& function<Scalar, domain>::operator=(const function<Scalar2, domain>& other) {
   if (size() != other.size()) {
-    throw(std::logic_error("Function size does not match."));
+    // Domain had not been initialized when the functions were created.
+    // Reset this function and check again.
+    reset();
+  }
+
+  if (size() != other.size()) {
+    // The other function has not been resetted after the domain was initialized.
+    throw std::logic_error("Copy assignment from a not yet resetted function.");
   }
 
-  fnc_values_ = other.fnc_values_;
+  std::copy(other.begin(), other.end(), begin());
 
   return *this;
 }

include/dca/function/util/real_complex_conversion.hpp

@@ -18,6 +18,10 @@
 
 #include "dca/function/function.hpp"
 
+#ifdef DCA_HAVE_MPI
+#include <mpi.h>
+#endif
+
 namespace dca {
 namespace func {
 namespace util {
@@ -40,14 +44,35 @@ auto complex(const function<Scalartype, Dmn>& f) {
 template <typename Scalartype, typename Dmn>
 auto real(const function<std::complex<Scalartype>, Dmn>& f, const bool check_imaginary = false) {
   function<Scalartype, Dmn> f_real;
+  Scalartype max_img = 0, rel_max = 0;
 
   for (int i = 0; i < f_real.size(); ++i) {
-    if (check_imaginary && std::abs(f(i).imag()) > 500 * std::numeric_limits<Scalartype>::epsilon())
-      throw(std::logic_error("The function is not purely real."));
+    const auto val = std::abs(f(i).imag());
+    if (val > max_img) {
+      max_img = val;
+      rel_max = val / std::abs(f(i).real());
+    }
 
     f_real(i) = f(i).real();
   }
 
+  if (check_imaginary && max_img > 1e3 * std::numeric_limits<Scalartype>::epsilon()) {
+    // throw(std::logic_error("The function is not purely real."));
+
+    int id = 0;
+#ifdef DCA_HAVE_MPI
+    int initialized;
+    MPI_Initialized(&initialized);
+    if (initialized) {
+      MPI_Comm_rank(MPI_COMM_WORLD, &id);
+    }
+#endif  // DCA_HAVE_MPI
+    if (id == 0) {
+      std::cerr << "WARNING: The function" << f.get_name() << " is not purely real. Max img "
+                << max_img << " relative of max " << rel_max << std::endl;
+    }
+  }
+
   return f_real;
 }