46 #ifndef MUELU_SEMICOARSENPFACTORY_DEF_HPP 47 #define MUELU_SEMICOARSENPFACTORY_DEF_HPP 51 #include <Teuchos_LAPACK.hpp> 53 #include <Xpetra_CrsMatrixWrap.hpp> 54 #include <Xpetra_ImportFactory.hpp> 55 #include <Xpetra_Matrix.hpp> 56 #include <Xpetra_MultiVectorFactory.hpp> 57 #include <Xpetra_VectorFactory.hpp> 67 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
69 RCP<ParameterList> validParamList = rcp(
new ParameterList());
71 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 74 #undef SET_VALID_ENTRY 75 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A");
76 validParamList->set< RCP<const FactoryBase> >(
"Nullspace", Teuchos::null,
"Generating factory of the nullspace");
77 validParamList->set< RCP<const FactoryBase> >(
"Coordinates", Teuchos::null,
"Generating factory for coordinates");
79 validParamList->set< RCP<const FactoryBase> >(
"LineDetection_VertLineIds", Teuchos::null,
"Generating factory for LineDetection vertical line ids");
80 validParamList->set< RCP<const FactoryBase> >(
"LineDetection_Layers", Teuchos::null,
"Generating factory for LineDetection layer ids");
81 validParamList->set< RCP<const FactoryBase> >(
"CoarseNumZLayers", Teuchos::null,
"Generating factory for number of coarse z-layers");
83 return validParamList;
86 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
88 Input(fineLevel,
"A");
89 Input(fineLevel,
"Nullspace");
91 Input(fineLevel,
"LineDetection_VertLineIds");
92 Input(fineLevel,
"LineDetection_Layers");
93 Input(fineLevel,
"CoarseNumZLayers");
101 bTransferCoordinates_ =
true;
103 }
else if (bTransferCoordinates_ ==
true){
107 RCP<const FactoryBase> myCoordsFact = GetFactory(
"Coordinates");
108 if (myCoordsFact == Teuchos::null) { myCoordsFact = fineLevel.
GetFactoryManager()->GetFactory(
"Coordinates"); }
109 if (fineLevel.
IsAvailable(
"Coordinates", myCoordsFact.get())) {
110 fineLevel.
DeclareInput(
"Coordinates", myCoordsFact.get(),
this);
115 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
117 return BuildP(fineLevel, coarseLevel);
120 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
125 RCP<Matrix> A = Get< RCP<Matrix> > (fineLevel,
"A");
126 RCP<MultiVector> fineNullspace = Get< RCP<MultiVector> > (fineLevel,
"Nullspace");
129 const ParameterList& pL = GetParameterList();
130 LO CoarsenRate = as<LO>(pL.get<
int>(
"semicoarsen: coarsen rate"));
133 LO BlkSize = A->GetFixedBlockSize();
134 RCP<const Map> rowMap = A->getRowMap();
135 LO Ndofs = rowMap->getNodeNumElements();
136 LO Nnodes = Ndofs/BlkSize;
139 LO FineNumZLayers = Get< LO >(fineLevel,
"CoarseNumZLayers");
140 Teuchos::ArrayRCP<LO> TVertLineId = Get< Teuchos::ArrayRCP<LO> > (fineLevel,
"LineDetection_VertLineIds");
141 Teuchos::ArrayRCP<LO> TLayerId = Get< Teuchos::ArrayRCP<LO> > (fineLevel,
"LineDetection_Layers");
142 LO* VertLineId = TVertLineId.getRawPtr();
143 LO* LayerId = TLayerId.getRawPtr();
146 RCP<const Map> theCoarseMap;
148 RCP<MultiVector> coarseNullspace;
149 GO Ncoarse = MakeSemiCoarsenP(Nnodes,FineNumZLayers,CoarsenRate,LayerId,VertLineId,
150 BlkSize, A, P, theCoarseMap, fineNullspace,coarseNullspace);
153 if (A->IsView(
"stridedMaps") ==
true)
154 P->CreateView(
"stridedMaps", A->getRowMap(
"stridedMaps"), theCoarseMap);
156 P->CreateView(
"stridedMaps", P->getRangeMap(), theCoarseMap);
158 if (pL.get<
bool>(
"semicoarsen: piecewise constant"))
159 RevertToPieceWiseConstant(P, BlkSize);
165 LO CoarseNumZLayers = FineNumZLayers*Ncoarse;
166 CoarseNumZLayers /= Ndofs;
170 Set(coarseLevel,
"P", P);
172 Set(coarseLevel,
"Nullspace", coarseNullspace);
175 if(bTransferCoordinates_) {
176 typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO> xdMV;
177 RCP<xdMV> fineCoords = Teuchos::null;
178 if (fineLevel.GetLevelID() == 0 &&
180 fineCoords = fineLevel.Get< RCP<xdMV> >(
"Coordinates",
NoFactory::get());
182 RCP<const FactoryBase> myCoordsFact = GetFactory(
"Coordinates");
183 if (myCoordsFact == Teuchos::null) { myCoordsFact = fineLevel.GetFactoryManager()->GetFactory(
"Coordinates"); }
184 if (fineLevel.IsAvailable(
"Coordinates", myCoordsFact.get())) {
185 fineCoords = fineLevel.Get< RCP<xdMV> >(
"Coordinates", myCoordsFact.get());
189 TEUCHOS_TEST_FOR_EXCEPTION(fineCoords==Teuchos::null,
Exceptions::RuntimeError,
"No Coordinates found provided by the user.");
191 TEUCHOS_TEST_FOR_EXCEPTION(fineCoords->getNumVectors() != 3,
Exceptions::RuntimeError,
"Three coordinates arrays must be supplied if line detection orientation not given.");
192 ArrayRCP<typename Teuchos::ScalarTraits<Scalar>::coordinateType> x = fineCoords->getDataNonConst(0);
193 ArrayRCP<typename Teuchos::ScalarTraits<Scalar>::coordinateType> y = fineCoords->getDataNonConst(1);
194 ArrayRCP<typename Teuchos::ScalarTraits<Scalar>::coordinateType> z = fineCoords->getDataNonConst(2);
197 typename Teuchos::ScalarTraits<Scalar>::coordinateType zval_max = -Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<Scalar>::coordinateType>::one() / Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<Scalar>::coordinateType>::sfmin();
198 typename Teuchos::ScalarTraits<Scalar>::coordinateType zval_min = Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<Scalar>::coordinateType>::one() / Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<Scalar>::coordinateType>::sfmin();
199 for (
auto it = z.begin(); it != z.end(); ++it) {
200 if(*it > zval_max) zval_max = *it;
201 if(*it < zval_min) zval_min = *it;
204 LO myCoarseZLayers = Teuchos::as<LO>(CoarseNumZLayers);
206 ArrayRCP<typename Teuchos::ScalarTraits<Scalar>::coordinateType> myZLayerCoords = Teuchos::arcp<typename Teuchos::ScalarTraits<Scalar>::coordinateType>(myCoarseZLayers);
207 if(myCoarseZLayers == 1) {
208 myZLayerCoords[0] = zval_min;
210 typename Teuchos::ScalarTraits<Scalar>::coordinateType dz = (zval_max-zval_min)/(myCoarseZLayers-1);
211 for(LO k = 0; k<myCoarseZLayers; ++k) {
212 myZLayerCoords[k] = k*dz;
220 LO numVertLines = Nnodes / FineNumZLayers;
221 LO numLocalCoarseNodes = numVertLines * myCoarseZLayers;
223 RCP<const Map> coarseCoordMap =
224 MapFactory::Build (fineCoords->getMap()->lib(),
225 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
226 Teuchos::as<size_t>(numLocalCoarseNodes),
227 fineCoords->getMap()->getIndexBase(),
228 fineCoords->getMap()->getComm());
229 RCP<xdMV> coarseCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO>::Build(coarseCoordMap, fineCoords->getNumVectors());
230 coarseCoords->putScalar(-1.0);
231 ArrayRCP<typename Teuchos::ScalarTraits<Scalar>::coordinateType> cx = coarseCoords->getDataNonConst(0);
232 ArrayRCP<typename Teuchos::ScalarTraits<Scalar>::coordinateType> cy = coarseCoords->getDataNonConst(1);
233 ArrayRCP<typename Teuchos::ScalarTraits<Scalar>::coordinateType> cz = coarseCoords->getDataNonConst(2);
236 LO cntCoarseNodes = 0;
237 for( LO vt = 0; vt < TVertLineId.size(); ++vt) {
239 LO curVertLineId = TVertLineId[vt];
241 if(cx[curVertLineId * myCoarseZLayers] == -1.0 &&
242 cy[curVertLineId * myCoarseZLayers] == -1.0) {
244 for (LO n=0; n<myCoarseZLayers; ++n) {
245 cx[curVertLineId * myCoarseZLayers + n] = x[vt];
246 cy[curVertLineId * myCoarseZLayers + n] = y[vt];
247 cz[curVertLineId * myCoarseZLayers + n] = myZLayerCoords[n];
249 cntCoarseNodes += myCoarseZLayers;
252 TEUCHOS_TEST_FOR_EXCEPTION(cntCoarseNodes > numLocalCoarseNodes,
Exceptions::RuntimeError,
"number of coarse nodes is inconsistent.");
253 if(cntCoarseNodes == numLocalCoarseNodes)
break;
257 Set(coarseLevel,
"Coordinates", coarseCoords);
262 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
286 typename Teuchos::ScalarTraits<Scalar>::coordinateType temp, RestStride, di;
291 temp = ((
typename Teuchos::ScalarTraits<Scalar>::coordinateType) (PtsPerLine+1))/((
typename Teuchos::ScalarTraits<Scalar>::coordinateType) (CoarsenRate)) - 1.0;
292 if (Thin == 1) NCpts = (LO) ceil(temp);
293 else NCpts = (LO) floor(temp);
295 TEUCHOS_TEST_FOR_EXCEPTION(PtsPerLine == 1,
Exceptions::RuntimeError,
"SemiCoarsenPFactory::FindCpts: cannot coarsen further.");
297 if (NCpts < 1) NCpts = 1;
299 FirstStride= (LO) ceil( ((
typename Teuchos::ScalarTraits<Scalar>::coordinateType) PtsPerLine+1)/( (
typename Teuchos::ScalarTraits<Scalar>::coordinateType) (NCpts+1)));
300 RestStride = ((
typename Teuchos::ScalarTraits<Scalar>::coordinateType) (PtsPerLine-FirstStride+1))/((
typename Teuchos::ScalarTraits<Scalar>::coordinateType) NCpts);
302 NCLayers = (LO) floor((((
typename Teuchos::ScalarTraits<Scalar>::coordinateType) (PtsPerLine-FirstStride+1))/RestStride)+.00001);
306 di = (
typename Teuchos::ScalarTraits<Scalar>::coordinateType) FirstStride;
307 for (i = 1; i <= NCpts; i++) {
308 (*LayerCpts)[i] = (LO) floor(di);
315 #define MaxHorNeighborNodes 75 317 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
320 LO
const VertLineId[], LO
const DofsPerNode, RCP<Matrix> & Amat, RCP<Matrix>& P,
321 RCP<const Map>& coarseMap,
const RCP<MultiVector> fineNullspace,
322 RCP<MultiVector>& coarseNullspace )
const {
374 int NLayers, NVertLines, MaxNnz, NCLayers, MyLine, MyLayer;
375 int *InvLineLayer=NULL, *CptLayers=NULL, StartLayer, NStencilNodes;
376 int BlkRow, dof_j, node_k, *Sub2FullMap=NULL, RowLeng;
377 int i, j, iii, col, count, index, loc, PtRow, PtCol;
378 SC *BandSol=NULL, *BandMat=NULL, TheSum;
379 int *IPIV=NULL, KL, KU, KLU, N, NRHS, LDAB,INFO;
383 int MaxStencilSize, MaxNnzPerRow;
385 int CurRow, LastGuy = -1, NewPtr;
388 LO *Layerdofs = NULL, *Col2Dof = NULL;
390 Teuchos::LAPACK<LO,SC> lapack;
398 Teuchos::ArrayRCP<LO> TLayDiff = Teuchos::arcp<LO>(1+MaxNnzPerRow); LayDiff = TLayDiff.getRawPtr();
400 Nghost = Amat->getColMap()->getNodeNumElements() - Amat->getDomainMap()->getNodeNumElements();
401 if (Nghost < 0) Nghost = 0;
402 Teuchos::ArrayRCP<LO> TLayerdofs= Teuchos::arcp<LO>(Ntotal*DofsPerNode+Nghost+1); Layerdofs = TLayerdofs.getRawPtr();
403 Teuchos::ArrayRCP<LO> TCol2Dof= Teuchos::arcp<LO>(Ntotal*DofsPerNode+Nghost+1); Col2Dof= TCol2Dof.getRawPtr();
405 RCP<Xpetra::Vector<LO,LO,GO,NO> > localdtemp = Xpetra::VectorFactory<LO,LO,GO,NO>::Build(Amat->getDomainMap());
406 RCP<Xpetra::Vector<LO,LO,GO,NO> > dtemp = Xpetra::VectorFactory<LO,LO,GO,NO>::Build(Amat->getColMap());
407 ArrayRCP<LO> valptr= localdtemp->getDataNonConst(0);
409 for (i = 0; i < Ntotal*DofsPerNode; i++)
410 valptr[i]= LayerId[i/DofsPerNode];
411 valptr=ArrayRCP<LO>();
413 RCP< const Import> importer;
414 importer = Amat->getCrsGraph()->getImporter();
415 if (importer == Teuchos::null) {
416 importer = ImportFactory::Build(Amat->getDomainMap(), Amat->getColMap());
418 dtemp->doImport(*localdtemp, *(importer), Xpetra::INSERT);
420 valptr= dtemp->getDataNonConst(0);
421 for (i = 0; i < Ntotal*DofsPerNode+Nghost; i++) Layerdofs[i]= valptr[i];
422 valptr= localdtemp->getDataNonConst(0);
423 for (i = 0; i < Ntotal*DofsPerNode; i++) valptr[i]= i%DofsPerNode;
424 valptr=ArrayRCP<LO>();
425 dtemp->doImport(*localdtemp, *(importer), Xpetra::INSERT);
427 valptr= dtemp->getDataNonConst(0);
428 for (i = 0; i < Ntotal*DofsPerNode+Nghost; i++) Col2Dof[i]= valptr[i];
429 valptr=ArrayRCP<LO>();
432 NLayers = LayerId[0];
433 NVertLines= VertLineId[0];
435 else { NLayers = -1; NVertLines = -1; }
437 for (i = 1; i < Ntotal; i++) {
438 if ( VertLineId[i] > NVertLines ) NVertLines = VertLineId[i];
439 if ( LayerId[i] > NLayers ) NLayers = LayerId[i];
449 Teuchos::ArrayRCP<LO> TInvLineLayer= Teuchos::arcp<LO>(1+NVertLines*NLayers); InvLineLayer = TInvLineLayer.getRawPtr();
450 for (i=0; i < Ntotal; i++) {
451 InvLineLayer[ VertLineId[i]+1+LayerId[i]*NVertLines ] = i;
458 Teuchos::ArrayRCP<LO> TCptLayers = Teuchos::arcp<LO>(nz+1); CptLayers = TCptLayers.getRawPtr();
459 NCLayers = FindCpts(nz,CoarsenRate,0, &CptLayers);
468 if (NCLayers < 2) MaxStencilSize = nz;
469 else MaxStencilSize = CptLayers[2];
471 for (i = 3; i <= NCLayers; i++) {
472 if (MaxStencilSize < CptLayers[i]- CptLayers[i-2])
473 MaxStencilSize = CptLayers[i]- CptLayers[i-2];
476 if (MaxStencilSize < nz - CptLayers[NCLayers-1]+1)
477 MaxStencilSize = nz - CptLayers[NCLayers-1]+1;
487 Teuchos::ArrayRCP<LO> TSub2FullMap= Teuchos::arcp<LO>((MaxStencilSize+1)*DofsPerNode); Sub2FullMap= TSub2FullMap.getRawPtr();
488 Teuchos::ArrayRCP<SC> TBandSol= Teuchos::arcp<SC>((MaxStencilSize+1)*DofsPerNode*DofsPerNode); BandSol = TBandSol.getRawPtr();
492 KL = 2*DofsPerNode-1;
493 KU = 2*DofsPerNode-1;
497 Teuchos::ArrayRCP<SC> TBandMat= Teuchos::arcp<SC>(LDAB*MaxStencilSize*DofsPerNode+1); BandMat = TBandMat.getRawPtr();
498 Teuchos::ArrayRCP<LO> TIPIV= Teuchos::arcp<LO>((MaxStencilSize+1)*DofsPerNode); IPIV = TIPIV.getRawPtr();
507 Ndofs = DofsPerNode*Ntotal;
508 MaxNnz = 2*DofsPerNode*Ndofs;
510 RCP<const Map> rowMap = Amat->getRowMap();
511 Xpetra::global_size_t GNdofs= rowMap->getGlobalNumElements();
513 std::vector<size_t> stridingInfo_;
514 stridingInfo_.push_back(DofsPerNode);
516 Xpetra::global_size_t itemp = GNdofs/nz;
517 coarseMap = StridedMapFactory::Build(rowMap->lib(),
519 NCLayers*NVertLines*DofsPerNode,
530 P = rcp(
new CrsMatrixWrap(rowMap, coarseMap , 0));
531 coarseNullspace = MultiVectorFactory::Build(coarseMap, fineNullspace->getNumVectors());
534 Teuchos::ArrayRCP<SC> TPvals= Teuchos::arcp<SC>(1+MaxNnz); Pvals= TPvals.getRawPtr();
535 Teuchos::ArrayRCP<size_t> TPptr = Teuchos::arcp<size_t>(DofsPerNode*(2+Ntotal)); Pptr = TPptr.getRawPtr();
536 Teuchos::ArrayRCP<LO> TPcols= Teuchos::arcp<LO>(1+MaxNnz); Pcols= TPcols.getRawPtr();
538 Pptr[0] = 0; Pptr[1] = 0;
548 for (i = 1; i <= MaxNnz; i++) Pcols[i] = -1;
550 for (i = 1; i <= DofsPerNode*Ntotal+1; i++) {
552 count += (2*DofsPerNode);
564 for (MyLine=1; MyLine <= NVertLines; MyLine += 1) {
565 for (iii=1; iii <= NCLayers; iii+= 1) {
567 MyLayer = CptLayers[iii];
580 if (iii != 1 ) StartLayer = CptLayers[iii-1]+1;
583 if (iii != NCLayers) NStencilNodes= CptLayers[iii+1]-StartLayer;
584 else NStencilNodes= NLayers - StartLayer + 1;
587 N = NStencilNodes*DofsPerNode;
594 for (i = 0; i < NStencilNodes*DofsPerNode*DofsPerNode; i++)
596 for (i = 0; i < LDAB*N; i++) BandMat[ i] = 0.0;
603 for (node_k=1; node_k <= NStencilNodes ; node_k++) {
609 BlkRow = InvLineLayer[MyLine+(StartLayer+node_k-2)*NVertLines]+1;
610 Sub2FullMap[node_k] = BlkRow;
623 if (StartLayer+node_k-1 != MyLayer) {
624 for (
int dof_i=0; dof_i < DofsPerNode; dof_i++) {
626 j = (BlkRow-1)*DofsPerNode+dof_i;
627 ArrayView<const LO> AAcols;
628 ArrayView<const SC> AAvals;
629 Amat->getLocalRowView(j, AAcols, AAvals);
630 const int *Acols = AAcols.getRawPtr();
631 const SC *Avals = AAvals.getRawPtr();
632 RowLeng = AAvals.size();
634 TEUCHOS_TEST_FOR_EXCEPTION(RowLeng >= MaxNnzPerRow,
Exceptions::RuntimeError,
"MakeSemiCoarsenP: recompile with larger Max(HorNeighborNodes)\n");
636 for (i = 0; i < RowLeng; i++) {
637 LayDiff[i] = Layerdofs[Acols[i]]-StartLayer-node_k+2;
645 PtRow = (node_k-1)*DofsPerNode+dof_i+1;
646 for (dof_j=0; dof_j < DofsPerNode; dof_j++) {
647 PtCol = (node_k-1)*DofsPerNode+dof_j + 1;
651 for (i = 0; i < RowLeng; i++) {
652 if ((LayDiff[i] == 0) && (Col2Dof[Acols[i]] == dof_j))
655 index = LDAB*(PtCol-1)+KLU+PtRow-PtCol;
656 BandMat[index] = TheSum;
657 if (node_k != NStencilNodes) {
661 for (i = 0; i < RowLeng; i++) {
662 if ((LayDiff[i] == 1) &&(Col2Dof[Acols[i]]== dof_j))
665 j = PtCol+DofsPerNode;
666 index=LDAB*(j-1)+KLU+PtRow-j;
667 BandMat[index] = TheSum;
673 for (i = 0; i < RowLeng; i++) {
674 if ((LayDiff[i]== -1) &&(Col2Dof[Acols[i]]== dof_j))
677 j = PtCol-DofsPerNode;
678 index=LDAB*(j-1)+KLU+PtRow-j;
679 BandMat[index] = TheSum;
689 for (
int k = 0; k < static_cast<int>(fineNullspace->getNumVectors()); k++) {
690 Teuchos::ArrayRCP<SC> OneCNull = coarseNullspace->getDataNonConst(k);
691 Teuchos::ArrayRCP<SC> OneFNull = fineNullspace->getDataNonConst(k);
692 for (
int dof_i = 0; dof_i < DofsPerNode; dof_i++) {
693 OneCNull[( col-1)*DofsPerNode+dof_i] = OneFNull[ (BlkRow-1)*DofsPerNode+dof_i];
697 for (
int dof_i = 0; dof_i < DofsPerNode; dof_i++) {
700 PtRow = (node_k-1)*DofsPerNode+dof_i+1;
701 index = LDAB*(PtRow-1)+KLU;
702 BandMat[index] = 1.0;
703 BandSol[(dof_i)*DofsPerNode*NStencilNodes+PtRow-1] = 1.;
710 lapack.GBTRF( N, N, KL, KU, BandMat, LDAB, IPIV, &INFO);
714 lapack.GBTRS(notrans[0], N, KL, KU, NRHS, BandMat, LDAB, IPIV,
719 for (dof_j=0; dof_j < DofsPerNode; dof_j++) {
720 for (
int dof_i=0; dof_i < DofsPerNode; dof_i++) {
721 for (i =1; i <= NStencilNodes ; i++) {
722 index = (Sub2FullMap[i]-1)*DofsPerNode+dof_i+1;
724 Pcols[loc] = (col-1)*DofsPerNode+dof_j+1;
725 Pvals[loc] = BandSol[dof_j*DofsPerNode*NStencilNodes +
726 (i-1)*DofsPerNode + dof_i];
727 Pptr[index]= Pptr[index] + 1;
743 for (
size_t ii=1; ii <= Pptr[Ntotal*DofsPerNode]-1; ii++) {
744 if (ii == Pptr[CurRow]) {
745 Pptr[CurRow] = LastGuy;
747 while (ii > Pptr[CurRow]) {
748 Pptr[CurRow] = LastGuy;
752 if (Pcols[ii] != -1) {
753 Pcols[NewPtr-1] = Pcols[ii]-1;
754 Pvals[NewPtr-1] = Pvals[ii];
759 for (i = CurRow; i <= Ntotal*DofsPerNode; i++) Pptr[CurRow] = LastGuy;
764 for (i=-Ntotal*DofsPerNode+1; i>= 2 ; i--) {
765 Pptr[i-1] = Pptr[i-2];
769 ArrayRCP<size_t> rcpRowPtr;
770 ArrayRCP<LO> rcpColumns;
771 ArrayRCP<SC> rcpValues;
773 RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
774 LO nnz = Pptr[Ndofs];
775 PCrs->allocateAllValues(nnz, rcpRowPtr, rcpColumns, rcpValues);
777 ArrayView<size_t> rowPtr = rcpRowPtr();
778 ArrayView<LO> columns = rcpColumns();
779 ArrayView<SC> values = rcpValues();
783 for (LO ii = 0; ii <= Ndofs; ii++) rowPtr[ii] = Pptr[ii];
784 for (LO ii = 0; ii < nnz; ii++) columns[ii] = Pcols[ii];
785 for (LO ii = 0; ii < nnz; ii++) values[ii] = Pvals[ii];
786 PCrs->setAllValues(rcpRowPtr, rcpColumns, rcpValues);
787 PCrs->expertStaticFillComplete(coarseMap, Amat->getDomainMap());
790 return NCLayers*NVertLines*DofsPerNode;
792 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
800 ArrayView<const LocalOrdinal> inds;
801 ArrayView<const Scalar> vals1;
802 ArrayView< Scalar> vals;
803 Scalar ZERO = Teuchos::ScalarTraits<Scalar>::zero();
804 Scalar ONE = Teuchos::ScalarTraits<Scalar>::one();
806 for (
size_t i = 0; i < P->getRowMap()->getNodeNumElements(); i++) {
807 P->getLocalRowView(i, inds, vals1);
809 size_t nnz = inds.size();
810 if (nnz != 0) vals = ArrayView<Scalar>(
const_cast<Scalar*
>(vals1.getRawPtr()), nnz);
812 LO largestIndex = -1;
813 Scalar largestValue = ZERO;
816 LO rowDof = i%BlkSize;
817 for (
size_t j =0; j < nnz; j++) {
818 if (Teuchos::ScalarTraits<SC>::magnitude(vals[ j ]) >= Teuchos::ScalarTraits<SC>::magnitude(largestValue)) {
819 if ( inds[j]%BlkSize == rowDof ) {
820 largestValue = vals[j];
821 largestIndex = (int) j;
826 if (largestIndex != -1) vals[largestIndex] = ONE;
828 TEUCHOS_TEST_FOR_EXCEPTION(nnz > 0,
Exceptions::RuntimeError,
"no nonzero column associated with a proper dof within node.");
830 if (Teuchos::ScalarTraits<SC>::magnitude(largestValue) == Teuchos::ScalarTraits<SC>::magnitude(ZERO)) vals[largestIndex] = ZERO;
836 #define MUELU_SEMICOARSENPFACTORY_SHORT 837 #endif // MUELU_SEMICOARSENPFACTORY_DEF_HPP MueLu::DefaultLocalOrdinal LocalOrdinal
void RevertToPieceWiseConstant(RCP< Matrix > &P, LO BlkSize) const
LO FindCpts(LO const PtsPerLine, LO const CoarsenRate, LO const Thin, LO **LayerCpts) const
Timer to be used in factories. Similar to Monitor but with additional timers.
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
Namespace for MueLu classes and methods.
static const NoFactory * get()
int GetLevelID() const
Return level number.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
MueLu::DefaultScalar Scalar
Class that holds all level-specific information.
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
#define SET_VALID_ENTRY(name)
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
#define MaxHorNeighborNodes
const RCP< const FactoryManagerBase > GetFactoryManager()
returns the current factory manager
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.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
LO MakeSemiCoarsenP(LO const Ntotal, LO const nz, LO const CoarsenRate, LO const LayerId[], LO const VertLineId[], LO const DofsPerNode, RCP< Matrix > &Amat, RCP< Matrix > &P, RCP< const Map > &coarseMap, const RCP< MultiVector > fineNullspace, RCP< MultiVector > &coarseNullspace) const
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.