46 #ifndef MUELU_FILTEREDAFACTORY_DEF_HPP 47 #define MUELU_FILTEREDAFACTORY_DEF_HPP 49 #include <Xpetra_Matrix.hpp> 50 #include <Xpetra_MatrixFactory.hpp> 54 #include "MueLu_FactoryManager.hpp" 61 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
63 RCP<ParameterList> validParamList = rcp(
new ParameterList());
65 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 69 #undef SET_VALID_ENTRY 71 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A used for filtering");
72 validParamList->set< RCP<const FactoryBase> >(
"Graph", Teuchos::null,
"Generating factory for coalesced filtered graph");
73 validParamList->set< RCP<const FactoryBase> >(
"Filtering", Teuchos::null,
"Generating factory for filtering boolean");
75 return validParamList;
78 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
80 Input(currentLevel,
"A");
81 Input(currentLevel,
"Filtering");
82 Input(currentLevel,
"Graph");
85 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
89 RCP<Matrix> A = Get< RCP<Matrix> >(currentLevel,
"A");
90 if (Get<bool>(currentLevel,
"Filtering") ==
false) {
91 GetOStream(
Runtime0) <<
"Filtered matrix is not being constructed as no filtering is being done" << std::endl;
92 Set(currentLevel,
"A", A);
96 const ParameterList& pL = GetParameterList();
97 bool lumping = pL.get<
bool>(
"filtered matrix: use lumping");
99 GetOStream(
Runtime0) <<
"Lumping dropped entries" << std::endl;
101 RCP<GraphBase> G = Get< RCP<GraphBase> >(currentLevel,
"Graph");
103 RCP<ParameterList> fillCompleteParams(
new ParameterList);
104 fillCompleteParams->set(
"No Nonlocal Changes",
true);
106 RCP<Matrix> filteredA;
107 if (pL.get<
bool>(
"filtered matrix: reuse graph")) {
108 filteredA = MatrixFactory::Build(A->getCrsGraph());
109 filteredA->resumeFill();
111 BuildReuse(*A, *G, lumping, *filteredA);
113 filteredA->fillComplete(fillCompleteParams);
116 filteredA = MatrixFactory::Build(A->getRowMap(), A->getColMap(), A->getNodeMaxNumRowEntries(), Xpetra::StaticProfile);
118 BuildNew(*A, *G, lumping, *filteredA);
120 filteredA->fillComplete(A->getDomainMap(), A->getRangeMap(), fillCompleteParams);
123 filteredA->SetFixedBlockSize(A->GetFixedBlockSize());
125 if (pL.get<
bool>(
"filtered matrix: reuse eigenvalue")) {
130 filteredA->SetMaxEigenvalueEstimate(A->GetMaxEigenvalueEstimate());
133 Set(currentLevel,
"A", filteredA);
142 #define ASSUME_DIRECT_ACCESS_TO_ROW 151 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
154 SC zero = Teuchos::ScalarTraits<SC>::zero();
156 size_t blkSize = A.GetFixedBlockSize();
158 ArrayView<const LO> inds;
159 ArrayView<const SC> valsA;
160 #ifdef ASSUME_DIRECT_ACCESS_TO_ROW 166 Array<char> filter(blkSize * G.
GetImportMap()->getNodeNumElements(), 0);
169 for (
size_t i = 0; i < numGRows; i++) {
172 for (
size_t j = 0; j < as<size_t>(indsG.size()); j++)
173 for (
size_t k = 0; k < blkSize; k++)
174 filter[indsG[j]*blkSize+k] = 1;
176 for (
size_t k = 0; k < blkSize; k++) {
177 LO row = i*blkSize + k;
179 A.getLocalRowView(row, inds, valsA);
181 size_t nnz = inds.size();
185 #ifdef ASSUME_DIRECT_ACCESS_TO_ROW 187 ArrayView<const SC> vals1;
188 filteredA.getLocalRowView(row, inds, vals1);
189 vals = ArrayView<SC>(
const_cast<SC*
>(vals1.getRawPtr()), nnz);
191 memcpy(vals.getRawPtr(), valsA.getRawPtr(), nnz*
sizeof(SC));
193 vals = Array<SC>(valsA);
196 if (lumping ==
false) {
197 for (
size_t j = 0; j < nnz; j++)
198 if (!filter[inds[j]])
205 for (
size_t j = 0; j < nnz; j++) {
206 if (filter[inds[j]]) {
207 if (inds[j] == row) {
214 diagExtra += vals[j];
224 vals[diagIndex] += diagExtra;
227 #ifndef ASSUME_DIRECT_ACCESS_TO_ROW 230 filteredA.replaceLocalValues(row, inds, vals);
235 for (
size_t j = 0; j < as<size_t> (indsG.size()); j++)
236 for (
size_t k = 0; k < blkSize; k++)
237 filter[indsG[j]*blkSize+k] = 0;
241 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
244 SC zero = Teuchos::ScalarTraits<SC>::zero();
246 size_t blkSize = A.GetFixedBlockSize();
248 ArrayView<const LO> indsA;
249 ArrayView<const SC> valsA;
253 Array<char> filter(blkSize * G.
GetImportMap()->getNodeNumElements(), 0);
256 for (
size_t i = 0; i < numGRows; i++) {
259 for (
size_t j = 0; j < as<size_t>(indsG.size()); j++)
260 for (
size_t k = 0; k < blkSize; k++)
261 filter[indsG[j]*blkSize+k] = 1;
263 for (
size_t k = 0; k < blkSize; k++) {
264 LO row = i*blkSize + k;
266 A.getLocalRowView(row, indsA, valsA);
268 size_t nnz = indsA.size();
272 inds.resize(indsA.size());
273 vals.resize(valsA.size());
276 if (lumping ==
false) {
277 for (
size_t j = 0; j < nnz; j++)
278 if (filter[indsA[j]]) {
279 inds[numInds] = indsA[j];
280 vals[numInds] = valsA[j];
288 for (
size_t j = 0; j < nnz; j++) {
289 if (filter[indsA[j]]) {
290 inds[numInds] = indsA[j];
291 vals[numInds] = valsA[j];
294 if (inds[numInds] == row)
300 diagExtra += valsA[j];
309 vals[diagIndex] += diagExtra;
311 inds.resize(numInds);
312 vals.resize(numInds);
316 filteredA.insertLocalValues(row, inds, vals);
320 for (
size_t j = 0; j < as<size_t> (indsG.size()); j++)
321 for (
size_t k = 0; k < blkSize; k++)
322 filter[indsG[j]*blkSize+k] = 0;
328 #endif // MUELU_FILTEREDAFACTORY_DEF_HPP RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
virtual size_t GetNodeNumVertices() const =0
Return number of vertices owned by the calling node.
Timer to be used in factories. Similar to Monitor but with additional timers.
One-liner description of what is happening.
Namespace for MueLu classes and methods.
void Build(Level ¤tLevel) const
Build method.
void BuildNew(const Matrix &A, const GraphBase &G, const bool lumping, Matrix &filteredA) const
Class that holds all level-specific information.
MueLu representation of a graph.
#define SET_VALID_ENTRY(name)
virtual const RCP< const Map > GetImportMap() const =0
void BuildReuse(const Matrix &A, const GraphBase &G, const bool lumping, Matrix &filteredA) const
void DeclareInput(Level ¤tLevel) const
Input.
virtual Teuchos::ArrayView< const LocalOrdinal > getNeighborVertices(LocalOrdinal v) const =0
Return the list of vertices adjacent to the vertex 'v'.