MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 <Teuchos_DefaultComm.hpp>
50 #include <Teuchos_ParameterList.hpp>
51 
52 #include "MueLu_ConfigDefs.hpp"
53 
54 #ifdef HAVE_MUELU_EPETRA
55 # ifdef HAVE_MPI
56 # include "Epetra_MpiComm.h"
57 # endif
58 #endif
59 
60 #include <Kokkos_ArithTraits.hpp>
61 #include <Kokkos_Core.hpp>
62 #include <KokkosSparse_CrsMatrix.hpp>
63 
64 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_EPETRAEXT)
65 #include <EpetraExt_MatrixMatrix.h>
66 #include <EpetraExt_RowMatrixOut.h>
67 #include <EpetraExt_MultiVectorOut.h>
68 #include <EpetraExt_CrsMatrixIn.h>
69 #include <EpetraExt_MultiVectorIn.h>
70 #include <EpetraExt_BlockMapIn.h>
71 #include <Xpetra_EpetraUtils.hpp>
72 #include <Xpetra_EpetraMultiVector.hpp>
73 #include <EpetraExt_BlockMapOut.h>
74 #endif
75 
76 #ifdef HAVE_MUELU_TPETRA
77 #include <MatrixMarket_Tpetra.hpp>
78 #include <Tpetra_RowMatrixTransposer.hpp>
79 #include <TpetraExt_MatrixMatrix.hpp>
80 #include <Xpetra_TpetraMultiVector.hpp>
81 #include <Xpetra_TpetraCrsMatrix.hpp>
82 #include <Xpetra_TpetraBlockCrsMatrix.hpp>
83 #endif
84 
85 #ifdef HAVE_MUELU_EPETRA
86 #include <Xpetra_EpetraMap.hpp>
87 #endif
88 
89 #include <Xpetra_BlockedCrsMatrix.hpp>
90 #include <Xpetra_DefaultPlatform.hpp>
91 #include <Xpetra_Import.hpp>
92 #include <Xpetra_ImportFactory.hpp>
93 #include <Xpetra_Map.hpp>
94 #include <Xpetra_MapFactory.hpp>
95 #include <Xpetra_Matrix.hpp>
96 #include <Xpetra_MatrixMatrix.hpp>
97 #include <Xpetra_MatrixFactory.hpp>
98 #include <Xpetra_MultiVector.hpp>
99 #include <Xpetra_MultiVectorFactory.hpp>
100 #include <Xpetra_Operator.hpp>
101 #include <Xpetra_Vector.hpp>
102 #include <Xpetra_VectorFactory.hpp>
103 
104 #include <MueLu_Utilities_kokkos_decl.hpp>
105 
106 namespace MueLu {
107 
108  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
109  Teuchos::ArrayRCP<Scalar> Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetMatrixDiagonal(const Matrix& A) {
110  // FIXME_KOKKOS
111 
112  size_t numRows = A.getRowMap()->getNodeNumElements();
113  Teuchos::ArrayRCP<SC> diag(numRows);
114 
115  Teuchos::ArrayView<const LO> cols;
116  Teuchos::ArrayView<const SC> vals;
117  for (size_t i = 0; i < numRows; ++i) {
118  A.getLocalRowView(i, cols, vals);
119 
120  LO j = 0;
121  for (; j < cols.size(); ++j) {
122  if (Teuchos::as<size_t>(cols[j]) == i) {
123  diag[i] = vals[j];
124  break;
125  }
126  }
127  if (j == cols.size()) {
128  // Diagonal entry is absent
129  diag[i] = Teuchos::ScalarTraits<SC>::zero();
130  }
131  }
132 
133  return diag;
134  } //GetMatrixDiagonal
135 
136  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
137  Teuchos::RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetMatrixDiagonalInverse(const Matrix& A, Magnitude tol) {
138  // FIXME_KOKKOS
139  RCP<const Map> rowMap = A.getRowMap();
140  RCP<Vector> diag = VectorFactory::Build(rowMap);
141  ArrayRCP<SC> diagVals = diag->getDataNonConst(0);
142 
143  size_t numRows = rowMap->getNodeNumElements();
144 
145  Teuchos::ArrayView<const LO> cols;
146  Teuchos::ArrayView<const SC> vals;
147  for (size_t i = 0; i < numRows; ++i) {
148  A.getLocalRowView(i, cols, vals);
149 
150  LO j = 0;
151  for (; j < cols.size(); ++j) {
152  if (Teuchos::as<size_t>(cols[j]) == i) {
153  if(Teuchos::ScalarTraits<SC>::magnitude(vals[j]) > tol)
154  diagVals[i] = Teuchos::ScalarTraits<SC>::one() / vals[j];
155  else
156  diagVals[i]=Teuchos::ScalarTraits<SC>::zero();
157  break;
158  }
159  }
160  if (j == cols.size()) {
161  // Diagonal entry is absent
162  diagVals[i]=Teuchos::ScalarTraits<SC>::zero();
163  }
164  }
165  diagVals=null;
166 
167  return diag;
168  } //GetMatrixDiagonalInverse
169 
170  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
171  Teuchos::ArrayRCP<Scalar> Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetLumpedMatrixDiagonal(const Matrix &A) {
172  // FIXME_KOKKOS
173  size_t numRows = A.getRowMap()->getNodeNumElements();
174  Teuchos::ArrayRCP<SC> diag(numRows);
175 
176  Teuchos::ArrayView<const LO> cols;
177  Teuchos::ArrayView<const SC> vals;
178  for (size_t i = 0; i < numRows; ++i) {
179  A.getLocalRowView(i, cols, vals);
180 
181  diag[i] = Teuchos::ScalarTraits<Scalar>::zero();
182  for (LO j = 0; j < cols.size(); ++j) {
183  diag[i] += Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]);
184  }
185  }
186 
187  return diag;
188  } //GetMatrixDiagonal
189 
190  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
191  RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetMatrixOverlappedDiagonal(const Matrix& A) {
192  // FIXME_KOKKOS
193  RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
194  RCP<Vector> localDiag = VectorFactory::Build(rowMap);
195 
196  try {
197  const CrsMatrixWrap* crsOp = dynamic_cast<const CrsMatrixWrap*>(&A);
198  if (crsOp == NULL) {
199  throw Exceptions::RuntimeError("cast to CrsMatrixWrap failed");
200  }
201  Teuchos::ArrayRCP<size_t> offsets;
202  crsOp->getLocalDiagOffsets(offsets);
203  crsOp->getLocalDiagCopy(*localDiag,offsets());
204  }
205  catch (...) {
206  ArrayRCP<SC> localDiagVals = localDiag->getDataNonConst(0);
207  Teuchos::ArrayRCP<SC> diagVals = GetMatrixDiagonal(A);
208  for (LO i = 0; i < localDiagVals.size(); i++)
209  localDiagVals[i] = diagVals[i];
210  localDiagVals = diagVals = null;
211  }
212 
213  RCP<Vector> diagonal = VectorFactory::Build(colMap);
214  RCP< const Import> importer;
215  importer = A.getCrsGraph()->getImporter();
216  if (importer == Teuchos::null) {
217  importer = ImportFactory::Build(rowMap, colMap);
218  }
219  diagonal->doImport(*localDiag, *(importer), Xpetra::INSERT);
220 
221  return diagonal;
222  } //GetMatrixOverlappedDiagonal
223 
224  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
225  void Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix(Matrix& Op, const Teuchos::ArrayRCP<const SC>& scalingVector, bool doInverse,
226  bool doFillComplete,
227  bool doOptimizeStorage)
228  {
229  SC one = Teuchos::ScalarTraits<SC>::one();
230  Teuchos::ArrayRCP<SC> sv(scalingVector.size());
231  if (doInverse) {
232  for (int i = 0; i < scalingVector.size(); ++i)
233  sv[i] = one / scalingVector[i];
234  } else {
235  for (int i = 0; i < scalingVector.size(); ++i)
236  sv[i] = scalingVector[i];
237  }
238 
239  switch (Op.getRowMap()->lib()) {
240  case Xpetra::UseTpetra:
241  MyOldScaleMatrix_Tpetra(Op, sv, doFillComplete, doOptimizeStorage);
242  break;
243 
244  case Xpetra::UseEpetra:
245  // FIXME_KOKKOS
246  // Utils2_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix_Epetra(Op, sv, doFillComplete, doOptimizeStorage);
247  throw std::runtime_error("FIXME");
248 #ifndef __NVCC__ //prevent nvcc warning
249  break;
250 #endif
251 
252  default:
253  throw Exceptions::RuntimeError("Only Epetra and Tpetra matrices can be scaled.");
254 #ifndef __NVCC__ //prevent nvcc warning
255  break;
256 #endif
257  }
258  }
259 
260  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
261  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 */) {
262  throw Exceptions::RuntimeError("MyOldScaleMatrix_Epetra: Epetra needs SC=double and LO=GO=int.");
263  }
264 
265  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
266  void Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix_Tpetra(Matrix& Op, const Teuchos::ArrayRCP<SC>& scalingVector,
267  bool doFillComplete,
268  bool doOptimizeStorage)
269  {
270 #ifdef HAVE_MUELU_TPETRA
271  try {
272  Tpetra::CrsMatrix<SC,LO,GO,NO>& tpOp = Op2NonConstTpetraCrs(Op);
273 
274  const RCP<const Tpetra::Map<LO,GO,NO> > rowMap = tpOp.getRowMap();
275  const RCP<const Tpetra::Map<LO,GO,NO> > domainMap = tpOp.getDomainMap();
276  const RCP<const Tpetra::Map<LO,GO,NO> > rangeMap = tpOp.getRangeMap();
277 
278  size_t maxRowSize = tpOp.getNodeMaxNumRowEntries();
279  if (maxRowSize == Teuchos::as<size_t>(-1)) // hasn't been determined yet
280  maxRowSize = 20;
281 
282  std::vector<SC> scaledVals(maxRowSize);
283  if (tpOp.isFillComplete())
284  tpOp.resumeFill();
285 
286  if (Op.isLocallyIndexed() == true) {
287  Teuchos::ArrayView<const LO> cols;
288  Teuchos::ArrayView<const SC> vals;
289 
290  for (size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
291  tpOp.getLocalRowView(i, cols, vals);
292  size_t nnz = tpOp.getNumEntriesInLocalRow(i);
293  if (nnz > maxRowSize) {
294  maxRowSize = nnz;
295  scaledVals.resize(maxRowSize);
296  }
297  for (size_t j = 0; j < nnz; ++j)
298  scaledVals[j] = vals[j]*scalingVector[i];
299 
300  if (nnz > 0) {
301  Teuchos::ArrayView<const SC> valview(&scaledVals[0], nnz);
302  tpOp.replaceLocalValues(i, cols, valview);
303  }
304  } //for (size_t i=0; ...
305 
306  } else {
307  Teuchos::ArrayView<const GO> cols;
308  Teuchos::ArrayView<const SC> vals;
309 
310  for (size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
311  GO gid = rowMap->getGlobalElement(i);
312  tpOp.getGlobalRowView(gid, cols, vals);
313  size_t nnz = tpOp.getNumEntriesInGlobalRow(gid);
314  if (nnz > maxRowSize) {
315  maxRowSize = nnz;
316  scaledVals.resize(maxRowSize);
317  }
318  // FIXME FIXME FIXME FIXME FIXME FIXME
319  for (size_t j = 0; j < nnz; ++j)
320  scaledVals[j] = vals[j]*scalingVector[i]; //FIXME i or gid?
321 
322  if (nnz > 0) {
323  Teuchos::ArrayView<const SC> valview(&scaledVals[0], nnz);
324  tpOp.replaceGlobalValues(gid, cols, valview);
325  }
326  } //for (size_t i=0; ...
327  }
328 
329  if (doFillComplete) {
330  if (domainMap == Teuchos::null || rangeMap == Teuchos::null)
331  throw Exceptions::RuntimeError("In Utils_kokkos::Scaling: cannot fillComplete because the domain and/or range map hasn't been defined");
332 
333  RCP<Teuchos::ParameterList> params = rcp(new Teuchos::ParameterList());
334  params->set("Optimize Storage", doOptimizeStorage);
335  params->set("No Nonlocal Changes", true);
336  Op.fillComplete(Op.getDomainMap(), Op.getRangeMap(), params);
337  }
338  } catch(...) {
339  throw Exceptions::RuntimeError("Only Tpetra::CrsMatrix types can be scaled (Err.1)");
340  }
341 #else
342  throw Exceptions::RuntimeError("Matrix scaling is not possible because Tpetra has not been enabled.");
343 #endif
344  } //MyOldScaleMatrix_Tpetra()
345 
346 
347  template <class SC, class LO, class GO, class NO>
348  Kokkos::View<const bool*, typename NO::device_type>
349  DetectDirichletRows(const Xpetra::Matrix<SC,LO,GO,NO>& A,
350  const typename Teuchos::ScalarTraits<SC>::magnitudeType& tol,
351  const bool count_twos_as_dirichlet) {
352  using ATS = Kokkos::ArithTraits<SC>;
353  using range_type = Kokkos::RangePolicy<LO, typename NO::execution_space>;
354 
355  auto localMatrix = A.getLocalMatrix();
356  LO numRows = A.getNodeNumRows();
357 
358  Kokkos::View<bool*, typename NO::device_type> boundaryNodes(Kokkos::ViewAllocateWithoutInitializing("boundaryNodes"), numRows);
359  if (count_twos_as_dirichlet)
360  Kokkos::parallel_for("MueLu:Utils::DetectDirichletRows_Twos_As_Dirichlet", range_type(0,numRows),
361  KOKKOS_LAMBDA(const LO row) {
362  auto rowView = localMatrix.row(row);
363  auto length = rowView.length;
364 
365  boundaryNodes(row) = true;
366  if (length > 2) {
367  decltype(length) colID;
368  for (colID = 0; colID < length; colID++)
369  if ((rowView.colidx(colID) != row) &&
370  (ATS::magnitude(rowView.value(colID)) > tol)) {
371  if (!boundaryNodes(row))
372  break;
373  boundaryNodes(row) = false;
374  }
375  if (colID == length)
376  boundaryNodes(row) = true;
377  }
378  });
379  else
380  Kokkos::parallel_for("MueLu:Utils::DetectDirichletRows", range_type(0,numRows),
381  KOKKOS_LAMBDA(const LO row) {
382  auto rowView = localMatrix.row(row);
383  auto length = rowView.length;
384 
385  boundaryNodes(row) = true;
386  for (decltype(length) colID = 0; colID < length; colID++)
387  if ((rowView.colidx(colID) != row) &&
388  (ATS::magnitude(rowView.value(colID)) > tol)) {
389  boundaryNodes(row) = false;
390  break;
391  }
392  });
393 
394  return boundaryNodes;
395  }
396 
397  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
398  Kokkos::View<const bool*, typename Node::device_type>
399  Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
400  DetectDirichletRows(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, const typename Teuchos::ScalarTraits<Scalar>::magnitudeType& tol, const bool count_twos_as_dirichlet) {
401  return MueLu::DetectDirichletRows<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, tol, count_twos_as_dirichlet);
402  }
403 
404  template <class Node>
405  Kokkos::View<const bool*, typename Node::device_type>
406  Utilities_kokkos<double,int,int,Node>::
407  DetectDirichletRows(const Xpetra::Matrix<double,int,int,Node>& A, const typename Teuchos::ScalarTraits<double>::magnitudeType& tol, const bool count_twos_as_dirichlet) {
408  return MueLu::DetectDirichletRows<double,int,int,Node>(A, tol,count_twos_as_dirichlet);
409  }
410 
411 
412  template <class SC, class LO, class GO, class NO>
413  Kokkos::View<const bool*, typename NO::device_type>
414  DetectDirichletCols(const Xpetra::Matrix<SC,LO,GO,NO>& A,
415  const Kokkos::View<const bool*, typename NO::device_type>& dirichletRows) {
416  using ATS = Kokkos::ArithTraits<SC>;
417  using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
418  using range_type = Kokkos::RangePolicy<LO, typename NO::execution_space>;
419 
420  SC zero = ATS::zero();
421  SC one = ATS::one();
422 
423  auto localMatrix = A.getLocalMatrix();
424  LO numRows = A.getNodeNumRows();
425 
426  Teuchos::RCP<const Xpetra::Map<LO,GO,NO> > domMap = A.getDomainMap();
427  Teuchos::RCP<const Xpetra::Map<LO,GO,NO> > colMap = A.getColMap();
428  Teuchos::RCP<Xpetra::MultiVector<SC,LO,GO,NO> > myColsToZero = Xpetra::MultiVectorFactory<SC,LO,GO,NO>::Build(colMap,1);
429  myColsToZero->putScalar(zero);
430  auto myColsToZeroView = myColsToZero->template getLocalView<typename NO::device_type>();
431  // Find all local column indices that are in Dirichlet rows, record in myColsToZero as 1.0
432  Kokkos::parallel_for("MueLu:Utils::DetectDirichletCols1", range_type(0,numRows),
433  KOKKOS_LAMBDA(const LO row) {
434  if (dirichletRows(row)) {
435  auto rowView = localMatrix.row(row);
436  auto length = rowView.length;
437 
438  for (decltype(length) colID = 0; colID < length; colID++)
439  myColsToZeroView(rowView.colidx(colID),0) = one;
440  }
441  });
442 
443  Teuchos::RCP<Xpetra::MultiVector<SC,LO,GO,NO> > globalColsToZero = Xpetra::MultiVectorFactory<SC,LO,GO,NO>::Build(domMap,1);
444  globalColsToZero->putScalar(zero);
445  Teuchos::RCP<Xpetra::Export<LO,GO,NO> > exporter = Xpetra::ExportFactory<LO,GO,NO>::Build(colMap,domMap);
446  // export to domain map
447  globalColsToZero->doExport(*myColsToZero,*exporter,Xpetra::ADD);
448  // import to column map
449  myColsToZero->doImport(*globalColsToZero,*exporter,Xpetra::INSERT);
450 
451  auto myCols = myColsToZero->template getLocalView<typename NO::device_type>();
452  size_t numColEntries = colMap->getNodeNumElements();
453  Kokkos::View<bool*, typename NO::device_type> dirichletCols(Kokkos::ViewAllocateWithoutInitializing("dirichletCols"), numColEntries);
454  const typename ATS::magnitudeType eps = 2.0*ATS::eps();
455 
456  Kokkos::parallel_for("MueLu:Utils::DetectDirichletCols2", range_type(0,numColEntries),
457  KOKKOS_LAMBDA (const size_t i) {
458  dirichletCols(i) = impl_ATS::magnitude(myCols(i,0))>eps;
459  });
460  return dirichletCols;
461  }
462 
463 
464  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
465  Kokkos::View<const bool*, typename Node::device_type>
466  Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
467  DetectDirichletCols(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
468  const Kokkos::View<const bool*, typename NO::device_type>& dirichletRows) {
469  return MueLu::DetectDirichletCols<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, dirichletRows);
470  }
471 
472  template <class Node>
473  Kokkos::View<const bool*, typename Node::device_type>
474  Utilities_kokkos<double,int,int,Node>::
475  DetectDirichletCols(const Xpetra::Matrix<double,int,int,Node>& A,
476  const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows) {
477  return MueLu::DetectDirichletCols<double,int,int,Node>(A, dirichletRows);
478  }
479 
480 
481  // Zeros out rows
482  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
483  void
484  ZeroDirichletRows(RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& A,
485  const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows,
486  Scalar replaceWith) {
487  using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
488 
489  auto localMatrix = A->getLocalMatrix();
490  LocalOrdinal numRows = A->getNodeNumRows();
491 
492  Kokkos::parallel_for("MueLu:Utils::ZeroDirichletRows", range_type(0,numRows),
493  KOKKOS_LAMBDA(const LocalOrdinal row) {
494  if (dirichletRows(row)) {
495  auto rowView = localMatrix.row(row);
496  auto length = rowView.length;
497  for (decltype(length) colID = 0; colID < length; colID++)
498  rowView.value(colID) = replaceWith;
499  }
500  });
501  }
502 
503  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
504  void
505  Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
506  ZeroDirichletRows(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
507  const Kokkos::View<const bool*, typename NO::device_type>& dirichletRows,
508  Scalar replaceWith) {
509  MueLu::ZeroDirichletRows<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, dirichletRows, replaceWith);
510  }
511 
512  template <class Node>
513  void
514  Utilities_kokkos<double,int,int,Node>::
515  ZeroDirichletRows(RCP<Xpetra::Matrix<double, int, int, Node> >& A,
516  const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows,
517  double replaceWith) {
518  return MueLu::ZeroDirichletRows<double,int,int,Node>(A, dirichletRows, replaceWith);
519  }
520 
521 
522  // Zeros out rows
523  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
524  void
525  ZeroDirichletRows(RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& X,
526  const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows,
527  Scalar replaceWith) {
528  using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
529  auto myCols = X->template getLocalView<typename Node::device_type>();
530  size_t numVecs = X->getNumVectors();
531  Kokkos::parallel_for("MueLu:Utils::ZeroDirichletRows_MV", range_type(0,dirichletRows.size()),
532  KOKKOS_LAMBDA(const size_t i) {
533  if (dirichletRows(i)) {
534  for(size_t j=0; j<numVecs; j++)
535  myCols(i,j) = replaceWith;
536  }
537  });
538  }
539 
540  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
541  void
542  Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
543  ZeroDirichletRows(RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& X,
544  const Kokkos::View<const bool*, typename NO::device_type>& dirichletRows,
545  Scalar replaceWith) {
546  MueLu::ZeroDirichletRows<Scalar,LocalOrdinal,GlobalOrdinal,Node>(X, dirichletRows, replaceWith);
547  }
548 
549  template <class Node>
550  void
551  Utilities_kokkos<double,int,int,Node>::
552  ZeroDirichletRows(RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, int, int, Node> >& X,
553  const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows,
554  double replaceWith) {
555  return MueLu::ZeroDirichletRows<double,int,int,Node>(X, dirichletRows, replaceWith);
556  }
557 
558 
559  // Zeros out columns
560  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
561  void
562  ZeroDirichletCols(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
563  const Kokkos::View<const bool*, typename Node::device_type>& dirichletCols,
564  Scalar replaceWith) {
565  using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
566 
567  auto localMatrix = A->getLocalMatrix();
568  LocalOrdinal numRows = A->getNodeNumRows();
569 
570  Kokkos::parallel_for("MueLu:Utils::ZeroDirichletCols", range_type(0,numRows),
571  KOKKOS_LAMBDA(const LocalOrdinal row) {
572  auto rowView = localMatrix.row(row);
573  auto length = rowView.length;
574  for (decltype(length) colID = 0; colID < length; colID++)
575  if (dirichletCols(rowView.colidx(colID))) {
576  rowView.value(colID) = replaceWith;
577  }
578  });
579  }
580 
581  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
582  void
583  Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
584  ZeroDirichletCols(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
585  const Kokkos::View<const bool*, typename NO::device_type>& dirichletCols,
586  Scalar replaceWith) {
587  MueLu::ZeroDirichletCols<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, dirichletCols, replaceWith);
588  }
589 
590  template <class Node>
591  void
592  Utilities_kokkos<double,int,int,Node>::
593  ZeroDirichletCols(RCP<Xpetra::Matrix<double,int,int,Node> >& A,
594  const Kokkos::View<const bool*, typename Node::device_type>& dirichletCols,
595  double replaceWith) {
596  return MueLu::ZeroDirichletCols<double,int,int,Node>(A, dirichletCols, replaceWith);
597  }
598 
599 
600  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
601  RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
602  Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
603  RealValuedToScalarMultiVector(RCP<RealValuedMultiVector > X) {
604  RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Xscalar;
605 #if defined(HAVE_XPETRA_TPETRA) && (defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) || defined(HAVE_TPETRA_INST_COMPLEX_FLOAT))
606  using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
607 
608  // Need to cast the real-valued multivector to Scalar=complex
609  if ((typeid(Scalar).name() == typeid(std::complex<double>).name()) ||
610  (typeid(Scalar).name() == typeid(std::complex<float>).name())) {
611  size_t numVecs = X->getNumVectors();
612  Xscalar = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(X->getMap(),numVecs);
613  auto XVec = X->template getLocalView<typename Node::device_type>();
614  auto XVecScalar = Xscalar->template getLocalView<typename Node::device_type>();
615 
616  Kokkos::parallel_for("MueLu:Utils::RealValuedToScalarMultiVector", range_type(0,X->getLocalLength()),
617  KOKKOS_LAMBDA(const size_t i) {
618  for (size_t j=0; j<numVecs; j++)
619  XVecScalar(i,j) = XVec(i,j);
620  });
621  } else
622 #endif
623  Xscalar = rcp_dynamic_cast<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(X);
624  return Xscalar;
625  }
626 
627  template <class Node>
628  RCP<Xpetra::MultiVector<double,int,int,Node> >
629  Utilities_kokkos<double,int,int,Node>::
630  RealValuedToScalarMultiVector(RCP<Xpetra::MultiVector<Magnitude,int,int,Node> > X) {
631  return X;
632  }
633 
634 } //namespace MueLu
635 
636 #define MUELU_UTILITIES_KOKKOS_SHORT
637 #endif // MUELU_UTILITIES_KOKKOS_DEF_HPP
Xpetra::Matrix::getNodeNumRows
virtual size_t getNodeNumRows() const =0
Xpetra::ExportFactory::Build
static RCP< Export< LocalOrdinal, GlobalOrdinal, Node > > Build(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &source, const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &target)
Xpetra_Vector.hpp
LocalOrdinal
MueLu::DefaultLocalOrdinal LocalOrdinal
Definition: MueLu_UseDefaultTypes.hpp:50
Kokkos::parallel_for
void parallel_for(const ExecPolicy &policy, const FunctorType &functor, const std::string &str="", typename Impl::enable_if< Kokkos::Impl::is_execution_policy< ExecPolicy >::value >::type *=0)
EpetraExt_CrsMatrixIn.h
Xpetra::MultiVectorFactory::Build
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node >> &map, size_t NumVectors, bool zeroOut=true)
GO
GlobalOrdinal GO
Teuchos::ArrayView::size
size_type size() const
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:562
Xpetra_EpetraUtils.hpp
MueLu::ZeroDirichletRows
void ZeroDirichletRows(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Kokkos::View< const bool *, typename Node::device_type > &dirichletRows, Scalar replaceWith)
Definition: MueLu_Utilities_kokkos_def.hpp:484
LO
LocalOrdinal LO
EpetraExt_MatrixMatrix.h
Xpetra_Matrix.hpp
Teuchos::ArrayRCP::size
size_type size() const
Xpetra_MatrixMatrix.hpp
MueLu_ConfigDefs.hpp
MueLu_Utilities_kokkos_decl.hpp
EpetraExt_RowMatrixOut.h
Xpetra_MultiVector.hpp
rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos_ParameterList.hpp
Scalar
MueLu::DefaultScalar Scalar
Definition: MueLu_UseDefaultTypes.hpp:49
EpetraExt_BlockMapOut.h
EpetraExt_MultiVectorOut.h
EpetraExt_MultiVectorIn.h
EpetraExt_BlockMapIn.h
Xpetra_EpetraMap.hpp
Teuchos::ScalarTraits::magnitude
static magnitudeType magnitude(T a)
Xpetra_BlockedCrsMatrix.hpp
Xpetra_MatrixFactory.hpp
Xpetra_DefaultPlatform.hpp
Xpetra_EpetraMultiVector.hpp
SC
Scalar SC
MueLu::DetectDirichletCols
Kokkos::View< const 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:414
Xpetra_VectorFactory.hpp
Xpetra_Import.hpp
Xpetra_Operator.hpp
MueLu::DetectDirichletRows
Kokkos::View< const 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:349
Xpetra::Operator::getDomainMap
virtual Teuchos::RCP< const Map > getDomainMap() const =0
Xpetra_ImportFactory.hpp