96 FactoryMonitor m(*
this,
"Matrix filtering (segregation)", currentLevel);
99 RCP<Matrix> Ain = Get< RCP<Matrix> >(currentLevel,
"A");
101 const ParameterList& pL = GetParameterList();
102 std::string mapName = pL.get<std::string> (
"map: name");
103 std::string mapFact = pL.get<std::string> (
"map: factory");
105 Teuchos::RCP<const FactoryBase> mapFactory = Teuchos::null;
106 if(mapFact ==
"NoFactory") {
108 }
else if (mapFact !=
"null") {
112 Teuchos::RCP<Map> map = currentLevel.
Get< Teuchos::RCP<Map> >(mapName,mapFactory.get());
115 Teuchos::RCP<Matrix> Aout = MatrixFactory::Build(Ain->getRowMap(), Ain->getGlobalMaxNumRowEntries());
117 size_t numLocalRows = Ain->getLocalNumRows();
118 for(
size_t row=0; row<numLocalRows; row++) {
119 GlobalOrdinal grid = Ain->getRowMap()->getGlobalElement(row);
120 bool isInMap = map->isNodeGlobalElement(grid);
123 Teuchos::ArrayView<const LocalOrdinal> indices;
124 Teuchos::ArrayView<const Scalar> vals;
125 Ain->getLocalRowView(row, indices, vals);
128 Teuchos::ArrayRCP<GlobalOrdinal> indout(indices.size(),Teuchos::ScalarTraits<GlobalOrdinal>::zero());
129 Teuchos::ArrayRCP<Scalar> valout(indices.size(),Teuchos::ScalarTraits<Scalar>::zero());
131 size_t nNonzeros = 0;
132 for(
size_t i=0; i<(size_t)indices.size(); i++) {
133 GlobalOrdinal gcid = Ain->getColMap()->getGlobalElement(indices[i]);
134 bool isInMap2 = map->isNodeGlobalElement(gcid);
136 if (isInMap == isInMap2) {
137 indout [nNonzeros] = gcid;
138 valout [nNonzeros] = vals[i];
142 indout.resize(nNonzeros);
143 valout.resize(nNonzeros);
145 Aout->insertGlobalValues(Ain->getRowMap()->getGlobalElement(row), indout.view(0,indout.size()), valout.view(0,valout.size()));
148 Aout->fillComplete(Ain->getDomainMap(), Ain->getRangeMap());
151 Aout->SetFixedBlockSize(Ain->GetFixedBlockSize());
153 GetOStream(
Statistics0, 0) <<
"Nonzeros in A (input): " << Ain->getGlobalNumEntries() <<
", Nonzeros after filtering A: " << Aout->getGlobalNumEntries() << std::endl;
155 Set(currentLevel,
"A", Aout);
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....