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" 70 #include "MueLu_NullspaceFactory.hpp" 71 #include "MueLu_PerfUtils.hpp" 72 #include "MueLu_Utilities.hpp" 76 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
78 RCP<ParameterList> validParamList = rcp(
new ParameterList());
80 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 83 #undef SET_VALID_ENTRY 85 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A");
86 validParamList->set< RCP<const FactoryBase> >(
"Aggregates", Teuchos::null,
"Generating factory of the aggregates");
87 validParamList->set< RCP<const FactoryBase> >(
"Nullspace", Teuchos::null,
"Generating factory of the nullspace");
88 validParamList->set< RCP<const FactoryBase> >(
"UnAmalgamationInfo", Teuchos::null,
"Generating factory of UnAmalgamationInfo");
89 validParamList->set< RCP<const FactoryBase> >(
"CoarseMap", Teuchos::null,
"Generating factory of the coarse map");
90 validParamList->set< RCP<const FactoryBase> >(
"Coordinates", Teuchos::null,
"Generating factory of the coordinates");
93 ParameterList norecurse;
94 norecurse.disableRecursiveValidation();
95 validParamList->set<ParameterList> (
"matrixmatrix: kernel params", norecurse,
"MatrixMatrix kernel parameters");
97 return validParamList;
100 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
103 const ParameterList& pL = GetParameterList();
105 Input(fineLevel,
"A");
106 Input(fineLevel,
"Aggregates");
107 Input(fineLevel,
"Nullspace");
108 Input(fineLevel,
"UnAmalgamationInfo");
109 Input(fineLevel,
"CoarseMap");
112 pL.get<
bool>(
"tentative: build coarse coordinates") ) {
113 bTransferCoordinates_ =
true;
114 Input(fineLevel,
"Coordinates");
115 }
else if (bTransferCoordinates_) {
116 Input(fineLevel,
"Coordinates");
120 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
122 return BuildP(fineLevel, coarseLevel);
125 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
128 typedef typename Teuchos::ScalarTraits<SC>::magnitudeType real_type;
129 typedef typename Xpetra::MultiVector<real_type,LO,GO,NO> RealValuedMultiVector;
131 RCP<Matrix> A = Get< RCP<Matrix> > (fineLevel,
"A");
132 RCP<Aggregates> aggregates = Get< RCP<Aggregates> > (fineLevel,
"Aggregates");
133 RCP<AmalgamationInfo> amalgInfo = Get< RCP<AmalgamationInfo> > (fineLevel,
"UnAmalgamationInfo");
134 RCP<MultiVector> fineNullspace = Get< RCP<MultiVector> > (fineLevel,
"Nullspace");
135 RCP<const Map> coarseMap = Get< RCP<const Map> > (fineLevel,
"CoarseMap");
136 RCP<RealValuedMultiVector> fineCoords;
137 if(bTransferCoordinates_) {
138 fineCoords = Get< RCP<RealValuedMultiVector> >(fineLevel,
"Coordinates");
141 TEUCHOS_TEST_FOR_EXCEPTION(A->getRowMap()->getNodeNumElements() != fineNullspace->getMap()->getNodeNumElements(),
144 RCP<Matrix> Ptentative;
145 RCP<MultiVector> coarseNullspace;
146 RCP<RealValuedMultiVector> coarseCoords;
148 if(bTransferCoordinates_) {
151 ArrayView<const GO> elementAList = coarseMap->getNodeElementList();
153 if (rcp_dynamic_cast<const StridedMap>(coarseMap) != Teuchos::null)
154 blkSize = rcp_dynamic_cast<
const StridedMap>(coarseMap)->getFixedBlockSize();
155 GO indexBase = coarseMap->getIndexBase();
156 size_t numNodes = elementAList.size() / blkSize;
157 Array<GO> nodeList(numNodes);
158 const int numDimensions = fineCoords->getNumVectors();
160 for (LO i = 0; i < Teuchos::as<LO>(numNodes); i++) {
161 nodeList[i] = (elementAList[i*blkSize]-indexBase)/blkSize + indexBase;
163 RCP<const Map> coarseCoordsMap = MapFactory::Build(fineCoords->getMap()->lib(),
164 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
167 fineCoords->getMap()->getComm());
168 coarseCoords = RealValuedMultiVectorFactory::Build(coarseCoordsMap, numDimensions);
171 RCP<RealValuedMultiVector> ghostedCoords;
172 if (aggregates->AggregatesCrossProcessors()) {
173 RCP<const Map> aggMap = aggregates->GetMap();
174 RCP<const Import> importer = ImportFactory::Build(fineCoords->getMap(), aggMap);
176 ghostedCoords = RealValuedMultiVectorFactory::Build(aggMap, numDimensions);
177 ghostedCoords->doImport(*fineCoords, *importer, Xpetra::INSERT);
179 ghostedCoords = fineCoords;
183 int myPID = coarseCoordsMap->getComm()->getRank();
184 LO numAggs = aggregates->GetNumAggregates();
185 ArrayRCP<LO> aggSizes = aggregates->ComputeAggregateSizes();
186 const ArrayRCP<const LO> vertex2AggID = aggregates->GetVertex2AggId()->getData(0);
187 const ArrayRCP<const LO> procWinner = aggregates->GetProcWinner()->getData(0);
190 for (
int dim = 0; dim < numDimensions; ++dim) {
191 ArrayRCP<const real_type> fineCoordsData = ghostedCoords->getData(dim);
192 ArrayRCP<real_type> coarseCoordsData = coarseCoords->getDataNonConst(dim);
194 for (LO lnode = 0; lnode < Teuchos::as<LO>(numNodes); lnode++) {
195 if (procWinner[lnode] == myPID &&
196 lnode < vertex2AggID.size() &&
197 lnode < fineCoordsData.size() &&
198 vertex2AggID[lnode] < coarseCoordsData.size() &&
199 Teuchos::ScalarTraits<double>::isnaninf(fineCoordsData[lnode]) ==
false) {
200 coarseCoordsData[vertex2AggID[lnode]] += fineCoordsData[lnode];
203 for (LO agg = 0; agg < numAggs; agg++) {
204 coarseCoordsData[agg] /= aggSizes[agg];
209 if (!aggregates->AggregatesCrossProcessors())
210 BuildPuncoupled(A, aggregates, amalgInfo, fineNullspace, coarseMap, Ptentative, coarseNullspace,coarseLevel.
GetLevelID());
212 BuildPcoupled (A, aggregates, amalgInfo, fineNullspace, coarseMap, Ptentative, coarseNullspace);
225 if (A->IsView(
"stridedMaps") ==
true)
226 Ptentative->CreateView(
"stridedMaps", A->getRowMap(
"stridedMaps"), coarseMap);
228 Ptentative->CreateView(
"stridedMaps", Ptentative->getRangeMap(), coarseMap);
230 if(bTransferCoordinates_) {
231 Set(coarseLevel,
"Coordinates", coarseCoords);
233 Set(coarseLevel,
"Nullspace", coarseNullspace);
234 Set(coarseLevel,
"P", Ptentative);
237 RCP<ParameterList> params = rcp(
new ParameterList());
238 params->set(
"printLoadBalancingInfo",
true);
243 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
245 BuildPuncoupled(RCP<Matrix> A, RCP<Aggregates> aggregates, RCP<AmalgamationInfo> amalgInfo, RCP<MultiVector> fineNullspace,
246 RCP<const Map> coarseMap, RCP<Matrix>& Ptentative, RCP<MultiVector>& coarseNullspace,
const int levelID)
const {
247 RCP<const Map> rowMap = A->getRowMap();
248 RCP<const Map> colMap = A->getColMap();
250 const size_t numRows = rowMap->getNodeNumElements();
252 typedef Teuchos::ScalarTraits<SC> STS;
253 typedef typename STS::magnitudeType Magnitude;
254 const SC zero = STS::zero();
255 const SC one = STS::one();
256 const LO INVALID = Teuchos::OrdinalTraits<LO>::invalid();
258 const GO numAggs = aggregates->GetNumAggregates();
259 const size_t NSDim = fineNullspace->getNumVectors();
264 bool goodMap = isGoodMap(*rowMap, *colMap);
270 ArrayRCP<LO> aggStart;
271 ArrayRCP<LO> aggToRowMapLO;
272 ArrayRCP<GO> aggToRowMapGO;
274 amalgInfo->UnamalgamateAggregatesLO(*aggregates, aggStart, aggToRowMapLO);
275 GetOStream(
Runtime1) <<
"Column map is consistent with the row map, good." << std::endl;
278 amalgInfo->UnamalgamateAggregates(*aggregates, aggStart, aggToRowMapGO);
279 GetOStream(
Warnings0) <<
"Column map is not consistent with the row map\n" 280 <<
"using GO->LO conversion with performance penalty" << std::endl;
283 coarseNullspace = MultiVectorFactory::Build(coarseMap, NSDim);
285 const ParameterList& pL = GetParameterList();
286 const bool &doQRStep = pL.get<
bool>(
"tentative: calculate qr");
289 ArrayRCP<ArrayRCP<const SC> > fineNS (NSDim);
290 ArrayRCP<ArrayRCP<SC> > coarseNS(NSDim);
291 for (
size_t i = 0; i < NSDim; i++) {
292 fineNS[i] = fineNullspace->getData(i);
293 if (coarseMap->getNodeNumElements() > 0)
294 coarseNS[i] = coarseNullspace->getDataNonConst(i);
297 size_t nnzEstimate = numRows * NSDim;
300 Ptentative = rcp(
new CrsMatrixWrap(rowMap, coarseMap, 0, Xpetra::StaticProfile));
301 RCP<CrsMatrix> PtentCrs = rcp_dynamic_cast<CrsMatrixWrap>(Ptentative)->getCrsMatrix();
303 ArrayRCP<size_t> iaPtent;
304 ArrayRCP<LO> jaPtent;
305 ArrayRCP<SC> valPtent;
307 PtentCrs->allocateAllValues(nnzEstimate, iaPtent, jaPtent, valPtent);
309 ArrayView<size_t> ia = iaPtent();
310 ArrayView<LO> ja = jaPtent();
311 ArrayView<SC> val = valPtent();
314 for (
size_t i = 1; i <= numRows; i++)
315 ia[i] = ia[i-1] + NSDim;
317 for (
size_t j = 0; j < nnzEstimate; j++) {
327 for (GO agg = 0; agg < numAggs; agg++) {
328 LO aggSize = aggStart[agg+1] - aggStart[agg];
330 Xpetra::global_size_t offset = agg*NSDim;
335 Teuchos::SerialDenseMatrix<LO,SC> localQR(aggSize, NSDim);
337 for (
size_t j = 0; j < NSDim; j++)
338 for (LO k = 0; k < aggSize; k++)
339 localQR(k,j) = fineNS[j][aggToRowMapLO[aggStart[agg]+k]];
341 for (
size_t j = 0; j < NSDim; j++)
342 for (LO k = 0; k < aggSize; k++)
343 localQR(k,j) = fineNS[j][rowMap->getLocalElement(aggToRowMapGO[aggStart[agg]+k])];
347 for (
size_t j = 0; j < NSDim; j++) {
348 bool bIsZeroNSColumn =
true;
350 for (LO k = 0; k < aggSize; k++)
351 if (localQR(k,j) != zero)
352 bIsZeroNSColumn =
false;
355 "MueLu::TentativePFactory::MakeTentative: fine level NS part has a zero column");
360 if (aggSize >= Teuchos::as<LO>(NSDim)) {
364 Magnitude norm = STS::magnitude(zero);
365 for (
size_t k = 0; k < Teuchos::as<size_t>(aggSize); k++)
366 norm += STS::magnitude(localQR(k,0)*localQR(k,0));
367 norm = Teuchos::ScalarTraits<Magnitude>::squareroot(norm);
370 coarseNS[0][offset] = norm;
373 for (LO i = 0; i < aggSize; i++)
374 localQR(i,0) /= norm;
377 Teuchos::SerialQRDenseSolver<LO,SC> qrSolver;
378 qrSolver.setMatrix(Teuchos::rcp(&localQR,
false));
382 for (
size_t j = 0; j < NSDim; j++)
383 for (
size_t k = 0; k <= j; k++)
384 coarseNS[j][offset+k] = localQR(k,j);
389 Teuchos::RCP<Teuchos::SerialDenseMatrix<LO,SC> > qFactor = qrSolver.getQ();
390 for (
size_t j = 0; j < NSDim; j++)
391 for (
size_t i = 0; i < Teuchos::as<size_t>(aggSize); i++)
392 localQR(i,j) = (*qFactor)(i,j);
423 for (
size_t j = 0; j < NSDim; j++)
424 for (
size_t k = 0; k < NSDim; k++)
425 if (k < as<size_t>(aggSize))
426 coarseNS[j][offset+k] = localQR(k,j);
428 coarseNS[j][offset+k] = (k == j ? one : zero);
431 for (
size_t i = 0; i < as<size_t>(aggSize); i++)
432 for (
size_t j = 0; j < NSDim; j++)
433 localQR(i,j) = (j == i ? one : zero);
439 for (LO j = 0; j < aggSize; j++) {
440 LO localRow = (goodMap ? aggToRowMapLO[aggStart[agg]+j] : rowMap->getLocalElement(aggToRowMapGO[aggStart[agg]+j]));
442 size_t rowStart = ia[localRow];
443 for (
size_t k = 0, lnnz = 0; k < NSDim; k++) {
445 if (localQR(j,k) != zero) {
446 ja [rowStart+lnnz] = offset + k;
447 val[rowStart+lnnz] = localQR(j,k);
455 GetOStream(
Runtime1) <<
"TentativePFactory : bypassing local QR phase" << std::endl;
457 GetOStream(
Warnings0) <<
"TentativePFactor : for nontrivial nullspace, this may degrade performance" << std::endl;
467 for (GO agg = 0; agg < numAggs; agg++) {
468 const LO aggSize = aggStart[agg+1] - aggStart[agg];
469 Xpetra::global_size_t offset = agg*NSDim;
473 for (LO j = 0; j < aggSize; j++) {
478 const LO localRow = aggToRowMapLO[aggStart[agg]+j];
480 const size_t rowStart = ia[localRow];
482 for (
size_t k = 0, lnnz = 0; k < NSDim; k++) {
484 const SC qr_jk = fineNS[k][aggToRowMapLO[aggStart[agg]+j]];
486 ja [rowStart+lnnz] = offset + k;
487 val[rowStart+lnnz] = qr_jk;
492 for (
size_t j = 0; j < NSDim; j++)
493 coarseNS[j][offset+j] = one;
497 for (GO agg = 0; agg < numAggs; agg++) {
498 const LO aggSize = aggStart[agg+1] - aggStart[agg];
499 Xpetra::global_size_t offset = agg*NSDim;
500 for (LO j = 0; j < aggSize; j++) {
502 const LO localRow = rowMap->getLocalElement(aggToRowMapGO[aggStart[agg]+j]);
504 const size_t rowStart = ia[localRow];
506 for (
size_t k = 0, lnnz = 0; k < NSDim; k++) {
508 const SC qr_jk = fineNS[k][rowMap->getLocalElement(aggToRowMapGO[aggStart[agg]+j])];
510 ja [rowStart+lnnz] = offset + k;
511 val[rowStart+lnnz] = qr_jk;
516 for (
size_t j = 0; j < NSDim; j++)
517 coarseNS[j][offset+j] = one;
526 size_t ia_tmp = 0, nnz = 0;
527 for (
size_t i = 0; i < numRows; i++) {
528 for (
size_t j = ia_tmp; j < ia[i+1]; j++)
529 if (ja[j] != INVALID) {
537 if (rowMap->lib() == Xpetra::UseTpetra) {
541 jaPtent .resize(nnz);
542 valPtent.resize(nnz);
545 GetOStream(
Runtime1) <<
"TentativePFactory : aggregates do not cross process boundaries" << std::endl;
547 PtentCrs->setAllValues(iaPtent, jaPtent, valPtent);
551 RCP<ParameterList> FCparams;
552 if(pL.isSublist(
"matrixmatrix: kernel params"))
553 FCparams=rcp(
new ParameterList(pL.sublist(
"matrixmatrix: kernel params")));
555 FCparams= rcp(
new ParameterList);
557 FCparams->set(
"compute global constants",FCparams->get(
"compute global constants",
false));
558 std::string levelIDs =
toString(levelID);
559 FCparams->set(
"Timer Label",std::string(
"MueLu::TentativeP-")+levelIDs);
560 RCP<const Export> dummy_e;
561 RCP<const Import> dummy_i;
563 PtentCrs->expertStaticFillComplete(coarseMap, A->getDomainMap(),dummy_i,dummy_e,FCparams);
566 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
568 BuildPcoupled(RCP<Matrix> A, RCP<Aggregates> aggregates, RCP<AmalgamationInfo> amalgInfo, RCP<MultiVector> fineNullspace,
569 RCP<const Map> coarseMap, RCP<Matrix>& Ptentative, RCP<MultiVector>& coarseNullspace)
const {
570 typedef Teuchos::ScalarTraits<SC> STS;
571 typedef typename STS::magnitudeType Magnitude;
572 const SC zero = STS::zero();
573 const SC one = STS::one();
576 GO numAggs = aggregates->GetNumAggregates();
582 ArrayRCP<LO> aggStart;
583 ArrayRCP< GO > aggToRowMap;
584 amalgInfo->UnamalgamateAggregates(*aggregates, aggStart, aggToRowMap);
588 for (GO i=0; i<numAggs; ++i) {
589 LO sizeOfThisAgg = aggStart[i+1] - aggStart[i];
590 if (sizeOfThisAgg > maxAggSize) maxAggSize = sizeOfThisAgg;
594 const size_t NSDim = fineNullspace->getNumVectors();
597 GO indexBase=A->getRowMap()->getIndexBase();
599 const RCP<const Map> nonUniqueMap = amalgInfo->ComputeUnamalgamatedImportDofMap(*aggregates);
600 const RCP<const Map> uniqueMap = A->getDomainMap();
601 RCP<const Import> importer = ImportFactory::Build(uniqueMap, nonUniqueMap);
602 RCP<MultiVector> fineNullspaceWithOverlap = MultiVectorFactory::Build(nonUniqueMap,NSDim);
603 fineNullspaceWithOverlap->doImport(*fineNullspace,*importer,Xpetra::INSERT);
606 ArrayRCP< ArrayRCP<const SC> > fineNS(NSDim);
607 for (
size_t i=0; i<NSDim; ++i)
608 fineNS[i] = fineNullspaceWithOverlap->getData(i);
611 coarseNullspace = MultiVectorFactory::Build(coarseMap, NSDim);
613 ArrayRCP< ArrayRCP<SC> > coarseNS(NSDim);
614 for (
size_t i=0; i<NSDim; ++i)
615 if (coarseMap->getNodeNumElements() > 0) coarseNS[i] = coarseNullspace->getDataNonConst(i);
620 RCP<const Map > rowMapForPtent = A->getRowMap();
621 const Map& rowMapForPtentRef = *rowMapForPtent;
625 RCP<const Map> colMap = A->getColMap();
627 RCP<const Map > ghostQMap;
628 RCP<MultiVector> ghostQvalues;
629 Array<RCP<Xpetra::Vector<GO,LO,GO,Node> > > ghostQcolumns;
630 RCP<Xpetra::Vector<GO,LO,GO,Node> > ghostQrowNums;
631 ArrayRCP< ArrayRCP<SC> > ghostQvals;
632 ArrayRCP< ArrayRCP<GO> > ghostQcols;
633 ArrayRCP< GO > ghostQrows;
636 for (LO j=0; j<numAggs; ++j) {
637 for (LO k=aggStart[j]; k<aggStart[j+1]; ++k) {
638 if (rowMapForPtentRef.isNodeGlobalElement(aggToRowMap[k]) ==
false) {
639 ghostGIDs.push_back(aggToRowMap[k]);
643 ghostQMap = MapFactory::Build(A->getRowMap()->lib(),
644 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
646 indexBase, A->getRowMap()->getComm());
648 ghostQvalues = MultiVectorFactory::Build(ghostQMap,NSDim);
652 ghostQcolumns.resize(NSDim);
653 for (
size_t i=0; i<NSDim; ++i)
654 ghostQcolumns[i] = Xpetra::VectorFactory<GO,LO,GO,Node>::Build(ghostQMap);
655 ghostQrowNums = Xpetra::VectorFactory<GO,LO,GO,Node>::Build(ghostQMap);
656 if (ghostQvalues->getLocalLength() > 0) {
657 ghostQvals.resize(NSDim);
658 ghostQcols.resize(NSDim);
659 for (
size_t i=0; i<NSDim; ++i) {
660 ghostQvals[i] = ghostQvalues->getDataNonConst(i);
661 ghostQcols[i] = ghostQcolumns[i]->getDataNonConst(0);
663 ghostQrows = ghostQrowNums->getDataNonConst(0);
667 importer = ImportFactory::Build(ghostQMap, A->getRowMap());
670 Teuchos::SerialQRDenseSolver<LO,SC> qrSolver;
673 Array<GO> globalColPtr(maxAggSize*NSDim,0);
674 Array<LO> localColPtr(maxAggSize*NSDim,0);
675 Array<SC> valPtr(maxAggSize*NSDim,0.);
678 const Map& coarseMapRef = *coarseMap;
681 ArrayRCP<size_t> ptent_rowptr;
682 ArrayRCP<LO> ptent_colind;
683 ArrayRCP<Scalar> ptent_values;
686 ArrayView<size_t> rowptr_v;
687 ArrayView<LO> colind_v;
688 ArrayView<Scalar> values_v;
691 Array<size_t> rowptr_temp;
692 Array<LO> colind_temp;
693 Array<Scalar> values_temp;
695 RCP<CrsMatrix> PtentCrs;
697 RCP<CrsMatrixWrap> PtentCrsWrap = rcp(
new CrsMatrixWrap(rowMapForPtent, NSDim, Xpetra::StaticProfile));
698 PtentCrs = PtentCrsWrap->getCrsMatrix();
699 Ptentative = PtentCrsWrap;
705 const Map& nonUniqueMapRef = *nonUniqueMap;
707 size_t total_nnz_count=0;
709 for (GO agg=0; agg<numAggs; ++agg)
711 LO myAggSize = aggStart[agg+1]-aggStart[agg];
714 Teuchos::SerialDenseMatrix<LO,SC> localQR(myAggSize, NSDim);
715 for (
size_t j=0; j<NSDim; ++j) {
716 bool bIsZeroNSColumn =
true;
717 for (LO k=0; k<myAggSize; ++k)
722 SC nsVal = fineNS[j][ nonUniqueMapRef.getLocalElement(aggToRowMap[aggStart[agg]+k]) ];
723 localQR(k,j) = nsVal;
724 if (nsVal != zero) bIsZeroNSColumn =
false;
727 GetOStream(
Runtime1,-1) <<
"length of fine level nsp: " << fineNullspace->getGlobalLength() << std::endl;
728 GetOStream(
Runtime1,-1) <<
"length of fine level nsp w overlap: " << fineNullspaceWithOverlap->getGlobalLength() << std::endl;
729 GetOStream(
Runtime1,-1) <<
"(local?) aggId=" << agg << std::endl;
730 GetOStream(
Runtime1,-1) <<
"aggSize=" << myAggSize << std::endl;
731 GetOStream(
Runtime1,-1) <<
"agg DOF=" << k << std::endl;
732 GetOStream(
Runtime1,-1) <<
"NS vector j=" << j << std::endl;
733 GetOStream(
Runtime1,-1) <<
"j*myAggSize + k = " << j*myAggSize + k << std::endl;
734 GetOStream(
Runtime1,-1) <<
"aggToRowMap["<<agg<<
"][" << k <<
"] = " << aggToRowMap[aggStart[agg]+k] << std::endl;
735 GetOStream(
Runtime1,-1) <<
"id aggToRowMap[agg][k]=" << aggToRowMap[aggStart[agg]+k] <<
" is global element in nonUniqueMap = " <<
736 nonUniqueMapRef.isNodeGlobalElement(aggToRowMap[aggStart[agg]+k]) << std::endl;
737 GetOStream(
Runtime1,-1) <<
"colMap local id aggToRowMap[agg][k]=" << nonUniqueMapRef.getLocalElement(aggToRowMap[aggStart[agg]+k]) << std::endl;
738 GetOStream(
Runtime1,-1) <<
"fineNS...=" << fineNS[j][ nonUniqueMapRef.getLocalElement(aggToRowMap[aggStart[agg]+k]) ] << std::endl;
739 GetOStream(
Errors,-1) <<
"caught an error!" << std::endl;
742 TEUCHOS_TEST_FOR_EXCEPTION(bIsZeroNSColumn ==
true,
Exceptions::RuntimeError,
"MueLu::TentativePFactory::MakeTentative: fine level NS part has a zero column. Error.");
745 Xpetra::global_size_t offset=agg*NSDim;
747 if(myAggSize >= Teuchos::as<LocalOrdinal>(NSDim)) {
752 SC tau = localQR(0,0);
757 for (
size_t k = 0; k < Teuchos::as<size_t>(myAggSize); ++k) {
758 Magnitude tmag = STS::magnitude(localQR(k,0));
761 dtemp = Teuchos::ScalarTraits<Magnitude>::squareroot(dtemp);
763 localQR(0,0) = dtemp;
765 qrSolver.setMatrix( Teuchos::rcp(&localQR,
false) );
772 for (
size_t j=0; j<NSDim; ++j) {
773 for (
size_t k=0; k<=j; ++k) {
775 if (coarseMapRef.isNodeLocalElement(offset+k)) {
776 coarseNS[j][offset+k] = localQR(k, j);
780 GetOStream(
Errors,-1) <<
"caught error in coarseNS insert, j="<<j<<
", offset+k = "<<offset+k<<std::endl;
790 Magnitude dtemp = Teuchos::ScalarTraits<SC>::magnitude(localQR(0,0));
793 for (LocalOrdinal i=0; i<myAggSize; ++i) {
794 localQR(i,0) *= dtemp ;
798 Teuchos::RCP<Teuchos::SerialDenseMatrix<LO,SC> > qFactor = qrSolver.getQ();
799 for (
size_t j=0; j<NSDim; j++) {
800 for (
size_t i = 0; i < Teuchos::as<size_t>(myAggSize); i++) {
801 localQR(i,j) = (*qFactor)(i,j);
811 for (
size_t j = 0; j < NSDim; j++)
812 for (
size_t k = 0; k < NSDim; k++) {
814 "Caught error in coarseNS insert, j=" << j <<
", offset+k = " << offset+k);
816 if (k < as<size_t>(myAggSize))
817 coarseNS[j][offset+k] = localQR(k,j);
819 coarseNS[j][offset+k] = (k == j ? one : zero);
823 for (
size_t i = 0; i < as<size_t>(myAggSize); i++)
824 for (
size_t j = 0; j < NSDim; j++)
825 localQR(i,j) = (j == i ? one : zero);
832 for (GO j=0; j<myAggSize; ++j) {
836 GO globalRow = aggToRowMap[aggStart[agg]+j];
839 if (rowMapForPtentRef.isNodeGlobalElement(globalRow) == false ) {
840 ghostQrows[qctr] = globalRow;
841 for (
size_t k=0; k<NSDim; ++k) {
842 ghostQcols[k][qctr] = coarseMapRef.getGlobalElement(agg*NSDim+k);
843 ghostQvals[k][qctr] = localQR(j,k);
848 for (
size_t k=0; k<NSDim; ++k) {
850 if (localQR(j,k) != Teuchos::ScalarTraits<SC>::zero()) {
851 localColPtr[nnz] = agg * NSDim + k;
852 globalColPtr[nnz] = coarseMapRef.getGlobalElement(localColPtr[nnz]);
853 valPtr[nnz] = localQR(j,k);
859 GetOStream(
Errors,-1) <<
"caught error in colPtr/valPtr insert, current index="<<nnz<<std::endl;
864 Ptentative->insertGlobalValues(globalRow,globalColPtr.view(0,nnz),valPtr.view(0,nnz));
867 GetOStream(
Errors,-1) <<
"pid " << A->getRowMap()->getComm()->getRank()
868 <<
"caught error during Ptent row insertion, global row " 869 << globalRow << std::endl;
880 GetOStream(
Runtime1) <<
"TentativePFactory : aggregates may cross process boundaries" << std::endl;
883 RCP<Xpetra::Vector<GO,LO,GO,Node> > targetQrowNums = Xpetra::VectorFactory<GO,LO,GO,Node>::Build(rowMapForPtent);
884 targetQrowNums->putScalar(-1);
885 targetQrowNums->doImport(*ghostQrowNums,*importer,Xpetra::INSERT);
886 ArrayRCP< GO > targetQrows = targetQrowNums->getDataNonConst(0);
889 Array<GO> gidsToImport;
890 gidsToImport.reserve(targetQrows.size());
891 for (
typename ArrayRCP<GO>::iterator r=targetQrows.begin(); r!=targetQrows.end(); ++r) {
893 gidsToImport.push_back(*r);
896 RCP<const Map > reducedMap = MapFactory::Build( A->getRowMap()->lib(),
897 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
898 gidsToImport, indexBase, A->getRowMap()->getComm() );
901 importer = ImportFactory::Build(ghostQMap, reducedMap);
903 Array<RCP<Xpetra::Vector<GO,LO,GO,Node> > > targetQcolumns(NSDim);
904 for (
size_t i=0; i<NSDim; ++i) {
905 targetQcolumns[i] = Xpetra::VectorFactory<GO,LO,GO,Node>::Build(reducedMap);
906 targetQcolumns[i]->doImport(*(ghostQcolumns[i]),*importer,Xpetra::INSERT);
908 RCP<MultiVector> targetQvalues = MultiVectorFactory::Build(reducedMap,NSDim);
909 targetQvalues->doImport(*ghostQvalues,*importer,Xpetra::INSERT);
911 ArrayRCP< ArrayRCP<SC> > targetQvals;
912 ArrayRCP<ArrayRCP<GO> > targetQcols;
913 if (targetQvalues->getLocalLength() > 0) {
914 targetQvals.resize(NSDim);
915 targetQcols.resize(NSDim);
916 for (
size_t i=0; i<NSDim; ++i) {
917 targetQvals[i] = targetQvalues->getDataNonConst(i);
918 targetQcols[i] = targetQcolumns[i]->getDataNonConst(0);
922 valPtr = Array<SC>(NSDim,0.);
923 globalColPtr = Array<GO>(NSDim,0);
924 for (
typename Array<GO>::iterator r=gidsToImport.begin(); r!=gidsToImport.end(); ++r) {
925 if (targetQvalues->getLocalLength() > 0) {
926 for (
size_t j=0; j<NSDim; ++j) {
927 valPtr[j] = targetQvals[j][reducedMap->getLocalElement(*r)];
928 globalColPtr[j] = targetQcols[j][reducedMap->getLocalElement(*r)];
930 Ptentative->insertGlobalValues(*r, globalColPtr.view(0,NSDim), valPtr.view(0,NSDim));
934 Ptentative->fillComplete(coarseMap, A->getDomainMap());
937 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
939 ArrayView<const GO> rowElements = rowMap.getNodeElementList();
940 ArrayView<const GO> colElements = colMap.getNodeElementList();
942 const size_t numElements = rowElements.size();
945 for (
size_t i = 0; i < numElements; i++)
946 if (rowElements[i] != colElements[i]) {
958 #define MUELU_TENTATIVEPFACTORY_SHORT 959 #endif // MUELU_TENTATIVEPFACTORY_DEF_HPP Important warning messages (one line)
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Timer to be used in factories. Similar to Monitor but with additional timers.
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.
static const NoFactory * get()
Print even more statistics.
int GetLevelID() const
Return level number.
#define SET_VALID_ENTRY(name)
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
Class that holds all level-specific information.
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 int levelID) const
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
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.
bool isGoodMap(const Map &rowMap, const Map &colMap) const