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)) 84 #undef SET_VALID_ENTRY 85 validParamList->set< std::string >(
"Nullspace name",
"Nullspace",
"Name for the input nullspace");
87 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A");
88 validParamList->set< RCP<const FactoryBase> >(
"Aggregates", Teuchos::null,
"Generating factory of the aggregates");
89 validParamList->set< RCP<const FactoryBase> >(
"Nullspace", Teuchos::null,
"Generating factory of the nullspace");
90 validParamList->set< RCP<const FactoryBase> >(
"Scaled Nullspace", Teuchos::null,
"Generating factory of the scaled nullspace");
91 validParamList->set< RCP<const FactoryBase> >(
"UnAmalgamationInfo", Teuchos::null,
"Generating factory of UnAmalgamationInfo");
92 validParamList->set< RCP<const FactoryBase> >(
"CoarseMap", Teuchos::null,
"Generating factory of the coarse map");
93 validParamList->set< RCP<const FactoryBase> >(
"Coordinates", Teuchos::null,
"Generating factory of the coordinates");
94 validParamList->set< RCP<const FactoryBase> >(
"Node Comm", Teuchos::null,
"Generating factory of the node level communicator");
97 ParameterList norecurse;
98 norecurse.disableRecursiveValidation();
99 validParamList->set<ParameterList> (
"matrixmatrix: kernel params", norecurse,
"MatrixMatrix kernel parameters");
101 return validParamList;
104 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
107 const ParameterList& pL = GetParameterList();
109 std::string nspName =
"Nullspace";
110 if(pL.isParameter(
"Nullspace name")) nspName = pL.get<std::string>(
"Nullspace name");
112 Input(fineLevel,
"A");
113 Input(fineLevel,
"Aggregates");
114 Input(fineLevel, nspName);
115 Input(fineLevel,
"UnAmalgamationInfo");
116 Input(fineLevel,
"CoarseMap");
119 pL.get<
bool>(
"tentative: build coarse coordinates") ) {
120 bTransferCoordinates_ =
true;
121 Input(fineLevel,
"Coordinates");
122 }
else if (bTransferCoordinates_) {
123 Input(fineLevel,
"Coordinates");
127 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
129 return BuildP(fineLevel, coarseLevel);
132 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
135 typedef typename Teuchos::ScalarTraits<Scalar>::coordinateType coordinate_type;
136 typedef Xpetra::MultiVector<coordinate_type,LO,GO,NO> RealValuedMultiVector;
137 typedef Xpetra::MultiVectorFactory<coordinate_type,LO,GO,NO> RealValuedMultiVectorFactory;
139 const ParameterList& pL = GetParameterList();
140 std::string nspName =
"Nullspace";
141 if(pL.isParameter(
"Nullspace name")) nspName = pL.get<std::string>(
"Nullspace name");
144 RCP<Matrix> A = Get< RCP<Matrix> > (fineLevel,
"A");
145 RCP<Aggregates> aggregates = Get< RCP<Aggregates> > (fineLevel,
"Aggregates");
146 RCP<AmalgamationInfo> amalgInfo = Get< RCP<AmalgamationInfo> > (fineLevel,
"UnAmalgamationInfo");
147 RCP<MultiVector> fineNullspace = Get< RCP<MultiVector> > (fineLevel, nspName);
148 RCP<const Map> coarseMap = Get< RCP<const Map> > (fineLevel,
"CoarseMap");
149 RCP<RealValuedMultiVector> fineCoords;
150 if(bTransferCoordinates_) {
151 fineCoords = Get< RCP<RealValuedMultiVector> >(fineLevel,
"Coordinates");
155 if(fineLevel.IsAvailable(
"Node Comm")) {
156 RCP<const Teuchos::Comm<int> > nodeComm = Get<RCP<const Teuchos::Comm<int> > >(fineLevel,
"Node Comm");
157 Set<RCP<const Teuchos::Comm<int> > >(coarseLevel,
"Node Comm", nodeComm);
160 TEUCHOS_TEST_FOR_EXCEPTION(A->getRowMap()->getNodeNumElements() != fineNullspace->getMap()->getNodeNumElements(),
163 RCP<Matrix> Ptentative;
164 RCP<MultiVector> coarseNullspace;
165 RCP<RealValuedMultiVector> coarseCoords;
167 if(bTransferCoordinates_) {
170 ArrayView<const GO> elementAList = coarseMap->getNodeElementList();
172 if (rcp_dynamic_cast<const StridedMap>(coarseMap) != Teuchos::null) {
173 blkSize = rcp_dynamic_cast<
const StridedMap>(coarseMap)->getFixedBlockSize();
175 GO indexBase = coarseMap->getIndexBase();
176 LO numCoarseNodes = Teuchos::as<LO>(elementAList.size() / blkSize);
177 Array<GO> nodeList(numCoarseNodes);
178 const int numDimensions = fineCoords->getNumVectors();
180 for (LO i = 0; i < numCoarseNodes; i++) {
181 nodeList[i] = (elementAList[i*blkSize]-indexBase)/blkSize + indexBase;
183 RCP<const Map> coarseCoordsMap = MapFactory::Build(fineCoords->getMap()->lib(),
184 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
187 fineCoords->getMap()->getComm());
188 coarseCoords = RealValuedMultiVectorFactory::Build(coarseCoordsMap, numDimensions);
191 RCP<RealValuedMultiVector> ghostedCoords;
192 if (aggregates->AggregatesCrossProcessors()) {
193 RCP<const Map> aggMap = aggregates->GetMap();
194 RCP<const Import> importer = ImportFactory::Build(fineCoords->getMap(), aggMap);
196 ghostedCoords = RealValuedMultiVectorFactory::Build(aggMap, numDimensions);
197 ghostedCoords->doImport(*fineCoords, *importer, Xpetra::INSERT);
199 ghostedCoords = fineCoords;
203 int myPID = coarseCoordsMap->getComm()->getRank();
204 LO numAggs = aggregates->GetNumAggregates();
205 ArrayRCP<LO> aggSizes = aggregates->ComputeAggregateSizes();
206 const ArrayRCP<const LO> vertex2AggID = aggregates->GetVertex2AggId()->getData(0);
207 const ArrayRCP<const LO> procWinner = aggregates->GetProcWinner()->getData(0);
210 for (
int dim = 0; dim < numDimensions; ++dim) {
211 ArrayRCP<const coordinate_type> fineCoordsData = ghostedCoords->getData(dim);
212 ArrayRCP<coordinate_type> coarseCoordsData = coarseCoords->getDataNonConst(dim);
214 for (LO lnode = 0; lnode < Teuchos::as<LO>(vertex2AggID.size()); lnode++) {
215 if (procWinner[lnode] == myPID &&
216 lnode < fineCoordsData.size() &&
217 vertex2AggID[lnode] < coarseCoordsData.size() &&
218 Teuchos::ScalarTraits<coordinate_type>::isnaninf(fineCoordsData[lnode]) ==
false) {
219 coarseCoordsData[vertex2AggID[lnode]] += fineCoordsData[lnode];
222 for (LO agg = 0; agg < numAggs; agg++) {
223 coarseCoordsData[agg] /= aggSizes[agg];
228 if (!aggregates->AggregatesCrossProcessors())
229 BuildPuncoupled(A, aggregates, amalgInfo, fineNullspace, coarseMap, Ptentative, coarseNullspace,coarseLevel.
GetLevelID());
231 BuildPcoupled (A, aggregates, amalgInfo, fineNullspace, coarseMap, Ptentative, coarseNullspace);
241 if (A->IsView(
"stridedMaps") ==
true)
242 Ptentative->CreateView(
"stridedMaps", A->getRowMap(
"stridedMaps"), coarseMap);
244 Ptentative->CreateView(
"stridedMaps", Ptentative->getRangeMap(), coarseMap);
246 if(bTransferCoordinates_) {
247 Set(coarseLevel,
"Coordinates", coarseCoords);
249 Set(coarseLevel,
"Nullspace", coarseNullspace);
250 Set(coarseLevel,
"P", Ptentative);
253 RCP<ParameterList> params = rcp(
new ParameterList());
254 params->set(
"printLoadBalancingInfo",
true);
259 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
261 BuildPuncoupled(RCP<Matrix> A, RCP<Aggregates> aggregates, RCP<AmalgamationInfo> amalgInfo, RCP<MultiVector> fineNullspace,
262 RCP<const Map> coarseMap, RCP<Matrix>& Ptentative, RCP<MultiVector>& coarseNullspace,
const int levelID)
const {
263 RCP<const Map> rowMap = A->getRowMap();
264 RCP<const Map> colMap = A->getColMap();
266 const size_t numRows = rowMap->getNodeNumElements();
268 typedef Teuchos::ScalarTraits<SC> STS;
269 typedef typename STS::magnitudeType Magnitude;
270 const SC zero = STS::zero();
271 const SC one = STS::one();
272 const LO INVALID = Teuchos::OrdinalTraits<LO>::invalid();
274 const GO numAggs = aggregates->GetNumAggregates();
275 const size_t NSDim = fineNullspace->getNumVectors();
276 ArrayRCP<LO> aggSizes = aggregates->ComputeAggregateSizes();
280 const ParameterList& pL = GetParameterList();
281 const bool &doQRStep = pL.get<
bool>(
"tentative: calculate qr");
282 const bool &constantColSums = pL.get<
bool>(
"tentative: constant column sums");
285 "MueLu::TentativePFactory::MakeTentative: cannot use 'constant column sums' and 'calculate qr' at the same time");
296 ArrayRCP<LO> aggStart;
297 ArrayRCP<LO> aggToRowMapLO;
298 ArrayRCP<GO> aggToRowMapGO;
300 amalgInfo->UnamalgamateAggregatesLO(*aggregates, aggStart, aggToRowMapLO);
301 GetOStream(
Runtime1) <<
"Column map is consistent with the row map, good." << std::endl;
304 amalgInfo->UnamalgamateAggregates(*aggregates, aggStart, aggToRowMapGO);
305 GetOStream(
Warnings0) <<
"Column map is not consistent with the row map\n" 306 <<
"using GO->LO conversion with performance penalty" << std::endl;
309 coarseNullspace = MultiVectorFactory::Build(coarseMap, NSDim);
312 ArrayRCP<ArrayRCP<const SC> > fineNS (NSDim);
313 ArrayRCP<ArrayRCP<SC> > coarseNS(NSDim);
314 for (
size_t i = 0; i < NSDim; i++) {
315 fineNS[i] = fineNullspace->getData(i);
316 if (coarseMap->getNodeNumElements() > 0)
317 coarseNS[i] = coarseNullspace->getDataNonConst(i);
320 size_t nnzEstimate = numRows * NSDim;
323 Ptentative = rcp(
new CrsMatrixWrap(rowMap, coarseMap, 0));
324 RCP<CrsMatrix> PtentCrs = rcp_dynamic_cast<CrsMatrixWrap>(Ptentative)->getCrsMatrix();
326 ArrayRCP<size_t> iaPtent;
327 ArrayRCP<LO> jaPtent;
328 ArrayRCP<SC> valPtent;
330 PtentCrs->allocateAllValues(nnzEstimate, iaPtent, jaPtent, valPtent);
332 ArrayView<size_t> ia = iaPtent();
333 ArrayView<LO> ja = jaPtent();
334 ArrayView<SC> val = valPtent();
337 for (
size_t i = 1; i <= numRows; i++)
338 ia[i] = ia[i-1] + NSDim;
340 for (
size_t j = 0; j < nnzEstimate; j++) {
350 for (GO agg = 0; agg < numAggs; agg++) {
351 LO aggSize = aggStart[agg+1] - aggStart[agg];
353 Xpetra::global_size_t offset = agg*NSDim;
358 Teuchos::SerialDenseMatrix<LO,SC> localQR(aggSize, NSDim);
360 for (
size_t j = 0; j < NSDim; j++)
361 for (LO k = 0; k < aggSize; k++)
362 localQR(k,j) = fineNS[j][aggToRowMapLO[aggStart[agg]+k]];
364 for (
size_t j = 0; j < NSDim; j++)
365 for (LO k = 0; k < aggSize; k++)
366 localQR(k,j) = fineNS[j][rowMap->getLocalElement(aggToRowMapGO[aggStart[agg]+k])];
370 for (
size_t j = 0; j < NSDim; j++) {
371 bool bIsZeroNSColumn =
true;
373 for (LO k = 0; k < aggSize; k++)
374 if (localQR(k,j) != zero)
375 bIsZeroNSColumn =
false;
378 "MueLu::TentativePFactory::MakeTentative: fine level NS part has a zero column in NS column " << j);
383 if (aggSize >= Teuchos::as<LO>(NSDim)) {
387 Magnitude norm = STS::magnitude(zero);
388 for (
size_t k = 0; k < Teuchos::as<size_t>(aggSize); k++)
389 norm += STS::magnitude(localQR(k,0)*localQR(k,0));
390 norm = Teuchos::ScalarTraits<Magnitude>::squareroot(norm);
393 coarseNS[0][offset] = norm;
396 for (LO i = 0; i < aggSize; i++)
397 localQR(i,0) /= norm;
400 Teuchos::SerialQRDenseSolver<LO,SC> qrSolver;
401 qrSolver.setMatrix(Teuchos::rcp(&localQR,
false));
405 for (
size_t j = 0; j < NSDim; j++)
406 for (
size_t k = 0; k <= j; k++)
407 coarseNS[j][offset+k] = localQR(k,j);
412 Teuchos::RCP<Teuchos::SerialDenseMatrix<LO,SC> > qFactor = qrSolver.getQ();
413 for (
size_t j = 0; j < NSDim; j++)
414 for (
size_t i = 0; i < Teuchos::as<size_t>(aggSize); i++)
415 localQR(i,j) = (*qFactor)(i,j);
446 for (
size_t j = 0; j < NSDim; j++)
447 for (
size_t k = 0; k < NSDim; k++)
448 if (k < as<size_t>(aggSize))
449 coarseNS[j][offset+k] = localQR(k,j);
451 coarseNS[j][offset+k] = (k == j ? one : zero);
454 for (
size_t i = 0; i < as<size_t>(aggSize); i++)
455 for (
size_t j = 0; j < NSDim; j++)
456 localQR(i,j) = (j == i ? one : zero);
462 for (LO j = 0; j < aggSize; j++) {
463 LO localRow = (goodMap ? aggToRowMapLO[aggStart[agg]+j] : rowMap->getLocalElement(aggToRowMapGO[aggStart[agg]+j]));
465 size_t rowStart = ia[localRow];
466 for (
size_t k = 0, lnnz = 0; k < NSDim; k++) {
468 if (localQR(j,k) != zero) {
469 ja [rowStart+lnnz] = offset + k;
470 val[rowStart+lnnz] = localQR(j,k);
478 GetOStream(
Runtime1) <<
"TentativePFactory : bypassing local QR phase" << std::endl;
480 GetOStream(
Warnings0) <<
"TentativePFactory : for nontrivial nullspace, this may degrade performance" << std::endl;
490 for (GO agg = 0; agg < numAggs; agg++) {
491 const LO aggSize = aggStart[agg+1] - aggStart[agg];
492 Xpetra::global_size_t offset = agg*NSDim;
496 for (LO j = 0; j < aggSize; j++) {
501 const LO localRow = aggToRowMapLO[aggStart[agg]+j];
503 const size_t rowStart = ia[localRow];
505 for (
size_t k = 0, lnnz = 0; k < NSDim; k++) {
507 SC qr_jk = fineNS[k][aggToRowMapLO[aggStart[agg]+j]];
508 if(constantColSums) qr_jk = qr_jk / (double)aggSizes[agg];
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;
521 for (GO agg = 0; agg < numAggs; agg++) {
522 const LO aggSize = aggStart[agg+1] - aggStart[agg];
523 Xpetra::global_size_t offset = agg*NSDim;
524 for (LO j = 0; j < aggSize; j++) {
526 const LO localRow = rowMap->getLocalElement(aggToRowMapGO[aggStart[agg]+j]);
528 const size_t rowStart = ia[localRow];
530 for (
size_t k = 0, lnnz = 0; k < NSDim; ++k) {
532 SC qr_jk = fineNS[k][rowMap->getLocalElement(aggToRowMapGO[aggStart[agg]+j])];
533 if(constantColSums) qr_jk = qr_jk / (double)aggSizes[agg];
535 ja [rowStart+lnnz] = offset + k;
536 val[rowStart+lnnz] = qr_jk;
541 for (
size_t j = 0; j < NSDim; j++)
542 coarseNS[j][offset+j] = one;
551 size_t ia_tmp = 0, nnz = 0;
552 for (
size_t i = 0; i < numRows; i++) {
553 for (
size_t j = ia_tmp; j < ia[i+1]; j++)
554 if (ja[j] != INVALID) {
562 if (rowMap->lib() == Xpetra::UseTpetra) {
566 jaPtent .resize(nnz);
567 valPtent.resize(nnz);
570 GetOStream(
Runtime1) <<
"TentativePFactory : aggregates do not cross process boundaries" << std::endl;
572 PtentCrs->setAllValues(iaPtent, jaPtent, valPtent);
576 RCP<ParameterList> FCparams;
577 if(pL.isSublist(
"matrixmatrix: kernel params"))
578 FCparams=rcp(
new ParameterList(pL.sublist(
"matrixmatrix: kernel params")));
580 FCparams= rcp(
new ParameterList);
582 FCparams->set(
"compute global constants",FCparams->get(
"compute global constants",
false));
583 std::string levelIDs =
toString(levelID);
584 FCparams->set(
"Timer Label",std::string(
"MueLu::TentativeP-")+levelIDs);
585 RCP<const Export> dummy_e;
586 RCP<const Import> dummy_i;
588 PtentCrs->expertStaticFillComplete(coarseMap, A->getDomainMap(),dummy_i,dummy_e,FCparams);
591 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
593 BuildPcoupled(RCP<Matrix> A, RCP<Aggregates> aggregates, RCP<AmalgamationInfo> amalgInfo, RCP<MultiVector> fineNullspace,
594 RCP<const Map> coarseMap, RCP<Matrix>& Ptentative, RCP<MultiVector>& coarseNullspace)
const {
595 typedef Teuchos::ScalarTraits<SC> STS;
596 typedef typename STS::magnitudeType Magnitude;
597 const SC zero = STS::zero();
598 const SC one = STS::one();
601 GO numAggs = aggregates->GetNumAggregates();
607 ArrayRCP<LO> aggStart;
608 ArrayRCP< GO > aggToRowMap;
609 amalgInfo->UnamalgamateAggregates(*aggregates, aggStart, aggToRowMap);
613 for (GO i=0; i<numAggs; ++i) {
614 LO sizeOfThisAgg = aggStart[i+1] - aggStart[i];
615 if (sizeOfThisAgg > maxAggSize) maxAggSize = sizeOfThisAgg;
619 const size_t NSDim = fineNullspace->getNumVectors();
622 GO indexBase=A->getRowMap()->getIndexBase();
624 const RCP<const Map> nonUniqueMap = amalgInfo->ComputeUnamalgamatedImportDofMap(*aggregates);
625 const RCP<const Map> uniqueMap = A->getDomainMap();
626 RCP<const Import> importer = ImportFactory::Build(uniqueMap, nonUniqueMap);
627 RCP<MultiVector> fineNullspaceWithOverlap = MultiVectorFactory::Build(nonUniqueMap,NSDim);
628 fineNullspaceWithOverlap->doImport(*fineNullspace,*importer,Xpetra::INSERT);
631 ArrayRCP< ArrayRCP<const SC> > fineNS(NSDim);
632 for (
size_t i=0; i<NSDim; ++i)
633 fineNS[i] = fineNullspaceWithOverlap->getData(i);
636 coarseNullspace = MultiVectorFactory::Build(coarseMap, NSDim);
638 ArrayRCP< ArrayRCP<SC> > coarseNS(NSDim);
639 for (
size_t i=0; i<NSDim; ++i)
640 if (coarseMap->getNodeNumElements() > 0) coarseNS[i] = coarseNullspace->getDataNonConst(i);
645 RCP<const Map > rowMapForPtent = A->getRowMap();
646 const Map& rowMapForPtentRef = *rowMapForPtent;
650 RCP<const Map> colMap = A->getColMap();
652 RCP<const Map > ghostQMap;
653 RCP<MultiVector> ghostQvalues;
654 Array<RCP<Xpetra::Vector<GO,LO,GO,Node> > > ghostQcolumns;
655 RCP<Xpetra::Vector<GO,LO,GO,Node> > ghostQrowNums;
656 ArrayRCP< ArrayRCP<SC> > ghostQvals;
657 ArrayRCP< ArrayRCP<GO> > ghostQcols;
658 ArrayRCP< GO > ghostQrows;
661 for (LO j=0; j<numAggs; ++j) {
662 for (LO k=aggStart[j]; k<aggStart[j+1]; ++k) {
663 if (rowMapForPtentRef.isNodeGlobalElement(aggToRowMap[k]) ==
false) {
664 ghostGIDs.push_back(aggToRowMap[k]);
668 ghostQMap = MapFactory::Build(A->getRowMap()->lib(),
669 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
671 indexBase, A->getRowMap()->getComm());
673 ghostQvalues = MultiVectorFactory::Build(ghostQMap,NSDim);
677 ghostQcolumns.resize(NSDim);
678 for (
size_t i=0; i<NSDim; ++i)
679 ghostQcolumns[i] = Xpetra::VectorFactory<GO,LO,GO,Node>::Build(ghostQMap);
680 ghostQrowNums = Xpetra::VectorFactory<GO,LO,GO,Node>::Build(ghostQMap);
681 if (ghostQvalues->getLocalLength() > 0) {
682 ghostQvals.resize(NSDim);
683 ghostQcols.resize(NSDim);
684 for (
size_t i=0; i<NSDim; ++i) {
685 ghostQvals[i] = ghostQvalues->getDataNonConst(i);
686 ghostQcols[i] = ghostQcolumns[i]->getDataNonConst(0);
688 ghostQrows = ghostQrowNums->getDataNonConst(0);
692 importer = ImportFactory::Build(ghostQMap, A->getRowMap());
695 Teuchos::SerialQRDenseSolver<LO,SC> qrSolver;
698 Array<GO> globalColPtr(maxAggSize*NSDim,0);
699 Array<LO> localColPtr(maxAggSize*NSDim,0);
700 Array<SC> valPtr(maxAggSize*NSDim,0.);
703 const Map& coarseMapRef = *coarseMap;
706 ArrayRCP<size_t> ptent_rowptr;
707 ArrayRCP<LO> ptent_colind;
708 ArrayRCP<Scalar> ptent_values;
711 ArrayView<size_t> rowptr_v;
712 ArrayView<LO> colind_v;
713 ArrayView<Scalar> values_v;
716 Array<size_t> rowptr_temp;
717 Array<LO> colind_temp;
718 Array<Scalar> values_temp;
720 RCP<CrsMatrix> PtentCrs;
722 RCP<CrsMatrixWrap> PtentCrsWrap = rcp(
new CrsMatrixWrap(rowMapForPtent, NSDim));
723 PtentCrs = PtentCrsWrap->getCrsMatrix();
724 Ptentative = PtentCrsWrap;
730 const Map& nonUniqueMapRef = *nonUniqueMap;
732 size_t total_nnz_count=0;
734 for (GO agg=0; agg<numAggs; ++agg)
736 LO myAggSize = aggStart[agg+1]-aggStart[agg];
739 Teuchos::SerialDenseMatrix<LO,SC> localQR(myAggSize, NSDim);
740 for (
size_t j=0; j<NSDim; ++j) {
741 bool bIsZeroNSColumn =
true;
742 for (LO k=0; k<myAggSize; ++k)
747 SC nsVal = fineNS[j][ nonUniqueMapRef.getLocalElement(aggToRowMap[aggStart[agg]+k]) ];
748 localQR(k,j) = nsVal;
749 if (nsVal != zero) bIsZeroNSColumn =
false;
752 GetOStream(
Runtime1,-1) <<
"length of fine level nsp: " << fineNullspace->getGlobalLength() << std::endl;
753 GetOStream(
Runtime1,-1) <<
"length of fine level nsp w overlap: " << fineNullspaceWithOverlap->getGlobalLength() << std::endl;
754 GetOStream(
Runtime1,-1) <<
"(local?) aggId=" << agg << std::endl;
755 GetOStream(
Runtime1,-1) <<
"aggSize=" << myAggSize << std::endl;
756 GetOStream(
Runtime1,-1) <<
"agg DOF=" << k << std::endl;
757 GetOStream(
Runtime1,-1) <<
"NS vector j=" << j << std::endl;
758 GetOStream(
Runtime1,-1) <<
"j*myAggSize + k = " << j*myAggSize + k << std::endl;
759 GetOStream(
Runtime1,-1) <<
"aggToRowMap["<<agg<<
"][" << k <<
"] = " << aggToRowMap[aggStart[agg]+k] << std::endl;
760 GetOStream(
Runtime1,-1) <<
"id aggToRowMap[agg][k]=" << aggToRowMap[aggStart[agg]+k] <<
" is global element in nonUniqueMap = " <<
761 nonUniqueMapRef.isNodeGlobalElement(aggToRowMap[aggStart[agg]+k]) << std::endl;
762 GetOStream(
Runtime1,-1) <<
"colMap local id aggToRowMap[agg][k]=" << nonUniqueMapRef.getLocalElement(aggToRowMap[aggStart[agg]+k]) << std::endl;
763 GetOStream(
Runtime1,-1) <<
"fineNS...=" << fineNS[j][ nonUniqueMapRef.getLocalElement(aggToRowMap[aggStart[agg]+k]) ] << std::endl;
764 GetOStream(
Errors,-1) <<
"caught an error!" << std::endl;
767 TEUCHOS_TEST_FOR_EXCEPTION(bIsZeroNSColumn ==
true,
Exceptions::RuntimeError,
"MueLu::TentativePFactory::MakeTentative: fine level NS part has a zero column. Error.");
770 Xpetra::global_size_t offset=agg*NSDim;
772 if(myAggSize >= Teuchos::as<LocalOrdinal>(NSDim)) {
777 SC tau = localQR(0,0);
782 for (
size_t k = 0; k < Teuchos::as<size_t>(myAggSize); ++k) {
783 Magnitude tmag = STS::magnitude(localQR(k,0));
786 dtemp = Teuchos::ScalarTraits<Magnitude>::squareroot(dtemp);
788 localQR(0,0) = dtemp;
790 qrSolver.setMatrix( Teuchos::rcp(&localQR,
false) );
797 for (
size_t j=0; j<NSDim; ++j) {
798 for (
size_t k=0; k<=j; ++k) {
800 if (coarseMapRef.isNodeLocalElement(offset+k)) {
801 coarseNS[j][offset+k] = localQR(k, j);
805 GetOStream(
Errors,-1) <<
"caught error in coarseNS insert, j="<<j<<
", offset+k = "<<offset+k<<std::endl;
815 Magnitude dtemp = Teuchos::ScalarTraits<SC>::magnitude(localQR(0,0));
819 localQR(i,0) *= dtemp ;
823 Teuchos::RCP<Teuchos::SerialDenseMatrix<LO,SC> > qFactor = qrSolver.getQ();
824 for (
size_t j=0; j<NSDim; j++) {
825 for (
size_t i = 0; i < Teuchos::as<size_t>(myAggSize); i++) {
826 localQR(i,j) = (*qFactor)(i,j);
836 for (
size_t j = 0; j < NSDim; j++)
837 for (
size_t k = 0; k < NSDim; k++) {
839 "Caught error in coarseNS insert, j=" << j <<
", offset+k = " << offset+k);
841 if (k < as<size_t>(myAggSize))
842 coarseNS[j][offset+k] = localQR(k,j);
844 coarseNS[j][offset+k] = (k == j ? one : zero);
848 for (
size_t i = 0; i < as<size_t>(myAggSize); i++)
849 for (
size_t j = 0; j < NSDim; j++)
850 localQR(i,j) = (j == i ? one : zero);
857 for (GO j=0; j<myAggSize; ++j) {
861 GO globalRow = aggToRowMap[aggStart[agg]+j];
864 if (rowMapForPtentRef.isNodeGlobalElement(globalRow) == false ) {
865 ghostQrows[qctr] = globalRow;
866 for (
size_t k=0; k<NSDim; ++k) {
867 ghostQcols[k][qctr] = coarseMapRef.getGlobalElement(agg*NSDim+k);
868 ghostQvals[k][qctr] = localQR(j,k);
873 for (
size_t k=0; k<NSDim; ++k) {
875 if (localQR(j,k) != Teuchos::ScalarTraits<SC>::zero()) {
876 localColPtr[nnz] = agg * NSDim + k;
877 globalColPtr[nnz] = coarseMapRef.getGlobalElement(localColPtr[nnz]);
878 valPtr[nnz] = localQR(j,k);
884 GetOStream(
Errors,-1) <<
"caught error in colPtr/valPtr insert, current index="<<nnz<<std::endl;
889 Ptentative->insertGlobalValues(globalRow,globalColPtr.view(0,nnz),valPtr.view(0,nnz));
892 GetOStream(
Errors,-1) <<
"pid " << A->getRowMap()->getComm()->getRank()
893 <<
"caught error during Ptent row insertion, global row " 894 << globalRow << std::endl;
905 GetOStream(
Runtime1) <<
"TentativePFactory : aggregates may cross process boundaries" << std::endl;
908 RCP<Xpetra::Vector<GO,LO,GO,Node> > targetQrowNums = Xpetra::VectorFactory<GO,LO,GO,Node>::Build(rowMapForPtent);
909 targetQrowNums->putScalar(-1);
910 targetQrowNums->doImport(*ghostQrowNums,*importer,Xpetra::INSERT);
911 ArrayRCP< GO > targetQrows = targetQrowNums->getDataNonConst(0);
914 Array<GO> gidsToImport;
915 gidsToImport.reserve(targetQrows.size());
916 for (
typename ArrayRCP<GO>::iterator r=targetQrows.begin(); r!=targetQrows.end(); ++r) {
918 gidsToImport.push_back(*r);
921 RCP<const Map > reducedMap = MapFactory::Build( A->getRowMap()->lib(),
922 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
923 gidsToImport, indexBase, A->getRowMap()->getComm() );
926 importer = ImportFactory::Build(ghostQMap, reducedMap);
928 Array<RCP<Xpetra::Vector<GO,LO,GO,Node> > > targetQcolumns(NSDim);
929 for (
size_t i=0; i<NSDim; ++i) {
930 targetQcolumns[i] = Xpetra::VectorFactory<GO,LO,GO,Node>::Build(reducedMap);
931 targetQcolumns[i]->doImport(*(ghostQcolumns[i]),*importer,Xpetra::INSERT);
933 RCP<MultiVector> targetQvalues = MultiVectorFactory::Build(reducedMap,NSDim);
934 targetQvalues->doImport(*ghostQvalues,*importer,Xpetra::INSERT);
936 ArrayRCP< ArrayRCP<SC> > targetQvals;
937 ArrayRCP<ArrayRCP<GO> > targetQcols;
938 if (targetQvalues->getLocalLength() > 0) {
939 targetQvals.resize(NSDim);
940 targetQcols.resize(NSDim);
941 for (
size_t i=0; i<NSDim; ++i) {
942 targetQvals[i] = targetQvalues->getDataNonConst(i);
943 targetQcols[i] = targetQcolumns[i]->getDataNonConst(0);
947 valPtr = Array<SC>(NSDim,0.);
948 globalColPtr = Array<GO>(NSDim,0);
949 for (
typename Array<GO>::iterator r=gidsToImport.begin(); r!=gidsToImport.end(); ++r) {
950 if (targetQvalues->getLocalLength() > 0) {
951 for (
size_t j=0; j<NSDim; ++j) {
952 valPtr[j] = targetQvals[j][reducedMap->getLocalElement(*r)];
953 globalColPtr[j] = targetQcols[j][reducedMap->getLocalElement(*r)];
955 Ptentative->insertGlobalValues(*r, globalColPtr.view(0,NSDim), valPtr.view(0,NSDim));
959 Ptentative->fillComplete(coarseMap, A->getDomainMap());
968 #define MUELU_TENTATIVEPFACTORY_SHORT 969 #endif // MUELU_TENTATIVEPFACTORY_DEF_HPP Important warning messages (one line)
MueLu::DefaultLocalOrdinal LocalOrdinal
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)
static bool MapsAreNested(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &rowMap, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &colMap)
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.