Commit 814dbc8f authored by gbalduzz's avatar gbalduzz
Browse files

Store function elements in a vector for safer assignments.

parent e05746bd

include/dca/function/function.hpp

+75 −102
Original line number Diff line number Diff line
@@ -94,8 +94,6 @@ public: 
  //                 The other function is in a non-specified state.

  function<scalartype, domain>& operator=(function<scalartype, domain>&& other);



  ~function();



  // Resets the function by resetting the domain object and reallocating the memory for the function

  // elements.

  // Postcondition: All elements are set to zero.
@@ -115,11 +113,12 @@ public: 
    return Nb_sbdms;

  }

  std::size_t size() const {

    return nb_elements_;

    return fnc_values_.size();

  }



  // TODO: remove as it breaks class' invariant.

  void resize(std::size_t nb_elements_new) {

    nb_elements_ = nb_elements_new;

    fnc_values_.resize(nb_elements_new);

  }

  // Returns the size of the leaf domain with the given index.

  // Does not return function values!
@@ -128,31 +127,31 @@ public: 
  }



  // Begin and end methods for compatibility with range for loop.

  scalartype* begin() {

    return fnc_values;

  auto begin() {

    return fnc_values_.begin();

  }

  scalartype* end() {

    return fnc_values + nb_elements_;

  auto end() {

    return fnc_values_.end();

  }

  const scalartype* begin() const {

    return fnc_values;

  const auto begin() const {

    return fnc_values_.begin();

  }

  const scalartype* end() const {

    return fnc_values + nb_elements_;

  const auto end() const {

    return fnc_values_.end();

  }



  // Returns a pointer to the function's elements.

  scalartype* values() {

    return fnc_values;

    return fnc_values_.data();

  }

  const scalartype* values() const {

    return fnc_values;

    return fnc_values_.data();

  }

  scalartype* data() {

    return fnc_values;

    return fnc_values_.data();

  }

  const scalartype* data() const {

    return fnc_values;

    return fnc_values_.data();

  }



  //
@@ -194,7 +193,7 @@ public: 
  template <typename T>

  int subind_2_linind(const T ind) const {

    static_assert(std::is_integral<T>::value, "Index ind must be an integer.");

    assert(ind >= 0 && ind < nb_elements_);

    assert(ind >= 0 && ind < size());

    return ind;

  }
@@ -208,24 +207,24 @@ public: 
  template <typename T>

  scalartype& operator()(const T linind) {

    static_assert(std::is_integral<T>::value, "Index linind must be an integer.");

    assert(linind >= 0 && linind < nb_elements_);

    return fnc_values[linind];

    assert(linind >= 0 && linind < size());

    return fnc_values_[linind];

  }

  template <typename T>

  const scalartype& operator()(const T linind) const {

    static_assert(std::is_integral<T>::value, "Index linind must be an integer.");

    assert(linind >= 0 && linind < nb_elements_);

    return fnc_values[linind];

    assert(linind >= 0 && linind < size());

    return fnc_values_[linind];

  }



  template <typename... Ts>

  scalartype& operator()(const Ts... subindices) {

    // We need to cast all indices to the same type for dmn_variadic.

    return fnc_values[dmn(static_cast<int>(subindices)...)];

    return fnc_values_[dmn(static_cast<int>(subindices)...)];

  }

  template <typename... Ts>

  const scalartype& operator()(const Ts... subindices) const {

    return fnc_values[dmn(static_cast<int>(subindices)...)];

    return fnc_values_[dmn(static_cast<int>(subindices)...)];

  }



  void operator+=(const function<scalartype, domain>& other);
@@ -240,7 +239,7 @@ public: 
  void operator/=(scalartype c);



  // Equal-comparison opertor

  // Returns true if the function's elements (fnc_values) are equal to other's elements, false

  // Returns true if the function's elements (fnc_values_) are equal to other's elements, false

  // otherwise.

  // TODO: Make the equal-comparison operator a non-member function.

  bool operator==(const function<scalartype, domain>& other) const;
@@ -284,14 +283,12 @@ private: 


  domain dmn;  // TODO: Remove domain object?



  std::size_t nb_elements_;



  // The subdomains (sbdmn) represent the leaf domains, not the branch domains.

  int Nb_sbdms;

  const std::vector<std::size_t>& size_sbdm;  // TODO: Remove?

  const std::vector<std::size_t>& step_sbdm;  // TODO: Remove?



  scalartype* fnc_values;

  std::vector<scalartype> fnc_values_;

};



template <typename scalartype, class domain>
@@ -302,14 +299,12 @@ function<scalartype, domain>::function(const std::string& name) 
    : name_(name),

      function_type(__PRETTY_FUNCTION__),

      dmn(),

      nb_elements_(dmn.get_size()),

      Nb_sbdms(dmn.get_leaf_domain_sizes().size()),

      size_sbdm(dmn.get_leaf_domain_sizes()),

      step_sbdm(dmn.get_leaf_domain_steps()),

      fnc_values(nullptr) {

  fnc_values = new scalartype[nb_elements_];

  for (int linind = 0; linind < nb_elements_; ++linind)

    setToZero(fnc_values[linind]);

      fnc_values_(dmn.get_size()) {

  for (int linind = 0; linind < size(); ++linind)

    setToZero(fnc_values_[linind]);

}



template <typename scalartype, class domain>
@@ -318,19 +313,17 @@ function<scalartype, domain>::function(const std::string& name, const Concurrenc 
    : name_(name),

      function_type(__PRETTY_FUNCTION__),

      dmn(),

      nb_elements_(dmn.get_size()),

      Nb_sbdms(dmn.get_leaf_domain_sizes().size()),

      size_sbdm(dmn.get_leaf_domain_sizes()),

      step_sbdm(dmn.get_leaf_domain_steps()),

      fnc_values(nullptr) {

      step_sbdm(dmn.get_leaf_domain_steps()) {

  // TODO: multi-index access to partitioned function is not safe.

  const std::size_t mpi_size = concurrency.number_of_processors();



  nb_elements_ = dca::util::ceilDiv(nb_elements_, mpi_size);

  const std::size_t nb_elements = dca::util::ceilDiv(dmn.get_size(), mpi_size);

  fnc_values_.resize(nb_elements);



  fnc_values = new scalartype[nb_elements_];

  for (int linind = 0; linind < nb_elements_; ++linind)

    setToZero(fnc_values[linind]);

  for (int linind = 0; linind < nb_elements; ++linind)

    setToZero(fnc_values_[linind]);

}



template <typename scalartype, class domain>
@@ -338,17 +331,13 @@ function<scalartype, domain>::function(const function<scalartype, domain>& other 
    : name_(other.name_),

      function_type(__PRETTY_FUNCTION__),

      dmn(),

      nb_elements_(dmn.get_size()),

      Nb_sbdms(dmn.get_leaf_domain_sizes().size()),

      size_sbdm(dmn.get_leaf_domain_sizes()),

      step_sbdm(dmn.get_leaf_domain_steps()),

      fnc_values(nullptr) {

      fnc_values_(other.fnc_values_) {

  if (dmn.get_size() != other.dmn.get_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.");



  fnc_values = new scalartype[nb_elements_];

  std::copy_n(other.fnc_values, nb_elements_, fnc_values);

}



template <typename scalartype, class domain>
@@ -356,18 +345,13 @@ function<scalartype, domain>::function(function<scalartype, domain>&& other) 
    : name_(std::move(other.name_)),

      function_type(__PRETTY_FUNCTION__),

      dmn(),

      nb_elements_(dmn.get_size()),

      Nb_sbdms(dmn.get_leaf_domain_sizes().size()),

      size_sbdm(dmn.get_leaf_domain_sizes()),

      step_sbdm(dmn.get_leaf_domain_steps()),

      fnc_values(nullptr) {

      fnc_values_(std::move(other.fnc_values_)) {

  if (dmn.get_size() != other.dmn.get_size())

    // The other function has not been resetted after the domain was initialized.

    throw std::logic_error("Move construction from a not yet resetted function.");



  fnc_values = other.fnc_values;

  other.nb_elements_ = 0;

  other.fnc_values = nullptr;

}



template <typename scalartype, class domain>
@@ -384,7 +368,7 @@ function<scalartype, domain>& function<scalartype, domain>::operator=( 
        throw std::logic_error("Copy assignment from a not yet resetted function.");

    }



    std::copy_n(other.values(), nb_elements_, fnc_values);

    fnc_values_ = other.fnc_values_;

  }



  return *this;
@@ -397,7 +381,7 @@ function<Scalar, domain>& function<Scalar, domain>::operator=(const function<Sca 
    throw(std::logic_error("Function size does not match."));

  }



  std::copy_n(other.values(), nb_elements_, fnc_values);

  fnc_values_ = other.fnc_values_;



  return *this;

}
@@ -416,33 +400,21 @@ function<scalartype, domain>& function<scalartype, domain>::operator=( 
        throw std::logic_error("Move assignment from a not yet resetted function.");

    }



    delete[] fnc_values;

    fnc_values = other.fnc_values;



    other.nb_elements_ = 0;

    other.fnc_values = nullptr;

    fnc_values_ = std::move(other.fnc_values_);

  }



  return *this;

}



template <typename scalartype, class domain>

function<scalartype, domain>::~function() {

  delete[] fnc_values;

}



template <typename scalartype, class domain>

void function<scalartype, domain>::reset() {

  dmn.reset();



  nb_elements_ = dmn.get_size();

  fnc_values_.resize(dmn.get_size());

  Nb_sbdms = dmn.get_leaf_domain_sizes().size();



  delete[] fnc_values;

  fnc_values = new scalartype[nb_elements_];



  for (int linind = 0; linind < nb_elements_; ++linind)

    setToZero(fnc_values[linind]);

  for (int linind = 0; linind < size(); ++linind)

    setToZero(fnc_values_[linind]);

}



template <typename scalartype, class domain>
@@ -494,8 +466,8 @@ scalartype& function<scalartype, domain>::operator()(const int* const subind) { 
  int linind;

  subind_2_linind(subind, linind);



  assert(linind >= 0 && linind < nb_elements_);

  return fnc_values[linind];

  assert(linind >= 0 && linind < size());

  return fnc_values_[linind];

}



template <typename scalartype, class domain>
@@ -503,64 +475,64 @@ const scalartype& function<scalartype, domain>::operator()(const int* const subi 
  int linind;

  subind_2_linind(subind, linind);



  assert(linind >= 0 && linind < nb_elements_);

  return fnc_values[linind];

  assert(linind >= 0 && linind < size());

  return fnc_values_[linind];

}



template <typename scalartype, class domain>

void function<scalartype, domain>::operator+=(const function<scalartype, domain>& other) {

  for (int linind = 0; linind < nb_elements_; ++linind)

    fnc_values[linind] += other(linind);

  for (int linind = 0; linind < size(); ++linind)

    fnc_values_[linind] += other(linind);

}



template <typename scalartype, class domain>

void function<scalartype, domain>::operator-=(const function<scalartype, domain>& other) {

  for (int linind = 0; linind < nb_elements_; ++linind)

    fnc_values[linind] -= other(linind);

  for (int linind = 0; linind < size(); ++linind)

    fnc_values_[linind] -= other(linind);

}



template <typename scalartype, class domain>

void function<scalartype, domain>::operator*=(const function<scalartype, domain>& other) {

  for (int linind = 0; linind < nb_elements_; ++linind)

    fnc_values[linind] *= other(linind);

  for (int linind = 0; linind < size(); ++linind)

    fnc_values_[linind] *= other(linind);

}



template <typename scalartype, class domain>

void function<scalartype, domain>::operator/=(const function<scalartype, domain>& other) {

  for (int linind = 0; linind < nb_elements_; ++linind) {

  for (int linind = 0; linind < size(); ++linind) {

    assert(std::abs(other(linind)) > 1.e-16);

    fnc_values[linind] /= other(linind);

    fnc_values_[linind] /= other(linind);

  }

}



template <typename scalartype, class domain>

void function<scalartype, domain>::operator=(const scalartype c) {

  for (int linind = 0; linind < nb_elements_; linind++)

    fnc_values[linind] = c;

  for (int linind = 0; linind < size(); linind++)

    fnc_values_[linind] = c;

}



template <typename scalartype, class domain>

void function<scalartype, domain>::operator+=(const scalartype c) {

  for (int linind = 0; linind < nb_elements_; linind++)

    fnc_values[linind] += c;

  for (int linind = 0; linind < size(); linind++)

    fnc_values_[linind] += c;

}



template <typename scalartype, class domain>

void function<scalartype, domain>::operator-=(const scalartype c) {

  for (int linind = 0; linind < nb_elements_; linind++)

    fnc_values[linind] -= c;

  for (int linind = 0; linind < size(); linind++)

    fnc_values_[linind] -= c;

}



template <typename scalartype, class domain>

void function<scalartype, domain>::operator*=(const scalartype c) {

  for (int linind = 0; linind < nb_elements_; linind++)

    fnc_values[linind] *= c;

  for (int linind = 0; linind < size(); linind++)

    fnc_values_[linind] *= c;

}



template <typename scalartype, class domain>

void function<scalartype, domain>::operator/=(const scalartype c) {

  for (int linind = 0; linind < nb_elements_; linind++)

    fnc_values[linind] /= c;

  for (int linind = 0; linind < size(); linind++)

    fnc_values_[linind] /= c;

}



template <typename scalartype, class domain>
@@ -569,8 +541,8 @@ bool function<scalartype, domain>::operator==(const function<scalartype, domain> 
    // One of the function has not been resetted after the domain was initialized.

    throw std::logic_error("Comparing functions of different sizes.");



  for (int i = 0; i < nb_elements_; ++i)

    if (other(i) != fnc_values[i])

  for (int i = 0; i < size(); ++i)

    if (other(i) != fnc_values_[i])

      return false;



  return true;
@@ -588,7 +560,8 @@ void function<scalartype, domain>::slice(const int sbdm_index, int* subind, 
  subind_2_linind(subind, linind);



  for (int i = 0; i < size_sbdm[sbdm_index]; i++)

    fnc_vals[i] = ScalarCast<new_scalartype>::execute(fnc_values[linind + i * step_sbdm[sbdm_index]]);

    fnc_vals[i] =

        ScalarCast<new_scalartype>::execute(fnc_values_[linind + i * step_sbdm[sbdm_index]]);

}



template <typename scalartype, class domain>
@@ -616,12 +589,12 @@ void function<scalartype, domain>::slice(const int sbdm_index_1, const int sbdm_ 


  for (int j = 0; j < size_sbdm_2; j++) {

    fnc_ptr_left = &fnc_vals[0 + j * size_sbdm_1];

    fnc_ptr_right = &fnc_values[linind + j * step_sbdm_2];

    fnc_ptr_right = &fnc_values_[linind + j * step_sbdm_2];



    for (int i = 0; i < size_sbdm_1; i++)

      fnc_ptr_left[i] = fnc_ptr_right[i * step_sbdm_1];

    //       fnc_vals[i+j*size_sbdm[sbdm_index_1]] = fnc_values[linind + i*step_sbdm[sbdm_index_1] +

    //       j*step_sbdm[sbdm_index_2]];

    //       fnc_vals[i+j*size_sbdm[sbdm_index_1]] = fnc_values_[linind + i*step_sbdm[sbdm_index_1]

    //       + j*step_sbdm[sbdm_index_2]];

  }

}
@@ -637,7 +610,7 @@ void function<scalartype, domain>::distribute(const int sbdm_index, int* subind, 
  subind_2_linind(subind, linind);



  for (int i = 0; i < size_sbdm[sbdm_index]; i++)

    fnc_values[linind + i * step_sbdm[sbdm_index]] = ScalarCast<scalartype>::execute(fnc_vals[i]);

    fnc_values_[linind + i * step_sbdm[sbdm_index]] = ScalarCast<scalartype>::execute(fnc_vals[i]);

}



template <typename scalartype, class domain>
@@ -656,7 +629,7 @@ void function<scalartype, domain>::distribute(const int sbdm_index_1, const int 


  for (int i = 0; i < size_sbdm[sbdm_index_1]; i++)

    for (int j = 0; j < size_sbdm[sbdm_index_2]; j++)

      fnc_values[linind + i * step_sbdm[sbdm_index_1] + j * step_sbdm[sbdm_index_2]] =

      fnc_values_[linind + i * step_sbdm[sbdm_index_1] + j * step_sbdm[sbdm_index_2]] =

          fnc_vals[i + j * size_sbdm[sbdm_index_1]];

}
@@ -675,8 +648,8 @@ void function<scalartype, domain>::print_fingerprint(std::ostream& stream) const 
    stream << "  " << size_sbdm[i];

  stream << "\n";



  stream << "# elements: " << nb_elements_ << "\n";

  stream << "memory: " << nb_elements_ * sizeof(scalartype) / (1024. * 1024.) << " MiB\n";

  stream << "# elements: " << size() << "\n";

  stream << "memory: " << size() * sizeof(scalartype) / (1024. * 1024.) << " MiB\n";

  stream << "****************************************\n" << std::endl;

}
@@ -687,11 +660,11 @@ void function<scalartype, domain>::print_elements(std::ostream& stream) const { 
  stream << "****************************************\n";



  std::vector<int> subind(Nb_sbdms);

  for (int lindex = 0; lindex < nb_elements_; ++lindex) {

  for (int lindex = 0; lindex < size(); ++lindex) {

    linind_2_subind(lindex, subind);

    for (int index : subind)

      stream << index << "\t";

    stream << " \t" << fnc_values[lindex] << "\n";

    stream << " \t" << fnc_values_[lindex] << "\n";

  }



  stream << "****************************************\n" << std::endl;

include/dca/math/inference/richardson_lucy_deconvolution.hpp

+1 −1
Original line number Diff line number Diff line
@@ -78,7 +78,7 @@ private: 
  linalg::Matrix<double, linalg::CPU> u_t_plus_1_;



  func::function<double, OtherDmn> shift_;

  func::function<bool, OtherDmn> is_finished_;

  func::function<char, OtherDmn> is_finished_;

  func::function<double, OtherDmn> error_;

};

include/dca/phys/dca_step/cluster_solver/ctint/structs/interaction_vertices.hpp

+1 −1
Original line number Diff line number Diff line
@@ -147,7 +147,7 @@ void InteractionVertices::initializeFromHamiltonian( 
  // Assume the density-density interaction Hamiltonian function is double counted, i.e.

  // H(b1, b2, r1 - r2) == H(b2, b1, r -r2) and both terms describe the same addendum to the

  // physical Hamiltonian.

  func::function<bool, func::dmn_variadic<Nu, Nu, Rdmn>> already_inserted;

  func::function<char, func::dmn_variadic<Nu, Nu, Rdmn>> already_inserted;

  const int r0 = Rdmn::parameter_type::origin_index();

  for (unsigned short nu1 = 0; nu1 < Nu::dmn_size(); nu1++) {

    for (unsigned short nu2 = 0; nu2 < Nu::dmn_size(); nu2++)

include/dca/phys/dca_step/cluster_solver/ss_ct_hyb/structures/ss_ct_hyb_configuration.hpp

+2 −2
Original line number Diff line number Diff line
@@ -40,7 +40,7 @@ public: 
    return vertices(i);

  }



  bool& get_full_line(int i) {

  char& get_full_line(int i) {

    return has_full_line(i);

  }
@@ -50,7 +50,7 @@ public: 


private:

  func::function<orbital_configuration_type, nu> vertices;

  func::function<bool, nu> has_full_line;

  func::function<char, nu> has_full_line;



  int N_spin_orbitals;

};