MueLu  Version of the Day
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>
58 #include <Teuchos_TestForException.hpp>
59 #include <Teuchos_Assert.hpp>
60 #include <Teuchos_Time.hpp>
61 #include <Teuchos_FancyOStream.hpp>
62 #include <Teuchos_VerbosityLevel.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 // #define IMPLICIT_TRANSPOSE
104 
105 namespace Thyra {
106 
107 #define MUELU_GPD(name, type, defaultValue) \
108  (paramList.isParameter(name) ? paramList.get<type>(name) : defaultValue)
109 
110  using Teuchos::RCP;
111  using Teuchos::rcp;
112  using Teuchos::rcp_const_cast;
113  using Teuchos::rcp_dynamic_cast;
114  using Teuchos::ParameterList;
115  using Teuchos::ArrayView;
116  using Teuchos::ArrayRCP;
117  using Teuchos::as;
118  using Teuchos::Array;
119 
120  // Constructors/initializers/accessors
121  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
123 
124 
125  // Overridden from PreconditionerFactoryBase
126  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
128  typedef Thyra ::TpetraLinearOp<SC,LO,GO,NO> ThyraTpetraLinOp;
129  typedef Tpetra::Operator <SC,LO,GO,NO> TpetraLinOp;
130  typedef Tpetra::CrsMatrix <SC,LO,GO,NO> TpetraCrsMat;
131 
132  const RCP<const LinearOpBase<SC> > fwdOp = fwdOpSrc.getOp();
133  const RCP<const ThyraTpetraLinOp> thyraTpetraFwdOp = rcp_dynamic_cast<const ThyraTpetraLinOp>(fwdOp);
134  const RCP<const TpetraLinOp> tpetraFwdOp = Teuchos::nonnull(thyraTpetraFwdOp) ? thyraTpetraFwdOp->getConstTpetraOperator() : Teuchos::null;
135  const RCP<const TpetraCrsMat> tpetraFwdCrsMat = rcp_dynamic_cast<const TpetraCrsMat>(tpetraFwdOp);
136 
137  return Teuchos::nonnull(tpetraFwdCrsMat);
138  }
139 
140  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
141  RCP<PreconditionerBase<Scalar> >
143  return rcp(new DefaultPreconditioner<SC>);
144  }
145 
146  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
148  initializePrec(const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc, PreconditionerBase<Scalar> *prec, const ESupportSolveUse supportSolveUse) const {
149  // Check precondition
150  TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
151  TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
152  TEUCHOS_ASSERT(prec);
153 
154  // Retrieve wrapped concrete Tpetra matrix from FwdOp
155  const RCP<const LinearOpBase<SC> > fwdOp = fwdOpSrc->getOp();
156  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
157 
158  typedef Thyra::TpetraLinearOp<SC,LO,GO,NO> ThyraTpetraLinOp;
159  const RCP<const ThyraTpetraLinOp> thyraTpetraFwdOp = rcp_dynamic_cast<const ThyraTpetraLinOp>(fwdOp);
160  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraTpetraFwdOp));
161 
162  typedef Tpetra::Operator<SC,LO,GO,NO> TpetraLinOp;
163  const RCP<const TpetraLinOp> tpetraFwdOp = thyraTpetraFwdOp->getConstTpetraOperator();
164  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpetraFwdOp));
165 
166  typedef Tpetra::CrsMatrix<SC,LO,GO,NO> TpetraCrsMat;
167  const RCP<const TpetraCrsMat> tpetraFwdCrsMat = rcp_dynamic_cast<const TpetraCrsMat>(tpetraFwdOp);
168  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpetraFwdCrsMat));
169 
170  // Retrieve concrete preconditioner object
171  const Teuchos::Ptr<DefaultPreconditioner<SC> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<SC> *>(prec));
172  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
173 
174  // Workaround since MueLu interface does not accept const matrix as input
175  const RCP<TpetraCrsMat> tpetraFwdCrsMatNonConst = rcp_const_cast<TpetraCrsMat>(tpetraFwdCrsMat);
176 
177  // Create and compute the initial preconditioner
178 
179  // Create a copy, as we may remove some things from the list
180  ParameterList paramList = *paramList_;
181 
182  typedef Tpetra::MultiVector<SC,LO,GO,NO> MultiVector;
183  RCP<MultiVector> coords, nullspace, velCoords, presCoords;
184  ArrayRCP<LO> p2vMap;
185  Teko::LinearOp thA11, thA12, thA21, thA11_9Pt;
186  if (paramList.isType<RCP<MultiVector> >("Coordinates")) { coords = paramList.get<RCP<MultiVector> >("Coordinates"); paramList.remove("Coordinates"); }
187  if (paramList.isType<RCP<MultiVector> >("Nullspace")) { nullspace = paramList.get<RCP<MultiVector> >("Nullspace"); paramList.remove("Nullspace"); }
188  if (paramList.isType<RCP<MultiVector> >("Velcoords")) { velCoords = paramList.get<RCP<MultiVector> >("Velcoords"); paramList.remove("Velcoords"); }
189  if (paramList.isType<RCP<MultiVector> >("Prescoords")) { presCoords = paramList.get<RCP<MultiVector> >("Prescoords"); paramList.remove("Prescoords"); }
190  if (paramList.isType<ArrayRCP<LO> > ("p2vMap")) { p2vMap = paramList.get<ArrayRCP<LO> > ("p2vMap"); paramList.remove("p2vMap"); }
191  if (paramList.isType<Teko::LinearOp> ("A11")) { thA11 = paramList.get<Teko::LinearOp> ("A11"); paramList.remove("A11"); }
192  if (paramList.isType<Teko::LinearOp> ("A12")) { thA12 = paramList.get<Teko::LinearOp> ("A12"); paramList.remove("A12"); }
193  if (paramList.isType<Teko::LinearOp> ("A21")) { thA21 = paramList.get<Teko::LinearOp> ("A21"); paramList.remove("A21"); }
194  if (paramList.isType<Teko::LinearOp> ("A11_9Pt")) { thA11_9Pt = paramList.get<Teko::LinearOp> ("A11_9Pt"); paramList.remove("A11_9Pt"); }
195 
196  typedef MueLu::TpetraOperator<SC,LO,GO,NO> MueLuOperator;
197  const RCP<MueLuOperator> mueluPrecOp = Q2Q1MkPrecond(paramList, velCoords, presCoords, p2vMap, thA11, thA12, thA21, thA11_9Pt);
198 
199  const RCP<LinearOpBase<SC> > thyraPrecOp = Thyra::createLinearOp(RCP<TpetraLinOp>(mueluPrecOp));
200  defaultPrec->initializeUnspecified(thyraPrecOp);
201  }
202 
203  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
205  uninitializePrec(PreconditionerBase<Scalar> *prec, RCP<const LinearOpSourceBase<Scalar> >* fwdOp, ESupportSolveUse* supportSolveUse) const {
206  // Check precondition
207  TEUCHOS_ASSERT(prec);
208 
209  // Retrieve concrete preconditioner object
210  const Teuchos::Ptr<DefaultPreconditioner<SC> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<SC> *>(prec));
211  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
212 
213  if (fwdOp) {
214  // TODO: Implement properly instead of returning default value
215  *fwdOp = Teuchos::null;
216  }
217 
218  if (supportSolveUse) {
219  // TODO: Implement properly instead of returning default value
220  *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
221  }
222 
223  defaultPrec->uninitialize();
224  }
225 
226 
227  // Overridden from ParameterListAcceptor
228  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
230  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
231  paramList_ = paramList;
232  }
233 
234  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
235  RCP<ParameterList>
237  return paramList_;
238  }
239 
240 
241  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
242  RCP<ParameterList>
244  RCP<ParameterList> savedParamList = paramList_;
245  paramList_ = Teuchos::null;
246  return savedParamList;
247  }
248 
249  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
250  RCP<const ParameterList>
252  return paramList_;
253  }
254 
255  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
256  RCP<const ParameterList>
258  static RCP<const ParameterList> validPL;
259 
260  if (validPL.is_null())
261  validPL = rcp(new ParameterList());
262 
263  return validPL;
264  }
265 
266  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
267  RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
269  Q2Q1MkPrecond(const ParameterList& paramList,
270  const RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& velCoords,
271  const RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& presCoords,
272  const ArrayRCP<LocalOrdinal>& p2vMap,
273  const Teko::LinearOp& thA11, const Teko::LinearOp& thA12, const Teko::LinearOp& thA21, const Teko::LinearOp& thA11_9Pt) const
274  {
275  using Teuchos::null;
276 
277  typedef Tpetra::CrsMatrix <SC,LO,GO,NO> TP_Crs;
278  typedef Tpetra::Operator <SC,LO,GO,NO> TP_Op;
279 
280  typedef Xpetra::BlockedCrsMatrix <SC,LO,GO,NO> BlockedCrsMatrix;
281  typedef Xpetra::CrsMatrix <SC,LO,GO,NO> CrsMatrix;
282  typedef Xpetra::CrsMatrixWrap <SC,LO,GO,NO> CrsMatrixWrap;
283  typedef Xpetra::MapExtractorFactory <SC,LO,GO,NO> MapExtractorFactory;
284  typedef Xpetra::MapExtractor <SC,LO,GO,NO> MapExtractor;
285  typedef Xpetra::Map <LO,GO,NO> Map;
286  typedef Xpetra::MapExtractor <SC,LO,GO,NO> MapExtractor;
287  typedef Xpetra::MapExtractorFactory <SC,LO,GO,NO> MapExtractorFactory;
288  typedef Xpetra::MapFactory <LO,GO,NO> MapFactory;
289  typedef Xpetra::Matrix <SC,LO,GO,NO> Matrix;
290  typedef Xpetra::MatrixFactory <SC,LO,GO,NO> MatrixFactory;
291  typedef Xpetra::StridedMapFactory <LO,GO,NO> StridedMapFactory;
292 
293  typedef MueLu::Hierarchy <SC,LO,GO,NO> Hierarchy;
294 
295  const RCP<const Teuchos::Comm<int> > comm = velCoords->getMap()->getComm();
296 
297  // Pull out Tpetra matrices
298  RCP<Thyra::LinearOpBase<SC> > ThNonConstA11 = rcp_const_cast<Thyra::LinearOpBase<double> >(thA11);
299  RCP<Thyra::LinearOpBase<SC> > ThNonConstA21 = rcp_const_cast<Thyra::LinearOpBase<double> >(thA21);
300  RCP<Thyra::LinearOpBase<SC> > ThNonConstA12 = rcp_const_cast<Thyra::LinearOpBase<double> >(thA12);
301  RCP<Thyra::LinearOpBase<SC> > ThNonConstA11_9Pt = rcp_const_cast<Thyra::LinearOpBase<double> >(thA11_9Pt);
302 
303  RCP<TP_Op> TpetA11 = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA11);
304  RCP<TP_Op> TpetA21 = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA21);
305  RCP<TP_Op> TpetA12 = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA12);
306  RCP<TP_Op> TpetA11_9Pt = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA11_9Pt);
307 
308  RCP<TP_Crs> TpetCrsA11 = rcp_dynamic_cast<TP_Crs>(TpetA11);
309  RCP<TP_Crs> TpetCrsA21 = rcp_dynamic_cast<TP_Crs>(TpetA21);
310  RCP<TP_Crs> TpetCrsA12 = rcp_dynamic_cast<TP_Crs>(TpetA12);
311  RCP<TP_Crs> TpetCrsA11_9Pt = rcp_dynamic_cast<TP_Crs>(TpetA11_9Pt);
312 
313  RCP<Matrix> A_11 = MueLu::TpetraCrs_To_XpetraMatrix(TpetCrsA11);
314  RCP<Matrix> tmp_A_21 = MueLu::TpetraCrs_To_XpetraMatrix(TpetCrsA21); // needs map modification
315  RCP<Matrix> tmp_A_12 = MueLu::TpetraCrs_To_XpetraMatrix(TpetCrsA12); // needs map modification
316  RCP<Matrix> A_11_9Pt = MueLu::TpetraCrs_To_XpetraMatrix(TpetCrsA11_9Pt);
317 
318  Xpetra::global_size_t numVel = A_11->getRowMap()->getNodeNumElements();
319  Xpetra::global_size_t numPres = tmp_A_21->getRowMap()->getNodeNumElements();
320 
321  // Create new A21 with map so that the global indices of the row map starts
322  // from numVel+1 (where numVel is the number of rows in the A11 block)
323  RCP<const Map> domainMap2 = tmp_A_12->getDomainMap();
324  RCP<const Map> rangeMap2 = tmp_A_21->getRangeMap();
325  Xpetra::global_size_t numRows2 = rangeMap2->getNodeNumElements();
326  Xpetra::global_size_t numCols2 = domainMap2->getNodeNumElements();
327  ArrayView<const GO> rangeElem2 = rangeMap2->getNodeElementList();
328  ArrayView<const GO> domainElem2 = domainMap2->getNodeElementList();
329  ArrayView<const GO> rowElem1 = tmp_A_12->getRowMap()->getNodeElementList();
330  ArrayView<const GO> colElem1 = tmp_A_21->getColMap()->getNodeElementList();
331 
332  Xpetra::UnderlyingLib lib = domainMap2->lib();
333  GO indexBase = domainMap2->getIndexBase();
334 
335  Array<GO> newRowElem2(numRows2, 0);
336  for (Xpetra::global_size_t i = 0; i < numRows2; i++)
337  newRowElem2[i] = numVel + rangeElem2[i];
338 
339  RCP<const Map> newRangeMap2 = MapFactory::Build(lib, numRows2, newRowElem2, indexBase, comm);
340 
341  // maybe should be column map???
342  Array<GO> newColElem2(numCols2, 0);
343  for (Xpetra::global_size_t i = 0; i < numCols2; i++)
344  newColElem2[i] = numVel + domainElem2[i];
345 
346  RCP<const Map> newDomainMap2 = MapFactory::Build(lib, numCols2, newColElem2, indexBase, comm);
347 
348  RCP<Matrix> A_12 = MatrixFactory::Build(tmp_A_12->getRangeMap(), newDomainMap2, tmp_A_12->getNodeMaxNumRowEntries());
349  RCP<Matrix> A_21 = MatrixFactory::Build(newRangeMap2, tmp_A_21->getDomainMap(), tmp_A_21->getNodeMaxNumRowEntries());
350 
351  RCP<CrsMatrix> A_11_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_11) ->getCrsMatrix();
352  RCP<CrsMatrix> A_12_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_12) ->getCrsMatrix();
353  RCP<CrsMatrix> A_21_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_21) ->getCrsMatrix();
354  RCP<CrsMatrix> A_11_crs_9Pt = rcp_dynamic_cast<CrsMatrixWrap>(A_11_9Pt)->getCrsMatrix();
355 
356 #if 0
357  RCP<Matrix> A_22 = MatrixFactory::Build(newRangeMap2, newDomainMap2, 1);
358  RCP<CrsMatrix> A_22_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_22) ->getCrsMatrix();
359 
360  // FIXME: why do we need to perturb A_22?
361  Array<SC> smallVal(1, 1.0e-9);
362 
363  // FIXME: could this be sped up using expertStaticFillComplete?
364  // There was an attempt on doing it, but it did not do the proper thing
365  // with empty columns. See git history
366  ArrayView<const LO> inds;
367  ArrayView<const SC> vals;
368  for (LO row = 0; row < as<LO>(numRows2); ++row) {
369  tmp_A_21->getLocalRowView(row, inds, vals);
370 
371  size_t nnz = inds.size();
372  Array<GO> newInds(nnz, 0);
373  for (LO colID = 0; colID < as<LO>(nnz); colID++)
374  newInds[colID] = colElem1[inds[colID]];
375 
376  A_21_crs->insertGlobalValues(newRowElem2[row], newInds, vals);
377  A_22_crs->insertGlobalValues(newRowElem2[row], Array<LO>(1, newRowElem2[row]), smallVal);
378  }
379  A_21_crs->fillComplete(tmp_A_21->getDomainMap(), newRangeMap2);
380  A_22_crs->fillComplete(newDomainMap2, newRangeMap2);
381 #else
382  RCP<Matrix> A_22 = Teuchos::null;
383  RCP<CrsMatrix> A_22_crs = Teuchos::null;
384 
385  ArrayView<const LO> inds;
386  ArrayView<const SC> vals;
387  for (LO row = 0; row < as<LO>(numRows2); ++row) {
388  tmp_A_21->getLocalRowView(row, inds, vals);
389 
390  size_t nnz = inds.size();
391  Array<GO> newInds(nnz, 0);
392  for (LO colID = 0; colID < as<LO>(nnz); colID++)
393  newInds[colID] = colElem1[inds[colID]];
394 
395  A_21_crs->insertGlobalValues(newRowElem2[row], newInds, vals);
396  }
397  A_21_crs->fillComplete(tmp_A_21->getDomainMap(), newRangeMap2);
398 #endif
399 
400  // Create new A12 with map so that the global indices of the ColMap starts
401  // from numVel+1 (where numVel is the number of rows in the A11 block)
402  for (LO row = 0; row < as<LO>(tmp_A_12->getRowMap()->getNodeNumElements()); ++row) {
403  tmp_A_12->getLocalRowView(row, inds, vals);
404 
405  size_t nnz = inds.size();
406  Array<GO> newInds(nnz, 0);
407  for (LO colID = 0; colID < as<LO>(nnz); colID++)
408  newInds[colID] = newColElem2[inds[colID]];
409 
410  A_12_crs->insertGlobalValues(rowElem1[row], newInds, vals);
411  }
412  A_12_crs->fillComplete(newDomainMap2, tmp_A_12->getRangeMap());
413 
414  RCP<Matrix> A_12_abs = Absolute(*A_12);
415  RCP<Matrix> A_21_abs = Absolute(*A_21);
416 
417  // =========================================================================
418  // Preconditioner construction - I (block)
419  // =========================================================================
420  RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
421  Teuchos::FancyOStream& out = *fancy;
422  out.setOutputToRootOnly(0);
423  RCP<Matrix> BBt = Xpetra::MatrixMatrix<SC,LO,GO,NO>::Multiply(*A_21, false, *A_12, false, out);
424  RCP<Matrix> BBt_abs = Xpetra::MatrixMatrix<SC,LO,GO,NO>::Multiply(*A_21_abs, false, *A_12_abs, false, out);
425 
426  SC dropTol = (paramList.get<int>("useFilters") ? 0.06 : 0.00);
427  RCP<Matrix> filteredA = FilterMatrix(*A_11, *A_11, dropTol);
428  RCP<Matrix> filteredB = FilterMatrix(*BBt, *BBt_abs, dropTol);
429 
430  RCP<CrsMatrix> fA_11_crs = rcp_dynamic_cast<CrsMatrixWrap>(filteredA)->getCrsMatrix();
431  RCP<CrsMatrix> fA_12_crs = Teuchos::null;
432  RCP<CrsMatrix> fA_21_crs = Teuchos::null;
433  RCP<CrsMatrix> fA_22_crs = rcp_dynamic_cast<CrsMatrixWrap>(filteredB)->getCrsMatrix();
434 
435  // Build the large filtered matrix which requires strided maps
436  std::vector<size_t> stridingInfo(1, 1);
437  int stridedBlockId = -1;
438 
439  Array<GO> elementList(numVel+numPres); // Not RCP ... does this get cleared ?
440  Array<GO> velElem = A_12_crs->getRangeMap()->getNodeElementList();
441  Array<GO> presElem = A_21_crs->getRangeMap()->getNodeElementList();
442 
443  for (Xpetra::global_size_t i = 0 ; i < numVel; i++) elementList[i] = velElem[i];
444  for (Xpetra::global_size_t i = numVel; i < numVel+numPres; i++) elementList[i] = presElem[i-numVel];
445  RCP<const Map> fullMap = StridedMapFactory::Build(Xpetra::UseTpetra, numVel+numPres, elementList(), indexBase, stridingInfo, comm);
446 
447  std::vector<RCP<const Map> > partMaps(2);
448  partMaps[0] = StridedMapFactory::Build(Xpetra::UseTpetra, numVel, velElem, indexBase, stridingInfo, comm);
449  partMaps[1] = StridedMapFactory::Build(Xpetra::UseTpetra, numPres, presElem, indexBase, stridingInfo, comm, stridedBlockId, numVel);
450 
451  // Map extractors are necessary for Xpetra's block operators
452  RCP<const MapExtractor> mapExtractor = MapExtractorFactory::Build(fullMap, partMaps);
453  RCP<BlockedCrsMatrix> fA = rcp(new BlockedCrsMatrix(mapExtractor, mapExtractor, 10));
454  fA->setMatrix(0, 0, fA_11_crs);
455  fA->setMatrix(0, 1, fA_12_crs);
456  fA->setMatrix(1, 0, fA_21_crs);
457  fA->setMatrix(1, 1, fA_22_crs);
458  fA->fillComplete();
459 
460  // -------------------------------------------------------------------------
461  // Preconditioner construction - I.a (filtered hierarchy)
462  // -------------------------------------------------------------------------
464  SetDependencyTree(M, paramList);
465 
466  RCP<Hierarchy> H = rcp(new Hierarchy);
467  RCP<MueLu::Level> finestLevel = H->GetLevel(0);
468  finestLevel->Set("A", rcp_dynamic_cast<Matrix>(fA));
469  finestLevel->Set("p2vMap", p2vMap);
470  finestLevel->Set("CoordinatesVelocity", Xpetra::toXpetra(velCoords));
471  finestLevel->Set("CoordinatesPressure", Xpetra::toXpetra(presCoords));
472  finestLevel->Set("AForPat", A_11_9Pt);
473  H->SetMaxCoarseSize(MUELU_GPD("coarse: max size", int, 1));
474 
475 #ifdef IMPLICIT_TRANSPOSE
476  out << "Using implicit transpose" << std::endl;
477 
478  H->SetImplicitTranspose(true);
479 #endif
480 
481  // The first invocation of Setup() builds the hierarchy using the filtered
482  // matrix. This build includes the grid transfers but not the creation of the
483  // smoothers.
484  // NOTE: we need to indicate what should be kept from the first invocation
485  // for the second invocation, which then focuses on building the smoothers
486  // for the unfiltered matrix.
487  H->Keep("P", M.GetFactory("P") .get());
488  H->Keep("R", M.GetFactory("R") .get());
489  H->Keep("Ptent", M.GetFactory("Ptent").get());
490  H->Setup(M, 0, MUELU_GPD("max levels", int, 3));
491 
492 #if 0
493  for (int i = 1; i < H->GetNumLevels(); i++) {
494  RCP<Matrix> P = H->GetLevel(i)->template Get<RCP<Matrix> >("P");
495  RCP<BlockedCrsMatrix> Pcrs = rcp_dynamic_cast<BlockedCrsMatrix>(P);
496  RCP<CrsMatrix> Ppcrs = Pcrs->getMatrix(1,1);
497  RCP<Matrix> Pp = rcp(new CrsMatrixWrap(Ppcrs));
498  RCP<CrsMatrix> Pvcrs = Pcrs->getMatrix(0,0);
499  RCP<Matrix> Pv = rcp(new CrsMatrixWrap(Pvcrs));
500 
501  Xpetra::IO<SC,LO,GO,NO>::Write("Pp_l" + MueLu::toString(i) + ".mm", *Pp);
502  Xpetra::IO<SC,LO,GO,NO>::Write("Pv_l" + MueLu::toString(i) + ".mm", *Pv);
503  }
504 #endif
505 
506  // -------------------------------------------------------------------------
507  // Preconditioner construction - I.b (smoothers for unfiltered matrix)
508  // -------------------------------------------------------------------------
509  std::string smootherType = MUELU_GPD("smoother: type", std::string, "vanka");
510  ParameterList smootherParams;
511  if (paramList.isSublist("smoother: params"))
512  smootherParams = paramList.sublist("smoother: params");
513  M.SetFactory("Smoother", GetSmoother(smootherType, smootherParams, false/*coarseSolver?*/));
514 
515  std::string coarseType = MUELU_GPD("coarse: type", std::string, "direct");
516  ParameterList coarseParams;
517  if (paramList.isSublist("coarse: params"))
518  coarseParams = paramList.sublist("coarse: params");
519  M.SetFactory("CoarseSolver", GetSmoother(coarseType, coarseParams, true/*coarseSolver?*/));
520 
521 #ifdef HAVE_MUELU_DEBUG
522  M.ResetDebugData();
523 #endif
524 
525  RCP<BlockedCrsMatrix> A = rcp(new BlockedCrsMatrix(mapExtractor, mapExtractor, 10));
526  A->setMatrix(0, 0, A_11_crs);
527  A->setMatrix(0, 1, A_12_crs);
528  A->setMatrix(1, 0, A_21_crs);
529  A->setMatrix(1, 1, A_22_crs);
530  A->fillComplete();
531 
532  H->GetLevel(0)->Set("A", rcp_dynamic_cast<Matrix>(A));
533 
534  H->Setup(M, 0, H->GetNumLevels());
535 
536  return rcp(new MueLu::TpetraOperator<SC,LO,GO,NO>(H));
537  }
538 
539  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
540  RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
542  FilterMatrix(Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Pattern, Scalar dropTol) const {
543  typedef Xpetra::Matrix<SC,LO,GO,NO> Matrix;
544  typedef MueLu::AmalgamationFactory<SC,LO,GO,NO> AmalgamationFactory;
545  typedef MueLu::CoalesceDropFactory<SC,LO,GO,NO> CoalesceDropFactory;
546  typedef MueLu::FactoryManager<SC,LO,GO,NO> FactoryManager;
547  typedef MueLu::FilteredAFactory<SC,LO,GO,NO> FilteredAFactory;
548  typedef MueLu::GraphBase<LO,GO,NO> GraphBase;
549 
550  RCP<GraphBase> filteredGraph;
551  {
552  // Get graph pattern for the pattern matrix
553  MueLu::Level level;
554  level.SetLevelID(1);
555 
556  level.Set<RCP<Matrix> >("A", rcpFromRef(Pattern));
557 
558  RCP<AmalgamationFactory> amalgFactory = rcp(new AmalgamationFactory());
559 
560  RCP<CoalesceDropFactory> dropFactory = rcp(new CoalesceDropFactory());
561  ParameterList dropParams = *(dropFactory->GetValidParameterList());
562  dropParams.set("lightweight wrap", true);
563  dropParams.set("aggregation: drop scheme", "classical");
564  dropParams.set("aggregation: drop tol", dropTol);
565  // dropParams.set("Dirichlet detection threshold", <>);
566  dropFactory->SetParameterList(dropParams);
567  dropFactory->SetFactory("UnAmalgamationInfo", amalgFactory);
568 
569  // Build
570  level.Request("Graph", dropFactory.get());
571  dropFactory->Build(level);
572 
573  level.Get("Graph", filteredGraph, dropFactory.get());
574  }
575 
576  RCP<Matrix> filteredA;
577  {
578  // Filter the original matrix, not the pattern one
579  MueLu::Level level;
580  level.SetLevelID(1);
581 
582  level.Set("A", rcpFromRef(A));
583  level.Set("Graph", filteredGraph);
584  level.Set("Filtering", true);
585 
586  RCP<FilteredAFactory> filterFactory = rcp(new FilteredAFactory());
587  ParameterList filterParams = *(filterFactory->GetValidParameterList());
588  // We need a graph that has proper structure in it. Therefore, we need to
589  // drop older pattern, i.e. not to reuse it
590  filterParams.set("filtered matrix: reuse graph", false);
591  filterParams.set("filtered matrix: use lumping", false);
592  filterFactory->SetParameterList(filterParams);
593 
594  // Build
595  level.Request("A", filterFactory.get());
596  filterFactory->Build(level);
597 
598  level.Get("A", filteredA, filterFactory.get());
599  }
600 
601  // Zero out row sums by fixing the diagonal
602  filteredA->resumeFill();
603  size_t numRows = filteredA->getRowMap()->getNodeNumElements();
604  for (size_t i = 0; i < numRows; i++) {
605  ArrayView<const LO> inds;
606  ArrayView<const SC> vals;
607  filteredA->getLocalRowView(i, inds, vals);
608 
609  size_t nnz = inds.size();
610 
611  Array<SC> valsNew = vals;
612 
613  LO diagIndex = -1;
614  SC diag = Teuchos::ScalarTraits<SC>::zero();
615  for (size_t j = 0; j < inds.size(); j++) {
616  diag += vals[j];
617  if (inds[j] == i)
618  diagIndex = j;
619  }
620  TEUCHOS_TEST_FOR_EXCEPTION(diagIndex == -1, MueLu::Exceptions::RuntimeError,
621  "No diagonal found");
622  if (nnz <= 1)
623  continue;
624 
625  valsNew[diagIndex] -= diag;
626 
627  filteredA->replaceLocalValues(i, inds, valsNew);
628  }
629  filteredA->fillComplete();
630 
631  return filteredA;
632  }
633 
634  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
635  void
638  typedef MueLu::BlockedPFactory <SC,LO,GO,NO> BlockedPFactory;
639  typedef MueLu::GenericRFactory <SC,LO,GO,NO> GenericRFactory;
640  typedef MueLu::BlockedRAPFactory <SC,LO,GO,NO> BlockedRAPFactory;
641  typedef MueLu::SmootherFactory <SC,LO,GO,NO> SmootherFactory;
642  typedef MueLu::BlockedDirectSolver <SC,LO,GO,NO> BlockedDirectSolver;
643  typedef MueLu::FactoryManager <SC,LO,GO,NO> FactoryManager;
644 
645  // Pressure and velocity dependency trees are identical. The only
646  // difference is that pressure has to go first, so that velocity can use
647  // some of pressure data
648  RCP<FactoryManager> M11 = rcp(new FactoryManager()), M22 = rcp(new FactoryManager());
649  SetBlockDependencyTree(*M11, 0, 0, "velocity", paramList);
650  SetBlockDependencyTree(*M22, 1, 1, "pressure", paramList);
651 
652  RCP<BlockedPFactory> PFact = rcp(new BlockedPFactory());
653  ParameterList pParamList = *(PFact->GetValidParameterList());
654  pParamList.set("backwards", true); // do pressure first
655  PFact->SetParameterList(pParamList);
656  PFact->AddFactoryManager(M11);
657  PFact->AddFactoryManager(M22);
658  M.SetFactory("P", PFact);
659 
660  RCP<MueLu::Factory > AcFact = rcp(new BlockedRAPFactory());
661 #ifdef IMPLICIT_TRANSPOSE
662  M.SetFactory("R", Teuchos::null);
663 
664  ParameterList RAPparams;
665  RAPparams.set("transpose: use implicit", true);
666  AcFact->SetParameterList(RAPparams);
667 #else
668  RCP<GenericRFactory> RFact = rcp(new GenericRFactory());
669  RFact->SetFactory("P", PFact);
670  M.SetFactory("R", RFact);
671 
672  AcFact->SetFactory("R", RFact);
673 #endif
674  AcFact->SetFactory("P", PFact);
675  M.SetFactory("A", AcFact);
676 
677  // Smoothers will be set later
678  M.SetFactory("Smoother", Teuchos::null);
679 
680  RCP<MueLu::Factory> coarseFact = rcp(new SmootherFactory(rcp(new BlockedDirectSolver()), Teuchos::null));
681  // M.SetFactory("CoarseSolver", coarseFact);
682  M.SetFactory("CoarseSolver", Teuchos::null);
683  }
684 
685  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
686  void
688  SetBlockDependencyTree(MueLu::FactoryManager<Scalar,LocalOrdinal,GlobalOrdinal,Node>& M, LocalOrdinal row, LocalOrdinal col, const std::string& mode, const ParameterList& paramList) const {
689  typedef MueLu::ConstraintFactory <SC,LO,GO,NO> ConstraintFactory;
690  typedef MueLu::EminPFactory <SC,LO,GO,NO> EminPFactory;
691  typedef MueLu::PatternFactory <SC,LO,GO,NO> PatternFactory;
692  typedef MueLu::Q2Q1PFactory <SC,LO,GO,NO> Q2Q1PFactory;
693  typedef MueLu::Q2Q1uPFactory <SC,LO,GO,NO> Q2Q1uPFactory;
694  typedef MueLu::SubBlockAFactory <SC,LO,GO,NO> SubBlockAFactory;
695 
696  RCP<SubBlockAFactory> AFact = rcp(new SubBlockAFactory());
697  AFact->SetFactory ("A", MueLu::NoFactory::getRCP());
698  AFact->SetParameter("block row", Teuchos::ParameterEntry(row));
699  AFact->SetParameter("block col", Teuchos::ParameterEntry(col));
700  M.SetFactory("A", AFact);
701 
702  RCP<MueLu::Factory> Q2Q1Fact;
703 
704  const bool isStructured = false;
705 
706  if (isStructured) {
707  Q2Q1Fact = rcp(new Q2Q1PFactory);
708 
709  } else {
710  Q2Q1Fact = rcp(new Q2Q1uPFactory);
711  ParameterList q2q1ParamList = *(Q2Q1Fact->GetValidParameterList());
712  q2q1ParamList.set("mode", mode);
713  if (paramList.isParameter("dump status"))
714  q2q1ParamList.set("dump status", paramList.get<bool>("dump status"));
715  if (paramList.isParameter("phase2"))
716  q2q1ParamList.set("phase2", paramList.get<bool>("phase2"));
717  Q2Q1Fact->SetParameterList(q2q1ParamList);
718  }
719  Q2Q1Fact->SetFactory("A", AFact);
720  M.SetFactory("Ptent", Q2Q1Fact);
721 
722  RCP<PatternFactory> patternFact = rcp(new PatternFactory);
723  ParameterList patternParams = *(patternFact->GetValidParameterList());
724  // Our prolongator constructs the exact pattern we are going to use,
725  // therefore we do not expand it
726  patternParams.set("emin: pattern order", 0);
727  patternFact->SetParameterList(patternParams);
728  patternFact->SetFactory("A", AFact);
729  patternFact->SetFactory("P", Q2Q1Fact);
730  M.SetFactory("Ppattern", patternFact);
731 
732  RCP<ConstraintFactory> CFact = rcp(new ConstraintFactory);
733  CFact->SetFactory("Ppattern", patternFact);
734  M.SetFactory("Constraint", CFact);
735 
736  RCP<EminPFactory> EminPFact = rcp(new EminPFactory());
737  ParameterList eminParams = *(EminPFact->GetValidParameterList());
738  if (paramList.isParameter("emin: num iterations"))
739  eminParams.set("emin: num iterations", paramList.get<int>("emin: num iterations"));
740  EminPFact->SetParameterList(eminParams);
741  EminPFact->SetFactory("A", AFact);
742  EminPFact->SetFactory("Constraint", CFact);
743  EminPFact->SetFactory("P", Q2Q1Fact);
744  M.SetFactory("P", EminPFact);
745  }
746 
747  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
748  RCP<MueLu::FactoryBase>
749  MueLuTpetraQ2Q1PreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::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 == "direct") {
786  smootherPrototype = rcp(new BlockedDirectSolver());
787 
788  } else if (type == "braess-sarazin") {
789  // Define smoother/solver for BraessSarazin
790  SC omega = MUELU_GPD("bs: omega", double, 1.0);
791  bool lumping = MUELU_GPD("bs: lumping", bool, false);
792 
793  RCP<SchurComplementFactory> schurFact = rcp(new SchurComplementFactory());
794  schurFact->SetParameter("omega", ParameterEntry(omega));
795  schurFact->SetParameter("lumping", ParameterEntry(lumping));
796  schurFact->SetFactory ("A", MueLu::NoFactory::getRCP());
797 
798  // Schur complement solver
799  RCP<SmootherPrototype> schurSmootherPrototype;
800  std::string schurSmootherType = (paramList.isParameter("schur smoother: type") ? paramList.get<std::string>("schur smoother: type") : "RELAXATION");
801  if (schurSmootherType == "RELAXATION") {
802  ParameterList schurSmootherParams = paramList.sublist("schur smoother: params");
803  // schurSmootherParams.set("relaxation: damping factor", omega);
804  schurSmootherPrototype = rcp(new TrilinosSmoother(schurSmootherType, schurSmootherParams));
805  } else {
806  schurSmootherPrototype = rcp(new DirectSolver());
807  }
808  schurSmootherPrototype->SetFactory("A", schurFact);
809 
810  RCP<SmootherFactory> schurSmootherFact = rcp(new SmootherFactory(schurSmootherPrototype));
811 
812  // Define temporary FactoryManager that is used as input for BraessSarazin smoother
813  RCP<FactoryManager> braessManager = rcp(new FactoryManager());
814  braessManager->SetFactory("A", schurFact); // SchurComplement operator for correction step (defined as "A")
815  braessManager->SetFactory("Smoother", schurSmootherFact); // solver/smoother for correction step
816  braessManager->SetFactory("PreSmoother", schurSmootherFact);
817  braessManager->SetFactory("PostSmoother", schurSmootherFact);
818  braessManager->SetIgnoreUserData(true); // always use data from factories defined in factory manager
819 
820  smootherPrototype = rcp(new BraessSarazinSmoother());
821  smootherPrototype->SetParameter("Sweeps", ParameterEntry(MUELU_GPD("bs: sweeps", int, 1)));
822  smootherPrototype->SetParameter("lumping", ParameterEntry(lumping));
823  smootherPrototype->SetParameter("Damping factor", ParameterEntry(omega));
824  rcp_dynamic_cast<BraessSarazinSmoother>(smootherPrototype)->AddFactoryManager(braessManager, 0); // set temporary factory manager in BraessSarazin smoother
825  }
826 
827  return coarseSolver ? rcp(new SmootherFactory(smootherPrototype, Teuchos::null)) : rcp(new SmootherFactory(smootherPrototype));
828  }
829 
830  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
831  RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
832  MueLuTpetraQ2Q1PreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Absolute(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A) const {
833  typedef Xpetra::CrsMatrix <SC,LO,GO,NO> CrsMatrix;
834  typedef Xpetra::CrsMatrixWrap<SC,LO,GO,NO> CrsMatrixWrap;
835  typedef Xpetra::Matrix <SC,LO,GO,NO> Matrix;
836 
837  const CrsMatrixWrap& Awrap = dynamic_cast<const CrsMatrixWrap&>(A);
838 
839  ArrayRCP<const size_t> iaA;
840  ArrayRCP<const LO> jaA;
841  ArrayRCP<const SC> valA;
842  Awrap.getCrsMatrix()->getAllValues(iaA, jaA, valA);
843 
844  ArrayRCP<size_t> iaB (iaA .size());
845  ArrayRCP<LO> jaB (jaA .size());
846  ArrayRCP<SC> valB(valA.size());
847  for (int i = 0; i < iaA .size(); i++) iaB [i] = iaA[i];
848  for (int i = 0; i < jaA .size(); i++) jaB [i] = jaA[i];
849  for (int i = 0; i < valA.size(); i++) valB[i] = Teuchos::ScalarTraits<SC>::magnitude(valA[i]);
850 
851  RCP<Matrix> B = rcp(new CrsMatrixWrap(A.getRowMap(), A.getColMap(), 0, Xpetra::StaticProfile));
852  RCP<CrsMatrix> Bcrs = rcp_dynamic_cast<CrsMatrixWrap>(B)->getCrsMatrix();
853  Bcrs->setAllValues(iaB, jaB, valB);
854  Bcrs->expertStaticFillComplete(A.getDomainMap(), A.getRangeMap());
855 
856  return B;
857  }
858 
859  // Public functions overridden from Teuchos::Describable
860  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
862  return "Thyra::MueLuTpetraQ2Q1PreconditionerFactory";
863  }
864 
865 } // namespace Thyra
866 
867 #endif
868 #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...
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->Get< RCP<Matrix> >("A", factory) if factory == NULL => 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.
Class that encapsulates external library smoothers.
Factory for building coarse matrices.
Teuchos::RCP< const Teuchos::ParameterList > getParameterList() const
void uninitializePrec(PreconditionerBase< SC > *prec, Teuchos::RCP< const LinearOpSourceBase< SC > > *fwdOp, ESupportSolveUse *supportSolveUse) const
Base class for smoother prototypes.
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.
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
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.
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
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.
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)
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
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.