MueLu  Version of the Day
MueLu_Utilities_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_DEF_HPP
47 #define MUELU_UTILITIES_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>
64 #include <EpetraExt_CrsMatrixIn.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_MatrixFactory.hpp>
93 #include <Xpetra_MultiVector.hpp>
94 #include <Xpetra_MultiVectorFactory.hpp>
95 #include <Xpetra_Operator.hpp>
96 #include <Xpetra_Vector.hpp>
97 #include <Xpetra_VectorFactory.hpp>
98 
99 #include <Xpetra_MatrixMatrix.hpp>
100 
101 #include <MueLu_Utilities_decl.hpp>
102 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_ML)
103 #include <ml_operator.h>
104 #include <ml_epetra_utils.h>
105 #endif
106 
107 namespace MueLu {
108 
109 #ifdef HAVE_MUELU_EPETRA
110  //using Xpetra::EpetraCrsMatrix; // TODO: mv in Xpetra_UseShortNamesScalar
111  //using Xpetra::EpetraMultiVector;
112 #endif
113 
114 #ifdef HAVE_MUELU_EPETRA
115  //defined after Utilities class
116  template<typename SC,typename LO,typename GO,typename NO>
117  RCP<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> > Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap(RCP<Epetra_CrsMatrix> &epAB);
118 #endif
119 
120 #ifdef HAVE_MUELU_EPETRA
121  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
122  RCP<const Epetra_MultiVector> Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2EpetraMV(const RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > vec) {
123  RCP<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> > tmpVec = rcp_dynamic_cast<Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >(vec);
124  if (tmpVec == Teuchos::null)
125  throw Exceptions::BadCast("Cast from Xpetra::MultiVector to Xpetra::EpetraMultiVector failed");
126  return tmpVec->getEpetra_MultiVector();
127  }
128 
129  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
130  RCP<Epetra_MultiVector> Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2NonConstEpetraMV(RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > vec) {
131  RCP<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> > tmpVec = rcp_dynamic_cast<Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >(vec);
132  if (tmpVec == Teuchos::null)
133  throw Exceptions::BadCast("Cast from Xpetra::MultiVector to Xpetra::EpetraMultiVector failed");
134  return tmpVec->getEpetra_MultiVector();
135  }
136 
137  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
138  Epetra_MultiVector& Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2NonConstEpetraMV(Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &vec) {
139  const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> & tmpVec = dynamic_cast<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> &>(vec);
140  return *(tmpVec.getEpetra_MultiVector());
141  }
142 
143  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
144  const Epetra_MultiVector& Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2EpetraMV(const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& vec) {
145  const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> & tmpVec = dynamic_cast<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> &>(vec);
146  return *(tmpVec.getEpetra_MultiVector());
147  }
148 
149  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
150  RCP<const Epetra_CrsMatrix> Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2EpetraCrs(RCP<const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Op) {
151  RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op);
152  if (crsOp == Teuchos::null)
153  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
154  const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node>>& tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(crsOp->getCrsMatrix());
155  if (tmp_ECrsMtx == Teuchos::null)
156  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
157  return tmp_ECrsMtx->getEpetra_CrsMatrix();
158  }
159 
160  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
161  RCP<Epetra_CrsMatrix> Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2NonConstEpetraCrs(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Op) {
162  RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op);
163  if (crsOp == Teuchos::null)
164  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
165  const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> > &tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(crsOp->getCrsMatrix());
166  if (tmp_ECrsMtx == Teuchos::null)
167  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
168  return tmp_ECrsMtx->getEpetra_CrsMatrixNonConst();
169  }
170 
171  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
172  const Epetra_CrsMatrix& Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2EpetraCrs(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op) {
173  try {
174  const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>& crsOp = dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(Op);
175  try {
176  const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node>& tmp_ECrsMtx = dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node>&>(*crsOp.getCrsMatrix());
177  return *tmp_ECrsMtx.getEpetra_CrsMatrix();
178  } catch (std::bad_cast) {
179  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
180  }
181  } catch (std::bad_cast) {
182  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
183  }
184  }
185 
186  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
187  Epetra_CrsMatrix& Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2NonConstEpetraCrs(Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op) {
188  try {
189  Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>& crsOp = dynamic_cast<Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(Op);
190  try {
191  Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node>& tmp_ECrsMtx = dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node>&>(*crsOp.getCrsMatrix());
192  return *tmp_ECrsMtx.getEpetra_CrsMatrixNonConst();
193  } catch (std::bad_cast) {
194  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
195  }
196  } catch (std::bad_cast) {
197  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
198  }
199  }
200 
201  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
202  const Epetra_Map& Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Map2EpetraMap(const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& map) {
203  RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> > xeMap = rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(rcpFromRef(map));
204  if (xeMap == Teuchos::null)
205  throw Exceptions::BadCast("Utilities::Map2EpetraMap : Cast from Xpetra::Map to Xpetra::EpetraMap failed");
206  return xeMap->getEpetra_Map();
207  }
208 #endif
209 
210 #ifdef HAVE_MUELU_TPETRA
211  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
212  RCP<const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
213  Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2TpetraMV(RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > const vec) {
214  RCP<const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmpVec = rcp_dynamic_cast<Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(vec);
215  if (tmpVec == Teuchos::null)
216  throw Exceptions::BadCast("Cast from Xpetra::MultiVector to Xpetra::TpetraMultiVector failed");
217  return tmpVec->getTpetra_MultiVector();
218  }
219 
220  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
221  RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2NonConstTpetraMV(RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > vec) {
222  RCP<const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmpVec = rcp_dynamic_cast<Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(vec);
223  if (tmpVec == Teuchos::null)
224  throw Exceptions::BadCast("Cast from Xpetra::MultiVector to Xpetra::TpetraMultiVector failed");
225  return tmpVec->getTpetra_MultiVector();
226  }
227 
228  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
229  Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2NonConstTpetraMV(Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& vec) {
230  const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmpVec = dynamic_cast<const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(vec);
231  return *(tmpVec.getTpetra_MultiVector());
232  }
233 
234  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
235  RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2NonConstTpetraMV2(Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &vec) {
236  const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmpVec = dynamic_cast<const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(vec);
237  return tmpVec.getTpetra_MultiVector();
238  }
239 
240  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
241  const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>&
242  Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2TpetraMV(const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& vec) {
243  const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmpVec = dynamic_cast<const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(vec);
244  return *(tmpVec.getTpetra_MultiVector());
245  }
246 
247  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
248  RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2TpetraCrs(RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Op) {
249  // Get the underlying Tpetra Mtx
250  RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op);
251  if (crsOp == Teuchos::null)
252  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
253  const RCP<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsOp->getCrsMatrix());
254  if (tmp_ECrsMtx == Teuchos::null)
255  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
256  return tmp_ECrsMtx->getTpetra_CrsMatrix();
257  }
258 
259  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
260  RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2NonConstTpetraCrs(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Op) {
261  RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op);
262  if (crsOp == Teuchos::null)
263  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
264  const RCP<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsOp->getCrsMatrix());
265  if (tmp_ECrsMtx == Teuchos::null)
266  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
267  return tmp_ECrsMtx->getTpetra_CrsMatrixNonConst();
268  }
269 
270  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
271  const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2TpetraCrs(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op) {
272  try {
273  const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>& crsOp = dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(Op);
274  try {
275  const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmp_ECrsMtx = dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(*crsOp.getCrsMatrix());
276  return *tmp_ECrsMtx.getTpetra_CrsMatrix();
277  } catch (std::bad_cast) {
278  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
279  }
280  } catch (std::bad_cast) {
281  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
282  }
283  }
284 
285  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
286  Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2NonConstTpetraCrs(Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op) {
287  try {
288  Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>& crsOp = dynamic_cast<Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(Op);
289  try {
290  Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmp_ECrsMtx = dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(*crsOp.getCrsMatrix());
291  return *tmp_ECrsMtx.getTpetra_CrsMatrixNonConst();
292  } catch (std::bad_cast) {
293  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
294  }
295  } catch (std::bad_cast) {
296  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
297  }
298  }
299 
300  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
301  RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2TpetraRow(RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Op) {
302  RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op);
303  if (crsOp == Teuchos::null)
304  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
305 
306  RCP<const Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsMat = crsOp->getCrsMatrix();
307  const RCP<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmp_Crs = rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsMat);
308  RCP<const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmp_BlockCrs;
309  if(!tmp_Crs.is_null()) {
310  return tmp_Crs->getTpetra_CrsMatrixNonConst();
311  }
312  else {
313  tmp_BlockCrs= rcp_dynamic_cast<const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsMat);
314  if (tmp_BlockCrs.is_null())
315  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix and Xpetra::TpetraBlockCrsMatrix failed");
316  return tmp_BlockCrs->getTpetra_BlockCrsMatrixNonConst();
317  }
318  }
319 
320  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
321  RCP<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2NonConstTpetraRow(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Op) {
322  RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op);
323  if (crsOp == Teuchos::null)
324  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
325 
326  RCP<const Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsMat = crsOp->getCrsMatrix();
327  const RCP<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmp_Crs = rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsMat);
328  RCP<const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmp_BlockCrs;
329  if(!tmp_Crs.is_null()) {
330  return tmp_Crs->getTpetra_CrsMatrixNonConst();
331  }
332  else {
333  tmp_BlockCrs= rcp_dynamic_cast<const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsMat);
334  if (tmp_BlockCrs.is_null())
335  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix and Xpetra::TpetraBlockCrsMatrix failed");
336  return tmp_BlockCrs->getTpetra_BlockCrsMatrixNonConst();
337  }
338  }
339 
340 
341  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
342  const RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal,Node> > Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Map2TpetraMap(const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& map) {
343  const RCP<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node>>& tmp_TMap = rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >(rcpFromRef(map));
344  if (tmp_TMap == Teuchos::null)
345  throw Exceptions::BadCast("Utilities::Map2TpetraMap : Cast from Xpetra::Map to Xpetra::TpetraMap failed");
346  return tmp_TMap->getTpetra_Map();
347  }
348 #endif
349 
350  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
351  void Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix(Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op, const Teuchos::ArrayRCP<const Scalar>& scalingVector, bool doInverse,
352  bool doFillComplete,
353  bool doOptimizeStorage)
354  {
355  Scalar one = Teuchos::ScalarTraits<Scalar>::one();
356  Teuchos::ArrayRCP<Scalar> sv(scalingVector.size());
357  if (doInverse) {
358  for (int i = 0; i < scalingVector.size(); ++i)
359  sv[i] = one / scalingVector[i];
360  } else {
361  for (int i = 0; i < scalingVector.size(); ++i)
362  sv[i] = scalingVector[i];
363  }
364 
365  switch (Op.getRowMap()->lib()) {
366  case Xpetra::UseTpetra:
367  MyOldScaleMatrix_Tpetra(Op, sv, doFillComplete, doOptimizeStorage);
368  break;
369 
370  case Xpetra::UseEpetra:
371  MyOldScaleMatrix_Epetra(Op, sv, doFillComplete, doOptimizeStorage);
372  break;
373 
374  default:
375  throw Exceptions::RuntimeError("Only Epetra and Tpetra matrices can be scaled.");
376  break;
377  }
378  }
379 
380  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
381  void Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix_Epetra(Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op, const Teuchos::ArrayRCP<Scalar>& scalingVector, bool doFillComplete, bool doOptimizeStorage) {
382  throw Exceptions::RuntimeError("MyOldScaleMatrix_Epetra: Epetra needs SC=double and LO=GO=int.");
383  }
384 
385  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
386  void Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix_Tpetra(Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op, const Teuchos::ArrayRCP<Scalar>& scalingVector,
387  bool doFillComplete,
388  bool doOptimizeStorage)
389  {
390 #ifdef HAVE_MUELU_TPETRA
391  try {
392  Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tpOp = Op2NonConstTpetraCrs(Op);
393 
394  const RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = tpOp.getRowMap();
395  const RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domainMap = tpOp.getDomainMap();
396  const RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rangeMap = tpOp.getRangeMap();
397 
398  size_t maxRowSize = tpOp.getNodeMaxNumRowEntries();
399  if (maxRowSize == Teuchos::as<size_t>(-1)) // hasn't been determined yet
400  maxRowSize = 20;
401 
402  std::vector<Scalar> scaledVals(maxRowSize);
403  if (tpOp.isFillComplete())
404  tpOp.resumeFill();
405 
406  if (Op.isLocallyIndexed() == true) {
407  Teuchos::ArrayView<const LocalOrdinal> cols;
408  Teuchos::ArrayView<const Scalar> vals;
409 
410  for (size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
411  tpOp.getLocalRowView(i, cols, vals);
412  size_t nnz = tpOp.getNumEntriesInLocalRow(i);
413  if (nnz > maxRowSize) {
414  maxRowSize = nnz;
415  scaledVals.resize(maxRowSize);
416  }
417  for (size_t j = 0; j < nnz; ++j)
418  scaledVals[j] = vals[j]*scalingVector[i];
419 
420  if (nnz > 0) {
421  Teuchos::ArrayView<const Scalar> valview(&scaledVals[0], nnz);
422  tpOp.replaceLocalValues(i, cols, valview);
423  }
424  } //for (size_t i=0; ...
425 
426  } else {
427  Teuchos::ArrayView<const GlobalOrdinal> cols;
428  Teuchos::ArrayView<const Scalar> vals;
429 
430  for (size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
431  GlobalOrdinal gid = rowMap->getGlobalElement(i);
432  tpOp.getGlobalRowView(gid, cols, vals);
433  size_t nnz = tpOp.getNumEntriesInGlobalRow(gid);
434  if (nnz > maxRowSize) {
435  maxRowSize = nnz;
436  scaledVals.resize(maxRowSize);
437  }
438  // FIXME FIXME FIXME FIXME FIXME FIXME
439  for (size_t j = 0; j < nnz; ++j)
440  scaledVals[j] = vals[j]*scalingVector[i]; //FIXME i or gid?
441 
442  if (nnz > 0) {
443  Teuchos::ArrayView<const Scalar> valview(&scaledVals[0], nnz);
444  tpOp.replaceGlobalValues(gid, cols, valview);
445  }
446  } //for (size_t i=0; ...
447  }
448 
449  if (doFillComplete) {
450  if (domainMap == Teuchos::null || rangeMap == Teuchos::null)
451  throw Exceptions::RuntimeError("In Utilities::Scaling: cannot fillComplete because the domain and/or range map hasn't been defined");
452 
453  RCP<Teuchos::ParameterList> params = rcp(new Teuchos::ParameterList());
454  params->set("Optimize Storage", doOptimizeStorage);
455  params->set("No Nonlocal Changes", true);
456  Op.fillComplete(Op.getDomainMap(), Op.getRangeMap(), params);
457  }
458  } catch(...) {
459  throw Exceptions::RuntimeError("Only Tpetra::CrsMatrix types can be scaled (Err.1)");
460  }
461 #else
462  throw Exceptions::RuntimeError("Matrix scaling is not possible because Tpetra has not been enabled.");
463 #endif
464  } //MyOldScaleMatrix_Tpetra()
465 
466  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
467  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
469  Transpose (Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op, bool optimizeTranspose,const std::string & label) {
470 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_EPETRAEXT)
471  std::string TorE = "epetra";
472 #else
473  std::string TorE = "tpetra";
474 #endif
475 
476 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_EPETRAEXT)
477  try {
479  (void) epetraOp; // silence "unused variable" compiler warning
480  } catch (...) {
481  TorE = "tpetra";
482  }
483 #endif
484 
485 #ifdef HAVE_MUELU_TPETRA
486  if (TorE == "tpetra") {
487  try {
488  const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tpetraOp = Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2TpetraCrs(Op);
489 
490  RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A;
491  Tpetra::RowMatrixTransposer<Scalar, LocalOrdinal, GlobalOrdinal, Node> transposer(rcpFromRef(tpetraOp),label); //more than meets the eye
492  A = transposer.createTranspose();
493 
494  RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > AA = rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(A) );
495  RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > AAA = rcp_implicit_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(AA);
496  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > AAAA = rcp( new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(AAA) );
497  if (!AAAA->isFillComplete())
498  AAAA->fillComplete(Op.getRangeMap(),Op.getDomainMap());
499 
500  return AAAA;
501 
502  } catch (std::exception& e) {
503  std::cout << "threw exception '" << e.what() << "'" << std::endl;
504  throw Exceptions::RuntimeError("Utilities::Transpose failed, perhaps because matrix is not a Crs matrix");
505  }
506  } //if
507 #endif
508 
509  // Epetra case
510  std::cout << "Utilities::Transpose() not implemented for Epetra" << std::endl;
511  return Teuchos::null;
512 
513  } // Transpose
514 
515 
516 } //namespace MueLu
517 
518 #define MUELU_UTILITIES_SHORT
519 #endif // MUELU_UTILITIES_DEF_HPP
520 
521 // LocalWords: LocalOrdinal
Exception indicating invalid cast attempted.
static const RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > Map2TpetraMap(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &map)
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.
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)
static RCP< const Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2TpetraCrs(RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static const Epetra_Map & Map2EpetraMap(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &map)
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())
static RCP< Epetra_CrsMatrix > Op2NonConstEpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
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.
Namespace for MueLu classes and methods.
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
static RCP< Epetra_MultiVector > MV2NonConstEpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
static RCP< Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV2(Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &vec)
static RCP< const Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2TpetraRow(RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
RCP< Xpetra::CrsMatrixWrap< SC, LO, GO, NO > > Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap(RCP< Epetra_CrsMatrix > &epAB)
static void MyOldScaleMatrix_Epetra(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Teuchos::ArrayRCP< Scalar > &scalingVector, bool doFillComplete, bool doOptimizeStorage)
static RCP< const Epetra_CrsMatrix > Op2EpetraCrs(RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
Exception throws to report errors in the internal logical of the program.
static RCP< Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraRow(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static void MyOldScaleMatrix_Tpetra(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Teuchos::ArrayRCP< Scalar > &scalingVector, bool doFillComplete, bool doOptimizeStorage)