Commit 3aa8928e authored by gbalduzz's avatar gbalduzz
Browse files

Store tp channels as a vector.

parent a9f24894

include/dca/phys/dca_analysis/bse_solver/bse_lattice_solver.hpp

+2 −1
Original line number Diff line number Diff line
@@ -401,7 +401,8 @@ void BseLatticeSolver<ParametersType, DcaDataType, ScalarType>::diagonalizeGamma 
#ifndef DCA_ANALYSIS_TEST_WITH_FULL_DIAGONALIZATION

    // Diagonalize the symmetric matrix \sqrt{\chi_0}\Gamma\sqrt{\chi_0}.

    // The origin in momentum space has always index = 0.

    if (parameters.get_four_point_type() == PARTICLE_PARTICLE_UP_DOWN &&

    // TODO: loop over multiple channels.

    if (parameters.get_channel()[0] == PARTICLE_PARTICLE_UP_DOWN &&

        parameters.get_four_point_momentum_transfer_index() == 0 &&

        parameters.get_four_point_frequency_transfer() == 0) {

      diagonalizeGammaChi0Symmetric();

include/dca/phys/dca_data/dca_data.hpp

+10 −54
Original line number Diff line number Diff line
@@ -287,60 +287,11 @@ DcaData<Parameters>::DcaData(Parameters& parameters_ref) 
  // Reserve storage in advance such that we don't have to copy elements when we fill the vector.

  // We want to avoid copies because function's copy ctor does not copy the name (and because copies

  // are expensive).

  G4_.reserve(parameters_.numG4Channels());



  for (auto channel : parameters_.get_channels()) {

    // Allocate memory for G4.

  // Ensure backward compatibility.

  if (parameters_.get_four_point_type() != NONE)

    G4_.emplace_back("G4");



  // Check which four point types to accumulate.

  else {

    if (parameters_.accumulateG4ParticleHoleTransverse())

      G4_.emplace_back("G4_" + toString(PARTICLE_HOLE_TRANSVERSE));



    if (parameters_.accumulateG4ParticleHoleMagnetic())

      G4_.emplace_back("G4_" + toString(PARTICLE_HOLE_MAGNETIC));



    if (parameters_.accumulateG4ParticleHoleCharge())

      G4_.emplace_back("G4_" + toString(PARTICLE_HOLE_CHARGE));



    if (parameters_.accumulateG4ParticleHoleLongitudinalUpUp())

      G4_.emplace_back("G4_" + toString(PARTICLE_HOLE_LONGITUDINAL_UP_UP));



    if (parameters_.accumulateG4ParticleHoleLongitudinalUpDown())

      G4_.emplace_back("G4_" + toString(PARTICLE_HOLE_LONGITUDINAL_UP_DOWN));



    if (parameters_.accumulateG4ParticleParticleUpDown())

      G4_.emplace_back("G4_" + toString(PARTICLE_PARTICLE_UP_DOWN));

  }



    G4_.emplace_back("G4_" + toString(channel));

    // Allocate memory for error on G4.

  if (parameters_.get_error_computation_type() != ErrorComputationType::NONE) {

    G4_err_.reserve(parameters_.numG4Channels());



    if (parameters_.get_four_point_type() != NONE)

      G4_err_.emplace_back("G4-error");



    else {

      if (parameters_.accumulateG4ParticleHoleTransverse())

        G4_err_.emplace_back("G4_" + toString(PARTICLE_HOLE_TRANSVERSE) + "_err");



      if (parameters_.accumulateG4ParticleHoleMagnetic())

        G4_err_.emplace_back("G4_" + toString(PARTICLE_HOLE_MAGNETIC) + "_err");



      if (parameters_.accumulateG4ParticleHoleCharge())

        G4_err_.emplace_back("G4_" + toString(PARTICLE_HOLE_CHARGE) + "_err");



      if (parameters_.accumulateG4ParticleHoleLongitudinalUpUp())

        G4_err_.emplace_back("G4_" + toString(PARTICLE_HOLE_LONGITUDINAL_UP_UP) + "_err");



      if (parameters_.accumulateG4ParticleHoleLongitudinalUpDown())

        G4_err_.emplace_back("G4_" + toString(PARTICLE_HOLE_LONGITUDINAL_UP_DOWN) + "_err");



      if (parameters_.accumulateG4ParticleParticleUpDown())

        G4_err_.emplace_back("G4_" + toString(PARTICLE_PARTICLE_UP_DOWN) + "_err");

    }

    G4_err_.emplace_back("G4_" + toString(channel) + "_err");

  }

}
@@ -400,6 +351,11 @@ void DcaData<Parameters>::read(Reader& reader) { 
  if (parameters_.accumulateG4()) {

    reader.execute(G_k_w);



    // Try to read G4 with a legacy name.

    if (parameters_.get_channels().size() == 1) {

      reader.execute("G4", G4_[0]);

    }



    for (auto& G4_channel : G4_)

      reader.execute(G4_channel);

  }

include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_cluster_solver.hpp

+5 −5
Original line number Diff line number Diff line
@@ -281,7 +281,7 @@ double CtauxClusterSolver<device_t, Parameters, Data>::finalize(dca_info_struct_ 
  }



  if (compute_jack_knife_ && parameters_.accumulateG4()) {

    for (std::size_t channel = 0; channel < parameters_.numG4Channels(); ++channel)

    for (std::size_t channel = 0; channel < data_.get_G4_error().size(); ++channel)

      data_.get_G4_error()[channel] = concurrency_.jackknifeError(data_.get_G4()[channel], true);

  }
@@ -446,7 +446,7 @@ void CtauxClusterSolver<device_t, Parameters, Data>::collect_measurements() { 


    // sum G4

    if (accumulate_g4) {

      for (int channel = 0; channel < parameters_.numG4Channels(); ++channel) {

      for (int channel = 0; channel < data_.get_G4().size(); ++channel) {

        auto& G4 = data_.get_G4()[channel];

        G4 = accumulator_.get_sign_times_G4()[channel];

        if (compute_jack_knife_)

include/dca/phys/dca_step/cluster_solver/exact_diagonalization_advanced/ed_cluster_solver.hpp

+4 −4
Original line number Diff line number Diff line
@@ -172,7 +172,7 @@ void EDClusterSolver<device_t, parameters_type, MOMS_type>::execute() { 


    sp_Greens_function_obj.compute_all_sp_functions_slow(MOMS_imag, MOMS_real, false);



    if (parameters.get_four_point_type() != NONE) {

    if (parameters.accumulateG4()) {

      if (concurrency.id() == concurrency.first()) {

        std::cout << "\n" << dca::util::print_time() << "\n" << std::endl;

      }
@@ -202,7 +202,7 @@ void EDClusterSolver<device_t, parameters_type, MOMS_type>::execute() { 


  sp_Greens_function_obj.compute_all_sp_functions_slow(MOMS_imag, MOMS_real, true);



  if (parameters.get_four_point_type() != NONE) {

  if (parameters.accumulateG4()) {

    if (concurrency.id() == concurrency.first()) {

      std::cout << "\n" << dca::util::print_time() << "\n" << std::endl;

    }
@@ -274,7 +274,7 @@ void EDClusterSolver<device_t, parameters_type, MOMS_type>::write(std::string fi 
    MOMS_imag.write(writer);

    MOMS_real.write(writer);



    if (parameters.get_four_point_type() != NONE) {

    if (parameters.accumulateG4()) {

      std::cout << "\n\n\t\t start writing tp-Greens-function\n\n";

      tp_Greens_function_obj.write(writer);

    }
@@ -295,7 +295,7 @@ void EDClusterSolver<device_t, parameters_type, MOMS_type>::write(std::string fi 
    std::cout << "\n\n\t\t start writing MOMS_real\n\n";

    MOMS_real.write(writer);



    if (parameters.get_four_point_type() != NONE) {

    if (parameters.accumulateG4()) {

      std::cout << "\n\n\t\t start writing tp-Greens-function\n\n";

      tp_Greens_function_obj.write(writer);

    }

include/dca/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator.hpp

+11 −34
Original line number Diff line number Diff line
@@ -136,7 +136,7 @@ protected: 
  float computeM(const std::array<linalg::Matrix<double, linalg::CPU>, 2>& M_pair,

                 const std::array<Configuration, 2>& configs);



  double updateG4(TpGreensFunction& G4_channel);

  double updateG4(int channel_id);



  void inline updateG4Atomic(Complex* G4_ptr, int s_a, int k1_a, int k2_a, int w1_a, int w2_a,

                             int s_b, int k1_b, int k2_b, int w1_b, int w2_b, Real alpha,
@@ -163,6 +163,7 @@ protected: 
  SpGreenFunction G_;



  std::vector<TpGreensFunction> G4_;

  std::vector<FourPointType> channels_;



  func::function<Complex, func::dmn_variadic<BDmn, BDmn, SDmn, KDmn, WTpExtDmn>> G0_;
@@ -184,6 +185,7 @@ TpAccumulator<Parameters, linalg::CPU>::TpAccumulator( 
      thread_id_(thread_id),

      multiple_accumulators_(pars.get_accumulators() > 1),

      beta_(pars.get_beta()),

      channels_(pars.get_channels()),

      extension_index_offset_((WTpExtDmn::dmn_size() - WTpDmn::dmn_size()) / 2),

      n_pos_frqs_(WTpExtPosDmn::dmn_size()),

      G0_M_(n_bands_),
@@ -196,32 +198,8 @@ TpAccumulator<Parameters, linalg::CPU>::TpAccumulator( 
  // Reserve storage in advance such that we don't have to copy elements when we fill the vector.

  // We want to avoid copies because function's copy ctor does not copy the name (and because copies

  // are expensive).

  G4_.reserve(pars.numG4Channels());



  // Ensure backward compatibility.

  if (pars.get_four_point_type() != NONE) {

    G4_.emplace_back("G4_" + toString(pars.get_four_point_type()));

  }



  // Check which four point types to accumulate.

  else {

    if (pars.accumulateG4ParticleHoleTransverse())

      G4_.emplace_back("G4_" + toString(PARTICLE_HOLE_TRANSVERSE));



    if (pars.accumulateG4ParticleHoleMagnetic())

      G4_.emplace_back("G4_" + toString(PARTICLE_HOLE_MAGNETIC));



    if (pars.accumulateG4ParticleHoleCharge())

      G4_.emplace_back("G4_" + toString(PARTICLE_HOLE_CHARGE));



    if (pars.accumulateG4ParticleHoleLongitudinalUpUp())

      G4_.emplace_back("G4_" + toString(PARTICLE_HOLE_LONGITUDINAL_UP_UP));



    if (pars.accumulateG4ParticleHoleLongitudinalUpDown())

      G4_.emplace_back("G4_" + toString(PARTICLE_HOLE_LONGITUDINAL_UP_DOWN));



    if (pars.accumulateG4ParticleParticleUpDown())

      G4_.emplace_back("G4_" + toString(PARTICLE_PARTICLE_UP_DOWN));

  for (auto channel : channels_) {

    G4_.emplace_back("G4_" + toString(channel));

  }

}
@@ -261,8 +239,8 @@ double TpAccumulator<Parameters, linalg::CPU>::accumulate( 
  gflops += computeM(M_pair, configs);

  gflops += computeG();



  for (auto& G4_channel : G4_)

    gflops += updateG4(G4_channel);

  for (int channel_id = 0; channel_id < G4_.size(); ++channel_id)

    gflops += updateG4(channel_id);



  return gflops;

}
@@ -402,7 +380,7 @@ void TpAccumulator<Parameters, linalg::CPU>::getGMultiband(int s, int k1, int k2 
}



template <class Parameters>

double TpAccumulator<Parameters, linalg::CPU>::updateG4(TpGreensFunction& G4) {

double TpAccumulator<Parameters, linalg::CPU>::updateG4(const int channel_id) {

  // G4 is stored with the following band convention:

  // b1 ------------------------ b3

  //        |           |
@@ -429,9 +407,8 @@ double TpAccumulator<Parameters, linalg::CPU>::updateG4(TpGreensFunction& G4) { 
  const auto& exchange_frq = domains::FrequencyExchangeDomain::get_elements();

  const auto& exchange_mom = domains::MomentumExchangeDomain::get_elements();



  // Strip "G4_" prefix from G4 name and convert to FourPointType.

  const std::string channel_str = G4.get_name().substr(3, G4.get_name().size() - 3);

  const FourPointType channel = stringToFourPointType(channel_str);

  auto& G4 = G4_[channel_id];

  auto channel = channels_[channel_id];



  switch (channel) {

    case PARTICLE_HOLE_TRANSVERSE: