MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_Utilities_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // MueLu: A package for multigrid based preconditioning
4 //
5 // Copyright 2012 NTESS and the MueLu contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef MUELU_UTILITIES_DEF_HPP
11 #define MUELU_UTILITIES_DEF_HPP
12 
13 #include <Teuchos_DefaultComm.hpp>
14 #include <Teuchos_ParameterList.hpp>
15 #include <iostream>
16 
17 #include "MueLu_ConfigDefs.hpp"
18 #include "Xpetra_TpetraRowMatrix.hpp"
19 
20 #ifdef HAVE_MUELU_EPETRA
21 #ifdef HAVE_MPI
22 #include "Epetra_MpiComm.h"
23 #endif
24 #endif
25 
26 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_EPETRAEXT)
27 #include <EpetraExt_MatrixMatrix.h>
28 #include <EpetraExt_RowMatrixOut.h>
29 #include <EpetraExt_MultiVectorOut.h>
30 #include <EpetraExt_CrsMatrixIn.h>
31 #include <EpetraExt_MultiVectorIn.h>
32 #include <EpetraExt_BlockMapIn.h>
33 #include <Xpetra_EpetraUtils.hpp>
34 #include <Xpetra_EpetraMultiVector.hpp>
35 #include <EpetraExt_BlockMapOut.h>
36 #endif
37 
38 #include <MatrixMarket_Tpetra.hpp>
39 #include <Tpetra_RowMatrixTransposer.hpp>
40 #include <TpetraExt_MatrixMatrix.hpp>
41 #include <Xpetra_TpetraMultiVector.hpp>
42 #include <Xpetra_TpetraOperator.hpp>
43 #include <Xpetra_TpetraCrsMatrix.hpp>
44 #include <Xpetra_TpetraBlockCrsMatrix.hpp>
45 
46 #ifdef HAVE_MUELU_EPETRA
47 #include <Xpetra_EpetraMap.hpp>
48 #endif
49 
50 #include <Xpetra_BlockedCrsMatrix.hpp>
51 //#include <Xpetra_DefaultPlatform.hpp>
52 #include <Xpetra_IO.hpp>
53 #include <Xpetra_Map.hpp>
54 #include <Xpetra_MapFactory.hpp>
55 #include <Xpetra_Matrix.hpp>
56 #include <Xpetra_MatrixFactory.hpp>
57 #include <Xpetra_MultiVector.hpp>
58 #include <Xpetra_MultiVectorFactory.hpp>
59 #include <Xpetra_Operator.hpp>
60 #include <Xpetra_Vector.hpp>
61 #include <Xpetra_VectorFactory.hpp>
62 
63 #include <Xpetra_MatrixMatrix.hpp>
64 
65 #include <MueLu_Utilities_decl.hpp>
66 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_ML)
67 #include <ml_operator.h>
68 #include <ml_epetra_utils.h>
69 #endif
70 
71 namespace MueLu {
72 
73 #ifdef HAVE_MUELU_EPETRA
74 // using Xpetra::EpetraCrsMatrix; // TODO: mv in Xpetra_UseShortNamesScalar
75 // using Xpetra::EpetraMultiVector;
76 #endif
77 
78 #ifdef HAVE_MUELU_EPETRA
79 template <typename SC, typename LO, typename GO, typename NO>
80 RCP<Xpetra::CrsMatrixWrap<SC, LO, GO, NO>> Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap(RCP<Epetra_CrsMatrix>& epAB) {
81  return Xpetra::Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<SC, LO, GO, NO>(epAB);
82 }
83 #endif
84 
85 #ifdef HAVE_MUELU_EPETRA
86 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
87 RCP<const Epetra_MultiVector> Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2EpetraMV(const RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> vec) {
88  RCP<const Xpetra::EpetraMultiVectorT<GlobalOrdinal, Node>> tmpVec = rcp_dynamic_cast<Xpetra::EpetraMultiVectorT<GlobalOrdinal, Node>>(vec);
89  if (tmpVec == Teuchos::null)
90  throw Exceptions::BadCast("Cast from Xpetra::MultiVector to Xpetra::EpetraMultiVector failed");
91  return tmpVec->getEpetra_MultiVector();
92 }
93 
94 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
95 RCP<Epetra_MultiVector> Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2NonConstEpetraMV(RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> vec) {
96  RCP<const Xpetra::EpetraMultiVectorT<GlobalOrdinal, Node>> tmpVec = rcp_dynamic_cast<Xpetra::EpetraMultiVectorT<GlobalOrdinal, Node>>(vec);
97  if (tmpVec == Teuchos::null)
98  throw Exceptions::BadCast("Cast from Xpetra::MultiVector to Xpetra::EpetraMultiVector failed");
99  return tmpVec->getEpetra_MultiVector();
100 }
101 
102 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
103 Epetra_MultiVector& Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2NonConstEpetraMV(Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& vec) {
104  const Xpetra::EpetraMultiVectorT<GlobalOrdinal, Node>& tmpVec = dynamic_cast<const Xpetra::EpetraMultiVectorT<GlobalOrdinal, Node>&>(vec);
105  return *(tmpVec.getEpetra_MultiVector());
106 }
107 
108 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
109 const Epetra_MultiVector& Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2EpetraMV(const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& vec) {
110  const Xpetra::EpetraMultiVectorT<GlobalOrdinal, Node>& tmpVec = dynamic_cast<const Xpetra::EpetraMultiVectorT<GlobalOrdinal, Node>&>(vec);
111  return *(tmpVec.getEpetra_MultiVector());
112 }
113 
114 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
115 RCP<const Epetra_CrsMatrix> Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2EpetraCrs(RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Op) {
116  RCP<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(Op);
117  if (crsOp == Teuchos::null)
118  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
119  const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>& tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(crsOp->getCrsMatrix());
120  if (tmp_ECrsMtx == Teuchos::null)
121  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
122  return tmp_ECrsMtx->getEpetra_CrsMatrix();
123 }
124 
125 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
126 RCP<Epetra_CrsMatrix> Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2NonConstEpetraCrs(RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Op) {
127  RCP<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(Op);
128  if (crsOp == Teuchos::null)
129  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
130  const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>& tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(crsOp->getCrsMatrix());
131  if (tmp_ECrsMtx == Teuchos::null)
132  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
133  return tmp_ECrsMtx->getEpetra_CrsMatrixNonConst();
134 }
135 
136 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
137 const Epetra_CrsMatrix& Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2EpetraCrs(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op) {
138  try {
139  const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>& crsOp = dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>&>(Op);
140  try {
141  const Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>& tmp_ECrsMtx = dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>&>(*crsOp.getCrsMatrix());
142  return *tmp_ECrsMtx.getEpetra_CrsMatrix();
143  } catch (std::bad_cast&) {
144  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
145  }
146  } catch (std::bad_cast&) {
147  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
148  }
149 }
150 
151 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
152 Epetra_CrsMatrix& Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2NonConstEpetraCrs(Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op) {
153  try {
154  Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>& crsOp = dynamic_cast<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>&>(Op);
155  try {
156  Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>& tmp_ECrsMtx = dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>&>(*crsOp.getCrsMatrix());
157  return *tmp_ECrsMtx.getEpetra_CrsMatrixNonConst();
158  } catch (std::bad_cast&) {
159  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
160  }
161  } catch (std::bad_cast&) {
162  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
163  }
164 }
165 
166 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
167 const Epetra_Map& Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Map2EpetraMap(const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>& map) {
168  RCP<const Xpetra::EpetraMapT<GlobalOrdinal, Node>> xeMap = rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal, Node>>(rcpFromRef(map));
169  if (xeMap == Teuchos::null)
170  throw Exceptions::BadCast("Utilities::Map2EpetraMap : Cast from Xpetra::Map to Xpetra::EpetraMap failed");
171  return xeMap->getEpetra_Map();
172 }
173 #endif
174 
175 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
176 RCP<const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
177 Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2TpetraMV(RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> const vec) {
178  RCP<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tmpVec = rcp_dynamic_cast<Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(vec);
179  if (tmpVec == Teuchos::null)
180  throw Exceptions::BadCast("Cast from Xpetra::MultiVector to Xpetra::TpetraMultiVector failed");
181  return tmpVec->getTpetra_MultiVector();
182 }
183 
184 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
185 RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2NonConstTpetraMV(RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> vec) {
186  RCP<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tmpVec = rcp_dynamic_cast<Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(vec);
187  if (tmpVec == Teuchos::null)
188  throw Exceptions::BadCast("Cast from Xpetra::MultiVector to Xpetra::TpetraMultiVector failed");
189  return tmpVec->getTpetra_MultiVector();
190 }
191 
192 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
193 Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2NonConstTpetraMV(Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& vec) {
194  const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tmpVec = dynamic_cast<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>&>(vec);
195  return *(tmpVec.getTpetra_MultiVector());
196 }
197 
198 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
199 RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2NonConstTpetraMV2(Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& vec) {
200  const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tmpVec = dynamic_cast<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>&>(vec);
201  return tmpVec.getTpetra_MultiVector();
202 }
203 
204 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
205 const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>&
206 Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2TpetraMV(const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& vec) {
207  const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tmpVec = dynamic_cast<const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>&>(vec);
208  return *(tmpVec.getTpetra_MultiVector());
209 }
210 
211 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
212 RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2TpetraCrs(RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Op) {
213  // Get the underlying Tpetra Mtx
214  RCP<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(Op);
215  if (crsOp == Teuchos::null)
216  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
217  const RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(crsOp->getCrsMatrix());
218  if (tmp_ECrsMtx == Teuchos::null)
219  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
220  return tmp_ECrsMtx->getTpetra_CrsMatrix();
221 }
222 
223 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
224 RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2NonConstTpetraCrs(RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Op) {
225  RCP<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(Op);
226  if (crsOp == Teuchos::null)
227  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
228  const RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(crsOp->getCrsMatrix());
229  if (tmp_ECrsMtx == Teuchos::null)
230  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
231  return tmp_ECrsMtx->getTpetra_CrsMatrixNonConst();
232 }
233 
234 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
235 const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2TpetraCrs(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op) {
236  try {
237  const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>& crsOp = dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>&>(Op);
238  try {
239  const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tmp_ECrsMtx = dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>&>(*crsOp.getCrsMatrix());
240  return *tmp_ECrsMtx.getTpetra_CrsMatrix();
241  } catch (std::bad_cast&) {
242  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
243  }
244  } catch (std::bad_cast&) {
245  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
246  }
247 }
248 
249 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
250 Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2NonConstTpetraCrs(Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op) {
251  try {
252  Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>& crsOp = dynamic_cast<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>&>(Op);
253  try {
254  Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tmp_ECrsMtx = dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>&>(*crsOp.getCrsMatrix());
255  return *tmp_ECrsMtx.getTpetra_CrsMatrixNonConst();
256  } catch (std::bad_cast&) {
257  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
258  }
259  } catch (std::bad_cast&) {
260  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
261  }
262 }
263 
264 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
265 RCP<const Tpetra::BlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2TpetraBlockCrs(RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Op) {
266  using XCrsMatrixWrap = Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
267  // Get the underlying Tpetra Mtx
268  RCP<const XCrsMatrixWrap> crsOp = rcp_dynamic_cast<const XCrsMatrixWrap>(Op);
269  if (crsOp == Teuchos::null)
270  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
271  const RCP<const Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(crsOp->getCrsMatrix());
272  if (tmp_ECrsMtx == Teuchos::null)
273  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraBlockCrsMatrix failed");
274  return tmp_ECrsMtx->getTpetra_BlockCrsMatrix();
275 }
276 
277 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
278 RCP<Tpetra::BlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2NonConstTpetraBlockCrs(RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Op) {
279  using XCrsMatrixWrap = Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
280  RCP<const XCrsMatrixWrap> crsOp = rcp_dynamic_cast<const XCrsMatrixWrap>(Op);
281  if (crsOp == Teuchos::null)
282  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
283  const RCP<const Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(crsOp->getCrsMatrix());
284  if (tmp_ECrsMtx == Teuchos::null)
285  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraBlockCrsMatrix failed");
286  return tmp_ECrsMtx->getTpetra_BlockCrsMatrixNonConst();
287 }
288 
289 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
290 const Tpetra::BlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2TpetraBlockCrs(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op) {
291  try {
292  using XCrsMatrixWrap = Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
293  const XCrsMatrixWrap& crsOp = dynamic_cast<const XCrsMatrixWrap&>(Op);
294  try {
295  const Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tmp_ECrsMtx = dynamic_cast<const Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>&>(*crsOp.getCrsMatrix());
296  return *tmp_ECrsMtx.getTpetra_BlockCrsMatrix();
297  } catch (std::bad_cast&) {
298  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraBlockCrsMatrix failed");
299  }
300  } catch (std::bad_cast&) {
301  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
302  }
303 }
304 
305 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
306 Tpetra::BlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2NonConstTpetraBlockCrs(Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op) {
307  try {
308  using XCrsMatrixWrap = Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
309  XCrsMatrixWrap& crsOp = dynamic_cast<XCrsMatrixWrap&>(Op);
310  try {
311  Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tmp_ECrsMtx = dynamic_cast<Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>&>(*crsOp.getCrsMatrix());
312  return *tmp_ECrsMtx.getTpetra_BlockCrsMatrixNonConst();
313  } catch (std::bad_cast&) {
314  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraBlockCrsMatrix failed");
315  }
316  } catch (std::bad_cast&) {
317  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
318  }
319 }
320 
321 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
322 RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2TpetraRow(RCP<const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Op) {
323  RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> mat = rcp_dynamic_cast<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(Op);
324  RCP<const Xpetra::TpetraRowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> rmat = rcp_dynamic_cast<const Xpetra::TpetraRowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(Op);
325  if (!mat.is_null()) {
326  RCP<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat);
327  if (crsOp == Teuchos::null)
328  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
329 
330  RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> crsMat = crsOp->getCrsMatrix();
331  const RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tmp_Crs = rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(crsMat);
332  RCP<const Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tmp_BlockCrs;
333  if (!tmp_Crs.is_null()) {
334  return tmp_Crs->getTpetra_CrsMatrixNonConst();
335  } else {
336  tmp_BlockCrs = rcp_dynamic_cast<const Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(crsMat);
337  if (tmp_BlockCrs.is_null())
338  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix and Xpetra::TpetraBlockCrsMatrix failed");
339  return tmp_BlockCrs->getTpetra_BlockCrsMatrixNonConst();
340  }
341  } else if (!rmat.is_null()) {
342  return rmat->getTpetra_RowMatrix();
343  } else {
344  RCP<const Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpOp = rcp_dynamic_cast<const Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(Op, true);
345  RCP<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tOp = tpOp->getOperatorConst();
346  RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tRow = rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tOp, true);
347  return tRow;
348  }
349 }
350 
351 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
352 RCP<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2NonConstTpetraRow(RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Op) {
353  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> mat = rcp_dynamic_cast<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(Op);
354  RCP<Xpetra::TpetraRowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> rmat = rcp_dynamic_cast<Xpetra::TpetraRowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(Op);
355  if (!mat.is_null()) {
356  RCP<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat);
357  if (crsOp == Teuchos::null)
358  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
359 
360  RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> crsMat = crsOp->getCrsMatrix();
361  const RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tmp_Crs = rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(crsMat);
362  RCP<const Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tmp_BlockCrs;
363  if (!tmp_Crs.is_null()) {
364  return tmp_Crs->getTpetra_CrsMatrixNonConst();
365  } else {
366  tmp_BlockCrs = rcp_dynamic_cast<const Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(crsMat);
367  if (tmp_BlockCrs.is_null())
368  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix and Xpetra::TpetraBlockCrsMatrix failed");
369  return tmp_BlockCrs->getTpetra_BlockCrsMatrixNonConst();
370  }
371  } else if (!rmat.is_null()) {
372  return rmat->getTpetra_RowMatrixNonConst();
373  } else {
374  RCP<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpOp = rcp_dynamic_cast<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(Op, true);
375  RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tOp = tpOp->getOperator();
376  RCP<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tRow = rcp_dynamic_cast<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tOp, true);
377  return tRow;
378  }
379 }
380 
381 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
382 const RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Map2TpetraMap(const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>& map) {
383  const RCP<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>>& tmp_TMap = rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>>(rcpFromRef(map));
384  if (tmp_TMap == Teuchos::null)
385  throw Exceptions::BadCast("Utilities::Map2TpetraMap : Cast from Xpetra::Map to Xpetra::TpetraMap failed");
386  return tmp_TMap->getTpetra_Map();
387 }
388 
389 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
390 void Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix(Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op, const Teuchos::ArrayRCP<const Scalar>& scalingVector, bool doInverse,
391  bool doFillComplete,
392  bool doOptimizeStorage) {
393  Scalar one = Teuchos::ScalarTraits<Scalar>::one();
394  Teuchos::ArrayRCP<Scalar> sv(scalingVector.size());
395  if (doInverse) {
396  for (int i = 0; i < scalingVector.size(); ++i)
397  sv[i] = one / scalingVector[i];
398  } else {
399  for (int i = 0; i < scalingVector.size(); ++i)
400  sv[i] = scalingVector[i];
401  }
402 
403  switch (Op.getRowMap()->lib()) {
404  case Xpetra::UseTpetra:
405  MyOldScaleMatrix_Tpetra(Op, sv, doFillComplete, doOptimizeStorage);
406  break;
407 
408  case Xpetra::UseEpetra:
409  MyOldScaleMatrix_Epetra(Op, sv, doFillComplete, doOptimizeStorage);
410  break;
411 
412  default:
413  throw Exceptions::RuntimeError("Only Epetra and Tpetra matrices can be scaled.");
414  }
415 }
416 
417 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
418 void Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix_Epetra(Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& /* Op */, const Teuchos::ArrayRCP<Scalar>& /* scalingVector */, bool /* doFillComplete */, bool /* doOptimizeStorage */) {
419  throw Exceptions::RuntimeError("MyOldScaleMatrix_Epetra: Epetra needs SC=double and LO=GO=int.");
420 }
421 
422 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
423 void Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix_Tpetra(Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op, const Teuchos::ArrayRCP<Scalar>& scalingVector,
424  bool doFillComplete,
425  bool doOptimizeStorage) {
426  try {
427  Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tpOp = Op2NonConstTpetraCrs(Op);
428 
429  const RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> rowMap = tpOp.getRowMap();
430  const RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> domainMap = tpOp.getDomainMap();
431  const RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> rangeMap = tpOp.getRangeMap();
432 
433  size_t maxRowSize = tpOp.getLocalMaxNumRowEntries();
434  if (maxRowSize == Teuchos::as<size_t>(-1)) // hasn't been determined yet
435  maxRowSize = 20;
436 
437  if (tpOp.isFillComplete())
438  tpOp.resumeFill();
439 
440  if (Op.isLocallyIndexed() == true) {
441  typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_inds_host_view_type cols;
442  typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::values_host_view_type vals;
443 
444  for (size_t i = 0; i < rowMap->getLocalNumElements(); ++i) {
445  tpOp.getLocalRowView(i, cols, vals);
446 
447  size_t nnz = tpOp.getNumEntriesInLocalRow(i);
448  typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::nonconst_values_host_view_type scaledVals("ScaledVals", nnz);
449 
450  for (size_t j = 0; j < nnz; ++j) {
451  scaledVals[j] = scalingVector[i] * vals[j];
452  }
453 
454  if (nnz > 0) {
455  tpOp.replaceLocalValues(i, cols, scaledVals);
456  }
457  } // for (size_t i=0; ...
458 
459  } else {
460  typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::global_inds_host_view_type cols;
461  typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::values_host_view_type vals;
462 
463  for (size_t i = 0; i < rowMap->getLocalNumElements(); ++i) {
464  GlobalOrdinal gid = rowMap->getGlobalElement(i);
465  tpOp.getGlobalRowView(gid, cols, vals);
466  size_t nnz = tpOp.getNumEntriesInGlobalRow(gid);
467  typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::nonconst_values_host_view_type scaledVals("ScaledVals", nnz);
468 
469  // FIXME FIXME FIXME FIXME FIXME FIXME
470  for (size_t j = 0; j < nnz; ++j)
471  scaledVals[j] = scalingVector[i] * vals[j]; // FIXME i or gid?
472 
473  if (nnz > 0) {
474  tpOp.replaceGlobalValues(gid, cols, scaledVals);
475  }
476  } // for (size_t i=0; ...
477  }
478 
479  if (doFillComplete) {
480  if (domainMap == Teuchos::null || rangeMap == Teuchos::null)
481  throw Exceptions::RuntimeError("In Utilities::Scaling: cannot fillComplete because the domain and/or range map hasn't been defined");
482 
483  RCP<Teuchos::ParameterList> params = rcp(new Teuchos::ParameterList());
484  params->set("Optimize Storage", doOptimizeStorage);
485  params->set("No Nonlocal Changes", true);
486  Op.fillComplete(Op.getDomainMap(), Op.getRangeMap(), params);
487  }
488  } catch (...) {
489  throw Exceptions::RuntimeError("Only Tpetra::CrsMatrix types can be scaled (Err.1)");
490  }
491 } // MyOldScaleMatrix_Tpetra()
492 
493 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
494 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
495 Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
496  Transpose(Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op, bool /* optimizeTranspose */, const std::string& label, const Teuchos::RCP<Teuchos::ParameterList>& params) {
497 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_EPETRAEXT)
498  std::string TorE = "epetra";
499 #else
500  std::string TorE = "tpetra";
501 #endif
502 
503 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_EPETRAEXT)
504  try {
505  const Epetra_CrsMatrix& epetraOp = Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2NonConstEpetraCrs(Op);
506  (void)epetraOp; // silence "unused variable" compiler warning
507  } catch (...) {
508  TorE = "tpetra";
509  }
510 #endif
511 
512  if (TorE == "tpetra") {
513  using Helpers = Xpetra::Helpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
514  /***************************************************************/
515  if (Helpers::isTpetraCrs(Op)) {
516  const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tpetraOp = Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2TpetraCrs(Op);
517 
518  RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> A;
519  Tpetra::RowMatrixTransposer<Scalar, LocalOrdinal, GlobalOrdinal, Node> transposer(rcpFromRef(tpetraOp), label); // more than meets the eye
520 
521  {
522  using Teuchos::ParameterList;
523  using Teuchos::rcp;
524  RCP<ParameterList> transposeParams = params.is_null() ? rcp(new ParameterList) : rcp(new ParameterList(*params));
525  transposeParams->set("sort", false);
526  A = transposer.createTranspose(transposeParams);
527  }
528 
529  RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> AA = rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(A));
530  RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> AAA = rcp_implicit_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(AA);
531  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> AAAA = rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(AAA));
532  if (!AAAA->isFillComplete())
533  AAAA->fillComplete(Op.getRangeMap(), Op.getDomainMap());
534 
535  if (Op.IsView("stridedMaps"))
536  AAAA->CreateView("stridedMaps", Teuchos::rcpFromRef(Op), true /*doTranspose*/);
537 
538  return AAAA;
539  } else if (Helpers::isTpetraBlockCrs(Op)) {
540  using XMatrix = Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
541  using XCrsMatrix = Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
542  using XCrsMatrixWrap = Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
543  using BCRS = Tpetra::BlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
544  // using CRS = Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
545  const BCRS& tpetraOp = Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2TpetraBlockCrs(Op);
546 
547  RCP<BCRS> At;
548  {
549  Tpetra::BlockCrsMatrixTransposer<Scalar, LocalOrdinal, GlobalOrdinal, Node> transposer(rcpFromRef(tpetraOp), label);
550 
551  using Teuchos::ParameterList;
552  using Teuchos::rcp;
553  RCP<ParameterList> transposeParams = params.is_null() ? rcp(new ParameterList) : rcp(new ParameterList(*params));
554  transposeParams->set("sort", false);
555  At = transposer.createTranspose(transposeParams);
556  }
557 
558  RCP<Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> AA = rcp(new Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(At));
559  RCP<XCrsMatrix> AAA = rcp_implicit_cast<XCrsMatrix>(AA);
560  RCP<XMatrix> AAAA = rcp(new XCrsMatrixWrap(AAA));
561 
562  if (Op.IsView("stridedMaps"))
563  AAAA->CreateView("stridedMaps", Teuchos::rcpFromRef(Op), true /*doTranspose*/);
564 
565  return AAAA;
566  } else {
567  throw Exceptions::RuntimeError("Utilities::Transpose failed, perhaps because matrix is not a Crs matrix");
568  }
569  } // if
570 
571  // Epetra case
572  std::cout << "Utilities::Transpose() not implemented for Epetra" << std::endl;
573  return Teuchos::null;
574 
575 } // Transpose
576 
577 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
578 RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
579 Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
580  RealValuedToScalarMultiVector(RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node>> X) {
581  RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Xscalar;
582 #if defined(HAVE_XPETRA_TPETRA) && (defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) || defined(HAVE_TPETRA_INST_COMPLEX_FLOAT))
583  using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
584 
585  // Need to cast the real-valued multivector to Scalar=complex
586  if ((typeid(Scalar).name() == typeid(std::complex<double>).name()) ||
587  (typeid(Scalar).name() == typeid(std::complex<float>).name())) {
588  size_t numVecs = X->getNumVectors();
589  Xscalar = Xpetra::MultiVectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(X->getMap(), numVecs);
590  auto XVec = X->getDeviceLocalView(Xpetra::Access::ReadOnly);
591  auto XVecScalar = Xscalar->getDeviceLocalView(Xpetra::Access::ReadWrite);
592 
593  Kokkos::parallel_for(
594  "MueLu:Utils::RealValuedToScalarMultiVector", range_type(0, X->getLocalLength()),
595  KOKKOS_LAMBDA(const size_t i) {
596  for (size_t j = 0; j < numVecs; j++)
597  XVecScalar(i, j) = XVec(i, j);
598  });
599  } else
600 #endif
601  Xscalar = rcp_dynamic_cast<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(X);
602  return Xscalar;
603 }
604 
605 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
606 RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LocalOrdinal, GlobalOrdinal, Node>>
607 Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
608  ExtractCoordinatesFromParameterList(ParameterList& paramList) {
609  RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LocalOrdinal, GlobalOrdinal, Node>> coordinates = Teuchos::null;
610 
611  // check whether coordinates are contained in parameter list
612  if (paramList.isParameter("Coordinates") == false)
613  return coordinates;
614 
615  // define Tpetra::MultiVector type with Scalar=float only if
616  // * ETI is turned off, since then the compiler will instantiate it automatically OR
617  // * Tpetra is instantiated on Scalar=float
618 #if !defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION) || defined(HAVE_TPETRA_INST_FLOAT)
619  typedef Tpetra::MultiVector<float, LocalOrdinal, GlobalOrdinal, Node> tfMV;
620  RCP<tfMV> floatCoords = Teuchos::null;
621 #endif
622 
623  // define Tpetra::MultiVector type with Scalar=double only if
624  // * ETI is turned off, since then the compiler will instantiate it automatically OR
625  // * Tpetra is instantiated on Scalar=double
626 #if !defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION) || defined(HAVE_TPETRA_INST_DOUBLE)
627  typedef Tpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LocalOrdinal, GlobalOrdinal, Node> tdMV;
628  RCP<tdMV> doubleCoords = Teuchos::null;
629  if (paramList.isType<RCP<tdMV>>("Coordinates")) {
630  // Coordinates are stored as a double vector
631  doubleCoords = paramList.get<RCP<tdMV>>("Coordinates");
632  paramList.remove("Coordinates");
633  }
634 #if !defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION) || defined(HAVE_TPETRA_INST_FLOAT)
635  else if (paramList.isType<RCP<tfMV>>("Coordinates")) {
636  // check if coordinates are stored as a float vector
637  floatCoords = paramList.get<RCP<tfMV>>("Coordinates");
638  paramList.remove("Coordinates");
639  doubleCoords = rcp(new tdMV(floatCoords->getMap(), floatCoords->getNumVectors()));
640  deep_copy(*doubleCoords, *floatCoords);
641  }
642 #endif
643  // We have the coordinates in a Tpetra double vector
644  if (doubleCoords != Teuchos::null) {
645  // rcp(new Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Vtpetra));
646  coordinates = Teuchos::rcp_dynamic_cast<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LocalOrdinal, GlobalOrdinal, Node>>(MueLu::TpetraMultiVector_To_XpetraMultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LocalOrdinal, GlobalOrdinal, Node>(doubleCoords));
647  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(coordinates));
648  TEUCHOS_TEST_FOR_EXCEPT(doubleCoords->getNumVectors() != coordinates->getNumVectors());
649  }
650 #else
651  // coordinates usually are stored as double vector
652  // Tpetra is not instantiated on scalar=double
653  throw Exceptions::RuntimeError("ExtractCoordinatesFromParameterList: The coordinates vector in parameter list is expected to be a Tpetra multivector with SC=double or float.");
654 #endif
655 
656  // check for Xpetra coordinates vector
657  if (paramList.isType<decltype(coordinates)>("Coordinates")) {
658  coordinates = paramList.get<decltype(coordinates)>("Coordinates");
659  }
660 
661  return coordinates;
662 } // ExtractCoordinatesFromParameterList
663 
664 } // namespace MueLu
665 
666 #define MUELU_UTILITIES_SHORT
667 #endif // MUELU_UTILITIES_DEF_HPP
668 
669 // LocalWords: LocalOrdinal
Xpetra::CrsMatrixWrap::getCrsMatrix
RCP< CrsMatrix > getCrsMatrix() const
Epetra_MultiVector
MueLu::Exceptions::BadCast
Exception indicating invalid cast attempted.
Definition: MueLu_Exceptions.hpp:21
Epetra_Map
Xpetra_Vector.hpp
MueLu::Utilities::Map2TpetraMap
static const RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > Map2TpetraMap(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &map)
Definition: MueLu_Utilities_def.hpp:382
LocalOrdinal
MueLu::DefaultLocalOrdinal LocalOrdinal
Definition: MueLu_UseDefaultTypes.hpp:13
Xpetra_TpetraOperator.hpp
MueLu::Utilities::MV2TpetraMV
static RCP< const Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2TpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node >> const vec)
Helper utility to pull out the underlying Tpetra objects from an Xpetra object.
Definition: MueLu_Utilities_def.hpp:177
EpetraExt_CrsMatrixIn.h
MueLu::Utilities::MyOldScaleMatrix
static void MyOldScaleMatrix(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Teuchos::ArrayRCP< const Scalar > &scalingVector, bool doInverse=true, bool doFillComplete=true, bool doOptimizeStorage=true)
Definition: MueLu_Utilities_def.hpp:390
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)
MueLu::Utilities::RealValuedToScalarMultiVector
static RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > RealValuedToScalarMultiVector(RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::coordinateType, LocalOrdinal, GlobalOrdinal, Node >> X)
Definition: MueLu_Utilities_def.hpp:580
MueLu::Utilities::Map2EpetraMap
static const Epetra_Map & Map2EpetraMap(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &map)
Definition: MueLu_Utilities_def.hpp:167
Tpetra::CrsMatrix::getRangeMap
Teuchos::RCP< const map_type > getRangeMap() const override
Xpetra::TpetraBlockCrsMatrix::getTpetra_BlockCrsMatrixNonConst
RCP< Tpetra::BlockCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetra_BlockCrsMatrixNonConst() const
Teuchos::ParameterList::get
T & get(const std::string &name, T def_value)
Xpetra::EpetraCrsMatrixT::getEpetra_CrsMatrixNonConst
RCP< Epetra_CrsMatrix > getEpetra_CrsMatrixNonConst() const
MueLu::Utilities::ExtractCoordinatesFromParameterList
static RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::magnitudeType, LocalOrdinal, GlobalOrdinal, Node > > ExtractCoordinatesFromParameterList(ParameterList &paramList)
Definition: MueLu_Utilities_def.hpp:608
MueLu::Utilities::Op2NonConstTpetraRow
static RCP< Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraRow(RCP< Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node >> Op)
Definition: MueLu_Utilities_def.hpp:352
Teuchos::is_null
bool is_null(const std::shared_ptr< T > &p)
MueLu::Utilities::Op2NonConstEpetraCrs
static RCP< Epetra_CrsMatrix > Op2NonConstEpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> Op)
Definition: MueLu_Utilities_def.hpp:126
Xpetra_TpetraRowMatrix.hpp
Tpetra::CrsMatrix::resumeFill
void resumeFill(const Teuchos::RCP< Teuchos::ParameterList > &params=Teuchos::null)
Xpetra_EpetraUtils.hpp
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
MueLu::Utilities::Op2TpetraCrs
static RCP< const Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2TpetraCrs(RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> Op)
Definition: MueLu_Utilities_def.hpp:212
EpetraExt_MatrixMatrix.h
Xpetra_Matrix.hpp
Teuchos::ArrayRCP::size
size_type size() const
Xpetra_MatrixMatrix.hpp
MueLu_ConfigDefs.hpp
Node
MueLu::DefaultNode Node
Definition: MueLu_UseDefaultTypes.hpp:15
MueLu::Utilities::Op2NonConstTpetraBlockCrs
static RCP< Tpetra::BlockCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraBlockCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> Op)
Definition: MueLu_Utilities_def.hpp:278
Xpetra_IO.hpp
MueLu_Utilities_decl.hpp
Teuchos::ParameterList::isParameter
bool isParameter(const std::string &name) const
EpetraExt_RowMatrixOut.h
deep_copy
void deep_copy(MultiVector< DS, DL, DG, DN > &dst, const MultiVector< SS, SL, SG, SN > &src)
Teuchos::ParameterList::remove
bool remove(std::string const &name, bool throwIfNotExists=true)
Xpetra::Matrix::fillComplete
virtual void fillComplete(const RCP< const Map > &domainMap, const RCP< const Map > &rangeMap, const RCP< ParameterList > &params=null)=0
Tpetra::CrsMatrix::getLocalMaxNumRowEntries
size_t getLocalMaxNumRowEntries() const override
Xpetra_MultiVector.hpp
Tpetra::CrsMatrix::isFillComplete
bool isFillComplete() const override
rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos_ParameterList.hpp
MueLu::Utilities::MV2NonConstEpetraMV
static RCP< Epetra_MultiVector > MV2NonConstEpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node >> vec)
Definition: MueLu_Utilities_def.hpp:95
Scalar
MueLu::DefaultScalar Scalar
Definition: MueLu_UseDefaultTypes.hpp:12
EpetraExt_BlockMapOut.h
GlobalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Definition: MueLu_UseDefaultTypes.hpp:14
Xpetra::TpetraBlockCrsMatrix::getTpetra_BlockCrsMatrix
RCP< const Tpetra::BlockCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetra_BlockCrsMatrix() const
EpetraExt_MultiVectorOut.h
MueLu::Utilities::MV2NonConstTpetraMV2
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV2(Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &vec)
Definition: MueLu_Utilities_def.hpp:199
MueLu::Utilities::Transpose
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Transpose(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, bool optimizeTranspose=false, const std::string &label=std::string(), const Teuchos::RCP< Teuchos::ParameterList > &params=Teuchos::null)
Definition: MueLu_Utilities_def.hpp:496
EpetraExt_MultiVectorIn.h
Xpetra::EpetraCrsMatrixT::getEpetra_CrsMatrix
RCP< const Epetra_CrsMatrix > getEpetra_CrsMatrix() const
Tpetra::CrsMatrix::replaceGlobalValues
local_ordinal_type replaceGlobalValues(const global_ordinal_type globalRow, const Kokkos::View< const global_ordinal_type *, Kokkos::AnonymousSpace > &inputInds, const Kokkos::View< const impl_scalar_type *, Kokkos::AnonymousSpace > &inputVals)
MueLu::Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap
RCP< Xpetra::CrsMatrixWrap< SC, LO, GO, NO > > Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap(RCP< Epetra_CrsMatrix > &epAB)
Definition: MueLu_Utilities_def.hpp:80
Tpetra::CrsMatrix::getNumEntriesInLocalRow
size_t getNumEntriesInLocalRow(local_ordinal_type localRow) const override
EpetraExt_BlockMapIn.h
Xpetra_EpetraMap.hpp
Tpetra::BlockCrsMatrixTransposer::createTranspose
Teuchos::RCP< bcrs_matrix_type > createTranspose(const Teuchos::RCP< Teuchos::ParameterList > &params=Teuchos::null)
Xpetra_BlockedCrsMatrix.hpp
Xpetra_MatrixFactory.hpp
Tpetra::RowMatrixTransposer::createTranspose
Teuchos::RCP< crs_matrix_type > createTranspose(const Teuchos::RCP< Teuchos::ParameterList > &params=Teuchos::null)
MueLu::Utilities::MyOldScaleMatrix_Epetra
static void MyOldScaleMatrix_Epetra(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Teuchos::ArrayRCP< Scalar > &scalingVector, bool doFillComplete, bool doOptimizeStorage)
Definition: MueLu_Utilities_def.hpp:418
Tpetra::CrsMatrix::getDomainMap
Teuchos::RCP< const map_type > getDomainMap() const override
Xpetra::Matrix::isLocallyIndexed
virtual bool isLocallyIndexed() const =0
Xpetra::TpetraCrsMatrix::getTpetra_CrsMatrix
RCP< const Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetra_CrsMatrix() const
Epetra_CrsMatrix
Xpetra_EpetraMultiVector.hpp
Xpetra::TpetraCrsMatrix::getTpetra_CrsMatrixNonConst
RCP< Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetra_CrsMatrixNonConst() const
Teuchos::ParameterList::isType
bool isType(const std::string &name) const
MueLu::Utilities::Op2TpetraRow
static RCP< const Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2TpetraRow(RCP< const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node >> Op)
Definition: MueLu_Utilities_def.hpp:322
Tpetra::CrsMatrix::getGlobalRowView
virtual void getGlobalRowView(GlobalOrdinal GlobalRow, global_inds_host_view_type &indices, values_host_view_type &values) const =0
MueLu::Utilities::MV2EpetraMV
static RCP< const Epetra_MultiVector > MV2EpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node >> const vec)
Helper utility to pull out the underlying Epetra objects from an Xpetra object.
Definition: MueLu_Utilities_def.hpp:87
Xpetra::EpetraMultiVectorT::getEpetra_MultiVector
RCP< Epetra_MultiVector > getEpetra_MultiVector() const
Xpetra::Operator::getRangeMap
virtual Teuchos::RCP< const Map > getRangeMap() const =0
Xpetra_VectorFactory.hpp
Xpetra::TpetraMultiVector::getTpetra_MultiVector
RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetra_MultiVector() const
Tpetra::CrsMatrix::getNumEntriesInGlobalRow
size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const override
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: MueLu_Exceptions.hpp:36
MatrixMarket_Tpetra.hpp
Tpetra::CrsMatrix::getLocalRowView
virtual void getLocalRowView(LocalOrdinal LocalRow, local_inds_host_view_type &indices, values_host_view_type &values) const =0
Tpetra::CrsMatrix::replaceLocalValues
local_ordinal_type replaceLocalValues(const local_ordinal_type localRow, const Kokkos::View< const local_ordinal_type *, Kokkos::AnonymousSpace > &inputInds, const Kokkos::View< const impl_scalar_type *, Kokkos::AnonymousSpace > &inputVals)
MueLu::Utilities::Op2TpetraBlockCrs
static RCP< const Tpetra::BlockCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2TpetraBlockCrs(RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> Op)
Definition: MueLu_Utilities_def.hpp:265
MueLu::Utilities::MyOldScaleMatrix_Tpetra
static void MyOldScaleMatrix_Tpetra(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Teuchos::ArrayRCP< Scalar > &scalingVector, bool doFillComplete, bool doOptimizeStorage)
Definition: MueLu_Utilities_def.hpp:423
Xpetra_Operator.hpp
MueLu::Utilities::Op2EpetraCrs
static RCP< const Epetra_CrsMatrix > Op2EpetraCrs(RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> Op)
Definition: MueLu_Utilities_def.hpp:115
TEUCHOS_TEST_FOR_EXCEPT
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
MueLu::Utilities::Op2NonConstTpetraCrs
static RCP< Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> Op)
Definition: MueLu_Utilities_def.hpp:224
Xpetra::Operator::getDomainMap
virtual Teuchos::RCP< const Map > getDomainMap() const =0
MueLu::Utilities::MV2NonConstTpetraMV
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node >> vec)
Definition: MueLu_Utilities_def.hpp:185
Tpetra::CrsMatrix::getRowMap
Teuchos::RCP< const map_type > getRowMap() const override
Teuchos::RCP::is_null
bool is_null() const