MueLu Version of the Day
MueLu_Utilities_kokkos_def.hpp
Go to the documentation of this file.
1// @HEADER
2//
3// ***********************************************************************
4//
5// MueLu: A package for multigrid based preconditioning
6// Copyright 2012 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact
39// Jonathan Hu (jhu@sandia.gov)
40// Andrey Prokopenko (aprokop@sandia.gov)
41// Ray Tuminaro (rstumin@sandia.gov)
42//
43// ***********************************************************************
44//
45// @HEADER
46#ifndef MUELU_UTILITIES_KOKKOS_DEF_HPP
47#define MUELU_UTILITIES_KOKKOS_DEF_HPP
48
49#include <algorithm>
50#include <Teuchos_DefaultComm.hpp>
51#include <Teuchos_ParameterList.hpp>
52
53#include "MueLu_ConfigDefs.hpp"
54
55#ifdef HAVE_MUELU_EPETRA
56# ifdef HAVE_MPI
57# include "Epetra_MpiComm.h"
58# endif
59#endif
60
61#include <Kokkos_Core.hpp>
62#include <KokkosSparse_CrsMatrix.hpp>
63#include <KokkosSparse_getDiagCopy.hpp>
64
65#if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_EPETRAEXT)
66#include <EpetraExt_MatrixMatrix.h>
67#include <EpetraExt_RowMatrixOut.h>
68#include <EpetraExt_MultiVectorOut.h>
69#include <EpetraExt_CrsMatrixIn.h>
70#include <EpetraExt_MultiVectorIn.h>
71#include <EpetraExt_BlockMapIn.h>
72#include <Xpetra_EpetraUtils.hpp>
73#include <Xpetra_EpetraMultiVector.hpp>
74#include <EpetraExt_BlockMapOut.h>
75#endif
76
77#ifdef HAVE_MUELU_TPETRA
78#include <MatrixMarket_Tpetra.hpp>
79#include <Tpetra_RowMatrixTransposer.hpp>
80#include <TpetraExt_MatrixMatrix.hpp>
81#include <Xpetra_TpetraMultiVector.hpp>
82#include <Xpetra_TpetraCrsMatrix.hpp>
83#include <Xpetra_TpetraBlockCrsMatrix.hpp>
84#endif
85
86#ifdef HAVE_MUELU_EPETRA
87#include <Xpetra_EpetraMap.hpp>
88#endif
89
90#include <Xpetra_BlockedCrsMatrix.hpp>
91#include <Xpetra_DefaultPlatform.hpp>
92#include <Xpetra_Import.hpp>
93#include <Xpetra_ImportFactory.hpp>
94#include <Xpetra_Map.hpp>
95#include <Xpetra_MapFactory.hpp>
96#include <Xpetra_Matrix.hpp>
97#include <Xpetra_MatrixMatrix.hpp>
98#include <Xpetra_MatrixFactory.hpp>
99#include <Xpetra_MultiVector.hpp>
100#include <Xpetra_MultiVectorFactory.hpp>
101#include <Xpetra_Operator.hpp>
102#include <Xpetra_Vector.hpp>
103#include <Xpetra_VectorFactory.hpp>
104
105#include <MueLu_Utilities_kokkos_decl.hpp>
106
107#include <KokkosKernels_Handle.hpp>
108#include <KokkosGraph_RCM.hpp>
109
110
111namespace MueLu {
112
113
114 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
115 Teuchos::RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
116 Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
117 GetMatrixDiagonal(const Matrix& A) {
118 const auto rowMap = A.getRowMap();
119 auto diag = Xpetra::VectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(rowMap, true);
120
121 A.getLocalDiagCopy(*diag);
122
123 return diag;
124 } //GetMatrixDiagonal
125
126 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
127 Teuchos::RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
128 Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
129 GetMatrixDiagonalInverse(const Matrix& A, Magnitude tol, const bool doLumped) {
130 Teuchos::TimeMonitor MM = *Teuchos::TimeMonitor::getNewTimer("Utilities_kokkos::GetMatrixDiagonalInverse");
131 // Some useful type definitions
132 using local_matrix_type = typename Matrix::local_matrix_type;
133 using local_graph_type = typename local_matrix_type::staticcrsgraph_type;
134 using value_type = typename local_matrix_type::value_type;
135 using ordinal_type = typename local_matrix_type::ordinal_type;
136 using execution_space = typename local_matrix_type::execution_space;
137 using memory_space = typename local_matrix_type::memory_space;
138 // Be careful with this one, if using Kokkos::ArithTraits<Scalar>
139 // you are likely to run into errors when handling std::complex<>
140 // a good way to work around that is to use the following:
141 // using KAT = Kokkos::ArithTraits<Kokkos::ArithTraits<Scalar>::val_type> >
142 // here we have: value_type = Kokkos::ArithTraits<Scalar>::val_type
143 using KAT = Kokkos::ArithTraits<value_type>;
144
145 // Get/Create distributed objects
146 RCP<const Map> rowMap = A.getRowMap();
147 RCP<Vector> diag = VectorFactory::Build(rowMap,false);
148
149 // Now generate local objects
150 local_matrix_type localMatrix = A.getLocalMatrixDevice();
151 auto diagVals = diag->getDeviceLocalView(Xpetra::Access::ReadWrite);
152
153 ordinal_type numRows = localMatrix.graph.numRows();
154
155 // Note: 2019-11-21, LBV
156 // This could be implemented with a TeamPolicy over the rows
157 // and a TeamVectorRange over the entries in a row if performance
158 // becomes more important here.
159 if (!doLumped)
160 Kokkos::parallel_for("Utilities_kokkos::GetMatrixDiagonalInverse",
161 Kokkos::RangePolicy<ordinal_type, execution_space>(0, numRows),
162 KOKKOS_LAMBDA(const ordinal_type rowIdx) {
163 bool foundDiagEntry = false;
164 auto myRow = localMatrix.rowConst(rowIdx);
165 for(ordinal_type entryIdx = 0; entryIdx < myRow.length; ++entryIdx) {
166 if(myRow.colidx(entryIdx) == rowIdx) {
167 foundDiagEntry = true;
168 if(KAT::magnitude(myRow.value(entryIdx)) > KAT::magnitude(tol)) {
169 diagVals(rowIdx, 0) = KAT::one() / myRow.value(entryIdx);
170 } else {
171 diagVals(rowIdx, 0) = KAT::zero();
172 }
173 break;
174 }
175 }
176
177 if(!foundDiagEntry) {diagVals(rowIdx, 0) = KAT::zero();}
178 });
179 else
180 Kokkos::parallel_for("Utilities_kokkos::GetMatrixDiagonalInverse",
181 Kokkos::RangePolicy<ordinal_type, execution_space>(0, numRows),
182 KOKKOS_LAMBDA(const ordinal_type rowIdx) {
183 auto myRow = localMatrix.rowConst(rowIdx);
184 for(ordinal_type entryIdx = 0; entryIdx < myRow.length; ++entryIdx) {
185 diagVals(rowIdx, 0) += KAT::magnitude(myRow.value(entryIdx));
186 }
187 if(KAT::magnitude(diagVals(rowIdx, 0)) > KAT::magnitude(tol))
188 diagVals(rowIdx, 0) = KAT::one() / diagVals(rowIdx, 0);
189 else
190 diagVals(rowIdx, 0) = KAT::zero();
191
192 });
193
194 return diag;
195 } //GetMatrixDiagonalInverse
196
197 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
198 RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
199 Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
200 GetMatrixOverlappedDiagonal(const Matrix& A) {
201 // FIXME_KOKKOS
202 RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
203 RCP<Vector> localDiag = VectorFactory::Build(rowMap);
204
205 const CrsMatrixWrap* crsOp = dynamic_cast<const CrsMatrixWrap*>(&A);
206 if (crsOp != NULL) {
207 Teuchos::ArrayRCP<size_t> offsets;
208 crsOp->getLocalDiagOffsets(offsets);
209 crsOp->getLocalDiagCopy(*localDiag,offsets());
210 }
211 else {
212 auto localDiagVals = localDiag->getDeviceLocalView(Xpetra::Access::ReadWrite);
213 const auto diagVals = GetMatrixDiagonal(A)->getDeviceLocalView(Xpetra::Access::ReadOnly);
214 Kokkos::deep_copy(localDiagVals, diagVals);
215 }
216
217 RCP<Vector> diagonal = VectorFactory::Build(colMap);
218 RCP< const Import> importer;
219 importer = A.getCrsGraph()->getImporter();
220 if (importer == Teuchos::null) {
221 importer = ImportFactory::Build(rowMap, colMap);
222 }
223 diagonal->doImport(*localDiag, *(importer), Xpetra::INSERT);
224
225 return diagonal;
226 } //GetMatrixOverlappedDiagonal
227
228 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
229 void Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
230 MyOldScaleMatrix(Matrix& Op, const Teuchos::ArrayRCP<const SC>& scalingVector,
231 bool doInverse, bool doFillComplete, bool doOptimizeStorage)
232 {
233 SC one = Teuchos::ScalarTraits<SC>::one();
234 Teuchos::ArrayRCP<SC> sv(scalingVector.size());
235 if (doInverse) {
236 for (int i = 0; i < scalingVector.size(); ++i)
237 sv[i] = one / scalingVector[i];
238 } else {
239 for (int i = 0; i < scalingVector.size(); ++i)
240 sv[i] = scalingVector[i];
241 }
242
243 switch (Op.getRowMap()->lib()) {
244 case Xpetra::UseTpetra:
245 MyOldScaleMatrix_Tpetra(Op, sv, doFillComplete, doOptimizeStorage);
246 break;
247
248 case Xpetra::UseEpetra:
249 // FIXME_KOKKOS
250 // Utils2_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix_Epetra(Op, sv, doFillComplete, doOptimizeStorage);
251 throw std::runtime_error("FIXME");
252#ifndef __NVCC__ //prevent nvcc warning
253 break;
254#endif
255
256 default:
257 throw Exceptions::RuntimeError("Only Epetra and Tpetra matrices can be scaled.");
258#ifndef __NVCC__ //prevent nvcc warning
259 break;
260#endif
261 }
262 }
263
264 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
265 void Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix_Epetra(Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& /* Op */, const Teuchos::ArrayRCP<Scalar>& /* scalingVector */, bool /* doFillComplete */, bool /* doOptimizeStorage */) {
266 throw Exceptions::RuntimeError("MyOldScaleMatrix_Epetra: Epetra needs SC=double and LO=GO=int.");
267 }
268
269 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
270 void Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix_Tpetra(Matrix& Op, const Teuchos::ArrayRCP<SC>& scalingVector,
271 bool doFillComplete,
272 bool doOptimizeStorage)
273 {
274#ifdef HAVE_MUELU_TPETRA
275 try {
276 Tpetra::CrsMatrix<SC,LO,GO,NO>& tpOp = Op2NonConstTpetraCrs(Op);
277
278 const RCP<const Tpetra::Map<LO,GO,NO> > rowMap = tpOp.getRowMap();
279 const RCP<const Tpetra::Map<LO,GO,NO> > domainMap = tpOp.getDomainMap();
280 const RCP<const Tpetra::Map<LO,GO,NO> > rangeMap = tpOp.getRangeMap();
281
282 size_t maxRowSize = tpOp.getNodeMaxNumRowEntries();
283 if (maxRowSize == Teuchos::as<size_t>(-1)) // hasn't been determined yet
284 maxRowSize = 20;
285
286 if (tpOp.isFillComplete())
287 tpOp.resumeFill();
288
289 if (Op.isLocallyIndexed() == true) {
290 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::local_inds_host_view_type cols;
291 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::values_host_view_type vals;
292
293 for (size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
294 tpOp.getLocalRowView(i, cols, vals);
295 size_t nnz = tpOp.getNumEntriesInLocalRow(i);
296 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::nonconst_values_host_view_type scaledVals("scaledVals", nnz);
297 for (size_t j = 0; j < nnz; ++j)
298 scaledVals[j] = scalingVector[i]*vals[j];
299
300 if (nnz > 0) {
301 tpOp.replaceLocalValues(i, cols, scaledVals);
302 }
303 } //for (size_t i=0; ...
304
305 } else {
306 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::global_inds_host_view_type cols;
307 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::values_host_view_type vals;
308
309 for (size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
310 GO gid = rowMap->getGlobalElement(i);
311 tpOp.getGlobalRowView(gid, cols, vals);
312 size_t nnz = tpOp.getNumEntriesInGlobalRow(gid);
313 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::nonconst_values_host_view_type scaledVals("scaledVals", nnz);
314
315 // FIXME FIXME FIXME FIXME FIXME FIXME
316 for (size_t j = 0; j < nnz; ++j)
317 scaledVals[j] = scalingVector[i]*vals[j]; //FIXME i or gid?
318
319 if (nnz > 0) {
320 tpOp.replaceGlobalValues(gid, cols, scaledVals);
321 }
322 } //for (size_t i=0; ...
323 }
324
325 if (doFillComplete) {
326 if (domainMap == Teuchos::null || rangeMap == Teuchos::null)
327 throw Exceptions::RuntimeError("In Utils_kokkos::Scaling: cannot fillComplete because the domain and/or range map hasn't been defined");
328
329 RCP<Teuchos::ParameterList> params = rcp(new Teuchos::ParameterList());
330 params->set("Optimize Storage", doOptimizeStorage);
331 params->set("No Nonlocal Changes", true);
332 Op.fillComplete(Op.getDomainMap(), Op.getRangeMap(), params);
333 }
334 } catch(...) {
335 throw Exceptions::RuntimeError("Only Tpetra::CrsMatrix types can be scaled (Err.1)");
336 }
337#else
338 throw Exceptions::RuntimeError("Matrix scaling is not possible because Tpetra has not been enabled.");
339#endif
340 } //MyOldScaleMatrix_Tpetra()
341
342
343 template <class SC, class LO, class GO, class NO>
344 Kokkos::View<bool*, typename NO::device_type>
345 DetectDirichletRows(const Xpetra::Matrix<SC,LO,GO,NO>& A,
346 const typename Teuchos::ScalarTraits<SC>::magnitudeType& tol,
347 const bool count_twos_as_dirichlet) {
348 using impl_scalar_type = typename Kokkos::ArithTraits<SC>::val_type;
349 using ATS = Kokkos::ArithTraits<impl_scalar_type>;
350 using range_type = Kokkos::RangePolicy<LO, typename NO::execution_space>;
351
352 auto localMatrix = A.getLocalMatrixDevice();
353 LO numRows = A.getNodeNumRows();
354
355 Kokkos::View<bool*, typename NO::device_type> boundaryNodes(Kokkos::ViewAllocateWithoutInitializing("boundaryNodes"), numRows);
356 if (count_twos_as_dirichlet)
357 Kokkos::parallel_for("MueLu:Utils::DetectDirichletRows_Twos_As_Dirichlet", range_type(0,numRows),
358 KOKKOS_LAMBDA(const LO row) {
359 auto rowView = localMatrix.row(row);
360 auto length = rowView.length;
361
362 boundaryNodes(row) = true;
363 if (length > 2) {
364 decltype(length) colID;
365 for (colID = 0; colID < length; colID++)
366 if ((rowView.colidx(colID) != row) &&
367 (ATS::magnitude(rowView.value(colID)) > tol)) {
368 if (!boundaryNodes(row))
369 break;
370 boundaryNodes(row) = false;
371 }
372 if (colID == length)
373 boundaryNodes(row) = true;
374 }
375 });
376 else
377 Kokkos::parallel_for("MueLu:Utils::DetectDirichletRows", range_type(0,numRows),
378 KOKKOS_LAMBDA(const LO row) {
379 auto rowView = localMatrix.row(row);
380 auto length = rowView.length;
381
382 boundaryNodes(row) = true;
383 for (decltype(length) colID = 0; colID < length; colID++)
384 if ((rowView.colidx(colID) != row) &&
385 (ATS::magnitude(rowView.value(colID)) > tol)) {
386 boundaryNodes(row) = false;
387 break;
388 }
389 });
390
391 return boundaryNodes;
392 }
393
394 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
395 Kokkos::View<bool*, typename Node::device_type>
396 Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
397 DetectDirichletRows(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, const typename Teuchos::ScalarTraits<Scalar>::magnitudeType& tol, const bool count_twos_as_dirichlet) {
398 return MueLu::DetectDirichletRows<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, tol, count_twos_as_dirichlet);
399 }
400
401 template <class Node>
402 Kokkos::View<bool*, typename Node::device_type>
403 Utilities_kokkos<double,int,int,Node>::
404 DetectDirichletRows(const Xpetra::Matrix<double,int,int,Node>& A, const typename Teuchos::ScalarTraits<double>::magnitudeType& tol, const bool count_twos_as_dirichlet) {
405 return MueLu::DetectDirichletRows<double,int,int,Node>(A, tol,count_twos_as_dirichlet);
406 }
407
408
409 template <class SC, class LO, class GO, class NO>
410 Kokkos::View<bool*, typename NO::device_type>
411 DetectDirichletCols(const Xpetra::Matrix<SC,LO,GO,NO>& A,
412 const Kokkos::View<const bool*, typename NO::device_type>& dirichletRows) {
413 using ATS = Kokkos::ArithTraits<SC>;
414 using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
415 using range_type = Kokkos::RangePolicy<LO, typename NO::execution_space>;
416
417 SC zero = ATS::zero();
418 SC one = ATS::one();
419
420 auto localMatrix = A.getLocalMatrixDevice();
421 LO numRows = A.getNodeNumRows();
422
423 Teuchos::RCP<const Xpetra::Map<LO,GO,NO> > domMap = A.getDomainMap();
424 Teuchos::RCP<const Xpetra::Map<LO,GO,NO> > colMap = A.getColMap();
425 Teuchos::RCP<Xpetra::MultiVector<SC,LO,GO,NO> > myColsToZero = Xpetra::MultiVectorFactory<SC,LO,GO,NO>::Build(colMap,1);
426 myColsToZero->putScalar(zero);
427 auto myColsToZeroView = myColsToZero->getDeviceLocalView(Xpetra::Access::ReadWrite);
428 // Find all local column indices that are in Dirichlet rows, record in myColsToZero as 1.0
429 Kokkos::parallel_for("MueLu:Utils::DetectDirichletCols1", range_type(0,numRows),
430 KOKKOS_LAMBDA(const LO row) {
431 if (dirichletRows(row)) {
432 auto rowView = localMatrix.row(row);
433 auto length = rowView.length;
434
435 for (decltype(length) colID = 0; colID < length; colID++)
436 myColsToZeroView(rowView.colidx(colID),0) = one;
437 }
438 });
439
440 Teuchos::RCP<Xpetra::MultiVector<SC,LO,GO,NO> > globalColsToZero = Xpetra::MultiVectorFactory<SC,LO,GO,NO>::Build(domMap,1);
441 globalColsToZero->putScalar(zero);
442 Teuchos::RCP<Xpetra::Export<LO,GO,NO> > exporter = Xpetra::ExportFactory<LO,GO,NO>::Build(colMap,domMap);
443 // export to domain map
444 globalColsToZero->doExport(*myColsToZero,*exporter,Xpetra::ADD);
445 // import to column map
446 myColsToZero->doImport(*globalColsToZero,*exporter,Xpetra::INSERT);
447
448 auto myCols = myColsToZero->getDeviceLocalView(Xpetra::Access::ReadWrite);
449 size_t numColEntries = colMap->getNodeNumElements();
450 Kokkos::View<bool*, typename NO::device_type> dirichletCols(Kokkos::ViewAllocateWithoutInitializing("dirichletCols"), numColEntries);
451 const typename ATS::magnitudeType eps = 2.0*ATS::eps();
452
453 Kokkos::parallel_for("MueLu:Utils::DetectDirichletCols2", range_type(0,numColEntries),
454 KOKKOS_LAMBDA (const size_t i) {
455 dirichletCols(i) = impl_ATS::magnitude(myCols(i,0))>eps;
456 });
457 return dirichletCols;
458 }
459
460
461 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
462 Kokkos::View<bool*, typename Node::device_type>
463 Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
464 DetectDirichletCols(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
465 const Kokkos::View<const bool*, typename NO::device_type>& dirichletRows) {
466 return MueLu::DetectDirichletCols<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, dirichletRows);
467 }
468
469 template <class Node>
470 Kokkos::View<bool*, typename Node::device_type>
471 Utilities_kokkos<double,int,int,Node>::
472 DetectDirichletCols(const Xpetra::Matrix<double,int,int,Node>& A,
473 const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows) {
474 return MueLu::DetectDirichletCols<double,int,int,Node>(A, dirichletRows);
475 }
476
477
478 // Find Nonzeros in a device view
479 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
480 void FindNonZeros(const typename Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::dual_view_type::t_dev_const_um vals,
481 Kokkos::View<bool*, typename Node::device_type> nonzeros) {
482 using ATS = Kokkos::ArithTraits<Scalar>;
483 using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
484 using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
485 TEUCHOS_ASSERT(vals.extent(0) == nonzeros.extent(0));
486 const typename ATS::magnitudeType eps = 2.0*impl_ATS::eps();
487
488 Kokkos::parallel_for("MueLu:Maxwell1::FindNonZeros", range_type(0,vals.extent(0)),
489 KOKKOS_LAMBDA (const size_t i) {
490 nonzeros(i) = (impl_ATS::magnitude(vals(i,0)) > eps);
491 });
492 }
493
494 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
495 void
496 Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
497 FindNonZeros(const typename Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::dual_view_type::t_dev_const_um vals,
498 Kokkos::View<bool*, typename Node::device_type> nonzeros) {
499 MueLu::FindNonZeros<Scalar,LocalOrdinal,GlobalOrdinal,Node>(vals,nonzeros);
500 }
501
502 template <class Node>
503 void
504 Utilities_kokkos<double,int,int,Node>::
505 FindNonZeros(const typename Xpetra::MultiVector<double,int,int,Node>::dual_view_type::t_dev_const_um vals,
506 Kokkos::View<bool*, typename Node::device_type> nonzeros) {
507 MueLu::FindNonZeros<double,int,int,Node>(vals,nonzeros);
508 }
509
510 // Detect Dirichlet cols and domains
511 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
512 void DetectDirichletColsAndDomains(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
513 const Kokkos::View<bool*, typename Node::device_type> & dirichletRows,
514 Kokkos::View<bool*, typename Node::device_type> dirichletCols,
515 Kokkos::View<bool*, typename Node::device_type> dirichletDomain) {
516 using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
517 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
518 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domMap = A.getDomainMap();
519 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = A.getRowMap();
520 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = A.getColMap();
521 TEUCHOS_ASSERT(dirichletRows.extent(0) == rowMap->getNodeNumElements());
522 TEUCHOS_ASSERT(dirichletCols.extent(0) == colMap->getNodeNumElements());
523 TEUCHOS_ASSERT(dirichletDomain.extent(0) == domMap->getNodeNumElements());
524 RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > myColsToZero = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap, /*zeroOut=*/true);
525 // Find all local column indices that are in Dirichlet rows, record in myColsToZero as 1.0
526 auto myColsToZeroView = myColsToZero->getDeviceLocalView(Xpetra::Access::ReadWrite);
527 auto localMatrix = A.getLocalMatrixDevice();
528 Kokkos::parallel_for("MueLu:Maxwell1::DetectDirichletCols", range_type(0,rowMap->getNodeNumElements()),
529 KOKKOS_LAMBDA(const LocalOrdinal row) {
530 if (dirichletRows(row)) {
531 auto rowView = localMatrix.row(row);
532 auto length = rowView.length;
533
534 for (decltype(length) colID = 0; colID < length; colID++)
535 myColsToZeroView(rowView.colidx(colID),0) = one;
536 }
537 });
538
539 RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > globalColsToZero;
540 RCP<const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer = A.getCrsGraph()->getImporter();
541 if (!importer.is_null()) {
542 globalColsToZero = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(domMap, /*zeroOut=*/true);
543 // export to domain map
544 globalColsToZero->doExport(*myColsToZero,*importer,Xpetra::ADD);
545 // import to column map
546 myColsToZero->doImport(*globalColsToZero,*importer,Xpetra::INSERT);
547 }
548 else
549 globalColsToZero = myColsToZero;
550 FindNonZeros<Scalar,LocalOrdinal,GlobalOrdinal,Node>(globalColsToZero->getDeviceLocalView(Xpetra::Access::ReadOnly),dirichletDomain);
551 FindNonZeros<Scalar,LocalOrdinal,GlobalOrdinal,Node>(myColsToZero->getDeviceLocalView(Xpetra::Access::ReadOnly),dirichletCols);
552 }
553
554
555 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
556 void
557 Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
558 DetectDirichletColsAndDomains(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
559 const Kokkos::View<bool*, typename Node::device_type> & dirichletRows,
560 Kokkos::View<bool*, typename Node::device_type> dirichletCols,
561 Kokkos::View<bool*, typename Node::device_type> dirichletDomain) {
562 MueLu::DetectDirichletColsAndDomains<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A,dirichletRows,dirichletCols,dirichletDomain);
563 }
564
565 template <class Node>
566 void
567 Utilities_kokkos<double,int,int,Node>::
568 DetectDirichletColsAndDomains(const Xpetra::Matrix<double,int,int,Node>& A,
569 const Kokkos::View<bool*, typename Node::device_type> & dirichletRows,
570 Kokkos::View<bool*, typename Node::device_type> dirichletCols,
571 Kokkos::View<bool*, typename Node::device_type> dirichletDomain) {
572 MueLu::DetectDirichletColsAndDomains<double,int,int,Node>(A,dirichletRows,dirichletCols,dirichletDomain);
573 }
574
575
576 // Zeros out rows
577 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
578 void
579 ZeroDirichletRows(RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& A,
580 const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows,
581 Scalar replaceWith) {
582 using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
583
584 auto localMatrix = A->getLocalMatrixDevice();
585 LocalOrdinal numRows = A->getNodeNumRows();
586
587 Kokkos::parallel_for("MueLu:Utils::ZeroDirichletRows", range_type(0,numRows),
588 KOKKOS_LAMBDA(const LocalOrdinal row) {
589 if (dirichletRows(row)) {
590 auto rowView = localMatrix.row(row);
591 auto length = rowView.length;
592 for (decltype(length) colID = 0; colID < length; colID++)
593 rowView.value(colID) = replaceWith;
594 }
595 });
596 }
597
598 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
599 void
600 Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
601 ZeroDirichletRows(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
602 const Kokkos::View<const bool*, typename NO::device_type>& dirichletRows,
603 Scalar replaceWith) {
604 MueLu::ZeroDirichletRows<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, dirichletRows, replaceWith);
605 }
606
607 template <class Node>
608 void
609 Utilities_kokkos<double,int,int,Node>::
610 ZeroDirichletRows(RCP<Xpetra::Matrix<double, int, int, Node> >& A,
611 const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows,
612 double replaceWith) {
613 return MueLu::ZeroDirichletRows<double,int,int,Node>(A, dirichletRows, replaceWith);
614 }
615
616
617 // Zeros out rows
618 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
619 void
620 ZeroDirichletRows(RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& X,
621 const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows,
622 Scalar replaceWith) {
623 using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
624 auto myCols = X->getDeviceLocalView(Xpetra::Access::ReadWrite);
625 size_t numVecs = X->getNumVectors();
626 Kokkos::parallel_for("MueLu:Utils::ZeroDirichletRows_MV", range_type(0,dirichletRows.size()),
627 KOKKOS_LAMBDA(const size_t i) {
628 if (dirichletRows(i)) {
629 for(size_t j=0; j<numVecs; j++)
630 myCols(i,j) = replaceWith;
631 }
632 });
633 }
634
635 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
636 void
637 Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
638 ZeroDirichletRows(RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& X,
639 const Kokkos::View<const bool*, typename NO::device_type>& dirichletRows,
640 Scalar replaceWith) {
641 MueLu::ZeroDirichletRows<Scalar,LocalOrdinal,GlobalOrdinal,Node>(X, dirichletRows, replaceWith);
642 }
643
644 template <class Node>
645 void
646 Utilities_kokkos<double,int,int,Node>::
647 ZeroDirichletRows(RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, int, int, Node> >& X,
648 const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows,
649 double replaceWith) {
650 return MueLu::ZeroDirichletRows<double,int,int,Node>(X, dirichletRows, replaceWith);
651 }
652
653
654 // Zeros out columns
655 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
656 void
657 ZeroDirichletCols(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
658 const Kokkos::View<const bool*, typename Node::device_type>& dirichletCols,
659 Scalar replaceWith) {
660 using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
661
662 auto localMatrix = A->getLocalMatrixDevice();
663 LocalOrdinal numRows = A->getNodeNumRows();
664
665 Kokkos::parallel_for("MueLu:Utils::ZeroDirichletCols", range_type(0,numRows),
666 KOKKOS_LAMBDA(const LocalOrdinal row) {
667 auto rowView = localMatrix.row(row);
668 auto length = rowView.length;
669 for (decltype(length) colID = 0; colID < length; colID++)
670 if (dirichletCols(rowView.colidx(colID))) {
671 rowView.value(colID) = replaceWith;
672 }
673 });
674 }
675
676 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
677 void
678 Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
679 ZeroDirichletCols(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
680 const Kokkos::View<const bool*, typename NO::device_type>& dirichletCols,
681 Scalar replaceWith) {
682 MueLu::ZeroDirichletCols<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, dirichletCols, replaceWith);
683 }
684
685 template <class Node>
686 void
687 Utilities_kokkos<double,int,int,Node>::
688 ZeroDirichletCols(RCP<Xpetra::Matrix<double,int,int,Node> >& A,
689 const Kokkos::View<const bool*, typename Node::device_type>& dirichletCols,
690 double replaceWith) {
691 return MueLu::ZeroDirichletCols<double,int,int,Node>(A, dirichletCols, replaceWith);
692 }
693
694 // Applies rowsum criterion
695 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
696 void ApplyRowSumCriterion(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
697 const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
698 Kokkos::View<bool*, typename Node::device_type> & dirichletRows)
699 {
700 typedef Teuchos::ScalarTraits<Scalar> STS;
701 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>> rowmap = A.getRowMap();
702 for (LocalOrdinal row = 0; row < Teuchos::as<LocalOrdinal>(rowmap->getNodeNumElements()); ++row) {
703 size_t nnz = A.getNumEntriesInLocalRow(row);
704 ArrayView<const LocalOrdinal> indices;
705 ArrayView<const Scalar> vals;
706 A.getLocalRowView(row, indices, vals);
707
708 Scalar rowsum = STS::zero();
709 Scalar diagval = STS::zero();
710 for (LocalOrdinal colID = 0; colID < Teuchos::as<LocalOrdinal>(nnz); colID++) {
711 LocalOrdinal col = indices[colID];
712 if (row == col)
713 diagval = vals[colID];
714 rowsum += vals[colID];
715 }
716 if (STS::real(rowsum) > STS::magnitude(diagval) * rowSumTol)
717 dirichletRows(row) = true;
718 }
719 }
720
721 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
722 void
723 Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
724 ApplyRowSumCriterion(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
725 const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
726 Kokkos::View<bool*, typename Node::device_type> & dirichletRows)
727 {
728 MueLu::ApplyRowSumCriterion<Scalar, LocalOrdinal, GlobalOrdinal, Node>(A,rowSumTol,dirichletRows);
729 }
730
731
732 template <class Node>
733 void
734 Utilities_kokkos<double,int,int,Node>::
735 ApplyRowSumCriterion(const Xpetra::Matrix<double,int,int,Node>& A,
736 const typename Teuchos::ScalarTraits<double>::magnitudeType rowSumTol,
737 Kokkos::View<bool*, typename Node::device_type> & dirichletRows)
738 {
739 MueLu::ApplyRowSumCriterion<double, int, int, Node>(A,rowSumTol,dirichletRows);
740 }
741
742
743 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
744 RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
745 Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
746 RealValuedToScalarMultiVector(RCP<RealValuedMultiVector > X) {
747 RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Xscalar;
748#if defined(HAVE_XPETRA_TPETRA) && (defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) || defined(HAVE_TPETRA_INST_COMPLEX_FLOAT))
749 using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
750
751 // Need to cast the real-valued multivector to Scalar=complex
752 if ((typeid(Scalar).name() == typeid(std::complex<double>).name()) ||
753 (typeid(Scalar).name() == typeid(std::complex<float>).name())) {
754 size_t numVecs = X->getNumVectors();
755 Xscalar = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(X->getMap(),numVecs);
756 auto XVec = X->getDeviceLocalView();
757 auto XVecScalar = Xscalar->getDeviceLocalView();
758
759 Kokkos::parallel_for("MueLu:Utils::RealValuedToScalarMultiVector", range_type(0,X->getLocalLength()),
760 KOKKOS_LAMBDA(const size_t i) {
761 for (size_t j=0; j<numVecs; j++)
762 XVecScalar(i,j) = XVec(i,j);
763 });
764 } else
765#endif
766 Xscalar = rcp_dynamic_cast<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(X);
767 return Xscalar;
768 }
769
770 template <class Node>
771 RCP<Xpetra::MultiVector<double,int,int,Node> >
772 Utilities_kokkos<double,int,int,Node>::
773 RealValuedToScalarMultiVector(RCP<Xpetra::MultiVector<Magnitude,int,int,Node> > X) {
774 return X;
775 }
776
777
778 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
779 Teuchos::RCP<Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> > ReverseCuthillMcKee(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Op) {
780 using local_matrix_type = typename Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::local_matrix_type;
781 using local_graph_type = typename local_matrix_type::staticcrsgraph_type;
782 using lno_nnz_view_t = typename local_graph_type::entries_type::non_const_type;
783 using device = typename local_graph_type::device_type;
784 using execution_space = typename local_matrix_type::execution_space;
785 using ordinal_type = typename local_matrix_type::ordinal_type;
786
787 local_graph_type localGraph = Op.getLocalMatrixDevice().graph;
788
789 lno_nnz_view_t rcmOrder = KokkosGraph::Experimental::graph_rcm
790 <device, typename local_graph_type::row_map_type, typename local_graph_type::entries_type, lno_nnz_view_t>
791 (localGraph.row_map, localGraph.entries);
792
793 RCP<Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> > retval =
794 Xpetra::VectorFactory<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node>::Build(Op.getRowMap());
795
796 // Copy out and reorder data
797 auto view1D = Kokkos::subview(retval->getDeviceLocalView(Xpetra::Access::ReadWrite),Kokkos::ALL (), 0);
798 Kokkos::parallel_for("Utilities_kokkos::ReverseCuthillMcKee",
799 Kokkos::RangePolicy<ordinal_type, execution_space>(0, localGraph.numRows()),
800 KOKKOS_LAMBDA(const ordinal_type rowIdx) {
801 view1D(rcmOrder(rowIdx)) = rowIdx;
802 });
803 return retval;
804 }
805
806 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
807 Teuchos::RCP<Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> > CuthillMcKee(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Op) {
808 using local_matrix_type = typename Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::local_matrix_type;
809 using local_graph_type = typename local_matrix_type::staticcrsgraph_type;
810 using lno_nnz_view_t = typename local_graph_type::entries_type::non_const_type;
811 using device = typename local_graph_type::device_type;
812 using execution_space = typename local_matrix_type::execution_space;
813 using ordinal_type = typename local_matrix_type::ordinal_type;
814
815 local_graph_type localGraph = Op.getLocalMatrixDevice().graph;
816 LocalOrdinal numRows = localGraph.numRows();
817
818 lno_nnz_view_t rcmOrder = KokkosGraph::Experimental::graph_rcm
819 <device, typename local_graph_type::row_map_type, typename local_graph_type::entries_type, lno_nnz_view_t>
820 (localGraph.row_map, localGraph.entries);
821
822 RCP<Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> > retval =
823 Xpetra::VectorFactory<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node>::Build(Op.getRowMap());
824
825 // Copy out data
826 auto view1D = Kokkos::subview(retval->getDeviceLocalView(Xpetra::Access::ReadWrite),Kokkos::ALL (), 0);
827 // Since KokkosKernels produced RCM, also reverse the order of the view to get CM
828 Kokkos::parallel_for("Utilities_kokkos::ReverseCuthillMcKee",
829 Kokkos::RangePolicy<ordinal_type, execution_space>(0, numRows),
830 KOKKOS_LAMBDA(const ordinal_type rowIdx) {
831 view1D(rcmOrder(numRows - 1 - rowIdx)) = rowIdx;
832 });
833 return retval;
834 }
835
836 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
837 Teuchos::RCP<Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> >
838 Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ReverseCuthillMcKee(const Matrix &Op) {
839 return MueLu::ReverseCuthillMcKee<Scalar,LocalOrdinal,GlobalOrdinal,Node>(Op);
840 }
841
842 template <class Node>
843 Teuchos::RCP<Xpetra::Vector<int,int,int,Node> >
844 Utilities_kokkos<double,int,int,Node>::ReverseCuthillMcKee(const Matrix &Op) {
845 return MueLu::ReverseCuthillMcKee<double,int,int,Node>(Op);
846 }
847
848 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
849 Teuchos::RCP<Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> >
850 Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::CuthillMcKee(const Matrix &Op) {
851 return MueLu::CuthillMcKee<Scalar,LocalOrdinal,GlobalOrdinal,Node>(Op);
852 }
853
854 template <class Node>
855 Teuchos::RCP<Xpetra::Vector<int,int,int,Node> >
856 Utilities_kokkos<double,int,int,Node>::CuthillMcKee(const Matrix &Op) {
857 return MueLu::CuthillMcKee<double,int,int,Node>(Op);
858 }
859
860 // Applies Ones-and-Zeros to matrix rows
861 // Takes a Boolean array.
862 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
863 void
864 ApplyOAZToMatrixRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
865 const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows) {
866 TEUCHOS_ASSERT(A->isFillComplete());
867 using ATS = Kokkos::ArithTraits<Scalar>;
868 using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
869 using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
870
871 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domMap = A->getDomainMap();
872 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > ranMap = A->getRangeMap();
873 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > Rmap = A->getRowMap();
874 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > Cmap = A->getColMap();
875
876 TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == Rmap->getNodeNumElements());
877
878 const Scalar one = impl_ATS::one();
879 const Scalar zero = impl_ATS::zero();
880
881 auto localMatrix = A->getLocalMatrixDevice();
882 auto localRmap = Rmap->getLocalMap();
883 auto localCmap = Cmap->getLocalMap();
884
885 Kokkos::parallel_for("MueLu::Utils::ApplyOAZ",range_type(0,dirichletRows.extent(0)),
886 KOKKOS_LAMBDA(const LocalOrdinal row) {
887 if (dirichletRows(row)){
888 auto rowView = localMatrix.row(row);
889 auto length = rowView.length;
890 auto row_gid = localRmap.getGlobalElement(row);
891 auto row_lid = localCmap.getLocalElement(row_gid);
892
893 for (decltype(length) colID = 0; colID < length; colID++)
894 if (rowView.colidx(colID) == row_lid)
895 rowView.value(colID) = one;
896 else
897 rowView.value(colID) = zero;
898 }
899 });
900 }
901
902 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
903 void
904 Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
905 ApplyOAZToMatrixRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
906 const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows) {
907 MueLu::ApplyOAZToMatrixRows<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, dirichletRows);
908 }
909
910 template <class Node>
911 void
912 Utilities_kokkos<double,int,int,Node>::
913 ApplyOAZToMatrixRows(Teuchos::RCP<Xpetra::Matrix<double,int,int,Node> >& A,
914 const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows) {
915 MueLu::ApplyOAZToMatrixRows<double,int,int,Node>(A, dirichletRows);
916 }
917
918} //namespace MueLu
919
920#define MUELU_UTILITIES_KOKKOS_SHORT
921#endif // MUELU_UTILITIES_KOKKOS_DEF_HPP
EpetraExt_BlockMapIn.h
EpetraExt_BlockMapOut.h
EpetraExt_CrsMatrixIn.h
EpetraExt_MatrixMatrix.h
EpetraExt_MultiVectorIn.h
EpetraExt_MultiVectorOut.h
EpetraExt_RowMatrixOut.h
MueLu_ConfigDefs.hpp
LocalOrdinal
MueLu::DefaultLocalOrdinal LocalOrdinal
Definition: MueLu_UseDefaultTypes.hpp:50
Scalar
MueLu::DefaultScalar Scalar
Definition: MueLu_UseDefaultTypes.hpp:49
Node
MueLu::DefaultNode Node
Definition: MueLu_UseDefaultTypes.hpp:52
MueLu_Utilities_kokkos_decl.hpp
Kokkos::View::extent
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< std::is_integral< iType >::value, size_t >::type extent(const iType &r) const noexcept
MueLu
Namespace for MueLu classes and methods.
Definition: MueLu_BrickAggregationFactory_decl.hpp:76
MueLu::ApplyOAZToMatrixRows
void ApplyOAZToMatrixRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Kokkos::View< const bool *, typename Node::device_type > &dirichletRows)
Definition: MueLu_Utilities_kokkos_def.hpp:864
MueLu::FindNonZeros
void FindNonZeros(const typename Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node >::dual_view_type::t_dev_const_um vals, Kokkos::View< bool *, typename Node::device_type > nonzeros)
Definition: MueLu_Utilities_kokkos_def.hpp:480
MueLu::ApplyRowSumCriterion
void ApplyRowSumCriterion(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType rowSumTol, Kokkos::View< bool *, typename Node::device_type > &dirichletRows)
Definition: MueLu_Utilities_kokkos_def.hpp:696
MueLu::ZeroDirichletCols
void ZeroDirichletCols(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Kokkos::View< const bool *, typename Node::device_type > &dirichletCols, Scalar replaceWith)
Definition: MueLu_Utilities_kokkos_def.hpp:657
MueLu::DetectDirichletCols
Kokkos::View< bool *, typename NO::device_type > DetectDirichletCols(const Xpetra::Matrix< SC, LO, GO, NO > &A, const Kokkos::View< const bool *, typename NO::device_type > &dirichletRows)
Definition: MueLu_Utilities_kokkos_def.hpp:411
MueLu::CuthillMcKee
Teuchos::RCP< Xpetra::Vector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > > CuthillMcKee(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op)
Definition: MueLu_Utilities_kokkos_def.hpp:807
MueLu::DetectDirichletRows
Kokkos::View< bool *, typename NO::device_type > DetectDirichletRows(const Xpetra::Matrix< SC, LO, GO, NO > &A, const typename Teuchos::ScalarTraits< SC >::magnitudeType &tol, const bool count_twos_as_dirichlet)
Definition: MueLu_Utilities_kokkos_def.hpp:345
MueLu::ZeroDirichletRows
void ZeroDirichletRows(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &X, const Kokkos::View< const bool *, typename Node::device_type > &dirichletRows, Scalar replaceWith)
Definition: MueLu_Utilities_kokkos_def.hpp:620
MueLu::DetectDirichletColsAndDomains
void DetectDirichletColsAndDomains(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Kokkos::View< bool *, typename Node::device_type > &dirichletRows, Kokkos::View< bool *, typename Node::device_type > dirichletCols, Kokkos::View< bool *, typename Node::device_type > dirichletDomain)
Definition: MueLu_Utilities_kokkos_def.hpp:512
MueLu::ReverseCuthillMcKee
Teuchos::RCP< Xpetra::Vector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > > ReverseCuthillMcKee(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op)
Definition: MueLu_Utilities_kokkos_def.hpp:779