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);
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;