46 #ifndef MUELU_SEMICOARSENPFACTORY_KOKKOS_DEF_HPP 47 #define MUELU_SEMICOARSENPFACTORY_KOKKOS_DEF_HPP 49 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 53 #include <Kokkos_Core.hpp> 55 #include <KokkosBatched_LU_Decl.hpp> 56 #include <KokkosBatched_LU_Serial_Impl.hpp> 57 #include <KokkosBatched_Trsm_Decl.hpp> 58 #include <KokkosBatched_Trsm_Serial_Impl.hpp> 59 #include <KokkosBatched_Util.hpp> 60 #include <KokkosSparse_CrsMatrix.hpp> 62 #include <Xpetra_CrsMatrixWrap.hpp> 63 #include <Xpetra_ImportFactory.hpp> 64 #include <Xpetra_Matrix.hpp> 65 #include <Xpetra_MultiVectorFactory.hpp> 66 #include <Xpetra_VectorFactory.hpp> 77 RCP<const ParameterList>
79 Kokkos::Compat::KokkosDeviceWrapperNode<
80 DeviceType>>::GetValidParameterList()
const {
81 RCP<ParameterList> validParamList = rcp(
new ParameterList());
83 std::string name =
"semicoarsen: coarsen rate";
85 validParamList->set<RCP<const FactoryBase>>(
86 "A", Teuchos::null,
"Generating factory of the matrix A");
87 validParamList->set<RCP<const FactoryBase>>(
88 "Nullspace", Teuchos::null,
"Generating factory of the nullspace");
89 validParamList->set<RCP<const FactoryBase>>(
90 "Coordinates", Teuchos::null,
"Generating factory for coordinates");
92 validParamList->set<RCP<const FactoryBase>>(
93 "LineDetection_VertLineIds", Teuchos::null,
94 "Generating factory for LineDetection vertical line ids");
95 validParamList->set<RCP<const FactoryBase>>(
96 "LineDetection_Layers", Teuchos::null,
97 "Generating factory for LineDetection layer ids");
98 validParamList->set<RCP<const FactoryBase>>(
99 "CoarseNumZLayers", Teuchos::null,
100 "Generating factory for number of coarse z-layers");
102 return validParamList;
107 void SemiCoarsenPFactory_kokkos<
109 Kokkos::Compat::KokkosDeviceWrapperNode<DeviceType>>::
110 DeclareInput(Level &fineLevel, Level & )
const {
111 Input(fineLevel,
"A");
112 Input(fineLevel,
"Nullspace");
114 Input(fineLevel,
"LineDetection_VertLineIds");
115 Input(fineLevel,
"LineDetection_Layers");
116 Input(fineLevel,
"CoarseNumZLayers");
121 if (fineLevel.GetLevelID() == 0) {
124 bTransferCoordinates_ =
true;
126 }
else if (bTransferCoordinates_ ==
true) {
130 RCP<const FactoryBase> myCoordsFact = GetFactory(
"Coordinates");
131 if (myCoordsFact == Teuchos::null) {
132 myCoordsFact = fineLevel.GetFactoryManager()->GetFactory(
"Coordinates");
134 if (fineLevel.IsAvailable(
"Coordinates", myCoordsFact.get())) {
135 fineLevel.DeclareInput(
"Coordinates", myCoordsFact.get(),
this);
143 Kokkos::Compat::KokkosDeviceWrapperNode<
144 DeviceType>>::Build(Level &fineLevel,
147 return BuildP(fineLevel, coarseLevel);
153 Kokkos::Compat::KokkosDeviceWrapperNode<
154 DeviceType>>::BuildP(Level &fineLevel,
157 FactoryMonitor m(*
this,
"Build", coarseLevel);
160 RCP<Matrix> A = Get<RCP<Matrix>>(fineLevel,
"A");
161 RCP<MultiVector> fineNullspace =
162 Get<RCP<MultiVector>>(fineLevel,
"Nullspace");
165 const ParameterList &pL = GetParameterList();
166 LO CoarsenRate = as<LO>(pL.get<
int>(
"semicoarsen: coarsen rate"));
167 TEUCHOS_TEST_FOR_EXCEPTION(
168 CoarsenRate < 2, Exceptions::RuntimeError,
169 "semicoarsen: coarsen rate must be greater than 1");
172 LO BlkSize = A->GetFixedBlockSize();
173 RCP<const Map> rowMap = A->getRowMap();
174 LO Ndofs = rowMap->getNodeNumElements();
175 LO Nnodes = Ndofs / BlkSize;
179 LO FineNumZLayers = Get<LO>(fineLevel,
"CoarseNumZLayers");
180 Teuchos::ArrayRCP<LO> TVertLineId =
181 Get<Teuchos::ArrayRCP<LO>>(fineLevel,
"LineDetection_VertLineIds");
182 Teuchos::ArrayRCP<LO> TLayerId =
183 Get<Teuchos::ArrayRCP<LO>>(fineLevel,
"LineDetection_Layers");
186 TEUCHOS_TEST_FOR_EXCEPTION(FineNumZLayers < 2, Exceptions::RuntimeError,
187 "Cannot coarsen further");
188 using coordT =
typename Teuchos::ScalarTraits<Scalar>::coordinateType;
189 LO CoarseNumZLayers =
190 (LO)floor(((coordT)(FineNumZLayers + 1)) / ((coordT)CoarsenRate) - 1.0);
191 if (CoarseNumZLayers < 1)
192 CoarseNumZLayers = 1;
196 RCP<MultiVector> coarseNullspace;
197 BuildSemiCoarsenP(coarseLevel, Ndofs, Nnodes, BlkSize, FineNumZLayers,
198 CoarseNumZLayers, TLayerId, TVertLineId, A, fineNullspace, P,
204 coarseLevel.Set(
"NumZLayers", Teuchos::as<LO>(CoarseNumZLayers),
208 Set(coarseLevel,
"P", P);
209 Set(coarseLevel,
"Nullspace", coarseNullspace);
212 if (bTransferCoordinates_) {
213 SubFactoryMonitor m2(*
this,
"TransferCoordinates", coarseLevel);
214 typedef Xpetra::MultiVector<
215 typename Teuchos::ScalarTraits<Scalar>::coordinateType, LO, GO, NO>
217 RCP<xdMV> fineCoords = Teuchos::null;
218 if (fineLevel.GetLevelID() == 0 &&
220 fineCoords = fineLevel.Get<RCP<xdMV>>(
"Coordinates",
NoFactory::get());
222 RCP<const FactoryBase> myCoordsFact = GetFactory(
"Coordinates");
223 if (myCoordsFact == Teuchos::null) {
224 myCoordsFact = fineLevel.GetFactoryManager()->GetFactory(
"Coordinates");
226 if (fineLevel.IsAvailable(
"Coordinates", myCoordsFact.get())) {
228 fineLevel.Get<RCP<xdMV>>(
"Coordinates", myCoordsFact.get());
232 TEUCHOS_TEST_FOR_EXCEPTION(fineCoords == Teuchos::null,
233 Exceptions::RuntimeError,
234 "No Coordinates found provided by the user.");
236 TEUCHOS_TEST_FOR_EXCEPTION(fineCoords->getNumVectors() != 3,
237 Exceptions::RuntimeError,
238 "Three coordinates arrays must be supplied if " 239 "line detection orientation not given.");
240 ArrayRCP<typename Teuchos::ScalarTraits<Scalar>::coordinateType> x =
241 fineCoords->getDataNonConst(0);
242 ArrayRCP<typename Teuchos::ScalarTraits<Scalar>::coordinateType> y =
243 fineCoords->getDataNonConst(1);
244 ArrayRCP<typename Teuchos::ScalarTraits<Scalar>::coordinateType> z =
245 fineCoords->getDataNonConst(2);
249 typename Teuchos::ScalarTraits<Scalar>::coordinateType zval_max =
250 -Teuchos::ScalarTraits<
251 typename Teuchos::ScalarTraits<Scalar>::coordinateType>::one() /
252 Teuchos::ScalarTraits<
253 typename Teuchos::ScalarTraits<Scalar>::coordinateType>::sfmin();
254 typename Teuchos::ScalarTraits<Scalar>::coordinateType zval_min =
255 Teuchos::ScalarTraits<
256 typename Teuchos::ScalarTraits<Scalar>::coordinateType>::one() /
257 Teuchos::ScalarTraits<
258 typename Teuchos::ScalarTraits<Scalar>::coordinateType>::sfmin();
259 for (
auto it = z.begin(); it != z.end(); ++it) {
266 LO myCoarseZLayers = Teuchos::as<LO>(CoarseNumZLayers);
268 ArrayRCP<typename Teuchos::ScalarTraits<Scalar>::coordinateType>
269 myZLayerCoords = Teuchos::arcp<
270 typename Teuchos::ScalarTraits<Scalar>::coordinateType>(
272 if (myCoarseZLayers == 1) {
273 myZLayerCoords[0] = zval_min;
275 typename Teuchos::ScalarTraits<Scalar>::coordinateType dz =
276 (zval_max - zval_min) / (myCoarseZLayers - 1);
277 for (LO k = 0; k < myCoarseZLayers; ++k) {
278 myZLayerCoords[k] = k * dz;
286 LO numVertLines = Nnodes / FineNumZLayers;
287 LO numLocalCoarseNodes = numVertLines * myCoarseZLayers;
289 RCP<const Map> coarseCoordMap = MapFactory::Build(
290 fineCoords->getMap()->lib(),
291 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
292 Teuchos::as<size_t>(numLocalCoarseNodes),
293 fineCoords->getMap()->getIndexBase(), fineCoords->getMap()->getComm());
294 RCP<xdMV> coarseCoords = Xpetra::MultiVectorFactory<
295 typename Teuchos::ScalarTraits<Scalar>::coordinateType, LO, GO,
296 NO>::Build(coarseCoordMap, fineCoords->getNumVectors());
297 coarseCoords->putScalar(-1.0);
298 ArrayRCP<typename Teuchos::ScalarTraits<Scalar>::coordinateType> cx =
299 coarseCoords->getDataNonConst(0);
300 ArrayRCP<typename Teuchos::ScalarTraits<Scalar>::coordinateType> cy =
301 coarseCoords->getDataNonConst(1);
302 ArrayRCP<typename Teuchos::ScalarTraits<Scalar>::coordinateType> cz =
303 coarseCoords->getDataNonConst(2);
306 LO cntCoarseNodes = 0;
307 for (LO vt = 0; vt < TVertLineId.size(); ++vt) {
309 LO curVertLineId = TVertLineId[vt];
311 if (cx[curVertLineId * myCoarseZLayers] == -1.0 &&
312 cy[curVertLineId * myCoarseZLayers] == -1.0) {
314 for (LO n = 0; n < myCoarseZLayers; ++n) {
315 cx[curVertLineId * myCoarseZLayers + n] = x[vt];
316 cy[curVertLineId * myCoarseZLayers + n] = y[vt];
317 cz[curVertLineId * myCoarseZLayers + n] = myZLayerCoords[n];
319 cntCoarseNodes += myCoarseZLayers;
322 TEUCHOS_TEST_FOR_EXCEPTION(cntCoarseNodes > numLocalCoarseNodes,
323 Exceptions::RuntimeError,
324 "number of coarse nodes is inconsistent.");
325 if (cntCoarseNodes == numLocalCoarseNodes)
330 Set(coarseLevel,
"Coordinates", coarseCoords);
336 void SemiCoarsenPFactory_kokkos<
338 Kokkos::Compat::KokkosDeviceWrapperNode<DeviceType>>::
339 BuildSemiCoarsenP(Level &coarseLevel,
const LO NFRows,
const LO NFNodes,
340 const LO DofsPerNode,
const LO NFLayers,
341 const LO NCLayers,
const ArrayRCP<LO> LayerId,
342 const ArrayRCP<LO> VertLineId,
const RCP<Matrix> &Amat,
343 const RCP<MultiVector> fineNullspace, RCP<Matrix> &P,
344 RCP<MultiVector> &coarseNullspace)
const {
345 SubFactoryMonitor m2(*
this,
"BuildSemiCoarsenP", coarseLevel);
346 using impl_SC =
typename Kokkos::ArithTraits<SC>::val_type;
347 using impl_ATS = Kokkos::ArithTraits<impl_SC>;
353 const auto FCol2LayerVector =
354 Xpetra::VectorFactory<LO, LO, GO, NO>::Build(Amat->getColMap());
355 const auto localTemp =
356 Xpetra::VectorFactory<LO, LO, GO, NO>::Build(Amat->getDomainMap());
357 RCP<const Import> importer = Amat->getCrsGraph()->getImporter();
358 if (importer == Teuchos::null)
359 importer = ImportFactory::Build(Amat->getDomainMap(), Amat->getColMap());
362 const auto localTempHost = localTemp->getHostLocalView(Xpetra::Access::ReadWrite);
363 for (
int row = 0; row < NFRows; row++)
364 localTempHost(row, 0) = LayerId[row / DofsPerNode];
365 const auto localTempView = localTemp->getDeviceLocalView(Xpetra::Access::ReadWrite);
366 Kokkos::deep_copy(localTempView, localTempHost);
367 FCol2LayerVector->doImport(*localTemp, *(importer), Xpetra::INSERT);
369 const auto FCol2LayerView = FCol2LayerVector->getDeviceLocalView(Xpetra::Access::ReadOnly);
370 const auto FCol2Layer = Kokkos::subview(FCol2LayerView, Kokkos::ALL(), 0);
374 const auto FCol2DofVector =
375 Xpetra::VectorFactory<LO, LO, GO, NO>::Build(Amat->getColMap());
378 const auto localTempHost = localTemp->getHostLocalView(Xpetra::Access::ReadWrite);
379 for (
int row = 0; row < NFRows; row++)
380 localTempHost(row, 0) = row % DofsPerNode;
381 const auto localTempView = localTemp->getDeviceLocalView(Xpetra::Access::ReadWrite);
382 Kokkos::deep_copy(localTempView, localTempHost);
383 FCol2DofVector->doImport(*localTemp, *(importer), Xpetra::INSERT);
385 const auto FCol2DofView = FCol2DofVector->getDeviceLocalView(Xpetra::Access::ReadOnly);
386 const auto FCol2Dof = Kokkos::subview(FCol2DofView, Kokkos::ALL(), 0);
392 NVertLines = VertLineId[0];
393 for (
int node = 1; node < NFNodes; ++node)
394 if (VertLineId[node] > NVertLines)
395 NVertLines = VertLineId[node];
399 LOView2D LineLayer2Node(
"LineLayer2Node", NVertLines, NFLayers);
400 typename LOView2D::HostMirror LineLayer2NodeHost =
401 Kokkos::create_mirror_view(LineLayer2Node);
402 for (
int node = 0; node < NFNodes; ++node)
403 LineLayer2NodeHost(VertLineId[node], LayerId[node]) = node;
404 Kokkos::deep_copy(LineLayer2Node, LineLayer2NodeHost);
407 LOView1D CLayer2FLayer(
"CLayer2FLayer", NCLayers);
408 typename LOView1D::HostMirror CLayer2FLayerHost =
409 Kokkos::create_mirror_view(CLayer2FLayer);
410 using coordT =
typename Teuchos::ScalarTraits<Scalar>::coordinateType;
411 const LO FirstStride =
412 (LO)ceil(((coordT)(NFLayers + 1)) / ((coordT)(NCLayers + 1)));
413 const coordT RestStride =
414 ((coordT)(NFLayers - FirstStride + 1)) / ((coordT)NCLayers);
416 (LO)floor((((coordT)(NFLayers - FirstStride + 1)) / RestStride) + .00001);
417 TEUCHOS_TEST_FOR_EXCEPTION(NCLayers != NCpts, Exceptions::RuntimeError,
418 "sizes do not match.");
419 coordT stride = (coordT)FirstStride;
420 for (
int clayer = 0; clayer < NCLayers; ++clayer) {
421 CLayer2FLayerHost(clayer) = (LO)floor(stride) - 1;
422 stride += RestStride;
424 Kokkos::deep_copy(CLayer2FLayer, CLayer2FLayerHost);
428 int MaxStencilSize = 1;
429 LOView1D CLayer2StartLayer(
"CLayer2StartLayer", NCLayers);
430 LOView1D CLayer2StencilSize(
"CLayer2StencilSize", NCLayers);
431 typename LOView1D::HostMirror CLayer2StartLayerHost =
432 Kokkos::create_mirror_view(CLayer2StartLayer);
433 typename LOView1D::HostMirror CLayer2StencilSizeHost =
434 Kokkos::create_mirror_view(CLayer2StencilSize);
435 for (
int clayer = 0; clayer < NCLayers; ++clayer) {
436 const int startLayer = (clayer > 0) ? CLayer2FLayerHost(clayer - 1) + 1 : 0;
437 const int stencilSize = (clayer < NCLayers - 1)
438 ? CLayer2FLayerHost(clayer + 1) - startLayer
439 : NFLayers - startLayer;
441 if (MaxStencilSize < stencilSize)
442 MaxStencilSize = stencilSize;
443 CLayer2StartLayerHost(clayer) = startLayer;
444 CLayer2StencilSizeHost(clayer) = stencilSize;
446 Kokkos::deep_copy(CLayer2StartLayer, CLayer2StartLayerHost);
447 Kokkos::deep_copy(CLayer2StencilSize, CLayer2StencilSizeHost);
454 int Nmax = MaxStencilSize * DofsPerNode;
456 "BandMat", NVertLines, Nmax, Nmax);
458 "BandSol", NVertLines, Nmax, DofsPerNode);
464 for (
int clayer = 0; clayer < NCLayers; ++clayer)
465 NnzP += CLayer2StencilSizeHost(clayer);
466 NnzP *= NVertLines * DofsPerNode * DofsPerNode;
475 Kokkos::create_mirror_view(Pptr);
476 Kokkos::deep_copy(PptrHost, 0);
477 for (
int line = 0; line < NVertLines; ++line) {
478 for (
int clayer = 0; clayer < NCLayers; ++clayer) {
479 const int stencilSize = CLayer2StencilSizeHost(clayer);
480 const int startLayer = CLayer2StartLayerHost(clayer);
481 for (
int snode = 0; snode < stencilSize; ++snode) {
482 for (
int dofi = 0; dofi < DofsPerNode; ++dofi) {
483 const int layer = startLayer + snode;
484 const int AmatBlkRow = LineLayer2NodeHost(line, layer);
485 const int AmatRow = AmatBlkRow * DofsPerNode + dofi;
486 PptrHost(AmatRow + 1) += DofsPerNode;
491 for (
int i = 2; i < NFRows + 1; ++i)
492 PptrHost(i) += PptrHost(i - 1);
493 TEUCHOS_TEST_FOR_EXCEPTION(NnzP != (
int)PptrHost(NFRows),
494 Exceptions::RuntimeError,
495 "Number of nonzeros in P does not match");
496 Kokkos::deep_copy(Pptr, PptrHost);
501 "layerBuckets", NFLayers);
502 Kokkos::deep_copy(layerBuckets, 0);
503 LOView2D CLayerSNode2PptrOffset(
"CLayerSNode2PptrOffset", NCLayers,
505 typename LOView2D::HostMirror CLayerSNode2PptrOffsetHost =
506 Kokkos::create_mirror_view(CLayerSNode2PptrOffset);
507 for (
int clayer = 0; clayer < NCLayers; ++clayer) {
508 const int stencilSize = CLayer2StencilSizeHost(clayer);
509 const int startLayer = CLayer2StartLayerHost(clayer);
510 for (
int snode = 0; snode < stencilSize; ++snode) {
511 const int layer = startLayer + snode;
512 CLayerSNode2PptrOffsetHost(clayer, snode) = layerBuckets(layer);
513 layerBuckets(layer)++;
516 Kokkos::deep_copy(CLayerSNode2PptrOffset, CLayerSNode2PptrOffsetHost);
519 SubFactoryMonitor m3(*
this,
"Fill P", coarseLevel);
521 const auto localAmat = Amat->getLocalMatrixDevice();
522 const auto zero = impl_ATS::zero();
523 const auto one = impl_ATS::one();
526 Kokkos::parallel_for(
527 "MueLu::SemiCoarsenPFactory_kokkos::BuildSemiCoarsenP Fill P",
528 range_policy(0, NVertLines), KOKKOS_LAMBDA(
const int line) {
529 for (
int clayer = 0; clayer < NCLayers; ++clayer) {
533 Kokkos::subview(BandSol, line, Kokkos::ALL(), Kokkos::ALL());
534 for (
int row = 0; row < Nmax; ++row)
535 for (
int dof = 0; dof < DofsPerNode; ++dof)
536 bandSol(row, dof) = zero;
539 const int stencilSize = CLayer2StencilSize(clayer);
540 const int N = stencilSize * DofsPerNode;
542 Kokkos::subview(BandMat, line, Kokkos::ALL(), Kokkos::ALL());
543 for (
int row = 0; row < Nmax; ++row)
544 for (
int col = 0; col < Nmax; ++col)
546 (row == col && row >= N) ? one : zero;
549 const int flayer = CLayer2FLayer(clayer);
550 const int startLayer = CLayer2StartLayer(clayer);
551 for (
int snode = 0; snode < stencilSize; ++snode) {
553 const int layer = startLayer + snode;
554 if (layer == flayer) {
555 for (
int dof = 0; dof < DofsPerNode; ++dof) {
556 const int row = snode * DofsPerNode + dof;
557 bandMat(row, row) = one;
558 bandSol(row, dof) = one;
561 const int AmatBlkRow = LineLayer2Node(line, layer);
562 for (
int dofi = 0; dofi < DofsPerNode; ++dofi) {
565 const int AmatRow = AmatBlkRow * DofsPerNode + dofi;
566 const auto localAmatRow = localAmat.rowConst(AmatRow);
567 const int AmatRowLeng = localAmatRow.length;
569 const int row = snode * DofsPerNode + dofi;
570 for (
int dofj = 0; dofj < DofsPerNode; ++dofj) {
571 const int col = snode * DofsPerNode + dofj;
576 for (
int i = 0; i < AmatRowLeng; ++i) {
577 const int colidx = localAmatRow.colidx(i);
578 if (FCol2Layer(colidx) == layer &&
579 FCol2Dof(colidx) == dofj)
580 val += localAmatRow.value(i);
582 bandMat(row, col) = val;
588 for (
int i = 0; i < AmatRowLeng; ++i) {
589 const int colidx = localAmatRow.colidx(i);
590 if (FCol2Layer(colidx) == layer - 1 &&
591 FCol2Dof(colidx) == dofj)
592 val += localAmatRow.value(i);
594 bandMat(row, col - DofsPerNode) = val;
597 if (snode < stencilSize - 1) {
601 for (
int i = 0; i < AmatRowLeng; ++i) {
602 const int colidx = localAmatRow.colidx(i);
603 if (FCol2Layer(colidx) == layer + 1 &&
604 FCol2Dof(colidx) == dofj)
605 val += localAmatRow.value(i);
607 bandMat(row, col + DofsPerNode) = val;
615 namespace KB = KokkosBatched;
616 using lu_type =
typename KB::SerialLU<KB::Algo::LU::Unblocked>;
617 lu_type::invoke(bandMat);
619 typename KB::SerialTrsm<KB::Side::Left, KB::Uplo::Lower,
620 KB::Trans::NoTranspose, KB::Diag::Unit,
621 KB::Algo::Trsm::Unblocked>;
622 trsv_l_type::invoke(one, bandMat, bandSol);
623 using trsv_u_type =
typename KB::SerialTrsm<
624 KB::Side::Left, KB::Uplo::Upper, KB::Trans::NoTranspose,
625 KB::Diag::NonUnit, KB::Algo::Trsm::Unblocked>;
626 trsv_u_type::invoke(one, bandMat, bandSol);
629 for (
int snode = 0; snode < stencilSize; ++snode) {
630 for (
int dofj = 0; dofj < DofsPerNode; ++dofj) {
631 for (
int dofi = 0; dofi < DofsPerNode; ++dofi) {
632 const int layer = startLayer + snode;
633 const int AmatBlkRow = LineLayer2Node(line, layer);
634 const int AmatRow = AmatBlkRow * DofsPerNode + dofi;
636 const int pptrOffset = CLayerSNode2PptrOffset(clayer, snode);
638 Pptr(AmatRow) + pptrOffset * DofsPerNode + dofj;
641 line * NCLayers + clayer;
642 Pcols(pptr) = col * DofsPerNode + dofj;
643 Pvals(pptr) = bandSol(snode * DofsPerNode + dofi, dofj);
652 RCP<const Map> rowMap = Amat->getRowMap();
653 Xpetra::global_size_t GNdofs = rowMap->getGlobalNumElements();
654 Xpetra::global_size_t itemp = GNdofs / NFLayers;
655 std::vector<size_t> stridingInfo_{(size_t)DofsPerNode};
656 RCP<const Map> coarseMap = StridedMapFactory::Build(
657 rowMap->lib(), NCLayers * itemp, NCLayers * NVertLines * DofsPerNode, 0,
658 stridingInfo_, rowMap->getComm(), -1, 0);
659 P = rcp(
new CrsMatrixWrap(rowMap, coarseMap, 0));
660 RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
661 PCrs->setAllValues(Pptr, Pcols, Pvals);
662 PCrs->expertStaticFillComplete(coarseMap, Amat->getDomainMap());
665 if (Amat->IsView(
"stridedMaps") ==
true)
666 P->CreateView(
"stridedMaps", Amat->getRowMap(
"stridedMaps"), coarseMap);
668 P->CreateView(
"stridedMaps", P->getRangeMap(), coarseMap);
672 MultiVectorFactory::Build(coarseMap, fineNullspace->getNumVectors());
673 const int numVectors = fineNullspace->getNumVectors();
674 const auto fineNullspaceView = fineNullspace->getDeviceLocalView(Xpetra::Access::ReadOnly);
675 const auto coarseNullspaceView = coarseNullspace->getDeviceLocalView(Xpetra::Access::ReadWrite);
677 Kokkos::parallel_for(
678 "MueLu::SemiCoarsenPFactory_kokkos::BuildSemiCoarsenP Inject Nullspace",
679 range_policy(0, NVertLines), KOKKOS_LAMBDA(
const int line) {
680 for (
int clayer = 0; clayer < NCLayers; ++clayer) {
681 const int layer = CLayer2FLayer(clayer);
682 const int AmatBlkRow =
683 LineLayer2Node(line, layer);
685 line * NCLayers + clayer;
686 for (
int k = 0; k < numVectors; ++k) {
687 for (
int dofi = 0; dofi < DofsPerNode; ++dofi) {
688 const int fRow = AmatBlkRow * DofsPerNode + dofi;
689 const int cRow = col * DofsPerNode + dofi;
690 coarseNullspaceView(cRow, k) = fineNullspaceView(fRow, k);
699 #define MUELU_SEMICOARSENPFACTORY_KOKKOS_SHORT 700 #endif // HAVE_MUELU_KOKKOS_REFACTOR 701 #endif // MUELU_SEMICOARSENPFACTORY_KOKKOS_DEF_HPP MueLu::DefaultLocalOrdinal LocalOrdinal
Namespace for MueLu classes and methods.
static const NoFactory * get()
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
static const Teuchos::ParameterEntry & getEntry(const std::string &name)
Returns default entry from the "master" list corresponding to the specified name. ...