Tpetra_Distribution1D.hpp

// @HEADER
// *****************************************************************************
//          Tpetra: Templated Linear Algebra Services Package
//
// Copyright 2008 NTESS and the Tpetra contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
 
//  1D Row-based Distribution class
//  Assumes square matrix
//  Karen Devine, SNL
//
 
#ifndef __TPETRA_DISTRIBUTION1D_HPP
#define __TPETRA_DISTRIBUTION1D_HPP
 
namespace Tpetra {
 
// Forward definition
template <typename gno_t, typename scalar_t>
class DistributionLowerTriangularBlock;
 
template <typename gno_t, typename scalar_t>
class Distribution1D : public Distribution<gno_t, scalar_t> {
  // 1D row-wise distribution of matrix and vector entries
  // Rows and vector entries may be linearly or randomly distributed, or
  // read from a file.
  // Row map and vector map are identical
 
 public:
  using Distribution<gno_t, scalar_t>::me;
  using Distribution<gno_t, scalar_t>::np;
  using Distribution<gno_t, scalar_t>::nrows;
  using Distribution<gno_t, scalar_t>::Mine;
 
  Distribution1D(size_t nrows_,
                 const Teuchos::RCP<const Teuchos::Comm<int> > &comm_,
                 const Teuchos::ParameterList &params)
    : Distribution<gno_t, scalar_t>(nrows_, comm_, params) {
    int npRow = -1;  // Number of processors among which to distribute rows;
                     // Will compute if not set by user
    const Teuchos::ParameterEntry *pe = params.getEntryPtr("nProcessorRows");
    if (pe != NULL) npRow = pe->getValue<int>(&npRow);
 
    TEUCHOS_TEST_FOR_EXCEPTION(npRow != -1 && npRow != np, std::logic_error,
                               " nProcessorRows " << npRow << " must equal"
                                                  << " nProcessors " << np << " for 1D distribution");
 
    if (me == 0) std::cout << "\n 1D Distribution: "
                           << "\n     np     = " << np << std::endl;
  }
 
  // Return whether this rank owns vector entry i.
  virtual bool VecMine(gno_t i) = 0;
 
  // Return whether this rank owns nonzero (i,j)
  // Vector map and row map are the same in 1D distribution.
  inline bool Mine(gno_t i, gno_t j) { return VecMine(i); }
  inline bool Mine(gno_t i, gno_t j, int p) { return VecMine(i); }
};
 
template <typename gno_t, typename scalar_t>
class Distribution1DLinear : public Distribution1D<gno_t, scalar_t> {
 public:
  using Distribution<gno_t, scalar_t>::me;
  using Distribution<gno_t, scalar_t>::np;
  using Distribution<gno_t, scalar_t>::nrows;
 
  Distribution1DLinear(size_t nrows_,
                       const Teuchos::RCP<const Teuchos::Comm<int> > &comm_,
                       const Teuchos::ParameterList &params)
    : Distribution1D<gno_t, scalar_t>(nrows_, comm_, params) {
    gno_t nMyRows = getNumRow(me);
    myFirstRow    = getFirstRow(me);
    myLastRow     = myFirstRow + nMyRows - 1;
  }
 
  inline enum DistributionType DistType() { return OneDLinear; }
 
  inline bool VecMine(gno_t i) { return (i >= myFirstRow && i <= myLastRow); }
 
 private:
  gno_t myFirstRow;
  gno_t myLastRow;
 
  inline size_t getNumRow(int p) {
    return (nrows / np + (int(nrows % np) > p));
  }
 
  inline gno_t getFirstRow(int p) {
    return (p * (nrows / np) + std::min<int>(int(nrows % np), p));
  }
 
  // DistributionLowerTriangularBlock class needs a 1DLinear distribution
  friend class DistributionLowerTriangularBlock<gno_t, scalar_t>;
};
 
template <typename gno_t, typename scalar_t>
class Distribution1DRandom : public Distribution1D<gno_t, scalar_t> {
 public:
  using Distribution<gno_t, scalar_t>::me;
 
  Distribution1DRandom(size_t nrows_,
                       const Teuchos::RCP<const Teuchos::Comm<int> > &comm_,
                       const Teuchos::ParameterList &params)
    : Distribution1D<gno_t, scalar_t>(nrows_, comm_, params) {
    if (me == 0) std::cout << "    randomize = true" << std::endl;
  }
 
  inline enum DistributionType DistType() { return OneDRandom; }
 
  inline bool VecMine(gno_t i) { return (this->HashToProc(i) == me); }
};
 
template <typename gno_t, typename scalar_t>
class Distribution1DVec : public Distribution1D<gno_t, scalar_t> {
  // Distribution of nonzeros is determined by the distribution of the
  // vector entries, as read from a file.
  //
  // Assumptions include:
  // -  Distribution file containing the vector part assignments (N lines)
  //    is provided.  This file is read during the constructor.
  //    Format for an NxN matrix:
  //        line 1 to N:  0-based part assignment of vector entry
 
 public:
  using Distribution<gno_t, scalar_t>::me;
  using Distribution<gno_t, scalar_t>::np;
  using Distribution<gno_t, scalar_t>::comm;
  using Distribution<gno_t, scalar_t>::nrows;
  using Distribution<gno_t, scalar_t>::Mine;
 
  Distribution1DVec(size_t nrows_,
                    const Teuchos::RCP<const Teuchos::Comm<int> > &comm_,
                    const Teuchos::ParameterList &params,
                    std::string &distributionfile)
    : Distribution1D<gno_t, scalar_t>(nrows_, comm_, params) {
    std::ifstream fpin;
    if (me == 0) {
      fpin.open(distributionfile.c_str(), std::ios::in);
      if (!fpin.is_open()) {
        std::cout << "Error:  distributionfile " << distributionfile
                  << " not found" << std::endl;
        exit(-1);
      }
    }
 
    // Read the vector part assignment and broadcast it to all processes.
    // Broadcast in chunks of bcastsize values.
    // TODO:  Make the vector part assignment more scalable instead of
    // TODO:  storing entire vector on every process.
 
    vecpart = new int[nrows];
 
    const int bcastsize = 1000000;
 
    gno_t start = 0;
    int cnt     = 0;
    for (size_t i = 0; i < nrows; i++) {
      if (me == 0) fpin >> vecpart[i];
      cnt++;
      if (cnt == bcastsize || i == nrows - 1) {
        Teuchos::broadcast<int, int>(*comm, 0, cnt, &(vecpart[start]));
        start += cnt;
        cnt = 0;
      }
    }
 
    if (me == 0) fpin.close();
  }
 
  ~Distribution1DVec() { delete[] vecpart; }
 
  inline enum DistributionType DistType() { return OneDVec; }
 
  // Vector distribution was read in.
  inline bool VecMine(gno_t i) { return (vecpart[i] == me); }
 
 protected:
  int *vecpart;  // part assignment of vector entries
};
 
}  // namespace Tpetra
 
#endif