MueLu Version of the Day
Loading...
Searching...
No Matches
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)
67#include <Xpetra_EpetraUtils.hpp>
68#include <Xpetra_EpetraMultiVector.hpp>
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_TpetraOperator.hpp>
78#include <Xpetra_TpetraCrsMatrix.hpp>
79#include <Xpetra_TpetraBlockCrsMatrix.hpp>
80#endif
81
82#ifdef HAVE_MUELU_EPETRA
83#include <Xpetra_EpetraMap.hpp>
84#endif
85
86#include <Xpetra_BlockedCrsMatrix.hpp>
87//#include <Xpetra_DefaultPlatform.hpp>
88#include <Xpetra_IO.hpp>
89#include <Xpetra_Import.hpp>
90#include <Xpetra_ImportFactory.hpp>
91#include <Xpetra_Map.hpp>
92#include <Xpetra_MapFactory.hpp>
93#include <Xpetra_Matrix.hpp>
94#include <Xpetra_MatrixFactory.hpp>
95#include <Xpetra_MultiVector.hpp>
96#include <Xpetra_MultiVectorFactory.hpp>
97#include <Xpetra_Operator.hpp>
98#include <Xpetra_Vector.hpp>
99#include <Xpetra_VectorFactory.hpp>
100
101#include <Xpetra_MatrixMatrix.hpp>
102
104#if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_ML)
105#include <ml_operator.h>
106#include <ml_epetra_utils.h>
107#endif
108
109namespace MueLu {
110
111#ifdef HAVE_MUELU_EPETRA
112 //using Xpetra::EpetraCrsMatrix; // TODO: mv in Xpetra_UseShortNamesScalar
113 //using Xpetra::EpetraMultiVector;
114#endif
115
116#ifdef HAVE_MUELU_EPETRA
117 template<typename SC,typename LO,typename GO,typename NO>
118 RCP<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> > Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap(RCP<Epetra_CrsMatrix> &epAB){
119 return Xpetra::Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap<SC,LO,GO,NO>(epAB);
120 }
121#endif
122
123#ifdef HAVE_MUELU_EPETRA
124 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
125 RCP<const Epetra_MultiVector> Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2EpetraMV(const RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > vec) {
126 RCP<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> > tmpVec = rcp_dynamic_cast<Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >(vec);
127 if (tmpVec == Teuchos::null)
128 throw Exceptions::BadCast("Cast from Xpetra::MultiVector to Xpetra::EpetraMultiVector failed");
129 return tmpVec->getEpetra_MultiVector();
130 }
131
132 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
133 RCP<Epetra_MultiVector> Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2NonConstEpetraMV(RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > vec) {
134 RCP<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> > tmpVec = rcp_dynamic_cast<Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >(vec);
135 if (tmpVec == Teuchos::null)
136 throw Exceptions::BadCast("Cast from Xpetra::MultiVector to Xpetra::EpetraMultiVector failed");
137 return tmpVec->getEpetra_MultiVector();
138 }
139
140 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
141 Epetra_MultiVector& Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2NonConstEpetraMV(Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &vec) {
142 const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> & tmpVec = dynamic_cast<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> &>(vec);
143 return *(tmpVec.getEpetra_MultiVector());
144 }
145
146 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
147 const Epetra_MultiVector& Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2EpetraMV(const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& vec) {
148 const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> & tmpVec = dynamic_cast<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> &>(vec);
149 return *(tmpVec.getEpetra_MultiVector());
150 }
151
152 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
153 RCP<const Epetra_CrsMatrix> Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2EpetraCrs(RCP<const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Op) {
154 RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op);
155 if (crsOp == Teuchos::null)
156 throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
157 const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node>>& tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(crsOp->getCrsMatrix());
158 if (tmp_ECrsMtx == Teuchos::null)
159 throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
160 return tmp_ECrsMtx->getEpetra_CrsMatrix();
161 }
162
163 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
164 RCP<Epetra_CrsMatrix> Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2NonConstEpetraCrs(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Op) {
165 RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op);
166 if (crsOp == Teuchos::null)
167 throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
168 const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> > &tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(crsOp->getCrsMatrix());
169 if (tmp_ECrsMtx == Teuchos::null)
170 throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
171 return tmp_ECrsMtx->getEpetra_CrsMatrixNonConst();
172 }
173
174 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
175 const Epetra_CrsMatrix& Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2EpetraCrs(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op) {
176 try {
177 const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>& crsOp = dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(Op);
178 try {
179 const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node>& tmp_ECrsMtx = dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node>&>(*crsOp.getCrsMatrix());
180 return *tmp_ECrsMtx.getEpetra_CrsMatrix();
181 } catch (std::bad_cast&) {
182 throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
183 }
184 } catch (std::bad_cast&) {
185 throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
186 }
187 }
188
189 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
190 Epetra_CrsMatrix& Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2NonConstEpetraCrs(Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op) {
191 try {
192 Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>& crsOp = dynamic_cast<Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(Op);
193 try {
194 Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node>& tmp_ECrsMtx = dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node>&>(*crsOp.getCrsMatrix());
195 return *tmp_ECrsMtx.getEpetra_CrsMatrixNonConst();
196 } catch (std::bad_cast&) {
197 throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
198 }
199 } catch (std::bad_cast&) {
200 throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
201 }
202 }
203
204 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
205 const Epetra_Map& Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Map2EpetraMap(const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& map) {
206 RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> > xeMap = rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(rcpFromRef(map));
207 if (xeMap == Teuchos::null)
208 throw Exceptions::BadCast("Utilities::Map2EpetraMap : Cast from Xpetra::Map to Xpetra::EpetraMap failed");
209 return xeMap->getEpetra_Map();
210 }
211#endif
212
213#ifdef HAVE_MUELU_TPETRA
214 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
215 RCP<const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
216 Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2TpetraMV(RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > const vec) {
217 RCP<const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmpVec = rcp_dynamic_cast<Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(vec);
218 if (tmpVec == Teuchos::null)
219 throw Exceptions::BadCast("Cast from Xpetra::MultiVector to Xpetra::TpetraMultiVector failed");
220 return tmpVec->getTpetra_MultiVector();
221 }
222
223 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
224 RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2NonConstTpetraMV(RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > vec) {
225 RCP<const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmpVec = rcp_dynamic_cast<Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(vec);
226 if (tmpVec == Teuchos::null)
227 throw Exceptions::BadCast("Cast from Xpetra::MultiVector to Xpetra::TpetraMultiVector failed");
228 return tmpVec->getTpetra_MultiVector();
229 }
230
231 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
232 Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2NonConstTpetraMV(Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& vec) {
233 const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmpVec = dynamic_cast<const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(vec);
234 return *(tmpVec.getTpetra_MultiVector());
235 }
236
237 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
238 RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2NonConstTpetraMV2(Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &vec) {
239 const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmpVec = dynamic_cast<const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(vec);
240 return tmpVec.getTpetra_MultiVector();
241 }
242
243 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
244 const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>&
245 Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MV2TpetraMV(const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& vec) {
246 const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmpVec = dynamic_cast<const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(vec);
247 return *(tmpVec.getTpetra_MultiVector());
248 }
249
250 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
251 RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2TpetraCrs(RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Op) {
252 // Get the underlying Tpetra Mtx
253 RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op);
254 if (crsOp == Teuchos::null)
255 throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
256 const RCP<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsOp->getCrsMatrix());
257 if (tmp_ECrsMtx == Teuchos::null)
258 throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
259 return tmp_ECrsMtx->getTpetra_CrsMatrix();
260 }
261
262 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
263 RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2NonConstTpetraCrs(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Op) {
264 RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op);
265 if (crsOp == Teuchos::null)
266 throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
267 const RCP<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsOp->getCrsMatrix());
268 if (tmp_ECrsMtx == Teuchos::null)
269 throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
270 return tmp_ECrsMtx->getTpetra_CrsMatrixNonConst();
271 }
272
273 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
274 const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2TpetraCrs(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op) {
275 try {
276 const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>& crsOp = dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(Op);
277 try {
278 const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmp_ECrsMtx = dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(*crsOp.getCrsMatrix());
279 return *tmp_ECrsMtx.getTpetra_CrsMatrix();
280 } catch (std::bad_cast&) {
281 throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
282 }
283 } catch (std::bad_cast&) {
284 throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
285 }
286 }
287
288 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
289 Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2NonConstTpetraCrs(Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op) {
290 try {
291 Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>& crsOp = dynamic_cast<Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(Op);
292 try {
293 Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmp_ECrsMtx = dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(*crsOp.getCrsMatrix());
294 return *tmp_ECrsMtx.getTpetra_CrsMatrixNonConst();
295 } catch (std::bad_cast&) {
296 throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
297 }
298 } catch (std::bad_cast&) {
299 throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
300 }
301 }
302
303
304 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
305 RCP<const Tpetra::BlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2TpetraBlockCrs(RCP<const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Op) {
306 using XCrsMatrixWrap = Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
307 // Get the underlying Tpetra Mtx
308 RCP<const XCrsMatrixWrap> crsOp = rcp_dynamic_cast<const XCrsMatrixWrap>(Op);
309 if (crsOp == Teuchos::null)
310 throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
311 const RCP<const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsOp->getCrsMatrix());
312 if (tmp_ECrsMtx == Teuchos::null)
313 throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraBlockCrsMatrix failed");
314 return tmp_ECrsMtx->getTpetra_BlockCrsMatrix();
315 }
316
317 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
318 RCP< Tpetra::BlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2NonConstTpetraBlockCrs(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Op){
319 using XCrsMatrixWrap = Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
320 RCP<const XCrsMatrixWrap> crsOp = rcp_dynamic_cast<const XCrsMatrixWrap>(Op);
321 if (crsOp == Teuchos::null)
322 throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
323 const RCP<const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsOp->getCrsMatrix());
324 if (tmp_ECrsMtx == Teuchos::null)
325 throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraBlockCrsMatrix failed");
326 return tmp_ECrsMtx->getTpetra_BlockCrsMatrixNonConst();
327 }
328
329 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
330 const Tpetra::BlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2TpetraBlockCrs(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op) {
331 try {
332 using XCrsMatrixWrap = Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
333 const XCrsMatrixWrap& crsOp = dynamic_cast<const XCrsMatrixWrap&>(Op);
334 try {
335 const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmp_ECrsMtx = dynamic_cast<const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(*crsOp.getCrsMatrix());
336 return *tmp_ECrsMtx.getTpetra_BlockCrsMatrix();
337 } catch (std::bad_cast&) {
338 throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraBlockCrsMatrix failed");
339 }
340 } catch (std::bad_cast&) {
341 throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
342 }
343 }
344
345 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
346 Tpetra::BlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2NonConstTpetraBlockCrs(Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op) {
347 try {
348 using XCrsMatrixWrap = Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
349 XCrsMatrixWrap& crsOp = dynamic_cast<XCrsMatrixWrap&>(Op);
350 try {
351 Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmp_ECrsMtx = dynamic_cast<Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(*crsOp.getCrsMatrix());
352 return *tmp_ECrsMtx.getTpetra_BlockCrsMatrixNonConst();
353 } catch (std::bad_cast&) {
354 throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraBlockCrsMatrix failed");
355 }
356 } catch (std::bad_cast&) {
357 throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
358 }
359 }
360
361
362 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
363 RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2TpetraRow(RCP<const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Op) {
364 RCP<const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > mat = rcp_dynamic_cast<const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op);
365 if (!mat.is_null()) {
366 RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(mat);
367 if (crsOp == Teuchos::null)
368 throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
369
370 RCP<const Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsMat = crsOp->getCrsMatrix();
371 const RCP<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmp_Crs = rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsMat);
372 RCP<const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmp_BlockCrs;
373 if(!tmp_Crs.is_null()) {
374 return tmp_Crs->getTpetra_CrsMatrixNonConst();
375 }
376 else {
377 tmp_BlockCrs= rcp_dynamic_cast<const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsMat);
378 if (tmp_BlockCrs.is_null())
379 throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix and Xpetra::TpetraBlockCrsMatrix failed");
380 return tmp_BlockCrs->getTpetra_BlockCrsMatrixNonConst();
381 }
382 } else {
383 RCP<const Xpetra::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tpOp = rcp_dynamic_cast<const Xpetra::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op, true);
384 RCP<const Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>> tOp = tpOp->getOperatorConst();
385 RCP<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tRow = rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(tOp, true);
386 return tRow;
387 }
388 }
389
390 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
391 RCP<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2NonConstTpetraRow(RCP<Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Op) {
392 RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > mat = rcp_dynamic_cast<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op);
393 if (!mat.is_null()) {
394 RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(mat);
395 if (crsOp == Teuchos::null)
396 throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
397
398 RCP<const Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsMat = crsOp->getCrsMatrix();
399 const RCP<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmp_Crs = rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsMat);
400 RCP<const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmp_BlockCrs;
401 if(!tmp_Crs.is_null()) {
402 return tmp_Crs->getTpetra_CrsMatrixNonConst();
403 }
404 else {
405 tmp_BlockCrs= rcp_dynamic_cast<const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsMat);
406 if (tmp_BlockCrs.is_null())
407 throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix and Xpetra::TpetraBlockCrsMatrix failed");
408 return tmp_BlockCrs->getTpetra_BlockCrsMatrixNonConst();
409 }
410 } else {
411 RCP<Xpetra::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tpOp = rcp_dynamic_cast<Xpetra::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op, true);
412 RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>> tOp = tpOp->getOperator();
413 RCP<Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tRow = rcp_dynamic_cast<Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(tOp, true);
414 return tRow;
415 }
416 }
417
418
419 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
420 const RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal,Node> > Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Map2TpetraMap(const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& map) {
421 const RCP<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node>>& tmp_TMap = rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >(rcpFromRef(map));
422 if (tmp_TMap == Teuchos::null)
423 throw Exceptions::BadCast("Utilities::Map2TpetraMap : Cast from Xpetra::Map to Xpetra::TpetraMap failed");
424 return tmp_TMap->getTpetra_Map();
425 }
426#endif
427
428 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
429 void Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix(Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op, const Teuchos::ArrayRCP<const Scalar>& scalingVector, bool doInverse,
430 bool doFillComplete,
431 bool doOptimizeStorage)
432 {
433 Scalar one = Teuchos::ScalarTraits<Scalar>::one();
434 Teuchos::ArrayRCP<Scalar> sv(scalingVector.size());
435 if (doInverse) {
436 for (int i = 0; i < scalingVector.size(); ++i)
437 sv[i] = one / scalingVector[i];
438 } else {
439 for (int i = 0; i < scalingVector.size(); ++i)
440 sv[i] = scalingVector[i];
441 }
442
443 switch (Op.getRowMap()->lib()) {
444 case Xpetra::UseTpetra:
445 MyOldScaleMatrix_Tpetra(Op, sv, doFillComplete, doOptimizeStorage);
446 break;
447
448 case Xpetra::UseEpetra:
449 MyOldScaleMatrix_Epetra(Op, sv, doFillComplete, doOptimizeStorage);
450 break;
451
452 default:
453 throw Exceptions::RuntimeError("Only Epetra and Tpetra matrices can be scaled.");
454 }
455 }
456
457 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
458 void Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix_Epetra(Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& /* Op */, const Teuchos::ArrayRCP<Scalar>& /* scalingVector */, bool /* doFillComplete */, bool /* doOptimizeStorage */) {
459 throw Exceptions::RuntimeError("MyOldScaleMatrix_Epetra: Epetra needs SC=double and LO=GO=int.");
460 }
461
462 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
463 void Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix_Tpetra(Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op, const Teuchos::ArrayRCP<Scalar>& scalingVector,
464 bool doFillComplete,
465 bool doOptimizeStorage)
466 {
467#ifdef HAVE_MUELU_TPETRA
468 try {
469 Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tpOp = Op2NonConstTpetraCrs(Op);
470
471 const RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = tpOp.getRowMap();
472 const RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domainMap = tpOp.getDomainMap();
473 const RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rangeMap = tpOp.getRangeMap();
474
475 size_t maxRowSize = tpOp.getLocalMaxNumRowEntries();
476 if (maxRowSize == Teuchos::as<size_t>(-1)) // hasn't been determined yet
477 maxRowSize = 20;
478
479
480 if (tpOp.isFillComplete())
481 tpOp.resumeFill();
482
483 if (Op.isLocallyIndexed() == true) {
484 typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_inds_host_view_type cols;
485 typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::values_host_view_type vals;
486
487 for (size_t i = 0; i < rowMap->getLocalNumElements(); ++i) {
488 tpOp.getLocalRowView(i, cols, vals);
489
490 size_t nnz = tpOp.getNumEntriesInLocalRow(i);
491 typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::nonconst_values_host_view_type scaledVals("ScaledVals",nnz);
492
493 for (size_t j = 0; j < nnz; ++j) {
494 scaledVals[j] = scalingVector[i]*vals[j];
495 }
496
497 if (nnz > 0) {
498 tpOp.replaceLocalValues(i, cols, scaledVals);
499 }
500 } //for (size_t i=0; ...
501
502 } else {
503 typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::global_inds_host_view_type cols;
504 typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::values_host_view_type vals;
505
506 for (size_t i = 0; i < rowMap->getLocalNumElements(); ++i) {
507 GlobalOrdinal gid = rowMap->getGlobalElement(i);
508 tpOp.getGlobalRowView(gid, cols, vals);
509 size_t nnz = tpOp.getNumEntriesInGlobalRow(gid);
510 typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::nonconst_values_host_view_type scaledVals("ScaledVals",nnz);
511
512 // FIXME FIXME FIXME FIXME FIXME FIXME
513 for (size_t j = 0; j < nnz; ++j)
514 scaledVals[j] = scalingVector[i]*vals[j]; //FIXME i or gid?
515
516 if (nnz > 0) {
517 tpOp.replaceGlobalValues(gid, cols, scaledVals);
518 }
519 } //for (size_t i=0; ...
520 }
521
522 if (doFillComplete) {
523 if (domainMap == Teuchos::null || rangeMap == Teuchos::null)
524 throw Exceptions::RuntimeError("In Utilities::Scaling: cannot fillComplete because the domain and/or range map hasn't been defined");
525
526 RCP<Teuchos::ParameterList> params = rcp(new Teuchos::ParameterList());
527 params->set("Optimize Storage", doOptimizeStorage);
528 params->set("No Nonlocal Changes", true);
529 Op.fillComplete(Op.getDomainMap(), Op.getRangeMap(), params);
530 }
531 } catch(...) {
532 throw Exceptions::RuntimeError("Only Tpetra::CrsMatrix types can be scaled (Err.1)");
533 }
534#else
535 throw Exceptions::RuntimeError("Matrix scaling is not possible because Tpetra has not been enabled.");
536#endif
537 } //MyOldScaleMatrix_Tpetra()
538
539 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
540 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
542 Transpose (Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op, bool /* optimizeTranspose */,const std::string & label,const Teuchos::RCP<Teuchos::ParameterList> &params) {
543#if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_EPETRAEXT)
544 std::string TorE = "epetra";
545#else
546 std::string TorE = "tpetra";
547#endif
548
549#if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_EPETRAEXT)
550 try {
552 (void) epetraOp; // silence "unused variable" compiler warning
553 } catch (...) {
554 TorE = "tpetra";
555 }
556#endif
557
558#ifdef HAVE_MUELU_TPETRA
559 if (TorE == "tpetra") {
560 using Helpers = Xpetra::Helpers<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
561 /***************************************************************/
562 if(Helpers::isTpetraCrs(Op)) {
563 const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tpetraOp = Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2TpetraCrs(Op);
564
565 RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A;
566 Tpetra::RowMatrixTransposer<Scalar, LocalOrdinal, GlobalOrdinal, Node> transposer(rcpFromRef(tpetraOp),label); //more than meets the eye
567
568 {
569 using Teuchos::ParameterList;
570 using Teuchos::rcp;
571 RCP<ParameterList> transposeParams = params.is_null () ?
572 rcp (new ParameterList) :
573 rcp (new ParameterList (*params));
574 transposeParams->set ("sort", false);
575 A = transposer.createTranspose (transposeParams);
576 }
577
578 RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > AA = rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(A) );
579 RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > AAA = rcp_implicit_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(AA);
580 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > AAAA = rcp( new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(AAA) );
581 if (!AAAA->isFillComplete())
582 AAAA->fillComplete(Op.getRangeMap(), Op.getDomainMap());
583
584 if (Op.IsView("stridedMaps"))
585 AAAA->CreateView("stridedMaps", Teuchos::rcpFromRef(Op), true/*doTranspose*/);
586
587 return AAAA;
588 }
589 else if(Helpers::isTpetraBlockCrs(Op)) {
590 using XMatrix = Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
591 using XCrsMatrix = Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
592 using XCrsMatrixWrap = Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
593 using BCRS = Tpetra::BlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
594 using CRS = Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
596
597 RCP<BCRS> At;
598 {
599 Tpetra::BlockCrsMatrixTransposer<Scalar, LocalOrdinal, GlobalOrdinal, Node> transposer(rcpFromRef(tpetraOp),label);
600
601 using Teuchos::ParameterList;
602 using Teuchos::rcp;
603 RCP<ParameterList> transposeParams = params.is_null () ?
604 rcp (new ParameterList) :
605 rcp (new ParameterList (*params));
606 transposeParams->set ("sort", false);
607 At = transposer.createTranspose(transposeParams);
608 }
609
610 RCP<Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > AA = rcp(new Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(At));
611 RCP<XCrsMatrix> AAA = rcp_implicit_cast<XCrsMatrix>(AA);
612 RCP<XMatrix> AAAA = rcp( new XCrsMatrixWrap(AAA));
613
614 if (Op.IsView("stridedMaps"))
615 AAAA->CreateView("stridedMaps", Teuchos::rcpFromRef(Op), true/*doTranspose*/);
616
617 return AAAA;
618 } else {
619 throw Exceptions::RuntimeError("Utilities::Transpose failed, perhaps because matrix is not a Crs matrix");
620 }
621 } //if
622#endif
623
624 // Epetra case
625 std::cout << "Utilities::Transpose() not implemented for Epetra" << std::endl;
626 return Teuchos::null;
627
628 } // Transpose
629
630
631 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
632 RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
634 RealValuedToScalarMultiVector(RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LocalOrdinal,GlobalOrdinal,Node> > X) {
635 RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Xscalar;
636#if defined(HAVE_XPETRA_TPETRA) && (defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) || defined(HAVE_TPETRA_INST_COMPLEX_FLOAT))
637 typedef typename Teuchos::ScalarTraits<Scalar>::coordinateType real_type;
638 // Need to cast the real-valued multivector to Scalar=complex
639 if ((typeid(Scalar).name() == typeid(std::complex<double>).name()) ||
640 (typeid(Scalar).name() == typeid(std::complex<float>).name())) {
641 Xscalar = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(X->getMap(),X->getNumVectors());
642 size_t numVecs = X->getNumVectors();
643 for (size_t j=0;j<numVecs;j++) {
644 Teuchos::ArrayRCP<const real_type> XVec = X->getData(j);
645 Teuchos::ArrayRCP<Scalar> XVecScalar = Xscalar->getDataNonConst(j);
646 for(size_t i = 0; i < static_cast<size_t>(XVec.size()); ++i)
647 XVecScalar[i]=XVec[i];
648 }
649 } else
650#endif
651 Xscalar = rcp_dynamic_cast<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(X);
652 return Xscalar;
653 }
654
655
656 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
657 RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LocalOrdinal,GlobalOrdinal,Node> >
659 ExtractCoordinatesFromParameterList (ParameterList& paramList) {
660
661 RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LocalOrdinal,GlobalOrdinal,Node> > coordinates = Teuchos::null;
662
663 // check whether coordinates are contained in parameter list
664 if(paramList.isParameter ("Coordinates") == false)
665 return coordinates;
666
667#if defined(HAVE_MUELU_TPETRA)
668 // only Tpetra code
669
670 // define Tpetra::MultiVector type with Scalar=float only if
671 // * ETI is turned off, since then the compiler will instantiate it automatically OR
672 // * Tpetra is instantiated on Scalar=float
673#if !defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION) || defined(HAVE_TPETRA_INST_FLOAT)
674 typedef Tpetra::MultiVector<float, LocalOrdinal, GlobalOrdinal, Node> tfMV;
675 RCP<tfMV> floatCoords = Teuchos::null;
676#endif
677
678 // define Tpetra::MultiVector type with Scalar=double only if
679 // * ETI is turned off, since then the compiler will instantiate it automatically OR
680 // * Tpetra is instantiated on Scalar=double
681#if !defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION) || defined(HAVE_TPETRA_INST_DOUBLE)
682 typedef Tpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LocalOrdinal, GlobalOrdinal, Node> tdMV;
683 RCP<tdMV> doubleCoords = Teuchos::null;
684 if (paramList.isType<RCP<tdMV> >("Coordinates")) {
685 // Coordinates are stored as a double vector
686 doubleCoords = paramList.get<RCP<tdMV> >("Coordinates");
687 paramList.remove("Coordinates");
688 }
689#if !defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION) || defined(HAVE_TPETRA_INST_FLOAT)
690 else if (paramList.isType<RCP<tfMV> >("Coordinates")) {
691 // check if coordinates are stored as a float vector
692 floatCoords = paramList.get<RCP<tfMV> >("Coordinates");
693 paramList.remove("Coordinates");
694 doubleCoords = rcp(new tdMV(floatCoords->getMap(), floatCoords->getNumVectors()));
695 deep_copy(*doubleCoords, *floatCoords);
696 }
697#endif
698 // We have the coordinates in a Tpetra double vector
699 if(doubleCoords != Teuchos::null) {
700 //rcp(new Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Vtpetra));
701 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));
702 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(coordinates));
703 TEUCHOS_TEST_FOR_EXCEPT(doubleCoords->getNumVectors() != coordinates->getNumVectors());
704 }
705#else
706 // coordinates usually are stored as double vector
707 // Tpetra is not instantiated on scalar=double
708 throw Exceptions::RuntimeError("ExtractCoordinatesFromParameterList: The coordinates vector in parameter list is expected to be a Tpetra multivector with SC=double or float.");
709#endif
710#endif // endif HAVE_TPETRA
711
712 // check for Xpetra coordinates vector
713 if(paramList.isType<decltype(coordinates)>("Coordinates")) {
714 coordinates = paramList.get<decltype(coordinates)>("Coordinates");
715 }
716
717 return coordinates;
718 } // ExtractCoordinatesFromParameterList
719
720} //namespace MueLu
721
722#define MUELU_UTILITIES_SHORT
723#endif // MUELU_UTILITIES_DEF_HPP
724
725// LocalWords: LocalOrdinal
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
MueLu::DefaultNode Node
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
static RCP< Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< Tpetra::BlockCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraBlockCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraRow(RCP< Xpetra::Operator< 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)
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)
static RCP< const Epetra_CrsMatrix > Op2EpetraCrs(RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV2(Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &vec)
static const Epetra_Map & Map2EpetraMap(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &map)
static RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > RealValuedToScalarMultiVector(RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::coordinateType, LocalOrdinal, GlobalOrdinal, Node > > X)
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 void MyOldScaleMatrix_Epetra(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Teuchos::ArrayRCP< Scalar > &scalingVector, bool doFillComplete, bool doOptimizeStorage)
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< Epetra_CrsMatrix > Op2NonConstEpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< const Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2TpetraCrs(RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< const Tpetra::BlockCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2TpetraBlockCrs(RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< const Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2TpetraRow(RCP< const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::magnitudeType, LocalOrdinal, GlobalOrdinal, Node > > ExtractCoordinatesFromParameterList(ParameterList &paramList)
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.
RCP< Xpetra::CrsMatrixWrap< SC, LO, GO, NO > > Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap(RCP< Epetra_CrsMatrix > &epAB)