46 #ifndef MUELU_TENTATIVEPFACTORY_DEF_HPP 47 #define MUELU_TENTATIVEPFACTORY_DEF_HPP 49 #include <Xpetra_MapFactory.hpp> 50 #include <Xpetra_Map.hpp> 51 #include <Xpetra_CrsMatrix.hpp> 52 #include <Xpetra_Matrix.hpp> 53 #include <Xpetra_MultiVector.hpp> 54 #include <Xpetra_MultiVectorFactory.hpp> 55 #include <Xpetra_VectorFactory.hpp> 56 #include <Xpetra_Import.hpp> 57 #include <Xpetra_ImportFactory.hpp> 58 #include <Xpetra_CrsMatrixWrap.hpp> 59 #include <Xpetra_StridedMap.hpp> 60 #include <Xpetra_StridedMapFactory.hpp> 64 #include "MueLu_Aggregates.hpp" 65 #include "MueLu_AmalgamationFactory.hpp" 66 #include "MueLu_AmalgamationInfo.hpp" 67 #include "MueLu_CoarseMapFactory.hpp" 68 #include "MueLu_NullspaceFactory.hpp" 69 #include "MueLu_PerfUtils.hpp" 71 #include "MueLu_Utilities.hpp" 75 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
82 validParamList->
set<
RCP<const FactoryBase> >(
"UnAmalgamationInfo", Teuchos::null,
"Generating factory of UnAmalgamationInfo");
84 return validParamList;
87 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
89 Input(fineLevel,
"A");
90 Input(fineLevel,
"Aggregates");
91 Input(fineLevel,
"Nullspace");
92 Input(fineLevel,
"UnAmalgamationInfo");
93 Input(fineLevel,
"CoarseMap");
96 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
98 return BuildP(fineLevel, coarseLevel);
101 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
105 RCP<Matrix> A = Get< RCP<Matrix> > (fineLevel,
"A");
106 RCP<Aggregates> aggregates = Get< RCP<Aggregates> > (fineLevel,
"Aggregates");
108 RCP<MultiVector> fineNullspace = Get< RCP<MultiVector> > (fineLevel,
"Nullspace");
109 RCP<const Map> coarseMap = Get< RCP<const Map> > (fineLevel,
"CoarseMap");
113 if (!aggregates->AggregatesCrossProcessors())
114 BuildPuncoupled(A, aggregates, amalgInfo, fineNullspace, coarseMap, Ptentative, coarseNullspace);
116 BuildPcoupled (A, aggregates, amalgInfo, fineNullspace, coarseMap, Ptentative, coarseNullspace);
126 if (A->IsView(
"stridedMaps") ==
true)
127 Ptentative->CreateView(
"stridedMaps", A->getRowMap(
"stridedMaps"), coarseMap);
129 Ptentative->CreateView(
"stridedMaps", Ptentative->getRangeMap(), coarseMap);
131 Set(coarseLevel,
"Nullspace", coarseNullspace);
132 Set(coarseLevel,
"P", Ptentative);
136 params->
set(
"printLoadBalancingInfo",
true);
141 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
148 const size_t numRows = rowMap->getNodeNumElements();
151 typedef typename STS::magnitudeType Magnitude;
152 const SC zero = STS::zero();
153 const SC one = STS::one();
157 const size_t NSDim = fineNullspace->getNumVectors();
162 bool goodMap = isGoodMap(*rowMap, *colMap);
173 GetOStream(
Runtime1) <<
"Column map is consistent with the row map, good." << std::endl;
177 GetOStream(
Warnings0) <<
"Column map is not consistent with the row map\n" 178 <<
"using GO->LO conversion with performance penalty" << std::endl;
181 coarseNullspace = MultiVectorFactory::Build(coarseMap, NSDim);
186 for (
size_t i = 0; i < NSDim; i++) {
187 fineNS[i] = fineNullspace->getData(i);
188 if (coarseMap->getNodeNumElements() > 0)
189 coarseNS[i] = coarseNullspace->getDataNonConst(i);
192 size_t nnzEstimate = numRows * NSDim;
195 Ptentative =
rcp(
new CrsMatrixWrap(rowMap, coarseMap, 0, Xpetra::StaticProfile));
196 RCP<CrsMatrix> PtentCrs = rcp_dynamic_cast<CrsMatrixWrap>(Ptentative)->getCrsMatrix();
202 PtentCrs->allocateAllValues(nnzEstimate, iaPtent, jaPtent, valPtent);
209 for (
size_t i = 1; i <= numRows; i++)
210 ia[i] = ia[i-1] + NSDim;
212 for (
size_t j = 0; j < nnzEstimate; j++) {
217 for (GO agg = 0; agg < numAggs; agg++) {
218 LO aggSize = aggStart[agg+1] - aggStart[agg];
220 Xpetra::global_size_t offset = agg*NSDim;
227 for (
size_t j = 0; j < NSDim; j++)
228 for (LO k = 0; k < aggSize; k++)
229 localQR(k,j) = fineNS[j][aggToRowMapLO[aggStart[agg]+k]];
231 for (
size_t j = 0; j < NSDim; j++)
232 for (LO k = 0; k < aggSize; k++)
233 localQR(k,j) = fineNS[j][rowMap->getLocalElement(aggToRowMapGO[aggStart[agg]+k])];
237 for (
size_t j = 0; j < NSDim; j++) {
238 bool bIsZeroNSColumn =
true;
240 for (LO k = 0; k < aggSize; k++)
241 if (localQR(k,j) != zero)
242 bIsZeroNSColumn =
false;
245 "MueLu::TentativePFactory::MakeTentative: fine level NS part has a zero column");
250 if (aggSize >= Teuchos::as<LO>(NSDim)) {
254 Magnitude norm = STS::magnitude(zero);
255 for (
size_t k = 0; k < Teuchos::as<size_t>(aggSize); k++)
256 norm += STS::magnitude(localQR(k,0)*localQR(k,0));
260 coarseNS[0][offset] = norm;
263 for (LO i = 0; i < aggSize; i++)
264 localQR(i,0) /= norm;
272 for (
size_t j = 0; j < NSDim; j++)
273 for (
size_t k = 0; k <= j; k++)
274 coarseNS[j][offset+k] = localQR(k,j);
280 for (
size_t j = 0; j < NSDim; j++)
281 for (
size_t i = 0; i < Teuchos::as<size_t>(aggSize); i++)
282 localQR(i,j) = (*qFactor)(i,j);
313 for (
size_t j = 0; j < NSDim; j++)
314 for (
size_t k = 0; k < NSDim; k++)
315 if (k < as<size_t>(aggSize))
316 coarseNS[j][offset+k] = localQR(k,j);
318 coarseNS[j][offset+k] = (k == j ? one : zero);
321 for (
size_t i = 0; i < as<size_t>(aggSize); i++)
322 for (
size_t j = 0; j < NSDim; j++)
323 localQR(i,j) = (j == i ? one : zero);
328 for (LO j = 0; j < aggSize; j++) {
329 LO localRow = (goodMap ? aggToRowMapLO[aggStart[agg]+j] : rowMap->getLocalElement(aggToRowMapGO[aggStart[agg]+j]));
331 size_t rowStart = ia[localRow];
332 for (
size_t k = 0, lnnz = 0; k < NSDim; k++) {
334 if (localQR(j,k) != zero) {
335 ja [rowStart+lnnz] = offset + k;
336 val[rowStart+lnnz] = localQR(j,k);
345 size_t ia_tmp = 0, nnz = 0;
346 for (
size_t i = 0; i < numRows; i++) {
347 for (
size_t j = ia_tmp; j < ia[i+1]; j++)
348 if (ja[j] != INVALID) {
356 if (rowMap->lib() == Xpetra::UseTpetra) {
364 GetOStream(
Runtime1) <<
"TentativePFactory : aggregates do not cross process boundaries" << std::endl;
366 PtentCrs->setAllValues(iaPtent, jaPtent, valPtent);
367 PtentCrs->expertStaticFillComplete(coarseMap, A->getDomainMap());
370 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
375 typedef typename STS::magnitudeType Magnitude;
376 const SC zero = STS::zero();
377 const SC one = STS::one();
392 for (GO i=0; i<numAggs; ++i) {
393 LO sizeOfThisAgg = aggStart[i+1] - aggStart[i];
394 if (sizeOfThisAgg > maxAggSize) maxAggSize = sizeOfThisAgg;
398 const size_t NSDim = fineNullspace->getNumVectors();
401 GO indexBase=A->getRowMap()->getIndexBase();
406 RCP<MultiVector> fineNullspaceWithOverlap = MultiVectorFactory::Build(nonUniqueMap,NSDim);
407 fineNullspaceWithOverlap->doImport(*fineNullspace,*importer,Xpetra::INSERT);
411 for (
size_t i=0; i<NSDim; ++i)
412 fineNS[i] = fineNullspaceWithOverlap->getData(i);
415 coarseNullspace = MultiVectorFactory::Build(coarseMap, NSDim);
418 for (
size_t i=0; i<NSDim; ++i)
419 if (coarseMap->getNodeNumElements() > 0) coarseNS[i] = coarseNullspace->getDataNonConst(i);
425 const Map& rowMapForPtentRef = *rowMapForPtent;
440 for (LO j=0; j<numAggs; ++j) {
441 for (LO k=aggStart[j]; k<aggStart[j+1]; ++k) {
442 if (rowMapForPtentRef.isNodeGlobalElement(aggToRowMap[k]) ==
false) {
447 ghostQMap = MapFactory::Build(A->getRowMap()->lib(),
450 indexBase, A->getRowMap()->getComm());
452 ghostQvalues = MultiVectorFactory::Build(ghostQMap,NSDim);
456 ghostQcolumns.
resize(NSDim);
457 for (
size_t i=0; i<NSDim; ++i)
458 ghostQcolumns[i] = Xpetra::VectorFactory<GO,LO,GO,Node>::Build(ghostQMap);
459 ghostQrowNums = Xpetra::VectorFactory<GO,LO,GO,Node>::Build(ghostQMap);
460 if (ghostQvalues->getLocalLength() > 0) {
463 for (
size_t i=0; i<NSDim; ++i) {
464 ghostQvals[i] = ghostQvalues->getDataNonConst(i);
465 ghostQcols[i] = ghostQcolumns[i]->getDataNonConst(0);
467 ghostQrows = ghostQrowNums->getDataNonConst(0);
471 importer = ImportFactory::Build(ghostQMap, A->getRowMap());
477 Array<GO> globalColPtr(maxAggSize*NSDim,0);
478 Array<LO> localColPtr(maxAggSize*NSDim,0);
482 const Map& coarseMapRef = *coarseMap;
501 RCP<CrsMatrixWrap> PtentCrsWrap =
rcp(
new CrsMatrixWrap(rowMapForPtent, NSDim, Xpetra::StaticProfile));
502 PtentCrs = PtentCrsWrap->getCrsMatrix();
503 Ptentative = PtentCrsWrap;
509 const Map& nonUniqueMapRef = *nonUniqueMap;
511 size_t total_nnz_count=0;
513 for (GO agg=0; agg<numAggs; ++agg)
515 LO myAggSize = aggStart[agg+1]-aggStart[agg];
519 for (
size_t j=0; j<NSDim; ++j) {
520 bool bIsZeroNSColumn =
true;
521 for (LO k=0; k<myAggSize; ++k)
526 SC nsVal = fineNS[j][ nonUniqueMapRef.getLocalElement(aggToRowMap[aggStart[agg]+k]) ];
527 localQR(k,j) = nsVal;
528 if (nsVal != zero) bIsZeroNSColumn =
false;
531 GetOStream(
Runtime1,-1) <<
"length of fine level nsp: " << fineNullspace->getGlobalLength() << std::endl;
532 GetOStream(
Runtime1,-1) <<
"length of fine level nsp w overlap: " << fineNullspaceWithOverlap->getGlobalLength() << std::endl;
533 GetOStream(
Runtime1,-1) <<
"(local?) aggId=" << agg << std::endl;
534 GetOStream(
Runtime1,-1) <<
"aggSize=" << myAggSize << std::endl;
535 GetOStream(
Runtime1,-1) <<
"agg DOF=" << k << std::endl;
536 GetOStream(
Runtime1,-1) <<
"NS vector j=" << j << std::endl;
537 GetOStream(
Runtime1,-1) <<
"j*myAggSize + k = " << j*myAggSize + k << std::endl;
538 GetOStream(
Runtime1,-1) <<
"aggToRowMap["<<agg<<
"][" << k <<
"] = " << aggToRowMap[aggStart[agg]+k] << std::endl;
539 GetOStream(
Runtime1,-1) <<
"id aggToRowMap[agg][k]=" << aggToRowMap[aggStart[agg]+k] <<
" is global element in nonUniqueMap = " <<
540 nonUniqueMapRef.isNodeGlobalElement(aggToRowMap[aggStart[agg]+k]) << std::endl;
541 GetOStream(
Runtime1,-1) <<
"colMap local id aggToRowMap[agg][k]=" << nonUniqueMapRef.getLocalElement(aggToRowMap[aggStart[agg]+k]) << std::endl;
542 GetOStream(
Runtime1,-1) <<
"fineNS...=" << fineNS[j][ nonUniqueMapRef.getLocalElement(aggToRowMap[aggStart[agg]+k]) ] << std::endl;
543 GetOStream(
Errors,-1) <<
"caught an error!" << std::endl;
549 Xpetra::global_size_t offset=agg*NSDim;
551 if(myAggSize >= Teuchos::as<LocalOrdinal>(NSDim)) {
556 SC tau = localQR(0,0);
561 for (
size_t k = 0; k < Teuchos::as<size_t>(myAggSize); ++k) {
562 Magnitude tmag = STS::magnitude(localQR(k,0));
567 localQR(0,0) = dtemp;
576 for (
size_t j=0; j<NSDim; ++j) {
577 for (
size_t k=0; k<=j; ++k) {
579 if (coarseMapRef.isNodeLocalElement(offset+k)) {
580 coarseNS[j][offset+k] = localQR(k, j);
584 GetOStream(
Errors,-1) <<
"caught error in coarseNS insert, j="<<j<<
", offset+k = "<<offset+k<<std::endl;
597 for (LocalOrdinal i=0; i<myAggSize; ++i) {
598 localQR(i,0) *= dtemp ;
603 for (
size_t j=0; j<NSDim; j++) {
604 for (
size_t i = 0; i < Teuchos::as<size_t>(myAggSize); i++) {
605 localQR(i,j) = (*qFactor)(i,j);
615 for (
size_t j = 0; j < NSDim; j++)
616 for (
size_t k = 0; k < NSDim; k++) {
618 "Caught error in coarseNS insert, j=" << j <<
", offset+k = " << offset+k);
620 if (k < as<size_t>(myAggSize))
621 coarseNS[j][offset+k] = localQR(k,j);
623 coarseNS[j][offset+k] = (k == j ? one : zero);
627 for (
size_t i = 0; i < as<size_t>(myAggSize); i++)
628 for (
size_t j = 0; j < NSDim; j++)
629 localQR(i,j) = (j == i ? one : zero);
636 for (GO j=0; j<myAggSize; ++j) {
640 GO globalRow = aggToRowMap[aggStart[agg]+j];
643 if (rowMapForPtentRef.isNodeGlobalElement(globalRow) == false ) {
644 ghostQrows[qctr] = globalRow;
645 for (
size_t k=0; k<NSDim; ++k) {
646 ghostQcols[k][qctr] = coarseMapRef.getGlobalElement(agg*NSDim+k);
647 ghostQvals[k][qctr] = localQR(j,k);
652 for (
size_t k=0; k<NSDim; ++k) {
655 localColPtr[nnz] = agg * NSDim + k;
656 globalColPtr[nnz] = coarseMapRef.getGlobalElement(localColPtr[nnz]);
657 valPtr[nnz] = localQR(j,k);
663 GetOStream(
Errors,-1) <<
"caught error in colPtr/valPtr insert, current index="<<nnz<<std::endl;
668 Ptentative->insertGlobalValues(globalRow,globalColPtr.
view(0,nnz),valPtr.
view(0,nnz));
671 GetOStream(
Errors,-1) <<
"pid " << A->getRowMap()->getComm()->getRank()
672 <<
"caught error during Ptent row insertion, global row " 673 << globalRow << std::endl;
684 GetOStream(
Runtime1) <<
"TentativePFactory : aggregates may cross process boundaries" << std::endl;
688 targetQrowNums->putScalar(-1);
689 targetQrowNums->doImport(*ghostQrowNums,*importer,Xpetra::INSERT);
702 gidsToImport, indexBase, A->getRowMap()->getComm() );
705 importer = ImportFactory::Build(ghostQMap, reducedMap);
708 for (
size_t i=0; i<NSDim; ++i) {
709 targetQcolumns[i] = Xpetra::VectorFactory<GO,LO,GO,Node>::Build(reducedMap);
710 targetQcolumns[i]->doImport(*(ghostQcolumns[i]),*importer,Xpetra::INSERT);
712 RCP<MultiVector> targetQvalues = MultiVectorFactory::Build(reducedMap,NSDim);
713 targetQvalues->doImport(*ghostQvalues,*importer,Xpetra::INSERT);
717 if (targetQvalues->getLocalLength() > 0) {
718 targetQvals.
resize(NSDim);
719 targetQcols.
resize(NSDim);
720 for (
size_t i=0; i<NSDim; ++i) {
721 targetQvals[i] = targetQvalues->getDataNonConst(i);
722 targetQcols[i] = targetQcolumns[i]->getDataNonConst(0);
729 if (targetQvalues->getLocalLength() > 0) {
730 for (
size_t j=0; j<NSDim; ++j) {
731 valPtr[j] = targetQvals[j][reducedMap->getLocalElement(*r)];
732 globalColPtr[j] = targetQcols[j][reducedMap->getLocalElement(*r)];
734 Ptentative->insertGlobalValues(*r, globalColPtr.
view(0,NSDim), valPtr.
view(0,NSDim));
738 Ptentative->fillComplete(coarseMap, A->getDomainMap());
741 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
746 const size_t numElements = rowElements.
size();
749 for (
size_t i = 0; i < numElements; i++)
750 if (rowElements[i] != colElements[i]) {
762 #define MUELU_TENTATIVEPFACTORY_SHORT 763 #endif // MUELU_TENTATIVEPFACTORY_DEF_HPP Important warning messages (one line)
void reserve(size_type n)
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getQ() const
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
ArrayView< T > view(size_type offset, size_type size)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Namespace for MueLu classes and methods.
Teuchos::RCP< Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > ComputeUnamalgamatedImportDofMap(const Aggregates &aggregates) const
ComputeUnamalgamatedImportDofMap build overlapping dof row map from aggregates needed for overlapping...
void BuildPuncoupled(RCP< Matrix > A, RCP< Aggregates > aggregates, RCP< AmalgamationInfo > amalgInfo, RCP< MultiVector > fineNullspace, RCP< const Map > coarseMap, RCP< Matrix > &Ptentative, RCP< MultiVector > &coarseNullspace) const
void resize(const size_type n, const T &val=T())
void BuildPcoupled(RCP< Matrix > A, RCP< Aggregates > aggregates, RCP< AmalgamationInfo > amalgInfo, RCP< MultiVector > fineNullspace, RCP< const Map > coarseMap, RCP< Matrix > &Ptentative, RCP< MultiVector > &coarseNullspace) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
void UnamalgamateAggregatesLO(const Aggregates &aggregates, Teuchos::ArrayRCP< LocalOrdinal > &aggStart, Teuchos::ArrayRCP< LO > &aggToRowMap) const
void resize(size_type new_size, const value_type &x=value_type())
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
static magnitudeType magnitude(T a)
void push_back(const value_type &x)
void UnamalgamateAggregates(const Aggregates &aggregates, Teuchos::ArrayRCP< LocalOrdinal > &aggStart, Teuchos::ArrayRCP< GlobalOrdinal > &aggToRowMap) const
UnamalgamateAggregates.
LO GetNumAggregates() const
returns the number of aggregates of the current processor. Note: could/should be renamed to GetNumLoc...
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Exception throws to report errors in the internal logical of the program.
Description of what is happening (more verbose)
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
std::vector< T >::iterator iterator
bool isGoodMap(const Map &rowMap, const Map &colMap) const