47 #ifndef MUELU_SMOOVECCOALESCEDROPFACTORY_DEF_HPP
48 #define MUELU_SMOOVECCOALESCEDROPFACTORY_DEF_HPP
52 #include <Xpetra_MapFactory.hpp>
53 #include <Xpetra_Map.hpp>
55 #include <Xpetra_MultiVectorFactory.hpp>
57 #include <Xpetra_StridedMap.hpp>
64 #include "MueLu_LWGraph.hpp"
67 #include "MueLu_LWGraph.hpp"
70 #include "MueLu_PreDropFunctionBaseClass.hpp"
86 #define poly0thOrderCoef 0
87 #define poly1stOrderCoef 1
88 #define poly2ndOrderCoef 2
89 #define poly3rdOrderCoef 3
90 #define poly4thOrderCoef 4
94 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
98 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
106 #undef SET_VALID_ENTRY
112 return validParamList;
115 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
117 : predrop_(Teuchos::null) {}
119 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
121 Input(currentLevel,
"A");
123 Input(currentLevel,
"PreSmoother");
125 else if (currentLevel.
IsAvailable(
"PostSmoother")) {
126 Input(currentLevel,
"PostSmoother");
130 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
136 if (predrop_ != Teuchos::null)
137 GetOStream(
Parameters0) << predrop_->description();
139 RCP<Matrix> A = Get<RCP<Matrix> >(currentLevel,
"A");
143 LO nPDEs = A->GetFixedBlockSize();
151 size_t numRandom = as<size_t>(pL.get<
int>(
"aggregation: number of random vectors"));
152 testVecs = MultiVectorFactory::Build(A->getRowMap(), numRandom,
true);
155 testVecs->randomize();
156 for (
size_t kk = 0; kk < testVecs->getNumVectors(); kk++) {
160 nearNull = MultiVectorFactory::Build(A->getRowMap(), nPDEs,
true);
163 for (
size_t kk = 0; kk < nearNull->getNumVectors(); kk++) {
171 zeroVec_TVecs = MultiVectorFactory::Build(A->getRowMap(), testVecs->getNumVectors(),
true);
172 zeroVec_Null = MultiVectorFactory::Build(A->getRowMap(), nPDEs,
true);
176 size_t nInvokeSmoother = as<size_t>(pL.get<
int>(
"aggregation: number of times to pre or post smooth"));
177 if (currentLevel.IsAvailable(
"PreSmoother")) {
179 for (
size_t ii = 0; ii < nInvokeSmoother; ii++) preSmoo->
Apply(*testVecs, *zeroVec_TVecs,
false);
180 for (
size_t ii = 0; ii < nInvokeSmoother; ii++) preSmoo->
Apply(*nearNull, *zeroVec_Null,
false);
181 }
else if (currentLevel.IsAvailable(
"PostSmoother")) {
183 for (
size_t ii = 0; ii < nInvokeSmoother; ii++) postSmoo->
Apply(*testVecs, *zeroVec_TVecs,
false);
184 for (
size_t ii = 0; ii < nInvokeSmoother; ii++) postSmoo->
Apply(*nearNull, *zeroVec_Null,
false);
197 if (pL.isParameter(
"aggregation: penalty parameters") && pL.get<
Teuchos::Array<double> >(
"aggregation: penalty parameters").size() > 0) {
202 for (
size_t i = 0; i < as<size_t>(inputPolyCoef.
size()); i++) penaltyPolyCoef[i] = as<Scalar>(inputPolyCoef[i]);
207 badGuysCoalesceDrop(*A, penaltyPolyCoef, nPDEs, *testVecs, *nearNull, filteredGraph);
212 FILE* fp = fopen(
"codeOutput",
"w");
217 for (
size_t j = 0; j < as<size_t>(inds.size()); j++) {
218 fprintf(fp,
"%d %d 1.00e+00\n", (
int)i + 1, (
int)inds[j] + 1);
225 Set<bool>(currentLevel,
"Filtering", (threshold != STS::zero()));
226 Set(currentLevel,
"Graph", filteredGraph);
227 Set(currentLevel,
"DofsPerNode", 1);
231 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
270 GO numMyNnz = Teuchos::as<GO>(Amat.getLocalNumEntries());
271 size_t nLoc = Amat.getRowMap()->getLocalNumElements();
273 size_t nBlks = nLoc / nPDEs;
274 if (nBlks * nPDEs != nLoc)
277 typename LWGraph::row_type::non_const_type newRowPtr(
"newRowPtr", nBlks + 1);
285 LO maxNzPerRow = 200;
306 Kokkos::deep_copy(boundaryNodes,
false);
308 for (
LO i = 0; i < maxNzPerRow; i++)
312 (penaltyPolyCoef[
poly3rdOrderCoef] * (as<Scalar>(i * i)) * (as<Scalar>(i))) +
313 (penaltyPolyCoef[
poly4thOrderCoef] * (as<Scalar>(i * i)) * (as<Scalar>(i * i)));
315 LO nzTotal = 0, numBCols = 0, row = -1, Nbcols, bcol;
319 for (
LO i = 0; i < as<LO>(nBlks); i++) {
320 newRowPtr[i + 1] = newRowPtr[i];
321 for (
LO j = 0; j < nPDEs; j++) {
327 Amat.getLocalRowView(row, indices, vals);
329 if (indices.
size() > maxNzPerRow) {
330 LO oldSize = maxNzPerRow;
331 maxNzPerRow = indices.
size() + 100;
332 penalties.
resize(as<size_t>(maxNzPerRow), 0.0);
333 for (
LO k = oldSize; k < maxNzPerRow; k++)
337 (penaltyPolyCoef[
poly3rdOrderCoef] * (as<Scalar>(i * i)) * (as<Scalar>(i))) +
338 (penaltyPolyCoef[
poly4thOrderCoef] * (as<Scalar>(i * i)) * (as<Scalar>(i * i)));
340 badGuysDropfunc(row, indices, vals, testVecs, nPDEs, penalties, nearNull, bcols, keepOrNot, Nbcols, nLoc);
341 for (
LO k = 0; k < Nbcols; k++) {
346 if (alreadyOnBColList[bcol] ==
false) {
347 bColList[numBCols++] = bcol;
348 alreadyOnBColList[bcol] =
true;
352 if (keepOrNot[k] ==
false) keepStatus[bcol] =
false;
360 if (numBCols < 2) boundaryNodes[i] =
true;
361 for (
LO j = 0; j < numBCols; j++) {
363 if (keepStatus[bcol] ==
true) {
364 newCols[nzTotal] = bColList[j];
366 nzTotal = nzTotal + 1;
368 keepStatus[bcol] =
true;
369 alreadyOnBColList[bcol] =
false;
377 typename LWGraph::entries_type::non_const_type finalCols(
"finalCols", nzTotal);
378 for (
LO i = 0; i < nzTotal; i++) finalCols(i) = newCols[i];
385 LO nAmalgNodesOnProc = rowMap->getLocalNumElements() / nPDEs;
388 for (
size_t i = 0; i < as<size_t>(nAmalgNodesOnProc); i++) {
389 GO gid = rowMap->getGlobalElement(i * nPDEs);
391 nodalGIDs[i] = as<GO>(floor(temp));
393 GO nAmalgNodesGlobal = rowMap->getGlobalNumElements();
394 GO nBlkGlobal = nAmalgNodesGlobal / nPDEs;
395 if (nBlkGlobal * nPDEs != nAmalgNodesGlobal)
399 nodalGIDs(), 0, rowMap->getComm());
401 filteredGraph =
rcp(
new LWGraph(newRowPtr, finalCols, AmalgRowMap, AmalgRowMap,
"thresholded graph of A"));
402 filteredGraph->SetBoundaryNodeMap(boundaryNodes);
405 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
406 void SmooVecCoalesceDropFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::badGuysDropfunc(
LO row,
const Teuchos::ArrayView<const LocalOrdinal>& cols,
const Teuchos::ArrayView<const Scalar>& vals,
const MultiVector& testVecs,
LO nPDEs,
Teuchos::ArrayRCP<Scalar>& penalties,
const MultiVector& nearNull,
Teuchos::ArrayRCP<LO>& Bcols,
Teuchos::ArrayRCP<bool>& keepOrNot,
LO& Nbcols,
LO nLoc)
const {
410 typename TST::coordinateType temp;
411 temp = ((
typename TST::coordinateType)(row)) / ((
typename TST::coordinateType)(nPDEs));
412 LO blkRow = as<LO>(floor(temp));
427 for (
LO i = 0; i < nLeng; i++) keepOrNot[i] =
false;
431 LO rowDof = row - blkRow * nPDEs;
434 for (
LO i = 0; i < nLeng; i++) {
435 if ((cols[i] < nLoc) && (TST::magnitude(vals[i]) != 0.0)) {
436 temp = ((
typename TST::coordinateType)(cols[i])) / ((
typename TST::coordinateType)(nPDEs));
437 LO colDof = cols[i] - (as<LO>(floor(temp))) * nPDEs;
438 if (colDof == rowDof) {
439 Bcols[Nbcols] = (cols[i] - colDof) / nPDEs;
440 subNull[Nbcols] = oneNull[cols[i]];
442 if (cols[i] != row) {
443 Scalar worstRatio = -TST::one();
444 Scalar targetRatio = subNull[Nbcols] / oneNull[row];
446 for (
size_t kk = 0; kk < testVecs.getNumVectors(); kk++) {
448 actualRatio = curVec[cols[i]] / curVec[row];
449 if (TST::magnitude(actualRatio - targetRatio) > TST::magnitude(worstRatio)) {
450 badGuy[Nbcols] = actualRatio;
456 keepOrNot[Nbcols] =
true;
467 Bcols[Nbcols] = (row - rowDof) / nPDEs;
468 subNull[Nbcols] = 1.;
470 keepOrNot[Nbcols] =
true;
475 Scalar currentRP = oneNull[row] * oneNull[row];
476 Scalar currentRTimesBadGuy = oneNull[row] * badGuy[diagInd];
477 Scalar currentScore = penalties[0];
488 LO nKeep = 1, flag = 1, minId;
489 Scalar minFit, minFitRP = 0., minFitRTimesBadGuy = 0.;
490 Scalar newRP, newRTimesBadGuy;
499 for (
LO i = 0; i < Nbcols; i++) {
500 if (keepOrNot[i] ==
false) {
502 newRP = currentRP + subNull[i] * subNull[i];
503 newRTimesBadGuy = currentRTimesBadGuy + subNull[i] * badGuy[i];
504 Scalar ratio = newRTimesBadGuy / newRP;
507 for (
LO k = 0; k < Nbcols; k++) {
508 if (keepOrNot[k] ==
true) {
509 Scalar diff = badGuy[k] - ratio * subNull[k];
510 newFit = newFit + diff * diff;
517 minFitRTimesBadGuy = newRTimesBadGuy;
519 keepOrNot[i] =
false;
525 minFit = sqrt(minFit);
526 Scalar newScore = penalties[nKeep] + minFit;
529 keepOrNot[minId] =
true;
530 currentScore = newScore;
531 currentRP = minFitRP;
532 currentRTimesBadGuy = minFitRTimesBadGuy;
541 #endif // MUELU_SMOOVECCOALESCEDROPFACTORY_DEF_HPP
void setValidator(RCP< const ParameterEntryValidator > const &validator)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
KOKKOS_INLINE_FUNCTION size_type GetNodeNumVertices() const
Return number of graph vertices.
KOKKOS_INLINE_FUNCTION size_type GetNodeNumEdges() const
Return number of graph edges.
void resize(const size_type n, const T &val=T())
void Build(Level ¤tLevel) const
Build an object with this factory.
Kokkos::View< bool *, memory_space > boundary_nodes_type
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultScalar Scalar
Class that holds all level-specific information.
#define SET_VALID_ENTRY(name)
static RCP< MultiVector > ReadMultiVector(const std::string &fileName, const RCP< const Map > &map)
SmooVecCoalesceDropFactory()
Constructor.
void DeclareInput(Level ¤tLevel) const
Input.
KOKKOS_INLINE_FUNCTION neighbor_vertices_type getNeighborVertices(LO i) const
Return the list of vertices adjacent to the vertex 'v'.
static magnitudeType magnitude(T a)
void badGuysDropfunc(LO row, const Teuchos::ArrayView< const LocalOrdinal > &indices, const Teuchos::ArrayView< const Scalar > &vals, const MultiVector &smoothedTVecs, LO nPDEs, Teuchos::ArrayRCP< Scalar > &penalties, const MultiVector &smoothedNull, Teuchos::ArrayRCP< LO > &Bcols, Teuchos::ArrayRCP< bool > &keepOrNot, LO &Nbcols, LO nLoc) const
Lightweight MueLu representation of a compressed row storage graph.
Exception throws to report errors in the internal logical of the program.
void badGuysCoalesceDrop(const Matrix &Amat, Teuchos::ArrayRCP< Scalar > &dropParams, LO nPDEs, const MultiVector &smoothedTVecs, const MultiVector &smoothedNull, RCP< LWGraph > &filteredGraph) const
Methods to support compatible-relaxation style dropping.
ParameterEntry & getEntry(const std::string &name)
virtual void Apply(MultiVector &x, const MultiVector &rhs, bool InitialGuessIsZero=false) const =0
Apply smoother.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.