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 <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 <Kokkos_CrsMatrix.hpp>
63 
64 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_EPETRAEXT)
65 #include <EpetraExt_MatrixMatrix.h>
66 #include <EpetraExt_RowMatrixOut.h>
68 #include <EpetraExt_CrsMatrixIn.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 
105 
106 namespace MueLu {
107 
108 #ifdef HAVE_MUELU_EPETRA
109  //using Xpetra::EpetraCrsMatrix; // TODO: mv in Xpetra_UseShortNamesScalar
110  //using Xpetra::EpetraMultiVector;
111 #endif
112 
113  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
114  Teuchos::ArrayRCP<Scalar> Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetMatrixDiagonal(const Matrix& A) {
115  // FIXME Kokkos
116 
117  size_t numRows = A.getRowMap()->getNodeNumElements();
118  Teuchos::ArrayRCP<SC> diag(numRows);
119 
120  Teuchos::ArrayView<const LO> cols;
121  Teuchos::ArrayView<const SC> vals;
122  for (size_t i = 0; i < numRows; ++i) {
123  A.getLocalRowView(i, cols, vals);
124 
125  LO j = 0;
126  for (; j < cols.size(); ++j) {
127  if (Teuchos::as<size_t>(cols[j]) == i) {
128  diag[i] = vals[j];
129  break;
130  }
131  }
132  if (j == cols.size()) {
133  // Diagonal entry is absent
134  diag[i] = Teuchos::ScalarTraits<SC>::zero();
135  }
136  }
137 
138  return diag;
139  } //GetMatrixDiagonal
140 
141  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
142  Teuchos::RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetMatrixDiagonalInverse(const Matrix& A, Magnitude tol) {
143  // FIXME Kokkos
144  RCP<const Map> rowMap = A.getRowMap();
145  RCP<Vector> diag = VectorFactory::Build(rowMap);
146  ArrayRCP<SC> diagVals = diag->getDataNonConst(0);
147 
148  size_t numRows = rowMap->getNodeNumElements();
149 
150  Teuchos::ArrayView<const LO> cols;
151  Teuchos::ArrayView<const SC> vals;
152  for (size_t i = 0; i < numRows; ++i) {
153  A.getLocalRowView(i, cols, vals);
154 
155  LO j = 0;
156  for (; j < cols.size(); ++j) {
157  if (Teuchos::as<size_t>(cols[j]) == i) {
158  if(Teuchos::ScalarTraits<SC>::magnitude(vals[j]) > tol)
159  diagVals[i] = Teuchos::ScalarTraits<SC>::one() / vals[j];
160  else
161  diagVals[i]=Teuchos::ScalarTraits<SC>::zero();
162  break;
163  }
164  }
165  if (j == cols.size()) {
166  // Diagonal entry is absent
167  diagVals[i]=Teuchos::ScalarTraits<SC>::zero();
168  }
169  }
170  diagVals=null;
171 
172  return diag;
173  } //GetMatrixDiagonalInverse
174 
175  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
176  Teuchos::ArrayRCP<Scalar> Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetLumpedMatrixDiagonal(const Matrix &A) {
177  // FIXME: Kokkos
178  size_t numRows = A.getRowMap()->getNodeNumElements();
179  Teuchos::ArrayRCP<SC> diag(numRows);
180 
181  Teuchos::ArrayView<const LO> cols;
182  Teuchos::ArrayView<const SC> vals;
183  for (size_t i = 0; i < numRows; ++i) {
184  A.getLocalRowView(i, cols, vals);
185 
186  diag[i] = Teuchos::ScalarTraits<Scalar>::zero();
187  for (LO j = 0; j < cols.size(); ++j) {
188  diag[i] += Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]);
189  }
190  }
191 
192  return diag;
193  } //GetMatrixDiagonal
194 
195  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
196  RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetMatrixOverlappedDiagonal(const Matrix& A) {
197  // FIXME: Kokkos
198  RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
199  RCP<Vector> localDiag = VectorFactory::Build(rowMap);
200 
201  try {
202  const CrsMatrixWrap* crsOp = dynamic_cast<const CrsMatrixWrap*>(&A);
203  if (crsOp == NULL) {
204  throw Exceptions::RuntimeError("cast to CrsMatrixWrap failed");
205  }
206  Teuchos::ArrayRCP<size_t> offsets;
207  crsOp->getLocalDiagOffsets(offsets);
208  crsOp->getLocalDiagCopy(*localDiag,offsets());
209  }
210  catch (...) {
211  ArrayRCP<SC> localDiagVals = localDiag->getDataNonConst(0);
212  Teuchos::ArrayRCP<SC> diagVals = GetMatrixDiagonal(A);
213  for (LO i = 0; i < localDiagVals.size(); i++)
214  localDiagVals[i] = diagVals[i];
215  localDiagVals = diagVals = null;
216  }
217 
218  RCP<Vector> diagonal = VectorFactory::Build(colMap);
219  RCP< const Import> importer;
220  importer = A.getCrsGraph()->getImporter();
221  if (importer == Teuchos::null) {
222  importer = ImportFactory::Build(rowMap, colMap);
223  }
224  diagonal->doImport(*localDiag, *(importer), Xpetra::INSERT);
225 
226  return diagonal;
227  } //GetMatrixOverlappedDiagonal
228 
229  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
230  void Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix(Matrix& Op, const Teuchos::ArrayRCP<const SC>& scalingVector, bool doInverse,
231  bool doFillComplete,
232  bool doOptimizeStorage)
233  {
234  SC one = Teuchos::ScalarTraits<SC>::one();
235  Teuchos::ArrayRCP<SC> sv(scalingVector.size());
236  if (doInverse) {
237  for (int i = 0; i < scalingVector.size(); ++i)
238  sv[i] = one / scalingVector[i];
239  } else {
240  for (int i = 0; i < scalingVector.size(); ++i)
241  sv[i] = scalingVector[i];
242  }
243 
244  switch (Op.getRowMap()->lib()) {
245  case Xpetra::UseTpetra:
246  MyOldScaleMatrix_Tpetra(Op, sv, doFillComplete, doOptimizeStorage);
247  break;
248 
249  case Xpetra::UseEpetra:
250  // FIXME?
251  // Utils2_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix_Epetra(Op, sv, doFillComplete, doOptimizeStorage);
252  throw std::runtime_error("FIXME");
253  break;
254 
255  default:
256  throw Exceptions::RuntimeError("Only Epetra and Tpetra matrices can be scaled.");
257  break;
258  }
259  }
260 
261  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
262  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) {
263  throw Exceptions::RuntimeError("MyOldScaleMatrix_Epetra: Epetra needs SC=double and LO=GO=int.");
264  }
265 
266  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
267  void Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix_Tpetra(Matrix& Op, const Teuchos::ArrayRCP<SC>& scalingVector,
268  bool doFillComplete,
269  bool doOptimizeStorage)
270  {
271 #ifdef HAVE_MUELU_TPETRA
272  try {
273  Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tpOp = Op2NonConstTpetraCrs(Op);
274 
275  const RCP<const Tpetra::Map<LO,GO,NO> > rowMap = tpOp.getRowMap();
276  const RCP<const Tpetra::Map<LO,GO,NO> > domainMap = tpOp.getDomainMap();
277  const RCP<const Tpetra::Map<LO,GO,NO> > rangeMap = tpOp.getRangeMap();
278 
279  size_t maxRowSize = tpOp.getNodeMaxNumRowEntries();
280  if (maxRowSize == Teuchos::as<size_t>(-1)) // hasn't been determined yet
281  maxRowSize = 20;
282 
283  std::vector<SC> scaledVals(maxRowSize);
284  if (tpOp.isFillComplete())
285  tpOp.resumeFill();
286 
287  if (Op.isLocallyIndexed() == true) {
288  Teuchos::ArrayView<const LO> cols;
289  Teuchos::ArrayView<const SC> vals;
290 
291  for (size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
292  tpOp.getLocalRowView(i, cols, vals);
293  size_t nnz = tpOp.getNumEntriesInLocalRow(i);
294  if (nnz > maxRowSize) {
295  maxRowSize = nnz;
296  scaledVals.resize(maxRowSize);
297  }
298  for (size_t j = 0; j < nnz; ++j)
299  scaledVals[j] = vals[j]*scalingVector[i];
300 
301  if (nnz > 0) {
302  Teuchos::ArrayView<const SC> valview(&scaledVals[0], nnz);
303  tpOp.replaceLocalValues(i, cols, valview);
304  }
305  } //for (size_t i=0; ...
306 
307  } else {
308  Teuchos::ArrayView<const GO> cols;
309  Teuchos::ArrayView<const SC> vals;
310 
311  for (size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
312  GO gid = rowMap->getGlobalElement(i);
313  tpOp.getGlobalRowView(gid, cols, vals);
314  size_t nnz = tpOp.getNumEntriesInGlobalRow(gid);
315  if (nnz > maxRowSize) {
316  maxRowSize = nnz;
317  scaledVals.resize(maxRowSize);
318  }
319  // FIXME FIXME FIXME FIXME FIXME FIXME
320  for (size_t j = 0; j < nnz; ++j)
321  scaledVals[j] = vals[j]*scalingVector[i]; //FIXME i or gid?
322 
323  if (nnz > 0) {
324  Teuchos::ArrayView<const SC> valview(&scaledVals[0], nnz);
325  tpOp.replaceGlobalValues(gid, cols, valview);
326  }
327  } //for (size_t i=0; ...
328  }
329 
330  if (doFillComplete) {
331  if (domainMap == Teuchos::null || rangeMap == Teuchos::null)
332  throw Exceptions::RuntimeError("In Utils_kokkos::Scaling: cannot fillComplete because the domain and/or range map hasn't been defined");
333 
334  RCP<Teuchos::ParameterList> params = rcp(new Teuchos::ParameterList());
335  params->set("Optimize Storage", doOptimizeStorage);
336  params->set("No Nonlocal Changes", true);
337  Op.fillComplete(Op.getDomainMap(), Op.getRangeMap(), params);
338  }
339  } catch(...) {
340  throw Exceptions::RuntimeError("Only Tpetra::CrsMatrix types can be scaled (Err.1)");
341  }
342 #else
343  throw Exceptions::RuntimeError("Matrix scaling is not possible because Tpetra has not been enabled.");
344 #endif
345  } //MyOldScaleMatrix_Tpetra()
346 
347  template <class SC, class LO, class GO, class NO>
348  Kokkos::View<const bool*, typename NO::device_type> DetectDirichletRows(const Xpetra::Matrix<SC,LO,GO,NO>& A, const typename Teuchos::ScalarTraits<SC>::magnitudeType& tol) {
349  typedef Kokkos::ArithTraits<SC> ATS;
350 
351  auto localMatrix = A.getLocalMatrix();
352  LO numRows = A.getNodeNumRows();
353 
354  Kokkos::View<bool*, typename NO::device_type> boundaryNodes("boundaryNodes", numRows);
355  Kokkos::parallel_for("Utils::DetectDirichletRows", numRows, KOKKOS_LAMBDA(const LO row) {
356  auto rowView = localMatrix.template row<LO>(row);
357  auto length = rowView.length;
358 
359  boundaryNodes(row) = true;
360  for (decltype(length) colID = 0; colID < length; colID++)
361  if ((rowView.colidx(colID) != row) && (ATS::magnitude(rowView.value(colID)) > tol)) {
362  boundaryNodes(row) = false;
363  break;
364  }
365  });
366 
367  return boundaryNodes;
368  }
369 
370  template <class SC, class LO, class GO, class NO>
371  Kokkos::View<const bool*, typename NO::device_type> Utilities_kokkos<SC, LO, GO, NO>::DetectDirichletRows(const Matrix& A, const typename Teuchos::ScalarTraits<SC>::magnitudeType& tol) {
372  return MueLu::DetectDirichletRows<SC,LO,GO,NO>(A, tol);
373  }
374 
375  template <class Node>
376  Kokkos::View<const bool*, typename Node::device_type> Utilities_kokkos<double,int,int,Node>::DetectDirichletRows(const Xpetra::Matrix<double,int,int,Node>& A, const typename Teuchos::ScalarTraits<double>::magnitudeType& tol) {
377  return MueLu::DetectDirichletRows<double,int,int,Node>(A, tol);
378  }
379 
380 } //namespace MueLu
381 
382 #define MUELU_UTILITIES_KOKKOS_SHORT
383 #endif // MUELU_UTILITIES_KOKKOS_DEF_HPP
Namespace for MueLu classes and methods.
Kokkos::View< const bool *, typename NO::device_type > DetectDirichletRows(const Xpetra::Matrix< SC, LO, GO, NO > &A, const typename Teuchos::ScalarTraits< SC >::magnitudeType &tol)
void parallel_for(const ExecPolicy &policy, const FunctorType &functor, const std::string &str="", typename Impl::enable_if< !Impl::is_integral< ExecPolicy >::value >::type *=0)