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 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_EPETRAEXT)
61 #include <EpetraExt_MatrixMatrix.h>
62 #include <EpetraExt_RowMatrixOut.h>
63 #include <EpetraExt_MultiVectorOut.h>
64 #include <EpetraExt_CrsMatrixIn.h>
65 #include <EpetraExt_MultiVectorIn.h>
66 #include <EpetraExt_BlockMapIn.h>
67 #include <Xpetra_EpetraUtils.hpp>
68 #include <Xpetra_EpetraMultiVector.hpp>
69 #include <EpetraExt_BlockMapOut.h>
70 #endif
71 
72 #ifdef HAVE_MUELU_TPETRA
73 #include <MatrixMarket_Tpetra.hpp>
74 #include <Tpetra_RowMatrixTransposer.hpp>
75 #include <TpetraExt_MatrixMatrix.hpp>
76 #include <Xpetra_TpetraMultiVector.hpp>
77 #include <Xpetra_TpetraCrsMatrix.hpp>
78 #include <Xpetra_TpetraBlockCrsMatrix.hpp>
79 #endif
80 
81 #ifdef HAVE_MUELU_EPETRA
82 #include <Xpetra_EpetraMap.hpp>
83 #endif
84 
85 #include <Xpetra_BlockedCrsMatrix.hpp>
86 #include <Xpetra_DefaultPlatform.hpp>
87 #include <Xpetra_Import.hpp>
88 #include <Xpetra_ImportFactory.hpp>
89 #include <Xpetra_Map.hpp>
90 #include <Xpetra_MapFactory.hpp>
91 #include <Xpetra_Matrix.hpp>
92 #include <Xpetra_MatrixMatrix.hpp>
93 #include <Xpetra_MatrixFactory.hpp>
94 #include <Xpetra_MultiVector.hpp>
95 #include <Xpetra_MultiVectorFactory.hpp>
96 #include <Xpetra_Operator.hpp>
97 #include <Xpetra_Vector.hpp>
98 #include <Xpetra_VectorFactory.hpp>
99 
100 #include <MueLu_Utilities_kokkos_decl.hpp>
101 
102 namespace MueLu {
103 
104 #ifdef HAVE_MUELU_EPETRA
105  using Xpetra::EpetraCrsMatrix; // TODO: mv in Xpetra_UseShortNamesScalar
106  using Xpetra::EpetraMultiVector;
107 #endif
108 
109  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
110  Teuchos::ArrayRCP<Scalar> Utils_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetMatrixDiagonal(const Matrix& A) {
111  // FIXME Kokkos
112 
113  size_t numRows = A.getRowMap()->getNodeNumElements();
114  Teuchos::ArrayRCP<SC> diag(numRows);
115 
116  Teuchos::ArrayView<const LO> cols;
117  Teuchos::ArrayView<const SC> vals;
118  for (size_t i = 0; i < numRows; ++i) {
119  A.getLocalRowView(i, cols, vals);
120 
121  LO j = 0;
122  for (; j < cols.size(); ++j) {
123  if (Teuchos::as<size_t>(cols[j]) == i) {
124  diag[i] = vals[j];
125  break;
126  }
127  }
128  if (j == cols.size()) {
129  // Diagonal entry is absent
130  diag[i] = Teuchos::ScalarTraits<SC>::zero();
131  }
132  }
133 
134  return diag;
135  } //GetMatrixDiagonal
136 
137  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
138  Teuchos::RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Utils_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetMatrixDiagonalInverse(const Matrix& A, Magnitude tol) {
139  // FIXME Kokkos
140  RCP<const Map> rowMap = A.getRowMap();
141  RCP<Vector> diag = VectorFactory::Build(rowMap);
142  ArrayRCP<SC> diagVals = diag->getDataNonConst(0);
143 
144  size_t numRows = rowMap->getNodeNumElements();
145 
146  Teuchos::ArrayView<const LO> cols;
147  Teuchos::ArrayView<const SC> vals;
148  for (size_t i = 0; i < numRows; ++i) {
149  A.getLocalRowView(i, cols, vals);
150 
151  LO j = 0;
152  for (; j < cols.size(); ++j) {
153  if (Teuchos::as<size_t>(cols[j]) == i) {
154  if(Teuchos::ScalarTraits<SC>::magnitude(vals[j]) > tol)
155  diagVals[i] = Teuchos::ScalarTraits<SC>::one() / vals[j];
156  else
157  diagVals[i]=Teuchos::ScalarTraits<SC>::zero();
158  break;
159  }
160  }
161  if (j == cols.size()) {
162  // Diagonal entry is absent
163  diagVals[i]=Teuchos::ScalarTraits<SC>::zero();
164  }
165  }
166  diagVals=null;
167 
168  return diag;
169  } //GetMatrixDiagonalInverse
170 
171  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
172  Teuchos::ArrayRCP<Scalar> Utils_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetLumpedMatrixDiagonal(const Matrix &A) {
173  // FIXME: Kokkos
174  size_t numRows = A.getRowMap()->getNodeNumElements();
175  Teuchos::ArrayRCP<SC> diag(numRows);
176 
177  Teuchos::ArrayView<const LO> cols;
178  Teuchos::ArrayView<const SC> vals;
179  for (size_t i = 0; i < numRows; ++i) {
180  A.getLocalRowView(i, cols, vals);
181 
182  diag[i] = Teuchos::ScalarTraits<Scalar>::zero();
183  for (LO j = 0; j < cols.size(); ++j) {
184  diag[i] += Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]);
185  }
186  }
187 
188  return diag;
189  } //GetMatrixDiagonal
190 
191  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
192  RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Utils_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetMatrixOverlappedDiagonal(const Matrix& A) {
193  // FIXME: Kokkos
194  RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
195  RCP<Vector> localDiag = VectorFactory::Build(rowMap);
196 
197  try {
198  const CrsMatrixWrap* crsOp = dynamic_cast<const CrsMatrixWrap*>(&A);
199  if (crsOp == NULL) {
200  throw Exceptions::RuntimeError("cast to CrsMatrixWrap failed");
201  }
202  Teuchos::ArrayRCP<size_t> offsets;
203  crsOp->getLocalDiagOffsets(offsets);
204  crsOp->getLocalDiagCopy(*localDiag,offsets());
205  }
206  catch (...) {
207  ArrayRCP<SC> localDiagVals = localDiag->getDataNonConst(0);
208  Teuchos::ArrayRCP<SC> diagVals = GetMatrixDiagonal(A);
209  for (LO i = 0; i < localDiagVals.size(); i++)
210  localDiagVals[i] = diagVals[i];
211  localDiagVals = diagVals = null;
212  }
213 
214  RCP<Vector> diagonal = VectorFactory::Build(colMap);
215  RCP< const Import> importer;
216  importer = A.getCrsGraph()->getImporter();
217  if (importer == Teuchos::null) {
218  importer = ImportFactory::Build(rowMap, colMap);
219  }
220  diagonal->doImport(*localDiag, *(importer), Xpetra::INSERT);
221 
222  return diagonal;
223  } //GetMatrixOverlappedDiagonal
224 
225  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
226  void Utils_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ScaleMatrix(Matrix& Op, const Teuchos::ArrayRCP<SC>& scalingVector, bool doInverse) {
227  // FIXME: Kokkos
228 #ifdef HAVE_MUELU_TPETRA
229  try {
230  Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tpOp = Op2NonConstTpetraCrs(Op);
231 
232  Tpetra::Vector<SC,LO,GO,NO> x(tpOp.getRowMap(), scalingVector());
233  if (doInverse){
234  Tpetra::Vector<SC,LO,GO,NO> xi(tpOp.getRowMap());
235  xi.reciprocal(x);
236  tpOp.leftScale(xi);
237 
238  } else {
239  tpOp.leftScale(x);
240  }
241  } catch(...) {
242  throw Exceptions::RuntimeError("Matrix scaling has not been implemented Epetra");
243  }
244 #else
245  throw Exceptions::RuntimeError("Matrix scaling has not been implemented Epetra");
246 #endif // HAVE_MUELU_TPETRA
247  } //ScaleMatrix()
248 
249  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
250  void Utils_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix(Matrix& Op, const Teuchos::ArrayRCP<const SC>& scalingVector, bool doInverse,
251  bool doFillComplete,
252  bool doOptimizeStorage)
253  {
254  SC one = Teuchos::ScalarTraits<SC>::one();
255  Teuchos::ArrayRCP<SC> sv(scalingVector.size());
256  if (doInverse) {
257  for (int i = 0; i < scalingVector.size(); ++i)
258  sv[i] = one / scalingVector[i];
259  } else {
260  for (int i = 0; i < scalingVector.size(); ++i)
261  sv[i] = scalingVector[i];
262  }
263 
264  switch (Op.getRowMap()->lib()) {
265  case Xpetra::UseTpetra:
266  MyOldScaleMatrix_Tpetra(Op, sv, doFillComplete, doOptimizeStorage);
267  break;
268 
269  case Xpetra::UseEpetra:
270  // FIXME?
271  // Utils2_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix_Epetra(Op, sv, doFillComplete, doOptimizeStorage);
272  throw std::runtime_error("FIXME");
273  break;
274 
275  default:
276  throw Exceptions::RuntimeError("Only Epetra and Tpetra matrices can be scaled.");
277  break;
278  }
279  }
280 
281  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
282  void Utils_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix_Tpetra(Matrix& Op, const Teuchos::ArrayRCP<SC>& scalingVector,
283  bool doFillComplete,
284  bool doOptimizeStorage)
285  {
286 #ifdef HAVE_MUELU_TPETRA
287  try {
288  Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tpOp = Op2NonConstTpetraCrs(Op);
289 
290  const RCP<const Tpetra::Map<LO,GO,NO> > rowMap = tpOp.getRowMap();
291  const RCP<const Tpetra::Map<LO,GO,NO> > domainMap = tpOp.getDomainMap();
292  const RCP<const Tpetra::Map<LO,GO,NO> > rangeMap = tpOp.getRangeMap();
293 
294  size_t maxRowSize = tpOp.getNodeMaxNumRowEntries();
295  if (maxRowSize == Teuchos::as<size_t>(-1)) // hasn't been determined yet
296  maxRowSize = 20;
297 
298  std::vector<SC> scaledVals(maxRowSize);
299  if (tpOp.isFillComplete())
300  tpOp.resumeFill();
301 
302  if (Op.isLocallyIndexed() == true) {
303  Teuchos::ArrayView<const LO> cols;
304  Teuchos::ArrayView<const SC> vals;
305 
306  for (size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
307  tpOp.getLocalRowView(i, cols, vals);
308  size_t nnz = tpOp.getNumEntriesInLocalRow(i);
309  if (nnz > maxRowSize) {
310  maxRowSize = nnz;
311  scaledVals.resize(maxRowSize);
312  }
313  for (size_t j = 0; j < nnz; ++j)
314  scaledVals[j] = vals[j]*scalingVector[i];
315 
316  if (nnz > 0) {
317  Teuchos::ArrayView<const SC> valview(&scaledVals[0], nnz);
318  tpOp.replaceLocalValues(i, cols, valview);
319  }
320  } //for (size_t i=0; ...
321 
322  } else {
323  Teuchos::ArrayView<const GO> cols;
324  Teuchos::ArrayView<const SC> vals;
325 
326  for (size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
327  GO gid = rowMap->getGlobalElement(i);
328  tpOp.getGlobalRowView(gid, cols, vals);
329  size_t nnz = tpOp.getNumEntriesInGlobalRow(gid);
330  if (nnz > maxRowSize) {
331  maxRowSize = nnz;
332  scaledVals.resize(maxRowSize);
333  }
334  // FIXME FIXME FIXME FIXME FIXME FIXME
335  for (size_t j = 0; j < nnz; ++j)
336  scaledVals[j] = vals[j]*scalingVector[i]; //FIXME i or gid?
337 
338  if (nnz > 0) {
339  Teuchos::ArrayView<const SC> valview(&scaledVals[0], nnz);
340  tpOp.replaceGlobalValues(gid, cols, valview);
341  }
342  } //for (size_t i=0; ...
343  }
344 
345  if (doFillComplete) {
346  if (domainMap == Teuchos::null || rangeMap == Teuchos::null)
347  throw Exceptions::RuntimeError("In Utils_kokkos::Scaling: cannot fillComplete because the domain and/or range map hasn't been defined");
348 
349  RCP<Teuchos::ParameterList> params = rcp(new Teuchos::ParameterList());
350  params->set("Optimize Storage", doOptimizeStorage);
351  params->set("No Nonlocal Changes", true);
352  Op.fillComplete(Op.getDomainMap(), Op.getRangeMap(), params);
353  }
354  } catch(...) {
355  throw Exceptions::RuntimeError("Only Tpetra::CrsMatrix types can be scaled (Err.1)");
356  }
357 #else
358  throw Exceptions::RuntimeError("Matrix scaling is not possible because Tpetra has not been enabled.");
359 #endif
360  } //MyOldScaleMatrix_Tpetra()
361 
362  template <class SC, class LO, class GO, class NO>
363  ArrayRCP<const bool> Utils_kokkos<SC, LO, GO, NO>::DetectDirichletRows(const Matrix& A, const typename Teuchos::ScalarTraits<SC>::magnitudeType& tol) {
364  LO numRows = A.getNodeNumRows();
365 
366  typedef Teuchos::ScalarTraits<SC> STS;
367 
368  ArrayRCP<bool> boundaryNodes(numRows, true);
369  for (LO row = 0; row < numRows; row++) {
370  ArrayView<const LO> indices;
371  ArrayView<const SC> vals;
372  A.getLocalRowView(row, indices, vals);
373 
374  size_t nnz = A.getNumEntriesInLocalRow(row);
375  if (nnz > 1)
376  for (size_t col = 0; col < nnz; col++)
377  if ( (indices[col] != row) && STS::magnitude(vals[col]) > tol) {
378  boundaryNodes[row] = false;
379  break;
380  }
381  }
382 
383  return boundaryNodes;
384  }
385 
386  template <class SC, class LO, class GO, class NO>
387  void Utils_kokkos<SC, LO, GO, NO>::findDirichletRows(Teuchos::RCP<Matrix> A,
388  std::vector<LO>& dirichletRows) {
389  dirichletRows.resize(0);
390  for(size_t i=0; i<A->getNodeNumRows(); i++) {
391  Teuchos::ArrayView<const LO> indices;
392  Teuchos::ArrayView<const SC> values;
393  A->getLocalRowView(i,indices,values);
394  int nnz=0;
395  for (int j=0; j<indices.size(); j++) {
396  if (abs(values[j]) > 1.0e-16) {
397  nnz++;
398  }
399  }
400  if (nnz == 1 || nnz == 2) {
401  dirichletRows.push_back(i);
402  }
403  }
404  }
405 
406  template<class SC, class LO, class GO, class NO>
407  void Utils_kokkos<SC, LO, GO, NO>::findDirichletCols(Teuchos::RCP<Matrix> A,
408  std::vector<LO>& dirichletRows,
409  std::vector<LO>& dirichletCols) {
410  Teuchos::RCP<const Map> domMap = A->getDomainMap();
411  Teuchos::RCP<const Map> colMap = A->getColMap();
412  Teuchos::RCP< Xpetra::Export<LO,GO,NO> > exporter
413  = Xpetra::ExportFactory<LO,GO,NO>::Build(colMap,domMap);
414  Teuchos::RCP<MultiVector> myColsToZero = MultiVectorFactory::Build(colMap,1);
415  Teuchos::RCP<MultiVector> globalColsToZero = MultiVectorFactory::Build(domMap,1);
416  myColsToZero->putScalar((SC)0.0);
417  globalColsToZero->putScalar((SC)0.0);
418  for(size_t i=0; i<dirichletRows.size(); i++) {
419  Teuchos::ArrayView<const LO> indices;
420  Teuchos::ArrayView<const SC> values;
421  A->getLocalRowView(dirichletRows[i],indices,values);
422  for(int j=0; j<indices.size(); j++)
423  myColsToZero->replaceLocalValue(indices[j],0,(SC)1.0);
424  }
425  globalColsToZero->doExport(*myColsToZero,*exporter,Xpetra::ADD);
426  myColsToZero->doImport(*globalColsToZero,*exporter,Xpetra::INSERT);
427  Teuchos::ArrayRCP<const SC> myCols = myColsToZero->getData(0);
428  dirichletCols.resize(colMap->getNodeNumElements());
429  for(size_t i=0; i<colMap->getNodeNumElements(); i++) {
430  if(abs(myCols[i])>0.0)
431  dirichletCols[i]=1;
432  else
433  dirichletCols[i]=0;
434  }
435  }
436 
437  template<class SC, class LO, class GO, class NO>
438  void Utils_kokkos<SC, LO, GO, NO>::Apply_BCsToMatrixRows(Teuchos::RCP<Matrix>& A,
439  std::vector<LO>& dirichletRows) {
440  for(size_t i=0; i<dirichletRows.size(); i++) {
441  Teuchos::ArrayView<const LO> indices;
442  Teuchos::ArrayView<const SC> values;
443  A->getLocalRowView(dirichletRows[i],indices,values);
444  std::vector<SC> vec;
445  vec.resize(indices.size());
446  Teuchos::ArrayView<SC> zerovalues(vec);
447  for(int j=0; j<indices.size(); j++)
448  zerovalues[j]=(SC)1.0e-32;
449  A->replaceLocalValues(dirichletRows[i],indices,zerovalues);
450  }
451  }
452 
453  template<class SC, class LO, class GO, class NO>
454  void Utils_kokkos<SC, LO, GO, NO>::Apply_BCsToMatrixCols(Teuchos::RCP<Matrix>& A,
455  std::vector<LO>& dirichletCols) {
456  for(size_t i=0; i<A->getNodeNumRows(); i++) {
457  Teuchos::ArrayView<const LO> indices;
458  Teuchos::ArrayView<const SC> values;
459  A->getLocalRowView(i,indices,values);
460  std::vector<SC> vec;
461  vec.resize(indices.size());
462  Teuchos::ArrayView<SC> zerovalues(vec);
463  for(int j=0; j<indices.size(); j++) {
464  if(dirichletCols[indices[j]]==1)
465  zerovalues[j]=(SC)1.0e-32;
466  else
467  zerovalues[j]=values[j];
468  }
469  A->replaceLocalValues(i,indices,zerovalues);
470  }
471  }
472 
473  template<class SC, class LO, class GO, class NO>
474  void Utils_kokkos<SC, LO, GO, NO>::Remove_Zeroed_Rows(Teuchos::RCP<Matrix>& A, double tol) {
475  Teuchos::RCP<const Map> rowMap = A->getRowMap();
476  RCP<Matrix> DiagMatrix = MatrixFactory::Build(rowMap,1);
477  RCP<Matrix> NewMatrix = MatrixFactory::Build(rowMap,1);
478  for(size_t i=0; i<A->getNodeNumRows(); i++) {
479  Teuchos::ArrayView<const LO> indices;
480  Teuchos::ArrayView<const SC> values;
481  A->getLocalRowView(i,indices,values);
482  int nnz=0;
483  for (int j=0; j<indices.size(); j++) {
484  if (abs(values[j]) > tol) {
485  nnz++;
486  }
487  }
488  SC one = (SC)1.0;
489  SC zero = (SC)0.0;
490  GO row = rowMap->getGlobalElement(i);
491  if (nnz == 0) {
492  DiagMatrix->insertGlobalValues(row,
493  Teuchos::ArrayView<GO>(&row,1),
494  Teuchos::ArrayView<SC>(&one,1));
495  }
496  else {
497  DiagMatrix->insertGlobalValues(row,
498  Teuchos::ArrayView<GO>(&row,1),
499  Teuchos::ArrayView<SC>(&zero,1));
500  }
501  }
502  DiagMatrix->fillComplete();
503  A->fillComplete();
504  // add matrices together
505  RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
506  Xpetra::MatrixMatrix<SC,LO,GO,NO>::TwoMatrixAdd(*DiagMatrix,false,(SC)1.0,*A,false,(SC)1.0,NewMatrix,*out);
507  NewMatrix->fillComplete();
508  A=NewMatrix;
509  }
510 
511 } //namespace MueLu
512 
513 #define MUELU_UTILITIES_KOKKOS_SHORT
514 #endif // MUELU_UTILITIES_KOKKOS_DEF_HPP
EpetraExt_CrsMatrixIn.h
MueLu
Namespace for MueLu classes and methods.
Definition: MueLu_BrickAggregationFactory_decl.hpp:76
EpetraExt_MatrixMatrix.h
MueLu_ConfigDefs.hpp
Kokkos::abs
KOKKOS_INLINE_FUNCTION RealType abs(const complex< RealType > &x)
MueLu_Utilities_kokkos_decl.hpp
EpetraExt_RowMatrixOut.h
EpetraExt_BlockMapOut.h
EpetraExt_MultiVectorOut.h
EpetraExt_MultiVectorIn.h
EpetraExt_BlockMapIn.h