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) 76 #ifdef HAVE_MUELU_TPETRA 77 #include <MatrixMarket_Tpetra.hpp> 78 #include <Tpetra_RowMatrixTransposer.hpp> 79 #include <TpetraExt_MatrixMatrix.hpp> 85 #ifdef HAVE_MUELU_EPETRA 108 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
109 Teuchos::ArrayRCP<Scalar> Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetMatrixDiagonal(
const Matrix& A) {
112 size_t numRows = A.getRowMap()->getNodeNumElements();
113 Teuchos::ArrayRCP<SC> diag(numRows);
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);
121 for (; j < cols.size(); ++j) {
122 if (Teuchos::as<size_t>(cols[j]) == i) {
127 if (j == cols.size()) {
129 diag[i] = Teuchos::ScalarTraits<SC>::zero();
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) {
139 RCP<const Map> rowMap = A.getRowMap();
140 RCP<Vector> diag = VectorFactory::Build(rowMap);
141 ArrayRCP<SC> diagVals = diag->getDataNonConst(0);
143 size_t numRows = rowMap->getNodeNumElements();
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);
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];
156 diagVals[i]=Teuchos::ScalarTraits<SC>::zero();
160 if (j == cols.size()) {
162 diagVals[i]=Teuchos::ScalarTraits<SC>::zero();
170 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
171 Teuchos::ArrayRCP<Scalar> Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetLumpedMatrixDiagonal(
const Matrix &A) {
173 size_t numRows = A.getRowMap()->getNodeNumElements();
174 Teuchos::ArrayRCP<SC> diag(numRows);
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);
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]);
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) {
193 RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
194 RCP<Vector> localDiag = VectorFactory::Build(rowMap);
197 const CrsMatrixWrap* crsOp =
dynamic_cast<const CrsMatrixWrap*
>(&A);
199 throw Exceptions::RuntimeError(
"cast to CrsMatrixWrap failed");
201 Teuchos::ArrayRCP<size_t> offsets;
202 crsOp->getLocalDiagOffsets(offsets);
203 crsOp->getLocalDiagCopy(*localDiag,offsets());
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;
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);
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,
227 bool doOptimizeStorage)
229 SC one = Teuchos::ScalarTraits<SC>::one();
230 Teuchos::ArrayRCP<SC> sv(scalingVector.size());
232 for (
int i = 0; i < scalingVector.size(); ++i)
233 sv[i] = one / scalingVector[i];
235 for (
int i = 0; i < scalingVector.size(); ++i)
236 sv[i] = scalingVector[i];
239 switch (Op.getRowMap()->lib()) {
241 MyOldScaleMatrix_Tpetra(Op, sv, doFillComplete, doOptimizeStorage);
247 throw std::runtime_error(
"FIXME");
251 throw Exceptions::RuntimeError(
"Only Epetra and Tpetra matrices can be scaled.");
256 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
258 throw Exceptions::RuntimeError(
"MyOldScaleMatrix_Epetra: Epetra needs SC=double and LO=GO=int.");
261 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
262 void Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix_Tpetra(Matrix& Op,
const Teuchos::ArrayRCP<SC>& scalingVector,
264 bool doOptimizeStorage)
266 #ifdef HAVE_MUELU_TPETRA 268 Tpetra::CrsMatrix<SC,LO,GO,NO>& tpOp = Op2NonConstTpetraCrs(Op);
270 const RCP<const Tpetra::Map<LO,GO,NO> > rowMap = tpOp.getRowMap();
271 const RCP<const Tpetra::Map<LO,GO,NO> > domainMap = tpOp.getDomainMap();
272 const RCP<const Tpetra::Map<LO,GO,NO> > rangeMap = tpOp.getRangeMap();
274 size_t maxRowSize = tpOp.getNodeMaxNumRowEntries();
275 if (maxRowSize == Teuchos::as<size_t>(-1))
278 std::vector<SC> scaledVals(maxRowSize);
279 if (tpOp.isFillComplete())
282 if (Op.isLocallyIndexed() ==
true) {
283 Teuchos::ArrayView<const LO> cols;
284 Teuchos::ArrayView<const SC> vals;
286 for (
size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
287 tpOp.getLocalRowView(i, cols, vals);
288 size_t nnz = tpOp.getNumEntriesInLocalRow(i);
289 if (nnz > maxRowSize) {
291 scaledVals.resize(maxRowSize);
293 for (
size_t j = 0; j < nnz; ++j)
294 scaledVals[j] = vals[j]*scalingVector[i];
297 Teuchos::ArrayView<const SC> valview(&scaledVals[0], nnz);
298 tpOp.replaceLocalValues(i, cols, valview);
303 Teuchos::ArrayView<const GO> cols;
304 Teuchos::ArrayView<const SC> vals;
306 for (
size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
307 GO gid = rowMap->getGlobalElement(i);
308 tpOp.getGlobalRowView(gid, cols, vals);
309 size_t nnz = tpOp.getNumEntriesInGlobalRow(gid);
310 if (nnz > maxRowSize) {
312 scaledVals.resize(maxRowSize);
315 for (
size_t j = 0; j < nnz; ++j)
316 scaledVals[j] = vals[j]*scalingVector[i];
319 Teuchos::ArrayView<const SC> valview(&scaledVals[0], nnz);
320 tpOp.replaceGlobalValues(gid, cols, valview);
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");
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);
335 throw Exceptions::RuntimeError(
"Only Tpetra::CrsMatrix types can be scaled (Err.1)");
338 throw Exceptions::RuntimeError(
"Matrix scaling is not possible because Tpetra has not been enabled.");
342 template<
class MatrixType,
class BNodesType>
345 typedef typename MatrixType::ordinal_type
LO;
346 typedef typename MatrixType::value_type
SC;
347 typedef Kokkos::ArithTraits<SC>
ATS;
351 typename ATS::mag_type
tol;
358 tol = ATS::magnitude(tol_);
361 KOKKOS_INLINE_FUNCTION
364 auto length = rowView.length;
367 for (decltype(length) colID = 0; colID < length; colID++)
368 if ((rowView.colidx(colID) != row) && (ATS::magnitude(rowView.value(colID)) >
tol)) {
375 template <
class SC,
class LO,
class GO,
class NO>
377 auto localMatrix = A.getLocalMatrix();
385 return boundaryNodes;
388 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
392 return MueLu::DetectDirichletRows<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, tol);
395 template <
class Node>
399 return MueLu::DetectDirichletRows<double,int,int,Node>(A, tol);
404 #define MUELU_UTILITIES_KOKKOS_SHORT 405 #endif // MUELU_UTILITIES_KOKKOS_DEF_HPP virtual size_t getNodeNumRows() const=0
DetectDirichletFunctor(MatrixType localMatrix_, BNodesType boundaryNodes_, SC tol_)
MatrixType::ordinal_type LO
Namespace for MueLu classes and methods.
KOKKOS_INLINE_FUNCTION void operator()(const LO row) const
Kokkos::ArithTraits< SC > ATS
Kokkos::View< const bool *, typename NO::device_type > DetectDirichletRows(const Xpetra::Matrix< SC, LO, GO, NO > &A, const typename Teuchos::ScalarTraits< SC >::magnitudeType &tol)
MatrixType::value_type SC
void parallel_for(const ExecPolicy &policy, const FunctorType &functor, const std::string &str="", typename Impl::enable_if< ! Impl::is_integral< ExecPolicy >::value >::type *=0)