46 #ifndef MUELU_UTILITIES_KOKKOS_DEF_HPP 47 #define MUELU_UTILITIES_KOKKOS_DEF_HPP 49 #include <Teuchos_DefaultComm.hpp> 50 #include <Teuchos_ParameterList.hpp> 54 #ifdef HAVE_MUELU_EPETRA 56 # include "Epetra_MpiComm.h" 60 #include <Kokkos_ArithTraits.hpp> 61 #include <Kokkos_Core.hpp> 62 #include <Kokkos_CrsMatrix.hpp> 64 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_EPETRAEXT) 71 #include <Xpetra_EpetraUtils.hpp> 72 #include <Xpetra_EpetraMultiVector.hpp> 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> 85 #ifdef HAVE_MUELU_EPETRA 86 #include <Xpetra_EpetraMap.hpp> 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> 108 #ifdef HAVE_MUELU_EPETRA 113 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
114 Teuchos::ArrayRCP<Scalar> Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetMatrixDiagonal(
const Matrix& A) {
117 size_t numRows = A.getRowMap()->getNodeNumElements();
118 Teuchos::ArrayRCP<SC> diag(numRows);
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);
126 for (; j < cols.size(); ++j) {
127 if (Teuchos::as<size_t>(cols[j]) == i) {
132 if (j == cols.size()) {
134 diag[i] = Teuchos::ScalarTraits<SC>::zero();
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) {
144 RCP<const Map> rowMap = A.getRowMap();
145 RCP<Vector> diag = VectorFactory::Build(rowMap);
146 ArrayRCP<SC> diagVals = diag->getDataNonConst(0);
148 size_t numRows = rowMap->getNodeNumElements();
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);
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];
161 diagVals[i]=Teuchos::ScalarTraits<SC>::zero();
165 if (j == cols.size()) {
167 diagVals[i]=Teuchos::ScalarTraits<SC>::zero();
175 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
176 Teuchos::ArrayRCP<Scalar> Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetLumpedMatrixDiagonal(
const Matrix &A) {
178 size_t numRows = A.getRowMap()->getNodeNumElements();
179 Teuchos::ArrayRCP<SC> diag(numRows);
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);
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]);
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) {
198 RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
199 RCP<Vector> localDiag = VectorFactory::Build(rowMap);
202 const CrsMatrixWrap* crsOp =
dynamic_cast<const CrsMatrixWrap*
>(&A);
204 throw Exceptions::RuntimeError(
"cast to CrsMatrixWrap failed");
206 Teuchos::ArrayRCP<size_t> offsets;
207 crsOp->getLocalDiagOffsets(offsets);
208 crsOp->getLocalDiagCopy(*localDiag,offsets());
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;
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);
224 diagonal->doImport(*localDiag, *(importer), Xpetra::INSERT);
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,
232 bool doOptimizeStorage)
234 SC one = Teuchos::ScalarTraits<SC>::one();
235 Teuchos::ArrayRCP<SC> sv(scalingVector.size());
237 for (
int i = 0; i < scalingVector.size(); ++i)
238 sv[i] = one / scalingVector[i];
240 for (
int i = 0; i < scalingVector.size(); ++i)
241 sv[i] = scalingVector[i];
244 switch (Op.getRowMap()->lib()) {
245 case Xpetra::UseTpetra:
246 MyOldScaleMatrix_Tpetra(Op, sv, doFillComplete, doOptimizeStorage);
249 case Xpetra::UseEpetra:
252 throw std::runtime_error(
"FIXME");
256 throw Exceptions::RuntimeError(
"Only Epetra and Tpetra matrices can be scaled.");
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.");
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,
269 bool doOptimizeStorage)
271 #ifdef HAVE_MUELU_TPETRA 273 Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tpOp = Op2NonConstTpetraCrs(Op);
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();
279 size_t maxRowSize = tpOp.getNodeMaxNumRowEntries();
280 if (maxRowSize == Teuchos::as<size_t>(-1))
283 std::vector<SC> scaledVals(maxRowSize);
284 if (tpOp.isFillComplete())
287 if (Op.isLocallyIndexed() ==
true) {
288 Teuchos::ArrayView<const LO> cols;
289 Teuchos::ArrayView<const SC> vals;
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) {
296 scaledVals.resize(maxRowSize);
298 for (
size_t j = 0; j < nnz; ++j)
299 scaledVals[j] = vals[j]*scalingVector[i];
302 Teuchos::ArrayView<const SC> valview(&scaledVals[0], nnz);
303 tpOp.replaceLocalValues(i, cols, valview);
308 Teuchos::ArrayView<const GO> cols;
309 Teuchos::ArrayView<const SC> vals;
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) {
317 scaledVals.resize(maxRowSize);
320 for (
size_t j = 0; j < nnz; ++j)
321 scaledVals[j] = vals[j]*scalingVector[i];
324 Teuchos::ArrayView<const SC> valview(&scaledVals[0], nnz);
325 tpOp.replaceGlobalValues(gid, cols, valview);
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");
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);
340 throw Exceptions::RuntimeError(
"Only Tpetra::CrsMatrix types can be scaled (Err.1)");
343 throw Exceptions::RuntimeError(
"Matrix scaling is not possible because Tpetra has not been enabled.");
347 template <
class SC,
class LO,
class GO,
class NO>
349 typedef Kokkos::ArithTraits<SC> ATS;
351 auto localMatrix = A.getLocalMatrix();
352 LO numRows = A.getNodeNumRows();
356 auto rowView = localMatrix.template row<LO>(row);
357 auto length = rowView.length;
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;
367 return boundaryNodes;
370 template <
class SC,
class LO,
class GO,
class NO>
372 return MueLu::DetectDirichletRows<SC,LO,GO,NO>(A, tol);
375 template <
class Node>
377 return MueLu::DetectDirichletRows<double,int,int,Node>(A, tol);
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)