HDFWalkerInput_0_4.cpp

Go to the documentation of this file.

//////////////////////////////////////////////////////////////////////////////////////
// This file is distributed under the University of Illinois/NCSA Open Source License.
// See LICENSE file in top directory for details.
//
// Copyright (c) 2016 Jeongnim Kim and QMCPACK developers.
//
// File developed by: Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
//                    Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
//                    Cynthia Gu, zg1@ornl.gov, Oak Ridge National Laboratory
//                    Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
//
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
//////////////////////////////////////////////////////////////////////////////////////


#include "HDFWalkerInput_0_4.h"
#include "hdf/hdf_archive.h"
#include "mpi/mpi_datatype.h"
#include "mpi/collectives.h"
#include "Utilities/FairDivide.h"

#include <array>

namespace qmcplusplus
{
HDFWalkerInput_0_4::HDFWalkerInput_0_4(WalkerConfigurations& wc_list,
                                       size_t num_ptcls,
                                       Communicate* c,
                                       const HDFVersion& v)
    : wc_list_(wc_list), num_ptcls_(num_ptcls), myComm(c), cur_version(0, 4)
{
  i_info.version = v;
}

HDFWalkerInput_0_4::~HDFWalkerInput_0_4()
{
  //if(h_plist != H5P_DEFAULT) H5Pclose(h_plist);
}

void HDFWalkerInput_0_4::checkOptions(xmlNodePtr cur)
{
  i_info.reset();
  std::string froot, cfile;
  std::string collected("no");
  OhmmsAttributeSet pAttrib;
  pAttrib.add(i_info.nprocs, "nprocs");
  pAttrib.add(i_info.rank, "node");
  pAttrib.add(cfile, "href");
  pAttrib.add(cfile, "file");
  pAttrib.add(froot, "fileroot");
  pAttrib.add(collected, "collected");
  pAttrib.put(cur);
  if (froot.empty())
    return;
  if (int ext = froot.find(hdf::config_ext); ext < froot.size())
  {
    //remove extension
    froot.erase(froot.begin() + ext, froot.end());
  }
  //file generated by a serial run is automatically collected
  i_info.collected = ((collected == "yes") || (i_info.nprocs == 1));
  if (i_info.collected)
    FileStack.push(froot);
  else
  {
    if (i_info.nprocs > 1) //need to process multiple files
    {
      int nprocs_now = myComm->size();
      if (i_info.nprocs > nprocs_now) //using less nodes
      {
        int np  = i_info.nprocs / nprocs_now;
        int pid = myComm->rank(), ip = 0;
        while (pid < i_info.nprocs && ip < np)
        {
          // 2 chars for '.p', 11 chars for pid, 1 char for \0
          std::array<char, 14> h5name;
          if (std::snprintf(h5name.data(), h5name.size(), ".p%03d", pid++) < 0)
            throw std::runtime_error("Error generating filename");
          auto fname = std::filesystem::path(froot).concat(h5name.data());
          FileStack.push(std::move(fname));
          pid += nprocs_now;
        }
      }
      else
      {
        int pid = myComm->rank() % i_info.nprocs;
        // 2 chars for '.p', 11 chars for pid, 1 char for \0
        std::array<char, 14> h5name;
        if (std::snprintf(h5name.data(), h5name.size(), ".p%03d", pid) < 0)
          throw std::runtime_error("Error generating filename");
        auto fname = std::filesystem::path(froot).concat(h5name.data());
        FileStack.push(std::move(fname));
      }
    }
    else
    {
      FileStack.push(froot);
    }
  }
}

bool HDFWalkerInput_0_4::put(xmlNodePtr cur)
{
  checkOptions(cur);
  if (FileStack.empty())
  {
    app_error() << "  No valid input hdf5 is found." << std::endl;
    return false;
  }
  bool success = true;
  while (FileStack.size())
  {
    FileName_noext = FileStack.top();
    FileStack.pop();
    auto h5name = std::filesystem::path(FileName_noext).concat(hdf::config_ext);
    //success |= read_hdf5_scatter(h5name);
#ifdef ENABLE_PHDF5
    success |= read_phdf5(h5name);
#else
    success |= read_hdf5(h5name);
#endif
  }
  return success;
}

bool HDFWalkerInput_0_4::read_hdf5(const std::filesystem::path& h5name)
{
  size_t nw_in = 0;

  hdf_archive hin(myComm, false); //everone reads this
  bool success = hin.open(h5name, H5F_ACC_RDONLY);
  //check if hdf and xml versions can work together
  HDFVersion aversion;

  hin.read(aversion, hdf::version);
  if (!(aversion < i_info.version))
  {
    int found_group = hin.is_group(hdf::main_state);
    hin.push(hdf::main_state);
    hin.read(nw_in, hdf::num_walkers);
  }
  else
  {
    app_error() << " Mismatched version. xml = " << i_info.version << " hdf = " << aversion << std::endl;
  }

  if (nw_in == 0)
  {
    app_error() << " No walkers in " << h5name << std::endl;
    return false;
  }

  using Buffer_t = std::vector<QMCTraits::RealType>;
  std::array<size_t, 3> dims{nw_in, num_ptcls_, OHMMS_DIM};
  Buffer_t posin(dims[0] * dims[1] * dims[2]);
  hin.readSlabReshaped(posin, dims, hdf::walkers);
  std::vector<QMCTraits::FullPrecRealType> weights_in(nw_in);
  const bool has_weights = hin.readEntry(weights_in, hdf::walker_weights);

  std::vector<int> woffsets;
  hin.read(woffsets, "walker_partition");

  int np1 = myComm->size() + 1;
  if (woffsets.size() != np1)
  {
    woffsets.resize(myComm->size() + 1, 0);
    FairDivideLow(nw_in, myComm->size(), woffsets);
  }

  app_log() << " HDFWalkerInput_0_4::put getting " << dims[0] << " walkers " << posin.size() << std::endl;
  nw_in = woffsets[myComm->rank() + 1] - woffsets[myComm->rank()];
  {
    const int nitems    = num_ptcls_ * OHMMS_DIM;
    const int curWalker = wc_list_.getActiveWalkers();
    wc_list_.createWalkers(nw_in, num_ptcls_);

    auto it = posin.begin() + woffsets[myComm->rank()] * nitems;
    for (int i = 0; i < nw_in; ++i, it += nitems)
      copy(it, it + nitems, get_first_address(wc_list_[i + curWalker]->R));
    if (has_weights)
    {
      const auto woffset = woffsets[myComm->rank()];
      for (int i = 0; i < nw_in; ++i)
        wc_list_[i + curWalker]->Weight = weights_in[i + woffset];
    }
  }

  return true;
}

bool HDFWalkerInput_0_4::read_hdf5_scatter(const std::filesystem::path& h5name)
{
  size_t nw_in = 0;

  if (myComm->rank() == 0)
  {
    hdf_archive hin(myComm); //everone reads this
    bool success = hin.open(h5name, H5F_ACC_RDONLY);
    //check if hdf and xml versions can work together
    HDFVersion aversion;

    hin.read(aversion, hdf::version);
    if (!(aversion < i_info.version))
    {
      int found_group = hin.is_group(hdf::main_state);
      hin.push(hdf::main_state);
      hin.read(nw_in, hdf::num_walkers);
    }
    else
    {
      app_error() << " Mismatched version. xml = " << i_info.version << " hdf = " << aversion << std::endl;
    }
  }

  myComm->barrier();
  mpi::bcast(*myComm, nw_in);

  if (nw_in == 0)
  {
    app_error() << " No walkers in " << h5name << std::endl;
    return false;
  }

  using Buffer_t = std::vector<QMCTraits::RealType>;

  const int np1 = myComm->size() + 1;
  std::vector<int> woffsets_weights(np1, 0);
  FairDivideLow(nw_in, myComm->size(), woffsets_weights);

  std::vector<int> counts_weights(myComm->size());
  for (int i = 0; i < counts_weights.size(); ++i)
    counts_weights[i] = woffsets_weights[i + 1] - woffsets_weights[i];

  // walker counts and offsets for electron coordinates
  std::vector<int> woffsets(np1, 0);
  const int nitems = num_ptcls_ * OHMMS_DIM;
  for (int i = 0; i < woffsets.size(); ++i)
    woffsets[i] = nitems * woffsets_weights[i];
  std::vector<int> counts(myComm->size());
  for (int i = 0; i < counts.size(); ++i)
    counts[i] = woffsets[i + 1] - woffsets[i];

  std::array<size_t, 3> dims{nw_in, num_ptcls_, OHMMS_DIM};
  Buffer_t posin(nw_in * nitems);
  std::vector<QMCTraits::FullPrecRealType> weights_in(nw_in);
  bool has_weights{false};
  if (myComm->rank() == 0)
  {
    hdf_archive hin(myComm);
    bool success = hin.open(h5name, H5F_ACC_RDONLY);
    hin.push(hdf::main_state);
    hin.readSlabReshaped(posin, dims, hdf::walkers);
    has_weights = hin.readEntry(weights_in, hdf::walker_weights);
  }

  Buffer_t posout(counts[myComm->rank()]);
  std::vector<QMCTraits::FullPrecRealType> weights_out(counts[myComm->rank()]);
  mpi::scatterv(*myComm, posin, posout, counts, woffsets);

  mpi::bcast(*myComm, has_weights);
  if (has_weights)
    mpi::scatterv(*myComm, weights_in, weights_out, counts_weights, woffsets_weights);

  const size_t nw_loc = woffsets[myComm->rank() + 1] - woffsets[myComm->rank()];
  const int curWalker = wc_list_.getActiveWalkers();
  wc_list_.createWalkers(nw_loc, num_ptcls_);

  auto it = posout.begin();
  for (int i = 0; i < nw_loc; ++i, it += nitems)
    std::copy(it, it + nitems, get_first_address(wc_list_[i + curWalker]->R));
  if (has_weights)
    for (int i = 0; i < nw_in; ++i)
      wc_list_[i + curWalker]->Weight = weights_out[i];
  return true;
}

bool HDFWalkerInput_0_4::read_phdf5(const std::filesystem::path& h5name)
{
  size_t nw_in = 0;
  std::vector<int> woffsets;
  int woffsets_size = 0;
  bool success      = false;

  { // handle small dataset with master rank
    hdf_archive hin(myComm, false);
    if (myComm->rank() == 0)
    {
      success = hin.open(h5name, H5F_ACC_RDONLY);
      //check if hdf and xml versions can work together
      HDFVersion aversion;

      hin.read(aversion, hdf::version);
      if (!(aversion < i_info.version))
      {
        int found_group = hin.is_group(hdf::main_state);
        hin.push(hdf::main_state);
        hin.read(nw_in, hdf::num_walkers);
        if (nw_in == 0)
        {
          app_error() << " No walkers in " << h5name << std::endl;
          success = false;
        }
      }
      else
      {
        app_error() << " Mismatched version. xml = " << i_info.version << " hdf = " << aversion << std::endl;
        success = false;
      }
    }
    mpi::bcast(*myComm, success);
    if (!success)
      return false;

    // load woffsets by master
    // can not read collectively since the size may differ from Nranks+1.
    if (myComm->rank() == 0)
    {
      hin.read(woffsets, "walker_partition");
      woffsets_size = woffsets.size();
      assert(woffsets[woffsets_size - 1] == nw_in);
    }

    mpi::bcast(*myComm, woffsets_size);
    woffsets.resize(woffsets_size);
    mpi::bcast(*myComm, woffsets.data(), woffsets_size);
    nw_in = woffsets[woffsets_size - 1];
  }

  hdf_archive hin(myComm, true); //everone reads this
  success         = hin.open(h5name, H5F_ACC_RDONLY);
  int found_group = hin.is_group(hdf::main_state);
  hin.push(hdf::main_state);

  using Buffer_t = std::vector<QMCTraits::RealType>;
  std::array<size_t, 3> dims{nw_in, num_ptcls_, OHMMS_DIM};
  std::array<size_t, 1> dims_w{nw_in};

  if (woffsets.size() != myComm->size() + 1)
  {
    woffsets.resize(myComm->size() + 1, 0);
    FairDivideLow(nw_in, myComm->size(), woffsets);
  }

  const size_t nw_loc = woffsets[myComm->rank() + 1] - woffsets[myComm->rank()];

  std::array<size_t, 3> counts{nw_loc, num_ptcls_, OHMMS_DIM};
  std::array<size_t, 1> counts_w{nw_loc};
  std::array<size_t, 3> offsets{static_cast<size_t>(woffsets[myComm->rank()]), 0, 0};
  std::array<size_t, 1> offsets_w{static_cast<size_t>(woffsets[myComm->rank()])};
  Buffer_t posin(nw_loc * dims[1] * dims[2]);
  std::vector<QMCTraits::FullPrecRealType> weights_in(nw_loc);

  hyperslab_proxy<Buffer_t, 3> slab(posin, dims, counts, offsets);
  hin.read(slab, hdf::walkers);

  hyperslab_proxy<std::vector<QMCTraits::FullPrecRealType>, 1> slab_w(weights_in, dims_w, counts_w, offsets_w);
  const bool has_weights = hin.readEntry(slab_w, hdf::walker_weights);

  app_log() << " HDFWalkerInput_0_4::put getting " << dims[0] << " walkers " << posin.size() << std::endl;
  nw_in = woffsets[myComm->rank() + 1] - woffsets[myComm->rank()];
  {
    const int nitems    = num_ptcls_ * OHMMS_DIM;
    const int curWalker = wc_list_.getActiveWalkers();
    wc_list_.createWalkers(nw_in, num_ptcls_);
    auto it = posin.begin();
    for (int i = 0; i < nw_in; ++i, it += nitems)
      copy(it, it + nitems, get_first_address(wc_list_[i + curWalker]->R));
    if (has_weights)
      for (int i = 0; i < nw_in; ++i)
        wc_list_[i + curWalker]->Weight = weights_in[i];
  }
  return true;
}

} // namespace qmcplusplus