Ifpack2_Details_UserPartitioner_def.hpp

// @HEADER
// *****************************************************************************
//       Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
//
// Copyright 2009 NTESS and the Ifpack2 contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER

#ifndef IFPACK2_USER_PARTITIONER_DEF_HPP
#define IFPACK2_USER_PARTITIONER_DEF_HPP

#include <Ifpack2_ConfigDefs.hpp>
#include <Ifpack2_Details_UserPartitioner_decl.hpp>
#include <Ifpack2_OverlappingPartitioner.hpp>

namespace Ifpack2 {
namespace Details {

template <class GraphType>
UserPartitioner<GraphType>::
    UserPartitioner(const Teuchos::RCP<const row_graph_type>& graph)
  : OverlappingPartitioner<GraphType>(graph)
  , userProvidedParts_(false)
  , userProvidedMap_(false) {}

template <class GraphType>
UserPartitioner<GraphType>::~UserPartitioner() {}

template <class GraphType>
void UserPartitioner<GraphType>::
    setPartitionParameters(Teuchos::ParameterList& List) {
  int partCount = 0;
  bool globalProvidedParts;
  userProvidedParts_  = List.isParameter("partitioner: parts");
  globalProvidedParts = List.isParameter("partitioner: global ID parts");
  if (userProvidedParts_) partCount++;
  if (globalProvidedParts) {
    partCount++;
    userProvidedParts_ = true;
  }

  userProvidedMap_ = List.isParameter("partitioner: map");
  if (userProvidedMap_) partCount++;
  if (partCount > 1) {
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
                               "Ifpack2::UserPartitioner::setPartitionParameters: "
                               "you may specify only one of \"partitioner: parts\", \"partitioner: global ID parts\" and \"partitioner: map\","
                               " not two of them");
  }
  if (userProvidedMap_) {
    map_ = List.get("partitioner: map", map_);
    TEUCHOS_TEST_FOR_EXCEPTION(
        map_.is_null(), std::runtime_error,
        "Ifpack2::UserPartitioner::"
        "setPartitionParameters: map_ is null.");
  }
  if (userProvidedParts_) {
    typedef Teuchos::Array<typename Teuchos::ArrayRCP<typename GraphType::local_ordinal_type>> parts_type;
    typedef Teuchos::Array<typename Teuchos::ArrayRCP<typename GraphType::global_ordinal_type>> gparts_type;
    // FIXME JJH 9-Dec-2015 - This is a potentially expensive deep copy.  Use an ArrayRCP<ArrayRCP<int>> instead.

    // convert glboal ids to local ids before storing in Parts_
    if (globalProvidedParts) {
      gparts_type gparts;
      gparts = List.get<gparts_type>("partitioner: global ID parts");
      typedef Teuchos::RCP<Tpetra::Map<typename GraphType::local_ordinal_type, typename GraphType::global_ordinal_type, typename GraphType::node_type> const> map_type;
      map_type OverlapMap = List.get<map_type>("OverlapRowMap");
      this->Parts_.resize(gparts.size());
      for (int ii = 0; ii < (int)gparts.size(); ii++) {
        this->Parts_[ii].resize(gparts[ii].size());
        for (int jj = 0; jj < (int)gparts[ii].size(); jj++) {
          local_ordinal_type itmp = (int)OverlapMap->getLocalElement(gparts[ii][jj]);
          if (itmp == Teuchos::OrdinalTraits<local_ordinal_type>::invalid()) printf("global id %d (%d,%d) not found\n", (int)gparts[ii][jj], (int)ii, (int)jj);
          TEUCHOS_TEST_FOR_EXCEPTION(itmp == Teuchos::OrdinalTraits<local_ordinal_type>::invalid(), std::runtime_error, " \"partitioner: global ID parts\" requires that all global IDs within a block reside on the MPI rank owning the block. This can be done using overlapping Schwarz with a BLOCK_RELAXATION subdomain solver making sure that \"schwarz: overlap level\" is sufficient. Note: zeroed out Dirichlet columns will never be included in overlap parts of domains.");
          this->Parts_[ii][jj] = itmp;
        }
      }
      List.remove("partitioner: global ID parts");  // JJH I observed that iterating through the ParameterList is much more
                                                    // expensive when "partitioner: collection" is a Parameter object.
      // Thus, I remove it once it's fetched from the list.
    } else {
      this->Parts_ = List.get<parts_type>("partitioner: parts");
      List.remove("partitioner: parts");  // JJH I observed that iterating through the ParameterList is much more
                                          // expensive when "partitioner: collection" is a Parameter object.
                                          // Thus, I remove it once it's fetched from the list.
    }
  }
}

template <class GraphType>
void UserPartitioner<GraphType>::computePartitions() {
  if (userProvidedParts_) {
    this->NumLocalParts_ = this->Parts_.size();
    // The user has explicitly defined Parts_, so there is no need to for Partition_.
    this->Partition_.resize(0);
  } else {
    TEUCHOS_TEST_FOR_EXCEPTION(
        map_.is_null(), std::logic_error,
        "Ifpack2::UserPartitioner::"
        "computePartitions: map_ is null.");
    const size_t localNumRows = this->Graph_->getLocalNumRows();
    for (size_t ii = 0; ii < localNumRows; ++ii) {
      this->Partition_[ii] = map_[ii];
    }
  }
  // IGNORING ALL THE SINGLETON STUFF IN IFPACK_USERPARTITIONER.CPP BY
  // ORDERS OF CHRIS SIEFERT
}

}  // namespace Details
}  // namespace Ifpack2

#endif  // IFPACK2_USERPARTITIONER_DEF_HPP