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>
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>
63 #include <EpetraExt_MultiVectorOut.h>
64 #include <EpetraExt_CrsMatrixIn.h>
65 #include <EpetraExt_MultiVectorIn.h>
66 #include <EpetraExt_BlockMapIn.h>
67 #include <Xpetra_EpetraUtils.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>
79 #endif
80 
81 #ifdef HAVE_MUELU_EPETRA
82 #include <Xpetra_EpetraMap.hpp>
83 #endif
84 
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>
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>
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>
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>
140  return *(tmpVec.getEpetra_MultiVector());
141  }
142 
143  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
146  return *(tmpVec.getEpetra_MultiVector());
147  }
148 
149  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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>
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>
173  try {
175  try {
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>
188  try {
190  try {
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>
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>
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>
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>
231  return *(tmpVec.getTpetra_MultiVector());
232  }
233 
234  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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>&
244  return *(tmpVec.getTpetra_MultiVector());
245  }
246 
247  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
249  // Get the underlying Tpetra Mtx
251  if (crsOp == Teuchos::null)
252  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
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>
262  if (crsOp == Teuchos::null)
263  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
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 {
274  try {
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>
287  try {
289  try {
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>
303  if (crsOp == Teuchos::null)
304  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
305 
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>
323  if (crsOp == Teuchos::null)
324  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
325 
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>
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>
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>
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>
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) {
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 {
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 
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>
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 
491  Tpetra::RowMatrixTransposer<Scalar, LocalOrdinal, GlobalOrdinal, Node> transposer(rcpFromRef(tpetraOp),label); //more than meets the eye
492  A = transposer.createTranspose();
493 
497  if (!AAAA->isFillComplete())
498  AAAA->fillComplete(Op.getRangeMap(), Op.getDomainMap());
499 
500  if (Op.IsView("stridedMaps"))
501  AAAA->CreateView("stridedMaps", Teuchos::rcpFromRef(Op), true/*doTranspose*/);
502 
503  return AAAA;
504 
505  } catch (std::exception& e) {
506  std::cout << "threw exception '" << e.what() << "'" << std::endl;
507  throw Exceptions::RuntimeError("Utilities::Transpose failed, perhaps because matrix is not a Crs matrix");
508  }
509  } //if
510 #endif
511 
512  // Epetra case
513  std::cout << "Utilities::Transpose() not implemented for Epetra" << std::endl;
514  return Teuchos::null;
515 
516  } // Transpose
517 
518 
519  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
523 
525 
526  // check whether coordinates are contained in parameter list
527  if(paramList.isParameter ("Coordinates") == false)
528  return coordinates;
529 
530 #if defined(HAVE_MUELU_TPETRA)
531  // only Tpetra code
532 
533  // define Tpetra::MultiVector type with Scalar=float only if
534  // * ETI is turned off, since then the compiler will instantiate it automatically OR
535  // * Tpetra is instantiated on Scalar=float
536 #if !defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION) || defined(HAVE_TPETRA_INST_FLOAT)
537  typedef Tpetra::MultiVector<float, LocalOrdinal, GlobalOrdinal, Node> tfMV;
538  RCP<tfMV> floatCoords = Teuchos::null;
539 #endif
540 
541  // define Tpetra::MultiVector type with Scalar=double only if
542  // * ETI is turned off, since then the compiler will instantiate it automatically OR
543  // * Tpetra is instantiated on Scalar=double
544 #if !defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION) || defined(HAVE_TPETRA_INST_DOUBLE)
545  typedef Tpetra::MultiVector<double, LocalOrdinal, GlobalOrdinal, Node> tdMV;
546  RCP<tdMV> doubleCoords = Teuchos::null;
547  if (paramList.isType<RCP<tdMV> >("Coordinates")) {
548  // Coordinates are stored as a double vector
549  doubleCoords = paramList.get<RCP<tdMV> >("Coordinates");
550  paramList.remove("Coordinates");
551  }
552 #if !defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION) || defined(HAVE_TPETRA_INST_FLOAT)
553  else if (paramList.isType<RCP<tfMV> >("Coordinates")) {
554  // check if coordinates are stored as a float vector
555  floatCoords = paramList.get<RCP<tfMV> >("Coordinates");
556  paramList.remove("Coordinates");
557  doubleCoords = rcp(new tdMV(floatCoords->getMap(), floatCoords->getNumVectors()));
558  deep_copy(*doubleCoords, *floatCoords);
559  }
560 #endif
561  // We have the coordinates in a Tpetra double vector
562  if(doubleCoords != Teuchos::null) {
563  //rcp(new Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Vtpetra));
564  coordinates = Teuchos::rcp_dynamic_cast<Xpetra::MultiVector<double,LocalOrdinal,GlobalOrdinal,Node> >(MueLu::TpetraMultiVector_To_XpetraMultiVector<double,LocalOrdinal,GlobalOrdinal,Node>(doubleCoords));
566  TEUCHOS_TEST_FOR_EXCEPT(doubleCoords->getNumVectors() != coordinates->getNumVectors());
567  }
568 #else
569  // coordinates usually are stored as double vector
570  // Tpetra is not instantiated on scalar=double
571  throw Exceptions::RuntimeError("ExtractCoordinatesFromParameterList: The coordinates vector in parameter list is expected to be a Tpetra multivector with SC=double or float.");
572 #endif
573 #endif // endif HAVE_TPETRA
574  return coordinates;
575  } // ExtractCoordinatesFromParameterList
576 
577 
578 } //namespace MueLu
579 
580 #define MUELU_UTILITIES_SHORT
581 #endif // MUELU_UTILITIES_DEF_HPP
582 
583 // LocalWords: LocalOrdinal
Exception indicating invalid cast attempted.
static const RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > Map2TpetraMap(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &map)
RCP< Epetra_CrsMatrix > getEpetra_CrsMatrixNonConst() const
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())
virtual const RCP< const Map > & getRowMap() const
static RCP< Epetra_CrsMatrix > Op2NonConstEpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
size_type size() const
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
bool is_null(const std::shared_ptr< T > &p)
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< Xpetra::MultiVector< double, LocalOrdinal, GlobalOrdinal, Node > > ExtractCoordinatesFromParameterList(ParameterList &paramList)
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
bool is_null() const
bool remove(std::string const &name, bool throwIfNotExists=true)
virtual void fillComplete(const RCP< const Map > &domainMap, const RCP< const Map > &rangeMap, const RCP< ParameterList > &params=null)=0
RCP< Epetra_MultiVector > getEpetra_MultiVector() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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)
RCP< Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetra_CrsMatrixNonConst() const
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV2(Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &vec)
RCP< CrsMatrix > getCrsMatrix() const
static RCP< const Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2TpetraRow(RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
bool isType(const std::string &name) const
RCP< Xpetra::CrsMatrixWrap< SC, LO, GO, NO > > Convert_Epetra_CrsMatrix_ToXpetra_CrsMatrixWrap(RCP< Epetra_CrsMatrix > &epAB)
RCP< const Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetra_CrsMatrix() const
static void MyOldScaleMatrix_Epetra(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Teuchos::ArrayRCP< Scalar > &scalingVector, bool doFillComplete, bool doOptimizeStorage)
RCP< const Epetra_CrsMatrix > getEpetra_CrsMatrix() const
static RCP< const Epetra_CrsMatrix > Op2EpetraCrs(RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
virtual bool isLocallyIndexed() const=0
bool isParameter(const std::string &name) const
virtual Teuchos::RCP< const Map > getRangeMap() const=0
bool IsView(const viewLabel_t viewLabel) const
Exception throws to report errors in the internal logical of the program.
RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetra_MultiVector() const
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)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
virtual Teuchos::RCP< const Map > getDomainMap() const=0