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))
77 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A");
78 validParamList->set< RCP<const FactoryBase> >(
"Nullspace", Teuchos::null,
"Generating factory of the nullspace");
79 validParamList->set< RCP<const FactoryBase> >(
"Coordinates", Teuchos::null,
"Generating factory for coordinates");
81 validParamList->set< RCP<const FactoryBase> >(
"LineDetection_VertLineIds", Teuchos::null,
"Generating factory for LineDetection vertical line ids");
82 validParamList->set< RCP<const FactoryBase> >(
"LineDetection_Layers", Teuchos::null,
"Generating factory for LineDetection layer ids");
83 validParamList->set< RCP<const FactoryBase> >(
"CoarseNumZLayers", Teuchos::null,
"Generating factory for number of coarse z-layers");
85 return validParamList;
88 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
90 Input(fineLevel,
"A");
91 Input(fineLevel,
"Nullspace");
93 Input(fineLevel,
"LineDetection_VertLineIds");
94 Input(fineLevel,
"LineDetection_Layers");
95 Input(fineLevel,
"CoarseNumZLayers");
103 bTransferCoordinates_ =
true;
105 }
else if (bTransferCoordinates_ ==
true){
109 RCP<const FactoryBase> myCoordsFact = GetFactory(
"Coordinates");
110 if (myCoordsFact == Teuchos::null) { myCoordsFact = fineLevel.
GetFactoryManager()->GetFactory(
"Coordinates"); }
111 if (fineLevel.
IsAvailable(
"Coordinates", myCoordsFact.get())) {
112 fineLevel.
DeclareInput(
"Coordinates", myCoordsFact.get(),
this);
117 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
119 return BuildP(fineLevel, coarseLevel);
122 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
127 RCP<Matrix> A = Get< RCP<Matrix> > (fineLevel,
"A");
128 RCP<MultiVector> fineNullspace = Get< RCP<MultiVector> > (fineLevel,
"Nullspace");
131 const ParameterList& pL = GetParameterList();
132 LO CoarsenRate = as<LO>(pL.get<
int>(
"semicoarsen: coarsen rate"));
133 bool buildRestriction = pL.get<
bool>(
"semicoarsen: calculate nonsym restriction");
134 bool doLinear = pL.get<
bool>(
"semicoarsen: piecewise linear");
137 LO BlkSize = A->GetFixedBlockSize();
138 RCP<const Map> rowMap = A->getRowMap();
139 LO Ndofs = rowMap->getLocalNumElements();
140 LO Nnodes = Ndofs/BlkSize;
143 LO FineNumZLayers = Get< LO >(fineLevel,
"CoarseNumZLayers");
144 Teuchos::ArrayRCP<LO> TVertLineId = Get< Teuchos::ArrayRCP<LO> > (fineLevel,
"LineDetection_VertLineIds");
145 Teuchos::ArrayRCP<LO> TLayerId = Get< Teuchos::ArrayRCP<LO> > (fineLevel,
"LineDetection_Layers");
146 LO* VertLineId = TVertLineId.getRawPtr();
147 LO* LayerId = TLayerId.getRawPtr();
150 RCP<const Map> theCoarseMap;
152 RCP<MultiVector> coarseNullspace;
153 GO Ncoarse = MakeSemiCoarsenP(Nnodes,FineNumZLayers,CoarsenRate,LayerId,VertLineId,
154 BlkSize, A, P, theCoarseMap, fineNullspace,coarseNullspace,R,buildRestriction,doLinear);
157 if (A->IsView(
"stridedMaps") ==
true)
158 P->CreateView(
"stridedMaps", A->getRowMap(
"stridedMaps"), theCoarseMap);
160 P->CreateView(
"stridedMaps", P->getRangeMap(), theCoarseMap);
162 if (buildRestriction) {
163 if (A->IsView(
"stridedMaps") ==
true)
164 R->CreateView(
"stridedMaps", theCoarseMap, A->getRowMap(
"stridedMaps"));
166 R->CreateView(
"stridedMaps", theCoarseMap, R->getDomainMap());
168 if (pL.get<
bool>(
"semicoarsen: piecewise constant")) {
169 TEUCHOS_TEST_FOR_EXCEPTION(buildRestriction,
Exceptions::RuntimeError,
"Cannot use calculate nonsym restriction with piecewise constant.");
170 RevertToPieceWiseConstant(P, BlkSize);
172 if (pL.get<
bool>(
"semicoarsen: piecewise linear")) {
173 TEUCHOS_TEST_FOR_EXCEPTION(buildRestriction,
Exceptions::RuntimeError,
"Cannot use calculate nonsym restriction with piecewise linear.");
174 TEUCHOS_TEST_FOR_EXCEPTION(pL.get<
bool>(
"semicoarsen: piecewise constant"),
Exceptions::RuntimeError,
"Cannot use piecewise constant with piecewise linear.");
181 LO CoarseNumZLayers = FineNumZLayers*Ncoarse;
182 CoarseNumZLayers /= Ndofs;
186 Set(coarseLevel,
"P", P);
187 if (buildRestriction) Set(coarseLevel,
"RfromPfactory", R);
189 Set(coarseLevel,
"Nullspace", coarseNullspace);
192 if(bTransferCoordinates_) {
193 typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO> xdMV;
194 RCP<xdMV> fineCoords = Teuchos::null;
199 RCP<const FactoryBase> myCoordsFact = GetFactory(
"Coordinates");
200 if (myCoordsFact == Teuchos::null) { myCoordsFact = fineLevel.
GetFactoryManager()->GetFactory(
"Coordinates"); }
201 if (fineLevel.
IsAvailable(
"Coordinates", myCoordsFact.get())) {
202 fineCoords = fineLevel.
Get< RCP<xdMV> >(
"Coordinates", myCoordsFact.get());
206 TEUCHOS_TEST_FOR_EXCEPTION(fineCoords==Teuchos::null,
Exceptions::RuntimeError,
"No Coordinates found provided by the user.");
208 TEUCHOS_TEST_FOR_EXCEPTION(fineCoords->getNumVectors() != 3,
Exceptions::RuntimeError,
"Three coordinates arrays must be supplied if line detection orientation not given.");
209 ArrayRCP<typename Teuchos::ScalarTraits<Scalar>::coordinateType> x = fineCoords->getDataNonConst(0);
210 ArrayRCP<typename Teuchos::ScalarTraits<Scalar>::coordinateType> y = fineCoords->getDataNonConst(1);
211 ArrayRCP<typename Teuchos::ScalarTraits<Scalar>::coordinateType> z = fineCoords->getDataNonConst(2);
214 typename Teuchos::ScalarTraits<Scalar>::coordinateType zval_max = -Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<Scalar>::coordinateType>::one() / Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<Scalar>::coordinateType>::sfmin();
215 typename Teuchos::ScalarTraits<Scalar>::coordinateType zval_min = Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<Scalar>::coordinateType>::one() / Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<Scalar>::coordinateType>::sfmin();
216 for (
auto it = z.begin(); it != z.end(); ++it) {
217 if(*it > zval_max) zval_max = *it;
218 if(*it < zval_min) zval_min = *it;
221 LO myCoarseZLayers = Teuchos::as<LO>(CoarseNumZLayers);
223 ArrayRCP<typename Teuchos::ScalarTraits<Scalar>::coordinateType> myZLayerCoords = Teuchos::arcp<typename Teuchos::ScalarTraits<Scalar>::coordinateType>(myCoarseZLayers);
224 if(myCoarseZLayers == 1) {
225 myZLayerCoords[0] = zval_min;
227 typename Teuchos::ScalarTraits<Scalar>::coordinateType dz = (zval_max-zval_min)/(myCoarseZLayers-1);
228 for(LO k = 0; k<myCoarseZLayers; ++k) {
229 myZLayerCoords[k] = k*dz;
237 LO numVertLines = Nnodes / FineNumZLayers;
238 LO numLocalCoarseNodes = numVertLines * myCoarseZLayers;
240 RCP<const Map> coarseCoordMap =
241 MapFactory::Build (fineCoords->getMap()->lib(),
242 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
243 Teuchos::as<size_t>(numLocalCoarseNodes),
244 fineCoords->getMap()->getIndexBase(),
245 fineCoords->getMap()->getComm());
246 RCP<xdMV> coarseCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO>::Build(coarseCoordMap, fineCoords->getNumVectors());
247 coarseCoords->putScalar(-1.0);
248 ArrayRCP<typename Teuchos::ScalarTraits<Scalar>::coordinateType> cx = coarseCoords->getDataNonConst(0);
249 ArrayRCP<typename Teuchos::ScalarTraits<Scalar>::coordinateType> cy = coarseCoords->getDataNonConst(1);
250 ArrayRCP<typename Teuchos::ScalarTraits<Scalar>::coordinateType> cz = coarseCoords->getDataNonConst(2);
253 LO cntCoarseNodes = 0;
254 for( LO vt = 0; vt < TVertLineId.size(); ++vt) {
256 LO curVertLineId = TVertLineId[vt];
258 if(cx[curVertLineId * myCoarseZLayers] == -1.0 &&
259 cy[curVertLineId * myCoarseZLayers] == -1.0) {
261 for (LO n=0; n<myCoarseZLayers; ++n) {
262 cx[curVertLineId * myCoarseZLayers + n] = x[vt];
263 cy[curVertLineId * myCoarseZLayers + n] = y[vt];
264 cz[curVertLineId * myCoarseZLayers + n] = myZLayerCoords[n];
266 cntCoarseNodes += myCoarseZLayers;
269 TEUCHOS_TEST_FOR_EXCEPTION(cntCoarseNodes > numLocalCoarseNodes,
Exceptions::RuntimeError,
"number of coarse nodes is inconsistent.");
270 if(cntCoarseNodes == numLocalCoarseNodes)
break;
274 Set(coarseLevel,
"Coordinates", coarseCoords);
279 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
303 typename Teuchos::ScalarTraits<Scalar>::coordinateType temp, RestStride, di;
308 temp = ((
typename Teuchos::ScalarTraits<Scalar>::coordinateType) (PtsPerLine+1))/((
typename Teuchos::ScalarTraits<Scalar>::coordinateType) (CoarsenRate)) - 1.0;
309 if (Thin == 1) NCpts = (LO) ceil(temp);
310 else NCpts = (LO) floor(temp);
312 TEUCHOS_TEST_FOR_EXCEPTION(PtsPerLine == 1,
Exceptions::RuntimeError,
"SemiCoarsenPFactory::FindCpts: cannot coarsen further.");
314 if (NCpts < 1) NCpts = 1;
316 FirstStride= (LO) ceil( ((
typename Teuchos::ScalarTraits<Scalar>::coordinateType) PtsPerLine+1)/( (
typename Teuchos::ScalarTraits<Scalar>::coordinateType) (NCpts+1)));
317 RestStride = ((
typename Teuchos::ScalarTraits<Scalar>::coordinateType) (PtsPerLine-FirstStride+1))/((
typename Teuchos::ScalarTraits<Scalar>::coordinateType) NCpts);
319 NCLayers = (LO) floor((((
typename Teuchos::ScalarTraits<Scalar>::coordinateType) (PtsPerLine-FirstStride+1))/RestStride)+.00001);
323 di = (
typename Teuchos::ScalarTraits<Scalar>::coordinateType) FirstStride;
324 for (i = 1; i <= NCpts; i++) {
325 (*LayerCpts)[i] = (LO) floor(di);
332#define MaxHorNeighborNodes 75
334 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
336 MakeSemiCoarsenP(LO
const Ntotal, LO
const nz, LO
const CoarsenRate, LO
const LayerId[],
337 LO
const VertLineId[], LO
const DofsPerNode, RCP<Matrix> & Amat, RCP<Matrix>& P,
338 RCP<const Map>& coarseMap,
const RCP<MultiVector> fineNullspace,
339 RCP<MultiVector>& coarseNullspace, RCP<Matrix>& R,
bool buildRestriction,
bool doLinear)
const {
391 int NLayers, NVertLines, MaxNnz, NCLayers, MyLine, MyLayer;
392 int *InvLineLayer=NULL, *CptLayers=NULL, StartLayer, NStencilNodes;
393 int BlkRow=-1, dof_j, node_k, *Sub2FullMap=NULL, RowLeng;
394 int i, j, iii, col, count, index, loc, PtRow, PtCol;
395 SC *BandSol=NULL, *BandMat=NULL, TheSum, *RestrictBandMat=NULL, *RestrictBandSol=NULL;
396 int *IPIV=NULL, KL, KU, KLU, N, NRHS, LDAB,INFO;
403 int MaxStencilSize, MaxNnzPerRow;
405 int CurRow, LastGuy = -1, NewPtr, RLastGuy = -1;
408 LO *Layerdofs = NULL, *Col2Dof = NULL;
410 Teuchos::LAPACK<LO,SC> lapack;
421 Teuchos::ArrayRCP<LO> TLayDiff = Teuchos::arcp<LO>(1+MaxNnzPerRow); LayDiff = TLayDiff.getRawPtr();
423 Nghost = Amat->getColMap()->getLocalNumElements() - Amat->getDomainMap()->getLocalNumElements();
424 if (Nghost < 0) Nghost = 0;
425 Teuchos::ArrayRCP<LO> TLayerdofs= Teuchos::arcp<LO>(Ntotal*DofsPerNode+Nghost+1); Layerdofs = TLayerdofs.getRawPtr();
426 Teuchos::ArrayRCP<LO> TCol2Dof= Teuchos::arcp<LO>(Ntotal*DofsPerNode+Nghost+1); Col2Dof= TCol2Dof.getRawPtr();
428 RCP<Xpetra::Vector<LO,LO,GO,NO> > localdtemp = Xpetra::VectorFactory<LO,LO,GO,NO>::Build(Amat->getDomainMap());
429 RCP<Xpetra::Vector<LO,LO,GO,NO> > dtemp = Xpetra::VectorFactory<LO,LO,GO,NO>::Build(Amat->getColMap());
430 ArrayRCP<LO> valptr= localdtemp->getDataNonConst(0);
432 for (i = 0; i < Ntotal*DofsPerNode; i++)
433 valptr[i]= LayerId[i/DofsPerNode];
434 valptr=ArrayRCP<LO>();
436 RCP< const Import> importer;
437 importer = Amat->getCrsGraph()->getImporter();
438 if (importer == Teuchos::null) {
439 importer = ImportFactory::Build(Amat->getDomainMap(), Amat->getColMap());
441 dtemp->doImport(*localdtemp, *(importer), Xpetra::INSERT);
443 valptr= dtemp->getDataNonConst(0);
444 for (i = 0; i < Ntotal*DofsPerNode+Nghost; i++) Layerdofs[i]= valptr[i];
445 valptr= localdtemp->getDataNonConst(0);
446 for (i = 0; i < Ntotal*DofsPerNode; i++) valptr[i]= i%DofsPerNode;
447 valptr=ArrayRCP<LO>();
448 dtemp->doImport(*localdtemp, *(importer), Xpetra::INSERT);
450 valptr= dtemp->getDataNonConst(0);
451 for (i = 0; i < Ntotal*DofsPerNode+Nghost; i++) Col2Dof[i]= valptr[i];
452 valptr=ArrayRCP<LO>();
455 NLayers = LayerId[0];
456 NVertLines= VertLineId[0];
458 else { NLayers = -1; NVertLines = -1; }
460 for (i = 1; i < Ntotal; i++) {
461 if ( VertLineId[i] > NVertLines ) NVertLines = VertLineId[i];
462 if ( LayerId[i] > NLayers ) NLayers = LayerId[i];
472 Teuchos::ArrayRCP<LO> TInvLineLayer= Teuchos::arcp<LO>(1+NVertLines*NLayers); InvLineLayer = TInvLineLayer.getRawPtr();
473 for (i=0; i < Ntotal; i++) {
474 InvLineLayer[ VertLineId[i]+1+LayerId[i]*NVertLines ] = i;
481 Teuchos::ArrayRCP<LO> TCptLayers = Teuchos::arcp<LO>(nz+1); CptLayers = TCptLayers.getRawPtr();
482 NCLayers = FindCpts(nz,CoarsenRate,0, &CptLayers);
491 if (NCLayers < 2) MaxStencilSize = nz;
492 else MaxStencilSize = CptLayers[2];
494 for (i = 3; i <= NCLayers; i++) {
495 if (MaxStencilSize < CptLayers[i]- CptLayers[i-2])
496 MaxStencilSize = CptLayers[i]- CptLayers[i-2];
499 if (MaxStencilSize < nz - CptLayers[NCLayers-1]+1)
500 MaxStencilSize = nz - CptLayers[NCLayers-1]+1;
510 Teuchos::ArrayRCP<LO> TSub2FullMap= Teuchos::arcp<LO>((MaxStencilSize+1)*DofsPerNode); Sub2FullMap= TSub2FullMap.getRawPtr();
511 Teuchos::ArrayRCP<SC> TBandSol= Teuchos::arcp<SC>((MaxStencilSize+1)*DofsPerNode*DofsPerNode); BandSol = TBandSol.getRawPtr();
512 Teuchos::ArrayRCP<SC> TResBandSol;
513 if (buildRestriction) {
514 TResBandSol= Teuchos::arcp<SC>((MaxStencilSize+1)*DofsPerNode*DofsPerNode); RestrictBandSol = TResBandSol.getRawPtr();
519 KL = 2*DofsPerNode-1;
520 KU = 2*DofsPerNode-1;
524 Teuchos::ArrayRCP<SC> TBandMat= Teuchos::arcp<SC>(LDAB*MaxStencilSize*DofsPerNode+1); BandMat = TBandMat.getRawPtr();
525 Teuchos::ArrayRCP<LO> TIPIV= Teuchos::arcp<LO>((MaxStencilSize+1)*DofsPerNode); IPIV = TIPIV.getRawPtr();
527 Teuchos::ArrayRCP<SC> TResBandMat;
528 if (buildRestriction) {
529 TResBandMat= Teuchos::arcp<SC>(LDAB*MaxStencilSize*DofsPerNode+1); RestrictBandMat = TResBandMat.getRawPtr();
539 Ndofs = DofsPerNode*Ntotal;
540 MaxNnz = 2*DofsPerNode*Ndofs;
542 RCP<const Map> rowMap = Amat->getRowMap();
543 Xpetra::global_size_t GNdofs= rowMap->getGlobalNumElements();
545 std::vector<size_t> stridingInfo_;
546 stridingInfo_.push_back(DofsPerNode);
548 Xpetra::global_size_t itemp = GNdofs/nz;
549 coarseMap = StridedMapFactory::Build(rowMap->lib(),
551 NCLayers*NVertLines*DofsPerNode,
562 P = rcp(
new CrsMatrixWrap(rowMap, coarseMap , 0));
563 coarseNullspace = MultiVectorFactory::Build(coarseMap, fineNullspace->getNumVectors());
566 Teuchos::ArrayRCP<SC> TPvals= Teuchos::arcp<SC>(1+MaxNnz); Pvals= TPvals.getRawPtr();
567 Teuchos::ArrayRCP<size_t> TPptr = Teuchos::arcp<size_t>(DofsPerNode*(2+Ntotal)); Pptr = TPptr.getRawPtr();
568 Teuchos::ArrayRCP<LO> TPcols= Teuchos::arcp<LO>(1+MaxNnz); Pcols= TPcols.getRawPtr();
571 Pptr[0] = 0; Pptr[1] = 0;
574 Teuchos::ArrayRCP<SC> TRvals;
575 Teuchos::ArrayRCP<size_t> TRptr;
576 Teuchos::ArrayRCP<LO> TRcols;
577 LO RmaxNnz=-1, RmaxPerRow=-1;
578 if (buildRestriction) {
579 RmaxPerRow = ((LO) ceil(((
double) (MaxNnz+7))/((double) (coarseMap->getLocalNumElements()))));
580 RmaxNnz = RmaxPerRow*coarseMap->getLocalNumElements();
581 TRvals= Teuchos::arcp<SC>(1+RmaxNnz); Rvals= TRvals.getRawPtr();
582 TRptr = Teuchos::arcp<size_t>((2+coarseMap->getLocalNumElements())); Rptr = TRptr.getRawPtr();
583 TRcols= Teuchos::arcp<LO>(1+RmaxNnz); Rcols= TRcols.getRawPtr();
585 Rptr[0] = 0; Rptr[1] = 0;
593 for (i = 1; i <= MaxNnz; i++) Pcols[i] = -1;
595 for (i = 1; i <= DofsPerNode*Ntotal+1; i++) {
597 count += (2*DofsPerNode);
599 if (buildRestriction) {
600 for (i = 1; i <= RmaxNnz; i++) Rcols[i] = -1;
602 for (i = 1; i <= (int) (coarseMap->getLocalNumElements()+1); i++) {
617 for (MyLine=1; MyLine <= NVertLines; MyLine += 1) {
618 for (iii=1; iii <= NCLayers; iii+= 1) {
620 MyLayer = CptLayers[iii];
633 if (iii != 1 ) StartLayer = CptLayers[iii-1]+1;
636 if (iii != NCLayers) NStencilNodes= CptLayers[iii+1]-StartLayer;
637 else NStencilNodes= NLayers - StartLayer + 1;
640 N = NStencilNodes*DofsPerNode;
647 for (i = 0; i < NStencilNodes*DofsPerNode*DofsPerNode; i++) BandSol[ i] = 0.0;
648 for (i = 0; i < LDAB*N; i++) BandMat[ i] = 0.0;
649 if (buildRestriction) {
650 for (i = 0; i < NStencilNodes*DofsPerNode*DofsPerNode; i++) RestrictBandSol[ i] = 0.0;
651 for (i = 0; i < LDAB*N; i++) RestrictBandMat[ i] = 0.0;
660 for (node_k=1; node_k <= NStencilNodes ; node_k++) {
666 BlkRow = InvLineLayer[MyLine+(StartLayer+node_k-2)*NVertLines]+1;
667 Sub2FullMap[node_k] = BlkRow;
680 for (
int dof_i=0; dof_i < DofsPerNode; dof_i++) {
682 j = (BlkRow-1)*DofsPerNode+dof_i;
683 ArrayView<const LO> AAcols;
684 ArrayView<const SC> AAvals;
685 Amat->getLocalRowView(j, AAcols, AAvals);
686 const int *Acols = AAcols.getRawPtr();
687 const SC *Avals = AAvals.getRawPtr();
688 RowLeng = AAvals.size();
690 TEUCHOS_TEST_FOR_EXCEPTION(RowLeng >= MaxNnzPerRow,
Exceptions::RuntimeError,
"MakeSemiCoarsenP: recompile with larger Max(HorNeighborNodes)\n");
692 for (i = 0; i < RowLeng; i++) {
693 LayDiff[i] = Layerdofs[Acols[i]]-StartLayer-node_k+2;
701 PtRow = (node_k-1)*DofsPerNode+dof_i+1;
702 for (dof_j=0; dof_j < DofsPerNode; dof_j++) {
703 PtCol = (node_k-1)*DofsPerNode+dof_j + 1;
707 for (i = 0; i < RowLeng; i++) {
708 if ((LayDiff[i] == 0) && (Col2Dof[Acols[i]] == dof_j))
711 index = LDAB*(PtCol-1)+KLU+PtRow-PtCol;
712 BandMat[index] = TheSum;
713 if (buildRestriction) RestrictBandMat[index] = TheSum;
714 if (node_k != NStencilNodes) {
718 for (i = 0; i < RowLeng; i++) {
719 if ((LayDiff[i] == 1) &&(Col2Dof[Acols[i]]== dof_j))
722 j = PtCol+DofsPerNode;
723 index=LDAB*(j-1)+KLU+PtRow-j;
724 BandMat[index] = TheSum;
725 if (buildRestriction) RestrictBandMat[index] = TheSum;
731 for (i = 0; i < RowLeng; i++) {
732 if ((LayDiff[i]== -1) &&(Col2Dof[Acols[i]]== dof_j))
735 j = PtCol-DofsPerNode;
736 index=LDAB*(j-1)+KLU+PtRow-j;
737 BandMat[index] = TheSum;
738 if (buildRestriction) RestrictBandMat[index] = TheSum;
743 node_k = MyLayer - StartLayer+1;
744 for (
int dof_i = 0; dof_i < DofsPerNode; dof_i++) {
747 PtRow = (node_k-1)*DofsPerNode+dof_i+1;
748 BandSol[(dof_i)*DofsPerNode*NStencilNodes+PtRow-1] = 1.;
749 if (buildRestriction) RestrictBandSol[(dof_i)*DofsPerNode*NStencilNodes+PtRow-1] = 1.;
751 for (dof_j = 0; dof_j < DofsPerNode; dof_j++) {
752 PtCol = (node_k-1)*DofsPerNode+dof_j+1;
753 index = LDAB*(PtCol-1)+KLU+PtRow-PtCol;
754 if (PtCol == PtRow) BandMat[index] = 1.0;
755 else BandMat[index] = 0.0;
756 if (buildRestriction) {
757 index = LDAB*(PtRow-1)+KLU+PtCol-PtRow;
758 if (PtCol == PtRow) RestrictBandMat[index] = 1.0;
759 else RestrictBandMat[index] = 0.0;
761 if (node_k != NStencilNodes) {
762 PtCol = (node_k )*DofsPerNode+dof_j+1;
763 index = LDAB*(PtCol-1)+KLU+PtRow-PtCol;
764 BandMat[index] = 0.0;
765 if (buildRestriction) {
766 index = LDAB*(PtRow-1)+KLU+PtCol-PtRow;
767 RestrictBandMat[index] = 0.0;
771 PtCol = (node_k-2)*DofsPerNode+dof_j+1;
772 index = LDAB*(PtCol-1)+KLU+PtRow-PtCol;
773 BandMat[index] = 0.0;
774 if (buildRestriction) {
775 index = LDAB*(PtRow-1)+KLU+PtCol-PtRow;
776 RestrictBandMat[index] = 0.0;
784 lapack.GBTRF( N, N, KL, KU, BandMat, LDAB, IPIV, &INFO);
788 lapack.GBTRS(notrans[0], N, KL, KU, NRHS, BandMat, LDAB, IPIV,
793 for (dof_j=0; dof_j < DofsPerNode; dof_j++) {
794 for (
int dof_i=0; dof_i < DofsPerNode; dof_i++) {
795 for (i =1; i <= NStencilNodes ; i++) {
796 index = (Sub2FullMap[i]-1)*DofsPerNode+dof_i+1;
798 Pcols[loc] = (col-1)*DofsPerNode+dof_j+1;
799 Pvals[loc] = BandSol[dof_j*DofsPerNode*NStencilNodes +
800 (i-1)*DofsPerNode + dof_i];
801 Pptr[index]= Pptr[index] + 1;
805 if (buildRestriction) {
806 lapack.GBTRF( N, N, KL, KU, RestrictBandMat, LDAB, IPIV, &INFO);
808 lapack.GBTRS(trans[0], N, KL, KU, NRHS, RestrictBandMat, LDAB, IPIV, RestrictBandSol, N, &INFO );
810 for (dof_j=0; dof_j < DofsPerNode; dof_j++) {
811 for (
int dof_i=0; dof_i < DofsPerNode; dof_i++) {
812 for (i =1; i <= NStencilNodes ; i++) {
813 index = (col-1)*DofsPerNode+dof_j+1;
815 Rcols[loc] = (Sub2FullMap[i]-1)*DofsPerNode+dof_i+1;
816 Rvals[loc] = RestrictBandSol[dof_j*DofsPerNode*NStencilNodes +
817 (i-1)*DofsPerNode + dof_i];
818 Rptr[index]= Rptr[index] + 1;
825 int denom1 = MyLayer-StartLayer+1;
826 int denom2 = StartLayer+NStencilNodes-MyLayer;
827 for (
int dof_i=0; dof_i < DofsPerNode; dof_i++) {
828 for (i =1; i <= NStencilNodes ; i++) {
829 index = (InvLineLayer[MyLine+(StartLayer+i-2)*NVertLines])*DofsPerNode+dof_i+1;
831 Pcols[loc] = (col-1)*DofsPerNode+dof_i+1;
832 if (i > denom1) Pvals[loc] = 1.0 + ((double)( denom1 - i))/((
double) denom2);
833 else Pvals[loc] = ((double)( i))/((
double) denom1);
834 Pptr[index]= Pptr[index] + 1;
842 BlkRow = InvLineLayer[MyLine+(MyLayer-1)*NVertLines]+1;
843 for (
int k = 0; k < static_cast<int>(fineNullspace->getNumVectors()); k++) {
844 Teuchos::ArrayRCP<SC> OneCNull = coarseNullspace->getDataNonConst(k);
845 Teuchos::ArrayRCP<SC> OneFNull = fineNullspace->getDataNonConst(k);
846 for (
int dof_i = 0; dof_i < DofsPerNode; dof_i++) {
847 OneCNull[( col-1)*DofsPerNode+dof_i] = OneFNull[ (BlkRow-1)*DofsPerNode+dof_i];
863 for (
size_t ii=1; ii <= Pptr[Ntotal*DofsPerNode]-1; ii++) {
864 if (ii == Pptr[CurRow]) {
865 Pptr[CurRow] = LastGuy;
867 while (ii > Pptr[CurRow]) {
868 Pptr[CurRow] = LastGuy;
872 if (Pcols[ii] != -1) {
873 Pcols[NewPtr-1] = Pcols[ii]-1;
874 Pvals[NewPtr-1] = Pvals[ii];
879 for (i = CurRow; i <= Ntotal*DofsPerNode; i++) Pptr[CurRow] = LastGuy;
884 for (i=-Ntotal*DofsPerNode+1; i>= 2 ; i--) {
885 Pptr[i-1] = Pptr[i-2];
890 if (buildRestriction) {
893 for (
size_t ii=1; ii <= Rptr[coarseMap->getLocalNumElements()]-1; ii++) {
894 if (ii == Rptr[CurRow]) {
895 Rptr[CurRow] = RLastGuy;
897 while (ii > Rptr[CurRow]) {
898 Rptr[CurRow] = RLastGuy;
902 if (Rcols[ii] != -1) {
903 Rcols[NewPtr-1] = Rcols[ii]-1;
904 Rvals[NewPtr-1] = Rvals[ii];
909 for (i = CurRow; i <= ((int) coarseMap->getLocalNumElements()); i++) Rptr[CurRow] = RLastGuy;
910 for (i=-coarseMap->getLocalNumElements()+1; i>= 2 ; i--) {
911 Rptr[i-1] = Rptr[i-2];
915 ArrayRCP<size_t> rcpRowPtr;
916 ArrayRCP<LO> rcpColumns;
917 ArrayRCP<SC> rcpValues;
919 RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
920 LO nnz = Pptr[Ndofs];
921 PCrs->allocateAllValues(nnz, rcpRowPtr, rcpColumns, rcpValues);
923 ArrayView<size_t> rowPtr = rcpRowPtr();
924 ArrayView<LO> columns = rcpColumns();
925 ArrayView<SC> values = rcpValues();
929 for (LO ii = 0; ii <= Ndofs; ii++) rowPtr[ii] = Pptr[ii];
930 for (LO ii = 0; ii < nnz; ii++) columns[ii] = Pcols[ii];
931 for (LO ii = 0; ii < nnz; ii++) values[ii] = Pvals[ii];
932 PCrs->setAllValues(rcpRowPtr, rcpColumns, rcpValues);
933 PCrs->expertStaticFillComplete(coarseMap, Amat->getDomainMap());
934 ArrayRCP<size_t> RrcpRowPtr;
935 ArrayRCP<LO> RrcpColumns;
936 ArrayRCP<SC> RrcpValues;
938 if (buildRestriction) {
939 R = rcp(
new CrsMatrixWrap(coarseMap, rowMap, 0));
940 RCrs = rcp_dynamic_cast<CrsMatrixWrap>(R)->getCrsMatrix();
941 nnz = Rptr[coarseMap->getLocalNumElements()];
942 RCrs->allocateAllValues(nnz, RrcpRowPtr, RrcpColumns, RrcpValues);
944 ArrayView<size_t> RrowPtr = RrcpRowPtr();
945 ArrayView<LO> Rcolumns = RrcpColumns();
946 ArrayView<SC> Rvalues = RrcpValues();
950 for (LO ii = 0; ii <= ((LO)coarseMap->getLocalNumElements()); ii++) RrowPtr[ii] = Rptr[ii];
951 for (LO ii = 0; ii < nnz; ii++) Rcolumns[ii] = Rcols[ii];
952 for (LO ii = 0; ii < nnz; ii++) Rvalues[ii] = Rvals[ii];
953 RCrs->setAllValues(RrcpRowPtr, RrcpColumns, RrcpValues);
954 RCrs->expertStaticFillComplete(Amat->getRangeMap(), coarseMap);
957 return NCLayers*NVertLines*DofsPerNode;
959 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
967 ArrayView<const LocalOrdinal> inds;
968 ArrayView<const Scalar> vals1;
969 ArrayView< Scalar> vals;
970 Scalar ZERO = Teuchos::ScalarTraits<Scalar>::zero();
971 Scalar ONE = Teuchos::ScalarTraits<Scalar>::one();
973 for (
size_t i = 0; i < P->getRowMap()->getLocalNumElements(); i++) {
974 P->getLocalRowView(i, inds, vals1);
976 size_t nnz = inds.size();
977 if (nnz != 0) vals = ArrayView<Scalar>(
const_cast<Scalar*
>(vals1.getRawPtr()), nnz);
979 LO largestIndex = -1;
980 Scalar largestValue = ZERO;
983 LO rowDof = i%BlkSize;
984 for (
size_t j =0; j < nnz; j++) {
985 if (Teuchos::ScalarTraits<SC>::magnitude(vals[ j ]) >= Teuchos::ScalarTraits<SC>::magnitude(largestValue)) {
986 if ( inds[j]%BlkSize == rowDof ) {
987 largestValue = vals[j];
988 largestIndex = (int) j;
993 if (largestIndex != -1) vals[largestIndex] = ONE;
995 TEUCHOS_TEST_FOR_EXCEPTION(nnz > 0,
Exceptions::RuntimeError,
"no nonzero column associated with a proper dof within node.");
997 if (Teuchos::ScalarTraits<SC>::magnitude(largestValue) == Teuchos::ScalarTraits<SC>::magnitude(ZERO)) vals[largestIndex] = ZERO;
1003#define MUELU_SEMICOARSENPFACTORY_SHORT
#define SET_VALID_ENTRY(name)
#define MaxHorNeighborNodes
MueLu::DefaultLocalOrdinal LocalOrdinal
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.
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.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
const RCP< const FactoryManagerBase > GetFactoryManager()
returns the current factory manager
int GetLevelID() const
Return level number.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
static const NoFactory * get()
LO FindCpts(LO const PtsPerLine, LO const CoarsenRate, LO const Thin, LO **LayerCpts) const
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
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, RCP< Matrix > &R, bool buildRestriction, bool doLinear) const
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
void RevertToPieceWiseConstant(RCP< Matrix > &P, LO BlkSize) const
Namespace for MueLu classes and methods.