46 #ifndef MUELU_SEMICOARSENPFACTORY_DEF_HPP 47 #define MUELU_SEMICOARSENPFACTORY_DEF_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>
71 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 73 #undef SET_VALID_ENTRY 78 validParamList->
set<
RCP<const FactoryBase> >(
"LineDetection_VertLineIds", Teuchos::null,
"Generating factory for LineDetection information");
79 validParamList->
set<
RCP<const FactoryBase> >(
"LineDetection_Layers", Teuchos::null,
"Generating factory for LineDetection information");
80 validParamList->
set<
RCP<const FactoryBase> >(
"CoarseNumZLayers", Teuchos::null,
"Generating factory for LineDetection information");
82 return validParamList;
85 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
87 Input(fineLevel,
"A");
88 Input(fineLevel,
"Nullspace");
90 Input(fineLevel,
"LineDetection_VertLineIds");
91 Input(fineLevel,
"LineDetection_Layers");
92 Input(fineLevel,
"CoarseNumZLayers");
100 bTransferCoordinates_ =
true;
102 }
else if (bTransferCoordinates_ ==
true){
107 if (myCoordsFact == Teuchos::null) { myCoordsFact = fineLevel.
GetFactoryManager()->GetFactory(
"Coordinates"); }
110 bTransferCoordinates_ =
true;
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");
130 LO CoarsenRate = as<LO>(pL.get<
int>(
"semicoarsen: coarsen rate"));
133 LO BlkSize = A->GetFixedBlockSize();
135 LO Ndofs = rowMap->getNodeNumElements();
136 LO Nnodes = Ndofs/BlkSize;
139 LO FineNumZLayers = Get< LO >(fineLevel,
"CoarseNumZLayers");
140 LO CoarseNumZLayers = FineNumZLayers;
141 Teuchos::ArrayRCP<LO> TVertLineId = Get< Teuchos::ArrayRCP<LO> > (fineLevel,
"LineDetection_VertLineIds");
143 LO* VertLineId = TVertLineId.
getRawPtr();
144 LO* LayerId = TLayerId.getRawPtr();
149 GO Ncoarse = MakeSemiCoarsenP(Nnodes,CoarseNumZLayers,CoarsenRate,LayerId,VertLineId,
150 BlkSize, A, P, theCoarseMap);
153 if (A->IsView(
"stridedMaps") ==
true)
154 P->CreateView(
"stridedMaps", A->getRowMap(
"stridedMaps"), theCoarseMap);
156 P->CreateView(
"stridedMaps", P->getRangeMap(), theCoarseMap);
164 Set(coarseLevel,
"P", P);
168 RCP<MultiVector> coarseNullspace = MultiVectorFactory::Build(P->getDomainMap(), fineNullspace->getNumVectors());
170 Set(coarseLevel,
"Nullspace", coarseNullspace);
173 if(bTransferCoordinates_) {
175 typedef Xpetra::MultiVector<double,LO,GO,NO> xdMV;
177 if (fineLevel.GetLevelID() == 0 &&
182 if (myCoordsFact == Teuchos::null) { myCoordsFact = fineLevel.GetFactoryManager()->GetFactory(
"Coordinates"); }
183 if (fineLevel.IsAvailable(
"Coordinates", myCoordsFact.
get())) {
184 fineCoords = fineLevel.Get<
RCP<xdMV> >(
"Coordinates", myCoordsFact.
get());
199 if(*it > zval_max) zval_max = *it;
200 if(*it < zval_min) zval_min = *it;
203 LO myCoarseZLayers = Teuchos::as<LO>(CoarseNumZLayers * Ncoarse/Ndofs);
206 if(myCoarseZLayers == 1) {
207 myZLayerCoords[0] = zval_min;
209 double dz = (zval_max-zval_min)/(myCoarseZLayers-1);
210 for(LO k = 0; k<myCoarseZLayers; ++k) {
211 myZLayerCoords[k] = k*dz;
219 LO numVertLines = Nnodes / FineNumZLayers;
220 LO numLocalCoarseNodes = numVertLines * myCoarseZLayers;
229 MapFactory::Build (fineCoords->getMap()->lib(),
231 Teuchos::as<size_t>(numLocalCoarseNodes),
232 fineCoords->getMap()->getIndexBase(),
233 fineCoords->getMap()->getComm());
234 RCP<xdMV> coarseCoords = Xpetra::MultiVectorFactory<double,LO,GO,NO>::Build(coarseCoordMap, fineCoords->getNumVectors());
235 coarseCoords->putScalar(-1.0);
241 LO cntCoarseNodes = 0;
242 for( LO vt = 0; vt < TVertLineId.
size(); ++vt) {
244 LO curVertLineId = TVertLineId[vt];
246 if(cx[curVertLineId * myCoarseZLayers] == -1.0 &&
247 cy[curVertLineId * myCoarseZLayers] == -1.0) {
249 for (LO n=0; n<myCoarseZLayers; ++n) {
250 cx[curVertLineId * myCoarseZLayers + n] = x[vt];
251 cy[curVertLineId * myCoarseZLayers + n] = y[vt];
252 cz[curVertLineId * myCoarseZLayers + n] = myZLayerCoords[n];
254 cntCoarseNodes += myCoarseZLayers;
258 if(cntCoarseNodes == numLocalCoarseNodes)
break;
262 Set(coarseLevel,
"Coordinates", coarseCoords);
267 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
291 double temp, RestStride, di;
296 temp = ((double) (PtsPerLine+1))/((
double) (CoarsenRate)) - 1.0;
297 if (Thin == 1) NCpts = (LO) ceil(temp);
298 else NCpts = (LO) floor(temp);
302 if (NCpts < 1) NCpts = 1;
304 FirstStride= (LO) ceil( ((
double) PtsPerLine+1)/( (double) (NCpts+1)));
305 RestStride = ((double) (PtsPerLine-FirstStride+1))/((
double) NCpts);
307 NCLayers = (LO) floor((((
double) (PtsPerLine-FirstStride+1))/RestStride)+.00001);
311 di = (double) FirstStride;
312 for (i = 1; i <= NCpts; i++) {
313 (*LayerCpts)[i] = (LO) floor(di);
320 #define MaxHorNeighborNodes 75 322 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
376 int NLayers, NVertLines, MaxNnz, NCLayers, MyLine, MyLayer;
377 int *InvLineLayer=NULL, *CptLayers=NULL, StartLayer, NStencilNodes;
378 int BlkRow, dof_i, dof_j, node_k, *Sub2FullMap=NULL, RowLeng;
379 int i, j, iii, col, count, index, loc, PtRow, PtCol;
380 SC *BandSol=NULL, *BandMat=NULL, TheSum;
381 int *IPIV=NULL, KL, KU, KLU, N, NRHS, LDAB,INFO;
385 int MaxStencilSize, MaxNnzPerRow;
387 int CurRow, LastGuy = -1, NewPtr;
390 LO *Layerdofs = NULL, *Col2Dof = NULL;
402 Nghost = Amat->getColMap()->getNodeNumElements() - Amat->getDomainMap()->getNodeNumElements();
403 if (Nghost < 0) Nghost = 0;
411 for (i = 0; i < Ntotal*DofsPerNode; i++)
412 valptr[i]= LayerId[i/DofsPerNode];
415 importer = Amat->getCrsGraph()->getImporter();
416 if (importer == Teuchos::null) {
417 importer = ImportFactory::Build(Amat->getDomainMap(), Amat->getColMap());
419 dtemp->doImport(*localdtemp, *(importer), Xpetra::INSERT);
421 valptr= dtemp->getDataNonConst(0);
422 for (i = 0; i < Ntotal*DofsPerNode+Nghost; i++) Layerdofs[i]= valptr[i];
423 valptr= localdtemp->getDataNonConst(0);
424 for (i = 0; i < Ntotal*DofsPerNode; i++) valptr[i]= i%DofsPerNode;
425 dtemp->doImport(*localdtemp, *(importer), Xpetra::INSERT);
426 valptr= dtemp->getDataNonConst(0);
427 for (i = 0; i < Ntotal*DofsPerNode+Nghost; i++) Col2Dof[i]= valptr[i];
430 NLayers = LayerId[0];
431 NVertLines= VertLineId[0];
433 else { NLayers = -1; NVertLines = -1; }
435 for (i = 1; i < Ntotal; i++) {
436 if ( VertLineId[i] > NVertLines ) NVertLines = VertLineId[i];
437 if ( LayerId[i] > NLayers ) NLayers = LayerId[i];
448 for (i=0; i < Ntotal; i++) {
449 InvLineLayer[ VertLineId[i]+1+LayerId[i]*NVertLines ] = i;
457 NCLayers = FindCpts(nz,CoarsenRate,0, &CptLayers);
466 if (NCLayers < 2) MaxStencilSize = nz;
467 else MaxStencilSize = CptLayers[2];
469 for (i = 3; i <= NCLayers; i++) {
470 if (MaxStencilSize < CptLayers[i]- CptLayers[i-2])
471 MaxStencilSize = CptLayers[i]- CptLayers[i-2];
474 if (MaxStencilSize < nz - CptLayers[NCLayers-1]+1)
475 MaxStencilSize = nz - CptLayers[NCLayers-1]+1;
490 KL = 2*DofsPerNode-1;
491 KU = 2*DofsPerNode-1;
505 Ndofs = DofsPerNode*Ntotal;
506 MaxNnz = 2*DofsPerNode*Ndofs;
509 int GNdofs= rowMap->getGlobalNumElements();
511 std::vector<size_t> stridingInfo_;
512 stridingInfo_.push_back(DofsPerNode);
514 coarseMap = StridedMapFactory::Build(rowMap->lib(),
515 (NCLayers*GNdofs)/nz,
516 NCLayers*NVertLines*DofsPerNode,
528 P =
rcp(
new CrsMatrixWrap(rowMap, coarseMap , 0, Xpetra::StaticProfile));
535 Pptr[0] = 0; Pptr[1] = 0;
545 for (i = 1; i <= MaxNnz; i++) Pcols[i] = -1;
547 for (i = 1; i <= DofsPerNode*Ntotal+1; i++) {
549 count += (2*DofsPerNode);
561 for (MyLine=1; MyLine <= NVertLines; MyLine += 1) {
562 for (iii=1; iii <= NCLayers; iii+= 1) {
564 MyLayer = CptLayers[iii];
577 if (iii != 1 ) StartLayer = CptLayers[iii-1]+1;
580 if (iii != NCLayers) NStencilNodes= CptLayers[iii+1]-StartLayer;
581 else NStencilNodes= NLayers - StartLayer + 1;
584 N = NStencilNodes*DofsPerNode;
591 for (i = 0; i < NStencilNodes*DofsPerNode*DofsPerNode; i++)
593 for (i = 0; i < LDAB*N; i++) BandMat[ i] = 0.0;
600 for (node_k=1; node_k <= NStencilNodes ; node_k++) {
606 BlkRow = InvLineLayer[MyLine+(StartLayer+node_k-2)*NVertLines]+1;
607 Sub2FullMap[node_k] = BlkRow;
620 if (StartLayer+node_k-1 != MyLayer) {
621 for (dof_i=0; dof_i < DofsPerNode; dof_i++) {
623 j = (BlkRow-1)*DofsPerNode+dof_i;
626 Amat->getLocalRowView(j, AAcols, AAvals);
629 RowLeng = AAvals.
size();
633 for (i = 0; i < RowLeng; i++) {
634 LayDiff[i] = Layerdofs[Acols[i]]-StartLayer-node_k+2;
642 PtRow = (node_k-1)*DofsPerNode+dof_i+1;
643 for (dof_j=0; dof_j < DofsPerNode; dof_j++) {
644 PtCol = (node_k-1)*DofsPerNode+dof_j + 1;
648 for (i = 0; i < RowLeng; i++) {
649 if ((LayDiff[i] == 0) && (Col2Dof[Acols[i]] == dof_j))
652 index = LDAB*(PtCol-1)+KLU+PtRow-PtCol;
653 BandMat[index] = TheSum;
654 if (node_k != NStencilNodes) {
658 for (i = 0; i < RowLeng; i++) {
659 if ((LayDiff[i] == 1) &&(Col2Dof[Acols[i]]== dof_j))
662 j = PtCol+DofsPerNode;
663 index=LDAB*(j-1)+KLU+PtRow-j;
664 BandMat[index] = TheSum;
670 for (i = 0; i < RowLeng; i++) {
671 if ((LayDiff[i]== -1) &&(Col2Dof[Acols[i]]== dof_j))
674 j = PtCol-DofsPerNode;
675 index=LDAB*(j-1)+KLU+PtRow-j;
676 BandMat[index] = TheSum;
682 for (dof_i = 0; dof_i < DofsPerNode; dof_i++) {
685 PtRow = (node_k-1)*DofsPerNode+dof_i+1;
686 index = LDAB*(PtRow-1)+KLU;
687 BandMat[index] = 1.0;
688 BandSol[(dof_i)*DofsPerNode*NStencilNodes+PtRow-1] = 1.;
695 lapack.
GBTRF( N, N, KL, KU, BandMat, LDAB, IPIV, &INFO);
699 lapack.
GBTRS(notrans[0], N, KL, KU, NRHS, BandMat, LDAB, IPIV,
704 for (dof_j=0; dof_j < DofsPerNode; dof_j++) {
705 for (dof_i=0; dof_i < DofsPerNode; dof_i++) {
706 for (i =1; i <= NStencilNodes ; i++) {
707 index = (Sub2FullMap[i]-1)*DofsPerNode+dof_i+1;
709 Pcols[loc] = (col-1)*DofsPerNode+dof_j+1;
710 Pvals[loc] = BandSol[dof_j*DofsPerNode*NStencilNodes +
711 (i-1)*DofsPerNode + dof_i];
712 Pptr[index]= Pptr[index] + 1;
728 for (
size_t ii=1; ii <= Pptr[Ntotal*DofsPerNode]-1; ii++) {
729 if (ii == Pptr[CurRow]) {
730 Pptr[CurRow] = LastGuy;
732 while (ii > Pptr[CurRow]) {
733 Pptr[CurRow] = LastGuy;
737 if (Pcols[ii] != -1) {
738 Pcols[NewPtr-1] = Pcols[ii]-1;
739 Pvals[NewPtr-1] = Pvals[ii];
744 for (i = CurRow; i <= Ntotal*DofsPerNode; i++) Pptr[CurRow] = LastGuy;
749 for (i=-Ntotal*DofsPerNode+1; i>= 2 ; i--) {
750 Pptr[i-1] = Pptr[i-2];
758 RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
759 LO nnz = Pptr[Ndofs];
760 PCrs->allocateAllValues(nnz, rcpRowPtr, rcpColumns, rcpValues);
768 for (LO ii = 0; ii <= Ndofs; ii++) rowPtr[ii] = Pptr[ii];
769 for (LO ii = 0; ii < nnz; ii++) columns[ii] = Pcols[ii];
770 for (LO ii = 0; ii < nnz; ii++) values[ii] = Pvals[ii];
771 PCrs->setAllValues(rcpRowPtr, rcpColumns, rcpValues);
772 PCrs->expertStaticFillComplete(coarseMap, Amat->getDomainMap());
775 return NCLayers*NVertLines*DofsPerNode;
779 #define MUELU_SEMICOARSENPFACTORY_SHORT 780 #endif // MUELU_SEMICOARSENPFACTORY_DEF_HPP LO FindCpts(LO const PtsPerLine, LO const CoarsenRate, LO const Thin, LO **LayerCpts) const
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
Namespace for MueLu classes and methods.
static magnitudeType sfmin()
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.
void GBTRS(const char TRANS, const OrdinalType n, const OrdinalType kl, const OrdinalType ku, const OrdinalType nrhs, const ScalarType *A, const OrdinalType lda, const OrdinalType *IPIV, ScalarType *B, const OrdinalType ldb, OrdinalType *info) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
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
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
void GBTRF(const OrdinalType m, const OrdinalType n, const OrdinalType kl, const OrdinalType ku, ScalarType *A, const OrdinalType lda, OrdinalType *IPIV, OrdinalType *info) const
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()
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.