8#ifndef MUELU_ISORROPIAINTERFACE_DEF_HPP_
9#define MUELU_ISORROPIAINTERFACE_DEF_HPP_
13#include <Teuchos_Utils.hpp>
17#include <Xpetra_MapFactory.hpp>
18#include <Xpetra_CrsGraphFactory.hpp>
19#include <Xpetra_BlockedMap.hpp>
20#include <Xpetra_BlockedCrsMatrix.hpp>
22#ifdef HAVE_MUELU_ISORROPIA
23#include <Isorropia_Exception.hpp>
26#ifdef HAVE_MUELU_EPETRA
27#include <Xpetra_EpetraCrsGraph.hpp>
28#include <Epetra_CrsGraph.h>
29#include <Isorropia_EpetraPartitioner.hpp>
32#ifdef HAVE_MUELU_TPETRA
33#include <Xpetra_TpetraCrsGraph.hpp>
40#include "MueLu_Graph.hpp"
41#include "MueLu_AmalgamationInfo.hpp"
42#include "MueLu_Utilities.hpp"
46 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
48 RCP<ParameterList> validParamList = rcp(
new ParameterList());
50 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Factory of the matrix A");
51 validParamList->set< RCP<const FactoryBase> >(
"number of partitions", Teuchos::null,
"Instance of RepartitionHeuristicFactory.");
52 validParamList->set< RCP<const FactoryBase> >(
"UnAmalgamationInfo", Teuchos::null,
"Generating factory of UnAmalgamationInfo");
54 return validParamList;
58 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
60 Input(currentLevel,
"A");
61 Input(currentLevel,
"number of partitions");
62 Input(currentLevel,
"UnAmalgamationInfo");
65 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
68 typedef Xpetra::BlockedMap<LocalOrdinal, GlobalOrdinal, Node> BlockMap;
70 RCP<Matrix> A = Get< RCP<Matrix> >(level,
"A");
71 RCP<AmalgamationInfo> amalInfo = Get< RCP<AmalgamationInfo> >(level,
"UnAmalgamationInfo");
72 GO numParts = Get< int >(level,
"number of partitions");
74 RCP<const Map> rowMap = A->getRowMap();
75 RCP<const Map> colMap = A->getColMap();
77 if (numParts == 1 || numParts == -1) {
79 RCP<Xpetra::Vector<GO, LO, GO, NO> > decomposition = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(rowMap,
true);
80 Set(level,
"AmalgamatedPartition", decomposition);
97 GO indexBase = rowMap->getIndexBase();
100 LO nStridedOffset = 0;
105 if(A->IsView(
"stridedMaps") &&
106 Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap(
"stridedMaps")) != Teuchos::null) {
107 Xpetra::viewLabel_t oldView = A->SwitchToView(
"stridedMaps");
108 RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap());
109 TEUCHOS_TEST_FOR_EXCEPTION(strMap == Teuchos::null,
Exceptions::BadCast,
"MueLu::IsorropiaInterface::Build: cast to strided row map failed.");
110 blockdim = strMap->getFixedBlockSize();
111 offset = strMap->getOffset();
112 blockid = strMap->getStridedBlockId();
114 std::vector<size_t> stridingInfo = strMap->getStridingData();
115 for (
size_t j = 0; j < Teuchos::as<size_t>(blockid); j++)
116 nStridedOffset += stridingInfo[j];
122 oldView = A->SwitchToView(oldView);
124 }
else GetOStream(
Statistics0) <<
"IsorropiaInterface::Build(): no striding information available. Use blockdim=1 with offset=0" << std::endl;
128 RCP<const Map> nodeMap= amalInfo->getNodeRowMap();
129 RCP<const BlockedMap> bnodeMap = Teuchos::rcp_dynamic_cast<const BlockedMap>(nodeMap);
130 if(!bnodeMap.is_null()) nodeMap=bnodeMap->getMap();
132 GetOStream(
Statistics0) <<
"IsorropiaInterface:Build(): nodeMap " << nodeMap->getLocalNumElements() <<
"/" << nodeMap->getGlobalNumElements() <<
" elements" << std::endl;
136 RCP<CrsGraph> crsGraph = CrsGraphFactory::Build(nodeMap, A->getLocalMaxNumRowEntries()*blockdim);
139 for(LO row=0; row<Teuchos::as<LO>(A->getRowMap()->getLocalNumElements()); row++) {
141 GO grid = rowMap->getGlobalElement(row);
146 GO nodeId = (grid - offset - indexBase) / blockdim + indexBase;
148 size_t nnz = A->getNumEntriesInLocalRow(row);
149 Teuchos::ArrayView<const LO> indices;
150 Teuchos::ArrayView<const SC> vals;
151 A->getLocalRowView(row, indices, vals);
153 RCP<std::vector<GO> > cnodeIds = Teuchos::rcp(
new std::vector<GO>);
155 for(LO col=0; col<Teuchos::as<LO>(nnz); col++) {
156 GO gcid = colMap->getGlobalElement(indices[col]);
161 GO cnodeId = (gcid - offset - indexBase) / blockdim + indexBase;
162 cnodeIds->push_back(cnodeId);
167 Teuchos::ArrayRCP<GO> arr_cnodeIds = Teuchos::arcp( cnodeIds );
169 if(arr_cnodeIds.size() > 0 )
170 crsGraph->insertGlobalIndices(nodeId, arr_cnodeIds());
173 crsGraph->fillComplete(nodeMap,nodeMap);
176#ifdef HAVE_MUELU_ISORROPIA
179 Teuchos::ParameterList paramlist;
180 paramlist.set(
"NUM PARTS",
toString(numParts));
188 Teuchos::ParameterList& sublist = paramlist.sublist(
"Zoltan");
189 sublist.set(
"LB_APPROACH",
"PARTITION");
193#ifdef HAVE_MUELU_EPETRA
194 RCP< Xpetra::EpetraCrsGraphT<GO, Node> > epCrsGraph = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsGraphT<GO, Node> >(crsGraph);
195 if(epCrsGraph != Teuchos::null) {
196 RCP< const Epetra_CrsGraph> epetraCrsGraph = epCrsGraph->getEpetra_CrsGraph();
198 RCP<Isorropia::Epetra::Partitioner> isoPart = Teuchos::rcp(
new Isorropia::Epetra::Partitioner(epetraCrsGraph,paramlist));
201 const int* array = NULL;
202 isoPart->extractPartsView(size,array);
204 TEUCHOS_TEST_FOR_EXCEPTION(size != Teuchos::as<int>(nodeMap->getLocalNumElements()),
Exceptions::RuntimeError,
"length of array returned from extractPartsView does not match local length of rowMap");
206 RCP<Xpetra::Vector<GO, LO, GO, NO> > decomposition = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(nodeMap,
false);
207 ArrayRCP<GO> decompEntries = decomposition->getDataNonConst(0);
210 for(
int i = 0; i<size; i++) {
211 decompEntries[i] = Teuchos::as<GO>(array[i]);
214 Set(level,
"AmalgamatedPartition", decomposition);
219#ifdef HAVE_MUELU_TPETRA
220#ifdef HAVE_MUELU_INST_DOUBLE_INT_INT
221 RCP< Xpetra::TpetraCrsGraph<LO, GO, Node> > tpCrsGraph = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsGraph<LO, GO, Node> >(crsGraph);
222 TEUCHOS_TEST_FOR_EXCEPTION(tpCrsGraph != Teuchos::null,
Exceptions::RuntimeError,
"Tpetra is not supported with Isorropia.");
224 TEUCHOS_TEST_FOR_EXCEPTION(
false,
Exceptions::RuntimeError,
"Isorropia is an interface to Zoltan which only has support for LO=GO=int and SC=double.");
232 RCP<Xpetra::Vector<GO, LO, GO, NO> > decomposition = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(rowMap,
true);
233 Set(level,
"AmalgamatedPartition", decomposition);
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
void DeclareInput(Level &level) const
Specifies the data that this class needs, and the factories that generate that data.
void Build(Level &level) const
Build an object with this factory.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Class that holds all level-specific information.
Namespace for MueLu classes and methods.
@ Statistics0
Print statistics that do not involve significant additional computation.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.