MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Thyra_MueLuTpetraQ2Q1PreconditionerFactory_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 THYRA_MUELU_TPETRA_Q2Q1PRECONDITIONER_FACTORY_DEF_HPP
47 #define THYRA_MUELU_TPETRA_Q2Q1PRECONDITIONER_FACTORY_DEF_HPP
48 
49 #ifdef HAVE_MUELU_EXPERIMENTAL
50 
52 
53 #include <Thyra_DefaultPreconditioner.hpp>
54 #include <Thyra_TpetraLinearOp.hpp>
55 #include <Thyra_TpetraThyraWrappers.hpp>
56 
57 #include <Teuchos_Ptr.hpp>
59 #include <Teuchos_Assert.hpp>
60 #include <Teuchos_Time.hpp>
61 #include <Teuchos_FancyOStream.hpp>
63 
64 #include <Teko_Utilities.hpp>
65 
66 #include <Xpetra_BlockedCrsMatrix.hpp>
67 #include <Xpetra_CrsMatrixWrap.hpp>
68 #include <Xpetra_IO.hpp>
69 #include <Xpetra_MapExtractorFactory.hpp>
70 #include <Xpetra_Matrix.hpp>
71 #include <Xpetra_MatrixMatrix.hpp>
72 
73 #include "MueLu.hpp"
74 
75 #include "../../research/q2q1/MueLu_Q2Q1PFactory.hpp"
76 #include "../../research/q2q1/MueLu_Q2Q1uPFactory.hpp"
77 
78 #include "MueLu_AmalgamationFactory.hpp"
79 #include "MueLu_BaseClass.hpp"
80 #include "MueLu_BlockedDirectSolver.hpp"
81 #include "MueLu_BlockedPFactory.hpp"
82 #include "MueLu_BlockedRAPFactory.hpp"
83 #include "MueLu_BraessSarazinSmoother.hpp"
84 #include "MueLu_CoalesceDropFactory.hpp"
85 #include "MueLu_ConstraintFactory.hpp"
87 #include "MueLu_DirectSolver.hpp"
88 #include "MueLu_EminPFactory.hpp"
89 #include "MueLu_FactoryManager.hpp"
90 #include "MueLu_FilteredAFactory.hpp"
91 #include "MueLu_GenericRFactory.hpp"
92 #include "MueLu_Level.hpp"
93 #include "MueLu_PatternFactory.hpp"
94 #include "MueLu_SchurComplementFactory.hpp"
95 #include "MueLu_SmootherFactory.hpp"
96 #include "MueLu_SmootherPrototype.hpp"
97 #include "MueLu_SubBlockAFactory.hpp"
98 #include "MueLu_TpetraOperator.hpp"
99 #include "MueLu_TrilinosSmoother.hpp"
100 
101 #include <string>
102 
103 namespace Thyra {
104 
105 #define MUELU_GPD(name, type, defaultValue) \
106  (paramList.isParameter(name) ? paramList.get<type>(name) : defaultValue)
107 
108  using Teuchos::RCP;
109  using Teuchos::rcp;
110  using Teuchos::rcp_const_cast;
111  using Teuchos::rcp_dynamic_cast;
113  using Teuchos::ArrayView;
114  using Teuchos::ArrayRCP;
115  using Teuchos::as;
116  using Teuchos::Array;
117 
118  // Constructors/initializers/accessors
119  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
121 
122 
123  // Overridden from PreconditionerFactoryBase
124  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
126  typedef Thyra ::TpetraLinearOp<SC,LO,GO,NO> ThyraTpetraLinOp;
127  typedef Tpetra::Operator <SC,LO,GO,NO> TpetraLinOp;
128  typedef Tpetra::CrsMatrix <SC,LO,GO,NO> TpetraCrsMat;
129 
130  const RCP<const LinearOpBase<SC> > fwdOp = fwdOpSrc.getOp();
131  const RCP<const ThyraTpetraLinOp> thyraTpetraFwdOp = rcp_dynamic_cast<const ThyraTpetraLinOp>(fwdOp);
132  const RCP<const TpetraLinOp> tpetraFwdOp = Teuchos::nonnull(thyraTpetraFwdOp) ? thyraTpetraFwdOp->getConstTpetraOperator() : Teuchos::null;
133  const RCP<const TpetraCrsMat> tpetraFwdCrsMat = rcp_dynamic_cast<const TpetraCrsMat>(tpetraFwdOp);
134 
135  return Teuchos::nonnull(tpetraFwdCrsMat);
136  }
137 
138  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
141  return rcp(new DefaultPreconditioner<SC>);
142  }
143 
144  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
146  initializePrec(const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc, PreconditionerBase<Scalar> *prec, const ESupportSolveUse supportSolveUse) const {
147  // Check precondition
148  TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
149  TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
150  TEUCHOS_ASSERT(prec);
151 
152  // Retrieve wrapped concrete Tpetra matrix from FwdOp
153  const RCP<const LinearOpBase<SC> > fwdOp = fwdOpSrc->getOp();
155 
156  typedef Thyra::TpetraLinearOp<SC,LO,GO,NO> ThyraTpetraLinOp;
157  const RCP<const ThyraTpetraLinOp> thyraTpetraFwdOp = rcp_dynamic_cast<const ThyraTpetraLinOp>(fwdOp);
158  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraTpetraFwdOp));
159 
160  typedef Tpetra::Operator<SC,LO,GO,NO> TpetraLinOp;
161  const RCP<const TpetraLinOp> tpetraFwdOp = thyraTpetraFwdOp->getConstTpetraOperator();
162  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpetraFwdOp));
163 
164  typedef Tpetra::CrsMatrix<SC,LO,GO,NO> TpetraCrsMat;
165  const RCP<const TpetraCrsMat> tpetraFwdCrsMat = rcp_dynamic_cast<const TpetraCrsMat>(tpetraFwdOp);
166  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpetraFwdCrsMat));
167 
168  // Retrieve concrete preconditioner object
169  const Teuchos::Ptr<DefaultPreconditioner<SC> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<SC> *>(prec));
170  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
171 
172  // Workaround since MueLu interface does not accept const matrix as input
173  const RCP<TpetraCrsMat> tpetraFwdCrsMatNonConst = rcp_const_cast<TpetraCrsMat>(tpetraFwdCrsMat);
174 
175  // Create and compute the initial preconditioner
176 
177  // Create a copy, as we may remove some things from the list
178  ParameterList paramList = *paramList_;
179 
180  typedef Tpetra::MultiVector<SC,LO,GO,NO> MultiVector;
181  RCP<MultiVector> coords, nullspace, velCoords, presCoords;
182  ArrayRCP<LO> p2vMap;
183  Teko::LinearOp thA11, thA12, thA21, thA11_9Pt;
184  if (paramList.isType<RCP<MultiVector> >("Coordinates")) { coords = paramList.get<RCP<MultiVector> >("Coordinates"); paramList.remove("Coordinates"); }
185  if (paramList.isType<RCP<MultiVector> >("Nullspace")) { nullspace = paramList.get<RCP<MultiVector> >("Nullspace"); paramList.remove("Nullspace"); }
186  if (paramList.isType<RCP<MultiVector> >("Velcoords")) { velCoords = paramList.get<RCP<MultiVector> >("Velcoords"); paramList.remove("Velcoords"); }
187  if (paramList.isType<RCP<MultiVector> >("Prescoords")) { presCoords = paramList.get<RCP<MultiVector> >("Prescoords"); paramList.remove("Prescoords"); }
188  if (paramList.isType<ArrayRCP<LO> > ("p2vMap")) { p2vMap = paramList.get<ArrayRCP<LO> > ("p2vMap"); paramList.remove("p2vMap"); }
189  if (paramList.isType<Teko::LinearOp> ("A11")) { thA11 = paramList.get<Teko::LinearOp> ("A11"); paramList.remove("A11"); }
190  if (paramList.isType<Teko::LinearOp> ("A12")) { thA12 = paramList.get<Teko::LinearOp> ("A12"); paramList.remove("A12"); }
191  if (paramList.isType<Teko::LinearOp> ("A21")) { thA21 = paramList.get<Teko::LinearOp> ("A21"); paramList.remove("A21"); }
192  if (paramList.isType<Teko::LinearOp> ("A11_9Pt")) { thA11_9Pt = paramList.get<Teko::LinearOp> ("A11_9Pt"); paramList.remove("A11_9Pt"); }
193 
194  typedef MueLu::TpetraOperator<SC,LO,GO,NO> MueLuOperator;
195  const RCP<MueLuOperator> mueluPrecOp = Q2Q1MkPrecond(paramList, velCoords, presCoords, p2vMap, thA11, thA12, thA21, thA11_9Pt);
196 
197  const RCP<LinearOpBase<SC> > thyraPrecOp = Thyra::createLinearOp(RCP<TpetraLinOp>(mueluPrecOp));
198  defaultPrec->initializeUnspecified(thyraPrecOp);
199  }
200 
201  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
203  uninitializePrec(PreconditionerBase<Scalar> *prec, RCP<const LinearOpSourceBase<Scalar> >* fwdOp, ESupportSolveUse* supportSolveUse) const {
204  // Check precondition
205  TEUCHOS_ASSERT(prec);
206 
207  // Retrieve concrete preconditioner object
208  const Teuchos::Ptr<DefaultPreconditioner<SC> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<SC> *>(prec));
210 
211  if (fwdOp) {
212  // TODO: Implement properly instead of returning default value
213  *fwdOp = Teuchos::null;
214  }
215 
216  if (supportSolveUse) {
217  // TODO: Implement properly instead of returning default value
218  *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
219  }
220 
221  defaultPrec->uninitialize();
222  }
223 
224 
225  // Overridden from ParameterListAcceptor
226  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
229  paramList_ = paramList;
230  }
231 
232  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
235  return paramList_;
236  }
237 
238 
239  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
242  RCP<ParameterList> savedParamList = paramList_;
243  paramList_ = Teuchos::null;
244  return savedParamList;
245  }
246 
247  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
250  return paramList_;
251  }
252 
253  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
256  static RCP<const ParameterList> validPL;
257 
258  if (validPL.is_null())
259  validPL = rcp(new ParameterList());
260 
261  return validPL;
262  }
263 
264  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
267  Q2Q1MkPrecond(const ParameterList& paramList,
268  const RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& velCoords,
269  const RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& presCoords,
270  const ArrayRCP<LocalOrdinal>& p2vMap,
271  const Teko::LinearOp& thA11, const Teko::LinearOp& thA12, const Teko::LinearOp& thA21, const Teko::LinearOp& thA11_9Pt) const
272  {
273  using Teuchos::null;
274 
275  typedef Tpetra::CrsMatrix <SC,LO,GO,NO> TP_Crs;
276  typedef Tpetra::Operator <SC,LO,GO,NO> TP_Op;
277 
278  typedef Xpetra::BlockedCrsMatrix <SC,LO,GO,NO> BlockedCrsMatrix;
279  typedef Xpetra::CrsMatrix <SC,LO,GO,NO> CrsMatrix;
280  typedef Xpetra::CrsMatrixWrap <SC,LO,GO,NO> CrsMatrixWrap;
281  typedef Xpetra::MapExtractorFactory <SC,LO,GO,NO> MapExtractorFactory;
282  typedef Xpetra::MapExtractor <SC,LO,GO,NO> MapExtractor;
283  typedef Xpetra::Map <LO,GO,NO> Map;
284  typedef Xpetra::MapFactory <LO,GO,NO> MapFactory;
285  typedef Xpetra::Matrix <SC,LO,GO,NO> Matrix;
286  typedef Xpetra::MatrixFactory <SC,LO,GO,NO> MatrixFactory;
287  typedef Xpetra::StridedMapFactory <LO,GO,NO> StridedMapFactory;
288 
289  typedef MueLu::Hierarchy <SC,LO,GO,NO> Hierarchy;
290 
291  const RCP<const Teuchos::Comm<int> > comm = velCoords->getMap()->getComm();
292 
293  // Pull out Tpetra matrices
294  RCP<Thyra::LinearOpBase<SC> > ThNonConstA11 = rcp_const_cast<Thyra::LinearOpBase<double> >(thA11);
295  RCP<Thyra::LinearOpBase<SC> > ThNonConstA21 = rcp_const_cast<Thyra::LinearOpBase<double> >(thA21);
296  RCP<Thyra::LinearOpBase<SC> > ThNonConstA12 = rcp_const_cast<Thyra::LinearOpBase<double> >(thA12);
297  RCP<Thyra::LinearOpBase<SC> > ThNonConstA11_9Pt = rcp_const_cast<Thyra::LinearOpBase<double> >(thA11_9Pt);
298 
299  RCP<TP_Op> TpetA11 = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA11);
300  RCP<TP_Op> TpetA21 = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA21);
301  RCP<TP_Op> TpetA12 = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA12);
302  RCP<TP_Op> TpetA11_9Pt = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA11_9Pt);
303 
304  RCP<TP_Crs> TpetCrsA11 = rcp_dynamic_cast<TP_Crs>(TpetA11);
305  RCP<TP_Crs> TpetCrsA21 = rcp_dynamic_cast<TP_Crs>(TpetA21);
306  RCP<TP_Crs> TpetCrsA12 = rcp_dynamic_cast<TP_Crs>(TpetA12);
307  RCP<TP_Crs> TpetCrsA11_9Pt = rcp_dynamic_cast<TP_Crs>(TpetA11_9Pt);
308 
310  RCP<Matrix> tmp_A_21 = MueLu::TpetraCrs_To_XpetraMatrix(TpetCrsA21); // needs map modification
311  RCP<Matrix> tmp_A_12 = MueLu::TpetraCrs_To_XpetraMatrix(TpetCrsA12); // needs map modification
312  RCP<Matrix> A_11_9Pt = MueLu::TpetraCrs_To_XpetraMatrix(TpetCrsA11_9Pt);
313 
314  Xpetra::global_size_t numVel = A_11->getRowMap()->getNodeNumElements();
315  Xpetra::global_size_t numPres = tmp_A_21->getRowMap()->getNodeNumElements();
316 
317  // Create new A21 with map so that the global indices of the row map starts
318  // from numVel+1 (where numVel is the number of rows in the A11 block)
319  RCP<const Map> domainMap2 = tmp_A_12->getDomainMap();
320  RCP<const Map> rangeMap2 = tmp_A_21->getRangeMap();
321  Xpetra::global_size_t numRows2 = rangeMap2->getNodeNumElements();
322  Xpetra::global_size_t numCols2 = domainMap2->getNodeNumElements();
323  ArrayView<const GO> rangeElem2 = rangeMap2->getNodeElementList();
324  ArrayView<const GO> domainElem2 = domainMap2->getNodeElementList();
325  ArrayView<const GO> rowElem1 = tmp_A_12->getRowMap()->getNodeElementList();
326  ArrayView<const GO> colElem1 = tmp_A_21->getColMap()->getNodeElementList();
327 
328  Xpetra::UnderlyingLib lib = domainMap2->lib();
329  GO indexBase = domainMap2->getIndexBase();
330 
331  Array<GO> newRowElem2(numRows2, 0);
332  for (Xpetra::global_size_t i = 0; i < numRows2; i++)
333  newRowElem2[i] = numVel + rangeElem2[i];
334 
335  RCP<const Map> newRangeMap2 = MapFactory::Build(lib, numRows2, newRowElem2, indexBase, comm);
336 
337  // maybe should be column map???
338  Array<GO> newColElem2(numCols2, 0);
339  for (Xpetra::global_size_t i = 0; i < numCols2; i++)
340  newColElem2[i] = numVel + domainElem2[i];
341 
342  RCP<const Map> newDomainMap2 = MapFactory::Build(lib, numCols2, newColElem2, indexBase, comm);
343 
344  RCP<Matrix> A_12 = MatrixFactory::Build(tmp_A_12->getRangeMap(), newDomainMap2, tmp_A_12->getNodeMaxNumRowEntries());
345  RCP<Matrix> A_21 = MatrixFactory::Build(newRangeMap2, tmp_A_21->getDomainMap(), tmp_A_21->getNodeMaxNumRowEntries());
346 
347  RCP<CrsMatrix> A_11_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_11) ->getCrsMatrix();
348  RCP<CrsMatrix> A_12_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_12) ->getCrsMatrix();
349  RCP<CrsMatrix> A_21_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_21) ->getCrsMatrix();
350  RCP<CrsMatrix> A_11_crs_9Pt = rcp_dynamic_cast<CrsMatrixWrap>(A_11_9Pt)->getCrsMatrix();
351 
352 #if 0
353  RCP<Matrix> A_22 = MatrixFactory::Build(newRangeMap2, newDomainMap2, 1);
354  RCP<CrsMatrix> A_22_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_22) ->getCrsMatrix();
355 
356  // FIXME: why do we need to perturb A_22?
357  Array<SC> smallVal(1, 1.0e-10);
358 
359  // FIXME: could this be sped up using expertStaticFillComplete?
360  // There was an attempt on doing it, but it did not do the proper thing
361  // with empty columns. See git history
362  ArrayView<const LO> inds;
363  ArrayView<const SC> vals;
364  for (LO row = 0; row < as<LO>(numRows2); ++row) {
365  tmp_A_21->getLocalRowView(row, inds, vals);
366 
367  size_t nnz = inds.size();
368  Array<GO> newInds(nnz, 0);
369  for (LO colID = 0; colID < as<LO>(nnz); colID++)
370  newInds[colID] = colElem1[inds[colID]];
371 
372  A_21_crs->insertGlobalValues(newRowElem2[row], newInds, vals);
373  A_22_crs->insertGlobalValues(newRowElem2[row], Array<LO>(1, newRowElem2[row]), smallVal);
374  }
375  A_21_crs->fillComplete(tmp_A_21->getDomainMap(), newRangeMap2);
376  A_22_crs->fillComplete(newDomainMap2, newRangeMap2);
377 #else
378  RCP<Matrix> A_22 = Teuchos::null;
379  RCP<CrsMatrix> A_22_crs = Teuchos::null;
380 
381  ArrayView<const LO> inds;
382  ArrayView<const SC> vals;
383  for (LO row = 0; row < as<LO>(numRows2); ++row) {
384  tmp_A_21->getLocalRowView(row, inds, vals);
385 
386  size_t nnz = inds.size();
387  Array<GO> newInds(nnz, 0);
388  for (LO colID = 0; colID < as<LO>(nnz); colID++)
389  newInds[colID] = colElem1[inds[colID]];
390 
391  A_21_crs->insertGlobalValues(newRowElem2[row], newInds, vals);
392  }
393  A_21_crs->fillComplete(tmp_A_21->getDomainMap(), newRangeMap2);
394 #endif
395 
396  // Create new A12 with map so that the global indices of the ColMap starts
397  // from numVel+1 (where numVel is the number of rows in the A11 block)
398  for (LO row = 0; row < as<LO>(tmp_A_12->getRowMap()->getNodeNumElements()); ++row) {
399  tmp_A_12->getLocalRowView(row, inds, vals);
400 
401  size_t nnz = inds.size();
402  Array<GO> newInds(nnz, 0);
403  for (LO colID = 0; colID < as<LO>(nnz); colID++)
404  newInds[colID] = newColElem2[inds[colID]];
405 
406  A_12_crs->insertGlobalValues(rowElem1[row], newInds, vals);
407  }
408  A_12_crs->fillComplete(newDomainMap2, tmp_A_12->getRangeMap());
409 
410  RCP<Matrix> A_12_abs = Absolute(*A_12);
411  RCP<Matrix> A_21_abs = Absolute(*A_21);
412 
413  // =========================================================================
414  // Preconditioner construction - I (block)
415  // =========================================================================
416  RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
417  Teuchos::FancyOStream& out = *fancy;
418  out.setOutputToRootOnly(0);
419  RCP<Matrix> BBt = Xpetra::MatrixMatrix<SC,LO,GO,NO>::Multiply(*A_21, false, *A_12, false, out);
420  RCP<Matrix> BBt_abs = Xpetra::MatrixMatrix<SC,LO,GO,NO>::Multiply(*A_21_abs, false, *A_12_abs, false, out);
421 
422  SC dropTol = (paramList.get<int>("useFilters") ? paramList.get<double>("tau_1") : 0.00);
423  RCP<Matrix> filteredA = FilterMatrix(*A_11, *A_11, dropTol);
424  RCP<Matrix> filteredB = FilterMatrix(*BBt, *BBt_abs, dropTol);
425 
426  RCP<Matrix> fA_11_crs = rcp_dynamic_cast<CrsMatrixWrap>(filteredA);
427  RCP<Matrix> fA_12_crs = Teuchos::null;
428  RCP<Matrix> fA_21_crs = Teuchos::null;
429  RCP<Matrix> fA_22_crs = rcp_dynamic_cast<CrsMatrixWrap>(filteredB);
430 
431  // Build the large filtered matrix which requires strided maps
432  std::vector<size_t> stridingInfo(1, 1);
433  int stridedBlockId = -1;
434 
435  Array<GO> elementList(numVel+numPres); // Not RCP ... does this get cleared ?
436  Array<GO> velElem = A_12_crs->getRangeMap()->getNodeElementList();
437  Array<GO> presElem = A_21_crs->getRangeMap()->getNodeElementList();
438 
439  for (Xpetra::global_size_t i = 0 ; i < numVel; i++) elementList[i] = velElem[i];
440  for (Xpetra::global_size_t i = numVel; i < numVel+numPres; i++) elementList[i] = presElem[i-numVel];
441  RCP<const Map> fullMap = StridedMapFactory::Build(Xpetra::UseTpetra, numVel+numPres, elementList(), indexBase, stridingInfo, comm);
442 
443  std::vector<RCP<const Map> > partMaps(2);
444  partMaps[0] = StridedMapFactory::Build(Xpetra::UseTpetra, numVel, velElem, indexBase, stridingInfo, comm);
445  partMaps[1] = StridedMapFactory::Build(Xpetra::UseTpetra, numPres, presElem, indexBase, stridingInfo, comm, stridedBlockId, numVel);
446 
447  // Map extractors are necessary for Xpetra's block operators
448  RCP<const MapExtractor> mapExtractor = MapExtractorFactory::Build(fullMap, partMaps);
449  RCP<BlockedCrsMatrix> fA = rcp(new BlockedCrsMatrix(mapExtractor, mapExtractor, 10));
450  fA->setMatrix(0, 0, fA_11_crs);
451  fA->setMatrix(0, 1, fA_12_crs);
452  fA->setMatrix(1, 0, fA_21_crs);
453  fA->setMatrix(1, 1, fA_22_crs);
454  fA->fillComplete();
455 
456  // -------------------------------------------------------------------------
457  // Preconditioner construction - I.a (filtered hierarchy)
458  // -------------------------------------------------------------------------
460  SetDependencyTree(M, paramList);
461 
462  RCP<Hierarchy> H = rcp(new Hierarchy);
463  RCP<MueLu::Level> finestLevel = H->GetLevel(0);
464  finestLevel->Set("A", rcp_dynamic_cast<Matrix>(fA));
465  finestLevel->Set("p2vMap", p2vMap);
466  finestLevel->Set("CoordinatesVelocity", Xpetra::toXpetra(velCoords));
467  finestLevel->Set("CoordinatesPressure", Xpetra::toXpetra(presCoords));
468  finestLevel->Set("AForPat", A_11_9Pt);
469  H->SetMaxCoarseSize(MUELU_GPD("coarse: max size", int, 1));
470 
471  // The first invocation of Setup() builds the hierarchy using the filtered
472  // matrix. This build includes the grid transfers but not the creation of the
473  // smoothers.
474  // NOTE: we need to indicate what should be kept from the first invocation
475  // for the second invocation, which then focuses on building the smoothers
476  // for the unfiltered matrix.
477  H->Keep("P", M.GetFactory("P") .get());
478  H->Keep("R", M.GetFactory("R") .get());
479  H->Keep("Ptent", M.GetFactory("Ptent").get());
480  H->Setup(M, 0, MUELU_GPD("max levels", int, 3));
481 
482 #if 0
483  for (int i = 1; i < H->GetNumLevels(); i++) {
484  RCP<Matrix> P = H->GetLevel(i)->template Get<RCP<Matrix> >("P");
485  RCP<BlockedCrsMatrix> Pcrs = rcp_dynamic_cast<BlockedCrsMatrix>(P);
486  RCP<Matrix> Pp = Pcrs->getMatrix(1,1);
487  RCP<Matrix> Pv = Pcrs->getMatrix(0,0);
488 
489  Xpetra::IO<SC,LO,GO,NO>::Write("Pp_l" + MueLu::toString(i) + ".mm", *Pp);
490  Xpetra::IO<SC,LO,GO,NO>::Write("Pv_l" + MueLu::toString(i) + ".mm", *Pv);
491  }
492 #endif
493 
494  // -------------------------------------------------------------------------
495  // Preconditioner construction - I.b (smoothers for unfiltered matrix)
496  // -------------------------------------------------------------------------
497  std::string smootherType = MUELU_GPD("smoother: type", std::string, "vanka");
498  ParameterList smootherParams;
499  if (paramList.isSublist("smoother: params"))
500  smootherParams = paramList.sublist("smoother: params");
501  M.SetFactory("Smoother", GetSmoother(smootherType, smootherParams, false/*coarseSolver?*/));
502 
503  std::string coarseType = MUELU_GPD("coarse: type", std::string, "direct");
504  ParameterList coarseParams;
505  if (paramList.isSublist("coarse: params"))
506  coarseParams = paramList.sublist("coarse: params");
507  M.SetFactory("CoarseSolver", GetSmoother(coarseType, coarseParams, true/*coarseSolver?*/));
508 
509 #ifdef HAVE_MUELU_DEBUG
510  M.ResetDebugData();
511 #endif
512 
513  RCP<BlockedCrsMatrix> A = rcp(new BlockedCrsMatrix(mapExtractor, mapExtractor, 10));
514  A->setMatrix(0, 0, A_11);
515  A->setMatrix(0, 1, A_12);
516  A->setMatrix(1, 0, A_21);
517  A->setMatrix(1, 1, A_22);
518  A->fillComplete();
519 
520  H->GetLevel(0)->Set("A", rcp_dynamic_cast<Matrix>(A));
521 
522  H->Setup(M, 0, H->GetNumLevels());
523 
525  }
526 
527  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
530  FilterMatrix(Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Pattern, Scalar dropTol) const {
531  typedef Xpetra::Matrix<SC,LO,GO,NO> Matrix;
532  typedef MueLu::AmalgamationFactory<SC,LO,GO,NO> AmalgamationFactory;
533  typedef MueLu::CoalesceDropFactory<SC,LO,GO,NO> CoalesceDropFactory;
534  typedef MueLu::FilteredAFactory<SC,LO,GO,NO> FilteredAFactory;
535  typedef MueLu::GraphBase<LO,GO,NO> GraphBase;
536 
537  RCP<GraphBase> filteredGraph;
538  {
539  // Get graph pattern for the pattern matrix
540  MueLu::Level level;
541  level.SetLevelID(1);
542 
543  level.Set<RCP<Matrix> >("A", rcpFromRef(Pattern));
544 
545  RCP<AmalgamationFactory> amalgFactory = rcp(new AmalgamationFactory());
546 
547  RCP<CoalesceDropFactory> dropFactory = rcp(new CoalesceDropFactory());
548  ParameterList dropParams = *(dropFactory->GetValidParameterList());
549  dropParams.set("lightweight wrap", true);
550  dropParams.set("aggregation: drop scheme", "classical");
551  dropParams.set("aggregation: drop tol", dropTol);
552  // dropParams.set("Dirichlet detection threshold", <>);
553  dropFactory->SetParameterList(dropParams);
554  dropFactory->SetFactory("UnAmalgamationInfo", amalgFactory);
555 
556  // Build
557  level.Request("Graph", dropFactory.get());
558  dropFactory->Build(level);
559 
560  level.Get("Graph", filteredGraph, dropFactory.get());
561  }
562 
563  RCP<Matrix> filteredA;
564  {
565  // Filter the original matrix, not the pattern one
566  MueLu::Level level;
567  level.SetLevelID(1);
568 
569  level.Set("A", rcpFromRef(A));
570  level.Set("Graph", filteredGraph);
571  level.Set("Filtering", true);
572 
573  RCP<FilteredAFactory> filterFactory = rcp(new FilteredAFactory());
574  ParameterList filterParams = *(filterFactory->GetValidParameterList());
575  // We need a graph that has proper structure in it. Therefore, we need to
576  // drop older pattern, i.e. not to reuse it
577  filterParams.set("filtered matrix: reuse graph", false);
578  filterParams.set("filtered matrix: use lumping", false);
579  filterFactory->SetParameterList(filterParams);
580 
581  // Build
582  level.Request("A", filterFactory.get());
583  filterFactory->Build(level);
584 
585  level.Get("A", filteredA, filterFactory.get());
586  }
587 
588  // Zero out row sums by fixing the diagonal
589  filteredA->resumeFill();
590  size_t numRows = filteredA->getRowMap()->getNodeNumElements();
591  for (size_t i = 0; i < numRows; i++) {
592  ArrayView<const LO> inds;
593  ArrayView<const SC> vals;
594  filteredA->getLocalRowView(i, inds, vals);
595 
596  size_t nnz = inds.size();
597 
598  Array<SC> valsNew = vals;
599 
600  LO diagIndex = -1;
602  for (size_t j = 0; j < nnz; j++) {
603  diag += vals[j];
604  if (inds[j] == Teuchos::as<int>(i))
605  diagIndex = j;
606  }
608  "No diagonal found");
609  if (nnz <= 1)
610  continue;
611 
612  valsNew[diagIndex] -= diag;
613 
614  filteredA->replaceLocalValues(i, inds, valsNew);
615  }
616  filteredA->fillComplete();
617 
618  return filteredA;
619  }
620 
621  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
622  void
625  typedef MueLu::BlockedPFactory <SC,LO,GO,NO> BlockedPFactory;
626  typedef MueLu::GenericRFactory <SC,LO,GO,NO> GenericRFactory;
627  typedef MueLu::BlockedRAPFactory <SC,LO,GO,NO> BlockedRAPFactory;
628  typedef MueLu::SmootherFactory <SC,LO,GO,NO> SmootherFactory;
629  typedef MueLu::BlockedDirectSolver <SC,LO,GO,NO> BlockedDirectSolver;
630  typedef MueLu::FactoryManager <SC,LO,GO,NO> FactoryManager;
631 
632  // Pressure and velocity dependency trees are identical. The only
633  // difference is that pressure has to go first, so that velocity can use
634  // some of pressure data
635  RCP<FactoryManager> M11 = rcp(new FactoryManager()), M22 = rcp(new FactoryManager());
636  M11->SetKokkosRefactor(paramList.get<bool>("use kokkos refactor"));
637  M22->SetKokkosRefactor(paramList.get<bool>("use kokkos refactor"));
638  SetBlockDependencyTree(*M11, 0, 0, "velocity", paramList);
639  SetBlockDependencyTree(*M22, 1, 1, "pressure", paramList);
640 
641  RCP<BlockedPFactory> PFact = rcp(new BlockedPFactory());
642  ParameterList pParamList = *(PFact->GetValidParameterList());
643  pParamList.set("backwards", true); // do pressure first
644  PFact->SetParameterList(pParamList);
645  PFact->AddFactoryManager(M11);
646  PFact->AddFactoryManager(M22);
647  M.SetFactory("P", PFact);
648 
649  RCP<GenericRFactory> RFact = rcp(new GenericRFactory());
650  RFact->SetFactory("P", PFact);
651  M.SetFactory("R", RFact);
652 
653  RCP<MueLu::Factory > AcFact = rcp(new BlockedRAPFactory());
654  AcFact->SetFactory("R", RFact);
655  AcFact->SetFactory("P", PFact);
656  M.SetFactory("A", AcFact);
657 
658  // Smoothers will be set later
659  M.SetFactory("Smoother", Teuchos::null);
660 
661  RCP<MueLu::Factory> coarseFact = rcp(new SmootherFactory(rcp(new BlockedDirectSolver()), Teuchos::null));
662  // M.SetFactory("CoarseSolver", coarseFact);
663  M.SetFactory("CoarseSolver", Teuchos::null);
664  }
665 
666  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
667  void
670  typedef MueLu::ConstraintFactory <SC,LO,GO,NO> ConstraintFactory;
671  typedef MueLu::EminPFactory <SC,LO,GO,NO> EminPFactory;
672  typedef MueLu::GenericRFactory <SC,LO,GO,NO> GenericRFactory;
673  typedef MueLu::PatternFactory <SC,LO,GO,NO> PatternFactory;
674  typedef MueLu::Q2Q1PFactory <SC,LO,GO,NO> Q2Q1PFactory;
675  typedef MueLu::Q2Q1uPFactory <SC,LO,GO,NO> Q2Q1uPFactory;
676  typedef MueLu::SubBlockAFactory <SC,LO,GO,NO> SubBlockAFactory;
677 
678  RCP<SubBlockAFactory> AFact = rcp(new SubBlockAFactory());
679  AFact->SetFactory ("A", MueLu::NoFactory::getRCP());
680  AFact->SetParameter("block row", Teuchos::ParameterEntry(row));
681  AFact->SetParameter("block col", Teuchos::ParameterEntry(col));
682  M.SetFactory("A", AFact);
683 
684  RCP<MueLu::Factory> Q2Q1Fact;
685 
686  const bool isStructured = false;
687 
688  if (isStructured) {
689  Q2Q1Fact = rcp(new Q2Q1PFactory);
690 
691  } else {
692  Q2Q1Fact = rcp(new Q2Q1uPFactory);
693  ParameterList q2q1ParamList = *(Q2Q1Fact->GetValidParameterList());
694  q2q1ParamList.set("mode", mode);
695  if (paramList.isParameter("dump status"))
696  q2q1ParamList.set("dump status", paramList.get<bool>("dump status"));
697  if (paramList.isParameter("phase2"))
698  q2q1ParamList.set("phase2", paramList.get<bool>("phase2"));
699  if (paramList.isParameter("tau_2"))
700  q2q1ParamList.set("tau_2", paramList.get<double>("tau_2"));
701  Q2Q1Fact->SetParameterList(q2q1ParamList);
702  }
703  Q2Q1Fact->SetFactory("A", AFact);
704  M.SetFactory("Ptent", Q2Q1Fact);
705 
706  RCP<PatternFactory> patternFact = rcp(new PatternFactory);
707  ParameterList patternParams = *(patternFact->GetValidParameterList());
708  // Our prolongator constructs the exact pattern we are going to use,
709  // therefore we do not expand it
710  patternParams.set("emin: pattern order", 0);
711  patternFact->SetParameterList(patternParams);
712  patternFact->SetFactory("A", AFact);
713  patternFact->SetFactory("P", Q2Q1Fact);
714  M.SetFactory("Ppattern", patternFact);
715 
716  RCP<ConstraintFactory> CFact = rcp(new ConstraintFactory);
717  CFact->SetFactory("Ppattern", patternFact);
718  M.SetFactory("Constraint", CFact);
719 
720  RCP<EminPFactory> EminPFact = rcp(new EminPFactory());
721  ParameterList eminParams = *(EminPFact->GetValidParameterList());
722  if (paramList.isParameter("emin: num iterations"))
723  eminParams.set("emin: num iterations", paramList.get<int>("emin: num iterations"));
724  if (mode == "pressure") {
725  eminParams.set("emin: iterative method", "cg");
726  } else {
727  eminParams.set("emin: iterative method", "gmres");
728  if (paramList.isParameter("emin: iterative method"))
729  eminParams.set("emin: iterative method", paramList.get<std::string>("emin: iterative method"));
730  }
731  EminPFact->SetParameterList(eminParams);
732  EminPFact->SetFactory("A", AFact);
733  EminPFact->SetFactory("Constraint", CFact);
734  EminPFact->SetFactory("P", Q2Q1Fact);
735  M.SetFactory("P", EminPFact);
736 
737  if (mode == "velocity" && (!paramList.isParameter("velocity: use transpose") || paramList.get<bool>("velocity: use transpose") == false)) {
738  // Pressure system is symmetric, so it does not matter
739  // Velocity system may benefit from running emin in restriction mode (with A^T)
740  RCP<GenericRFactory> RFact = rcp(new GenericRFactory());
741  RFact->SetFactory("P", EminPFact);
742  M.SetFactory("R", RFact);
743  }
744  }
745 
746  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
749  GetSmoother(const std::string& type, const ParameterList& paramList, bool coarseSolver) const {
750  typedef Teuchos::ParameterEntry ParameterEntry;
751 
752  typedef MueLu::BlockedDirectSolver <SC,LO,GO,NO> BlockedDirectSolver;
753  typedef MueLu::BraessSarazinSmoother <SC,LO,GO,NO> BraessSarazinSmoother;
754  typedef MueLu::DirectSolver <SC,LO,GO,NO> DirectSolver;
755  typedef MueLu::FactoryManager <SC,LO,GO,NO> FactoryManager;
756  typedef MueLu::SchurComplementFactory<SC,LO,GO,NO> SchurComplementFactory;
757  typedef MueLu::SmootherFactory <SC,LO,GO,NO> SmootherFactory;
758  typedef MueLu::SmootherPrototype <SC,LO,GO,NO> SmootherPrototype;
759  typedef MueLu::TrilinosSmoother <SC,LO,GO,NO> TrilinosSmoother;
760 
761  RCP<SmootherPrototype> smootherPrototype;
762  if (type == "none") {
763  return Teuchos::null;
764 
765  } else if (type == "vanka") {
766  // Set up Vanka smoothing via a combination of Schwarz and block relaxation.
767  ParameterList schwarzList;
768  schwarzList.set("schwarz: overlap level", as<int>(0));
769  schwarzList.set("schwarz: zero starting solution", false);
770  schwarzList.set("subdomain solver name", "Block_Relaxation");
771 
772  ParameterList& innerSolverList = schwarzList.sublist("subdomain solver parameters");
773  innerSolverList.set("partitioner: type", "user");
774  innerSolverList.set("partitioner: overlap", MUELU_GPD("partitioner: overlap", int, 1));
775  innerSolverList.set("relaxation: type", MUELU_GPD("relaxation: type", std::string, "Gauss-Seidel"));
776  innerSolverList.set("relaxation: sweeps", MUELU_GPD("relaxation: sweeps", int, 1));
777  innerSolverList.set("relaxation: damping factor", MUELU_GPD("relaxation: damping factor", double, 0.5));
778  innerSolverList.set("relaxation: zero starting solution", false);
779  // innerSolverList.set("relaxation: backward mode", MUELU_GPD("relaxation: backward mode", bool, true); NOT SUPPORTED YET
780 
781  std::string ifpackType = "SCHWARZ";
782 
783  smootherPrototype = rcp(new TrilinosSmoother(ifpackType, schwarzList));
784 
785  } else if (type == "schwarz") {
786 
787  std::string ifpackType = "SCHWARZ";
788 
789  smootherPrototype = rcp(new TrilinosSmoother(ifpackType, paramList));
790 
791  } else if (type == "braess-sarazin") {
792  // Define smoother/solver for BraessSarazin
793  SC omega = MUELU_GPD("bs: omega", double, 1.0);
794  bool lumping = MUELU_GPD("bs: lumping", bool, false);
795 
796  RCP<SchurComplementFactory> schurFact = rcp(new SchurComplementFactory());
797  schurFact->SetParameter("omega", ParameterEntry(omega));
798  schurFact->SetParameter("lumping", ParameterEntry(lumping));
799  schurFact->SetFactory ("A", MueLu::NoFactory::getRCP());
800 
801  // Schur complement solver
802  RCP<SmootherPrototype> schurSmootherPrototype;
803  std::string schurSmootherType = (paramList.isParameter("schur smoother: type") ? paramList.get<std::string>("schur smoother: type") : "RELAXATION");
804  if (schurSmootherType == "RELAXATION") {
805  ParameterList schurSmootherParams = paramList.sublist("schur smoother: params");
806  // schurSmootherParams.set("relaxation: damping factor", omega);
807  schurSmootherPrototype = rcp(new TrilinosSmoother(schurSmootherType, schurSmootherParams));
808  } else {
809  schurSmootherPrototype = rcp(new DirectSolver());
810  }
811  schurSmootherPrototype->SetFactory("A", schurFact);
812 
813  RCP<SmootherFactory> schurSmootherFact = rcp(new SmootherFactory(schurSmootherPrototype));
814 
815  // Define temporary FactoryManager that is used as input for BraessSarazin smoother
816  RCP<FactoryManager> braessManager = rcp(new FactoryManager());
817  braessManager->SetFactory("A", schurFact); // SchurComplement operator for correction step (defined as "A")
818  braessManager->SetFactory("Smoother", schurSmootherFact); // solver/smoother for correction step
819  braessManager->SetFactory("PreSmoother", schurSmootherFact);
820  braessManager->SetFactory("PostSmoother", schurSmootherFact);
821  braessManager->SetIgnoreUserData(true); // always use data from factories defined in factory manager
822 
823  smootherPrototype = rcp(new BraessSarazinSmoother());
824  smootherPrototype->SetParameter("Sweeps", ParameterEntry(MUELU_GPD("bs: sweeps", int, 1)));
825  smootherPrototype->SetParameter("lumping", ParameterEntry(lumping));
826  smootherPrototype->SetParameter("Damping factor", ParameterEntry(omega));
827  smootherPrototype->SetParameter("q2q1 mode", ParameterEntry(true));
828  rcp_dynamic_cast<BraessSarazinSmoother>(smootherPrototype)->AddFactoryManager(braessManager, 0); // set temporary factory manager in BraessSarazin smoother
829 
830  } else if (type == "ilu") {
831  std::string ifpackType = "RILUK";
832 
833  smootherPrototype = rcp(new TrilinosSmoother(ifpackType, paramList));
834 
835  } else if (type == "direct") {
836  smootherPrototype = rcp(new BlockedDirectSolver());
837 
838  } else {
839  throw MueLu::Exceptions::RuntimeError("Unknown smoother type: \"" + type + "\"");
840  }
841 
842  return coarseSolver ? rcp(new SmootherFactory(smootherPrototype, Teuchos::null)) : rcp(new SmootherFactory(smootherPrototype));
843  }
844 
845  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
847  MueLuTpetraQ2Q1PreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Absolute(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A) const {
848  typedef Xpetra::CrsMatrix <SC,LO,GO,NO> CrsMatrix;
849  typedef Xpetra::CrsMatrixWrap<SC,LO,GO,NO> CrsMatrixWrap;
850  typedef Xpetra::Matrix <SC,LO,GO,NO> Matrix;
851 
852  const CrsMatrixWrap& Awrap = dynamic_cast<const CrsMatrixWrap&>(A);
853 
855  ArrayRCP<const LO> jaA;
856  ArrayRCP<const SC> valA;
857  Awrap.getCrsMatrix()->getAllValues(iaA, jaA, valA);
858 
859  ArrayRCP<size_t> iaB (iaA .size());
860  ArrayRCP<LO> jaB (jaA .size());
861  ArrayRCP<SC> valB(valA.size());
862  for (int i = 0; i < iaA .size(); i++) iaB [i] = iaA[i];
863  for (int i = 0; i < jaA .size(); i++) jaB [i] = jaA[i];
864  for (int i = 0; i < valA.size(); i++) valB[i] = Teuchos::ScalarTraits<SC>::magnitude(valA[i]);
865 
866  RCP<Matrix> B = rcp(new CrsMatrixWrap(A.getRowMap(), A.getColMap(), 0));
867  RCP<CrsMatrix> Bcrs = rcp_dynamic_cast<CrsMatrixWrap>(B)->getCrsMatrix();
868  Bcrs->setAllValues(iaB, jaB, valB);
869  Bcrs->expertStaticFillComplete(A.getDomainMap(), A.getRangeMap());
870 
871  return B;
872  }
873 
874  // Public functions overridden from Teuchos::Describable
875  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
877  return "Thyra::MueLuTpetraQ2Q1PreconditionerFactory";
878  }
879 
880 } // namespace Thyra
881 
882 #endif
883 #endif // ifdef THYRA_MUELU_TPETRA_Q2Q1PRECONDITIONER_FACTORY_DEF_HPP
Generic Smoother Factory for generating the smoothers of the MG hierarchy.
void SetDependencyTree(MueLu::FactoryManager< SC, LO, GO, NO > &M, const ParameterList &paramList) const
This class specifies the default factory that should generate some data on a Level if the data does n...
MueLu::DefaultLocalOrdinal LocalOrdinal
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level-&gt;Get&lt; RCP&lt;Matrix&gt; &gt;(&quot;A&quot;, factory) if factory == NULL =&gt; use default factory.
BraessSarazin smoother for 2x2 block matrices.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &paramList)
Factory for building blocked, segregated prolongation operators.
T & get(const std::string &name, T def_value)
Class that encapsulates external library smoothers.
Factory for building coarse matrices.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
bool nonnull(const std::shared_ptr< T > &p)
bool is_null(const std::shared_ptr< T > &p)
Teuchos::RCP< const Teuchos::ParameterList > getParameterList() const
void uninitializePrec(PreconditionerBase< SC > *prec, Teuchos::RCP< const LinearOpSourceBase< SC > > *fwdOp, ESupportSolveUse *supportSolveUse) const
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
size_type size() const
Base class for smoother prototypes.
T * get() const
size_type size() const
const RCP< const FactoryBase > GetFactory(const std::string &varName) const
Get factory associated with a particular data name.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Class that encapsulates direct solvers. Autoselection of AmesosSmoother or Amesos2Smoother according ...
void initializePrec(const Teuchos::RCP< const LinearOpSourceBase< SC > > &fwdOp, PreconditionerBase< SC > *prec, const ESupportSolveUse supportSolveUse) const
Factory for building restriction operators using a prolongator factory.
bool isParameter(const std::string &name) const
virtual void SetParameterList(const Teuchos::ParameterList &paramList)
Set parameters from a parameter list and return with default values.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultScalar Scalar
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
Teuchos::RCP< Xpetra::Matrix< SC, LO, GO, NO > > FilterMatrix(Xpetra::Matrix< SC, LO, GO, NO > &A, Xpetra::Matrix< SC, LO, GO, NO > &Pattern, SC dropTol) const
bool isSublist(const std::string &name) const
Teuchos::RCP< MueLu::TpetraOperator< SC, LO, GO, NO > > Q2Q1MkPrecond(const ParameterList &paramList, const Teuchos::RCP< Tpetra::MultiVector< SC, LO, GO, NO > > &velCoords, const Teuchos::RCP< Tpetra::MultiVector< SC, LO, GO, NO > > &presCoords, const Teuchos::ArrayRCP< LO > &p2vMap, const Teko::LinearOp &thA11, const Teko::LinearOp &thA12, const Teko::LinearOp &thA21, const Teko::LinearOp &thA11_9Pt) const
Factory for building a thresholded operator.
AmalgamationFactory for subblocks of strided map based amalgamation data.
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
Factory for building the constraint operator.
direct solver for nxn blocked matrices
bool isCompatible(const LinearOpSourceBase< SC > &fwdOp) const
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
Various adapters that will create a MueLu preconditioner that is a Tpetra::Operator.
void SetFactory(const std::string &varName, const RCP< const FactoryBase > &factory)
Set Factory.
MueLu representation of a graph.
void SetBlockDependencyTree(MueLu::FactoryManager< SC, LO, GO, NO > &M, LO row, LO col, const std::string &mode, const ParameterList &paramList) const
Factory for building the Schur Complement for a 2x2 block matrix.
Factory for creating a graph base on a given matrix.
Teuchos::RCP< Xpetra::Matrix< SC, LO, GO, NO > > Absolute(const Xpetra::Matrix< SC, LO, GO, NO > &A) const
void SetLevelID(int levelID)
Set level number.
Definition: MueLu_Level.cpp:78
TypeTo as(const TypeFrom &t)
virtual void SetFactory(const std::string &varName, const RCP< const FactoryBase > &factory)=0
Configuration.
Factory for building nonzero patterns for energy minimization.
RCP< Xpetra::Matrix< SC, LO, GO, NO > > TpetraCrs_To_XpetraMatrix(const Teuchos::RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > &Atpetra)
Factory for building Energy Minimization prolongators.
Wraps an existing MueLu::Hierarchy as a Tpetra::Operator.
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Exception throws to report errors in the internal logical of the program.
Factory for building filtered matrices using filtered graphs.
#define MUELU_GPD(name, type, defaultValue)
#define TEUCHOS_ASSERT(assertion_test)
T * get() const
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
virtual Teuchos::RCP< const Teuchos::ParameterList > GetValidParameterList() const =0
Return a const parameter list of valid parameters that setParameterList() will accept.
void Request(const FactoryBase &factory)
Increment the storage counter for all the inputs of a factory.
RCP< MueLu::FactoryBase > GetSmoother(const std::string &type, const ParameterList &paramList, bool coarseSolver) const
static const RCP< const NoFactory > getRCP()
Static Get() functions.
bool is_null() const