47#ifndef MUELU_SMOOVECCOALESCEDROPFACTORY_DEF_HPP
48#define MUELU_SMOOVECCOALESCEDROPFACTORY_DEF_HPP
50#include <Xpetra_CrsGraphFactory.hpp>
51#include <Xpetra_CrsGraph.hpp>
52#include <Xpetra_ImportFactory.hpp>
53#include <Xpetra_MapFactory.hpp>
54#include <Xpetra_Map.hpp>
55#include <Xpetra_Matrix.hpp>
56#include <Xpetra_MultiVectorFactory.hpp>
57#include <Xpetra_MultiVector.hpp>
58#include <Xpetra_StridedMap.hpp>
59#include <Xpetra_VectorFactory.hpp>
60#include <Xpetra_Vector.hpp>
64#include "MueLu_AmalgamationFactory.hpp"
65#include "MueLu_AmalgamationInfo.hpp"
68#include "MueLu_Graph.hpp"
70#include "MueLu_LWGraph.hpp"
73#include "MueLu_PreDropFunctionBaseClass.hpp"
74#include "MueLu_PreDropFunctionConstVal.hpp"
75#include "MueLu_Utilities.hpp"
76#include "MueLu_TopSmootherFactory.hpp"
78#include "MueLu_SmootherFactory.hpp"
81#include <Xpetra_IO.hpp>
96#define poly0thOrderCoef 0
97#define poly1stOrderCoef 1
98#define poly2ndOrderCoef 2
99#define poly3rdOrderCoef 3
100#define poly4thOrderCoef 4
104 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
106 RCP<ParameterList> validParamList = rcp(
new ParameterList());
108#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
111 typedef Teuchos::StringToIntegralParameterEntryValidator<int> validatorType;
112 validParamList->getEntry(
"aggregation: drop scheme").setValidator(
113 rcp(
new validatorType(Teuchos::tuple<std::string>(
"unsupported vector smoothing"),
"aggregation: drop scheme")));
118#undef SET_VALID_ENTRY
120 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A");
121 validParamList->set< RCP<const FactoryBase> >(
"PreSmoother", Teuchos::null,
"Generating factory of the PreSmoother");
122 validParamList->set< RCP<const FactoryBase> >(
"PostSmoother", Teuchos::null,
"Generating factory of the PostSmoother");
124 return validParamList;
127 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
130 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
132 Input(currentLevel,
"A");
134 Input(currentLevel,
"PreSmoother");
136 else if (currentLevel.
IsAvailable(
"PostSmoother")) {
137 Input(currentLevel,
"PostSmoother");
141 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
146 typedef Teuchos::ScalarTraits<SC> STS;
148 if (predrop_ != Teuchos::null)
149 GetOStream(
Parameters0) << predrop_->description();
151 RCP<Matrix> A = Get< RCP<Matrix> >(currentLevel,
"A");
153 const ParameterList & pL = GetParameterList();
155 LO nPDEs = A->GetFixedBlockSize();
157 RCP< MultiVector > testVecs;
158 RCP< MultiVector > nearNull;
161 testVecs = Xpetra::IO<SC,LO,GO,Node>::ReadMultiVector(
"TpetraTVecs.mm", A->getRowMap());
163 size_t numRandom= as<size_t>(pL.get<
int>(
"aggregation: number of random vectors"));
164 testVecs = MultiVectorFactory::Build(A->getRowMap(), numRandom,
true);
167 testVecs->randomize();
168 for (
size_t kk = 0; kk < testVecs->getNumVectors(); kk++ ) {
169 Teuchos::ArrayRCP< Scalar > curVec = testVecs->getDataNonConst(kk);
170 for (
size_t ii = kk; ii < as<size_t>(A->getRowMap()->getLocalNumElements()); ii++ ) curVec[ii] = Teuchos::ScalarTraits<SC>::magnitude(curVec[ii]);
172 nearNull = MultiVectorFactory::Build(A->getRowMap(), nPDEs,
true);
175 for (
size_t kk = 0; kk < nearNull->getNumVectors(); kk++ ) {
176 Teuchos::ArrayRCP< Scalar > curVec = nearNull->getDataNonConst(kk);
177 for (
size_t ii = kk; ii < as<size_t>(A->getRowMap()->getLocalNumElements()); ii += nearNull->getNumVectors() ) curVec[ii] = Teuchos::ScalarTraits<Scalar>::one();
180 RCP< MultiVector > zeroVec_TVecs;
181 RCP< MultiVector > zeroVec_Null;
183 zeroVec_TVecs = MultiVectorFactory::Build(A->getRowMap(), testVecs->getNumVectors(),
true);
184 zeroVec_Null = MultiVectorFactory::Build(A->getRowMap(), nPDEs,
true);
185 zeroVec_TVecs->putScalar(Teuchos::ScalarTraits<Scalar>::zero());
186 zeroVec_Null->putScalar( Teuchos::ScalarTraits<Scalar>::zero());
188 size_t nInvokeSmoother=as<size_t>(pL.get<
int>(
"aggregation: number of times to pre or post smooth"));
190 RCP<SmootherBase> preSmoo = currentLevel.
Get< RCP<SmootherBase> >(
"PreSmoother");
191 for (
size_t ii = 0; ii < nInvokeSmoother; ii++) preSmoo->Apply(*testVecs,*zeroVec_TVecs,
false);
192 for (
size_t ii = 0; ii < nInvokeSmoother; ii++) preSmoo->Apply(*nearNull,*zeroVec_Null,
false);
194 else if (currentLevel.
IsAvailable(
"PostSmoother")) {
195 RCP<SmootherBase> postSmoo = currentLevel.
Get< RCP<SmootherBase> >(
"PostSmoother");
196 for (
size_t ii = 0; ii < nInvokeSmoother; ii++) postSmoo->Apply(*testVecs,*zeroVec_TVecs,
false);
197 for (
size_t ii = 0; ii < nInvokeSmoother; ii++) postSmoo->Apply(*nearNull, *zeroVec_Null,
false);
202 Teuchos::ArrayRCP<Scalar> penaltyPolyCoef(5);
203 Teuchos::ArrayView<const double> inputPolyCoef;
211 if(pL.isParameter(
"aggregation: penalty parameters") && pL.get<Teuchos::Array<double> >(
"aggregation: penalty parameters").size() > 0) {
212 if (pL.get<Teuchos::Array<double> >(
"aggregation: penalty parameters").size() > penaltyPolyCoef.size())
213 TEUCHOS_TEST_FOR_EXCEPTION(
true,
Exceptions::RuntimeError,
"Number of penalty parameters must be " << penaltyPolyCoef.size() <<
" or less");
214 inputPolyCoef = pL.get<Teuchos::Array<double> >(
"aggregation: penalty parameters")();
216 for (
size_t i = 0; i < as<size_t>(inputPolyCoef.size()) ; i++) penaltyPolyCoef[i] = as<Scalar>(inputPolyCoef[i]);
217 for (
size_t i = as<size_t>(inputPolyCoef.size()); i < as<size_t>(penaltyPolyCoef.size()); i++) penaltyPolyCoef[i] = Teuchos::ScalarTraits<Scalar>::zero();
221 RCP<GraphBase> filteredGraph;
222 badGuysCoalesceDrop(*A, penaltyPolyCoef, nPDEs, *testVecs, *nearNull, filteredGraph);
228 FILE* fp = fopen(
"codeOutput",
"w");
229 fprintf(fp,
"%d %d %d\n",(
int) filteredGraph->GetNodeNumVertices(),(
int) filteredGraph->GetNodeNumVertices(),
230 (
int) filteredGraph->GetNodeNumEdges());
231 for (
size_t i = 0; i < filteredGraph->GetNodeNumVertices(); i++) {
232 ArrayView<const LO> inds = filteredGraph->getNeighborVertices(as<LO>(i));
233 for (
size_t j = 0; j < as<size_t>(inds.size()); j++) {
234 fprintf(fp,
"%d %d 1.00e+00\n",(
int) i+1,(
int) inds[j]+1);
241 Set<bool>(currentLevel,
"Filtering", (threshold != STS::zero()));
242 Set(currentLevel,
"Graph", filteredGraph);
243 Set(currentLevel,
"DofsPerNode", 1);
247 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
286 GO numMyNnz = Teuchos::as<GO>(Amat.getLocalNumEntries());
287 size_t nLoc = Amat.getRowMap()->getLocalNumElements();
289 size_t nBlks = nLoc/nPDEs;
290 if (nBlks*nPDEs != nLoc )
293 Teuchos::ArrayRCP<LO> newRowPtr(nBlks+1);
294 Teuchos::ArrayRCP<LO> newCols(numMyNnz);
296 Teuchos::ArrayRCP<LO> bcols(nBlks);
297 Teuchos::ArrayRCP<bool> keepOrNot(nBlks);
301 LO maxNzPerRow = 200;
302 Teuchos::ArrayRCP<Scalar> penalties(maxNzPerRow);
305 Teuchos::ArrayRCP<bool> keepStatus(nBlks,
true);
306 Teuchos::ArrayRCP<LO> bColList(nBlks);
316 Teuchos::ArrayRCP<bool> alreadyOnBColList(nBlks,
false);
321 Teuchos::ArrayRCP<bool> boundaryNodes(nBlks,
false);
324 for (LO i = 0; i < maxNzPerRow; i++)
331 LO nzTotal = 0, numBCols = 0, row = -1, Nbcols, bcol;
335 for (LO i = 0; i < as<LO>(nBlks); i++) {
336 newRowPtr[i+1] = newRowPtr[i];
337 for (LO j = 0; j < nPDEs; j++) {
340 Teuchos::ArrayView<const LocalOrdinal> indices;
341 Teuchos::ArrayView<const Scalar> vals;
343 Amat.getLocalRowView(row, indices, vals);
345 if (indices.size() > maxNzPerRow) {
346 LO oldSize = maxNzPerRow;
347 maxNzPerRow = indices.size() + 100;
348 penalties.resize(as<size_t>(maxNzPerRow),0.0);
349 for (LO k = oldSize; k < maxNzPerRow; k++)
356 badGuysDropfunc(row, indices, vals, testVecs, nPDEs, penalties, nearNull, bcols,keepOrNot,Nbcols,nLoc);
357 for (LO k=0; k < Nbcols; k++) {
362 if (alreadyOnBColList[bcol] ==
false) {
363 bColList[numBCols++] = bcol;
364 alreadyOnBColList[bcol] =
true;
368 if (keepOrNot[k] ==
false) keepStatus[bcol] =
false;
376 if ( numBCols < 2) boundaryNodes[i] =
true;
377 for (LO j=0; j < numBCols; j++) {
379 if (keepStatus[bcol] ==
true) {
380 newCols[nzTotal] = bColList[j];
382 nzTotal = nzTotal + 1;
384 keepStatus[bcol] =
true;
385 alreadyOnBColList[bcol] =
false;
393 Teuchos::ArrayRCP<LO> finalCols(nzTotal);
394 for (LO i = 0; i < nzTotal; i++) finalCols[i] = newCols[i];
399 RCP<const Map> rowMap = Amat.getRowMap();
401 LO nAmalgNodesOnProc = rowMap->getLocalNumElements()/nPDEs;
402 Teuchos::Array<GO> nodalGIDs(nAmalgNodesOnProc);
403 typename Teuchos::ScalarTraits<Scalar>::coordinateType temp;
404 for (
size_t i = 0; i < as<size_t>(nAmalgNodesOnProc); i++ ) {
405 GO gid = rowMap->getGlobalElement(i*nPDEs);
406 temp = ((
typename Teuchos::ScalarTraits<Scalar>::coordinateType) (gid))/((
typename Teuchos::ScalarTraits<Scalar>::coordinateType) (nPDEs));
407 nodalGIDs[i] = as<GO>(floor(temp));
409 GO nAmalgNodesGlobal = rowMap->getGlobalNumElements();
410 GO nBlkGlobal = nAmalgNodesGlobal/nPDEs;
411 if (nBlkGlobal*nPDEs != nAmalgNodesGlobal)
414 Teuchos::RCP<Map> AmalgRowMap = MapFactory::Build(rowMap->lib(), nBlkGlobal,
415 nodalGIDs(),0,rowMap->getComm());
417 filteredGraph = rcp(
new LWGraph(newRowPtr, finalCols, AmalgRowMap, AmalgRowMap,
"thresholded graph of A"));
418 filteredGraph->SetBoundaryNodeMap(boundaryNodes);
422 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
423 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 {
424 using TST=Teuchos::ScalarTraits<Scalar>;
426 LO nLeng = cols.size();
427 typename TST::coordinateType temp;
428 temp = ((
typename TST::coordinateType) (row))/((
typename TST::coordinateType) (nPDEs));
429 LO blkRow = as<LO>(floor(temp));
430 Teuchos::ArrayRCP<Scalar> badGuy( nLeng, 0.0);
431 Teuchos::ArrayRCP<Scalar> subNull(nLeng, 0.0);
445 for (LO i = 0; i < nLeng; i++) keepOrNot[i] =
false;
449 LO rowDof = row - blkRow*nPDEs;
450 Teuchos::ArrayRCP< const Scalar > oneNull = nearNull.getData( as<size_t>(rowDof));
452 for (LO i = 0; i < nLeng; i++) {
453 if ((cols[i] < nLoc ) && (TST::magnitude(vals[i]) != 0.0)) {
454 temp = ((
typename TST::coordinateType) (cols[i]))/((
typename TST::coordinateType) (nPDEs));
455 LO colDof = cols[i] - (as<LO>(floor( temp )))*nPDEs;
456 if (colDof == rowDof) {
457 Bcols[ Nbcols] = (cols[i] - colDof)/nPDEs;
458 subNull[Nbcols] = oneNull[cols[i]];
460 if (cols[i] != row) {
461 Scalar worstRatio = -TST::one();
462 Scalar targetRatio = subNull[Nbcols]/oneNull[row];
464 for (
size_t kk = 0; kk < testVecs.getNumVectors(); kk++ ) {
465 Teuchos::ArrayRCP< const Scalar > curVec = testVecs.getData(kk);
466 actualRatio = curVec[cols[i]]/curVec[row];
467 if (TST::magnitude(actualRatio - targetRatio) > TST::magnitude(worstRatio)) {
468 badGuy[Nbcols] = actualRatio;
469 worstRatio = Teuchos::ScalarTraits<SC>::magnitude(actualRatio - targetRatio);
474 badGuy[ Nbcols] = 1.;
475 keepOrNot[Nbcols] =
true;
486 Bcols[ Nbcols] = (row - rowDof)/nPDEs;
487 subNull[ Nbcols] = 1.;
488 badGuy[ Nbcols] = 1.;
489 keepOrNot[Nbcols] =
true;
494 Scalar currentRP = oneNull[row]*oneNull[row];
495 Scalar currentRTimesBadGuy= oneNull[row]*badGuy[diagInd];
496 Scalar currentScore = penalties[0];
507 LO nKeep = 1, flag = 1, minId;
508 Scalar minFit, minFitRP = 0., minFitRTimesBadGuy = 0.;
509 Scalar newRP, newRTimesBadGuy;
518 for (LO i=0; i < Nbcols; i++) {
519 if (keepOrNot[i] ==
false) {
521 newRP = currentRP + subNull[i]*subNull[i];
522 newRTimesBadGuy= currentRTimesBadGuy + subNull[i]*badGuy[i];
523 Scalar ratio = newRTimesBadGuy/newRP;
526 for (LO k=0; k < Nbcols; k++) {
527 if (keepOrNot[k] ==
true) {
528 Scalar diff = badGuy[k] - ratio*subNull[k];
529 newFit = newFit + diff*diff;
532 if (Teuchos::ScalarTraits<SC>::magnitude(newFit) < Teuchos::ScalarTraits<SC>::magnitude(minFit)) {
536 minFitRTimesBadGuy= newRTimesBadGuy;
538 keepOrNot[i] =
false;
541 if (minId == -1) flag = 0;
543 minFit = sqrt(minFit);
544 Scalar newScore = penalties[nKeep] + minFit;
545 if (Teuchos::ScalarTraits<SC>::magnitude(newScore) < Teuchos::ScalarTraits<SC>::magnitude(currentScore)) {
547 keepOrNot[minId]=
true;
548 currentScore = newScore;
549 currentRP = minFitRP;
550 currentRTimesBadGuy= minFitRTimesBadGuy;
#define SET_VALID_ENTRY(name)
MueLu::DefaultScalar Scalar
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.
Lightweight MueLu representation of a compressed row storage graph.
Class that holds all level-specific information.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
void badGuysCoalesceDrop(const Matrix &Amat, Teuchos::ArrayRCP< Scalar > &dropParams, LO nPDEs, const MultiVector &smoothedTVecs, const MultiVector &smoothedNull, RCP< GraphBase > &filteredGraph) const
Methods to support compatible-relaxation style dropping.
void Build(Level ¤tLevel) const
Build an object with this factory.
SmooVecCoalesceDropFactory()
Constructor.
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
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void DeclareInput(Level ¤tLevel) const
Input.
Namespace for MueLu classes and methods.
@ Parameters0
Print class parameters.