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 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;
112  using Teuchos::ParameterList;
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>
139  RCP<PreconditionerBase<Scalar> >
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();
154  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
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));
209  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
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>
228  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
229  paramList_ = paramList;
230  }
231 
232  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
233  RCP<ParameterList>
235  return paramList_;
236  }
237 
238 
239  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
240  RCP<ParameterList>
242  RCP<ParameterList> savedParamList = paramList_;
243  paramList_ = Teuchos::null;
244  return savedParamList;
245  }
246 
247  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
248  RCP<const ParameterList>
250  return paramList_;
251  }
252 
253  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
254  RCP<const ParameterList>
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>
265  RCP<MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,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 
309  RCP<Matrix> A_11 = MueLu::TpetraCrs_To_XpetraMatrix(TpetCrsA11);
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 
524  return rcp(new MueLu::TpetraOperator<SC,LO,GO,NO>(H));
525  }
526 
527  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
528  RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,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;
601  SC diag = Teuchos::ScalarTraits<SC>::zero();
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  }
607  TEUCHOS_TEST_FOR_EXCEPTION(diagIndex == -1, MueLu::Exceptions::RuntimeError,
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
669  SetBlockDependencyTree(MueLu::FactoryManager<Scalar,LocalOrdinal,GlobalOrdinal,Node>& M, LocalOrdinal row, LocalOrdinal col, const std::string& mode, const ParameterList& paramList) const {
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>
747  RCP<MueLu::FactoryBase>
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>
846  RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,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 
854  ArrayRCP<const size_t> iaA;
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.
This class specifies the default factory that should generate some data on a Level if the data does n...
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
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->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.
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
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &paramList)
Factory for building blocked, segregated prolongation operators.
bool isCompatible(const LinearOpSourceBase< SC > &fwdOp) const
Class that encapsulates external library smoothers.
Factory for building coarse matrices.
Base class for smoother prototypes.
Class that encapsulates direct solvers. Autoselection of AmesosSmoother or Amesos2Smoother according ...
void SetBlockDependencyTree(MueLu::FactoryManager< SC, LO, GO, NO > &M, LO row, LO col, const std::string &mode, const ParameterList &paramList) const
Factory for building restriction operators using a prolongator factory.
MueLu::DefaultScalar Scalar
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
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
Teuchos::RCP< Xpetra::Matrix< SC, LO, GO, NO > > Absolute(const Xpetra::Matrix< SC, LO, GO, NO > &A) 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.
const RCP< const FactoryBase > GetFactory(const std::string &varName) const
Get factory associated with a particular data name.
MueLu representation of a graph.
Factory for building the Schur Complement for a 2x2 block matrix.
Factory for creating a graph base on a given matrix.
RCP< MueLu::FactoryBase > GetSmoother(const std::string &type, const ParameterList &paramList, bool coarseSolver) 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.
void SetDependencyTree(MueLu::FactoryManager< SC, LO, GO, NO > &M, const ParameterList &paramList) const
Teuchos::RCP< const Teuchos::ParameterList > getParameterList() const
Exception throws to report errors in the internal logical of the program.
Factory for building filtered matrices using filtered graphs.
void initializePrec(const Teuchos::RCP< const LinearOpSourceBase< SC > > &fwdOp, PreconditionerBase< SC > *prec, const ESupportSolveUse supportSolveUse) const
#define MUELU_GPD(name, type, defaultValue)
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
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.
void uninitializePrec(PreconditionerBase< SC > *prec, Teuchos::RCP< const LinearOpSourceBase< SC > > *fwdOp, ESupportSolveUse *supportSolveUse) const
static const RCP< const NoFactory > getRCP()
Static Get() functions.