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,
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::MapExtractor <SC,LO,GO,NO> MapExtractor;
285  typedef Xpetra::MapExtractorFactory <SC,LO,GO,NO> MapExtractorFactory;
286  typedef Xpetra::MapFactory <LO,GO,NO> MapFactory;
287  typedef Xpetra::Matrix <SC,LO,GO,NO> Matrix;
288  typedef Xpetra::MatrixFactory <SC,LO,GO,NO> MatrixFactory;
289  typedef Xpetra::StridedMapFactory <LO,GO,NO> StridedMapFactory;
290 
291  typedef MueLu::Hierarchy <SC,LO,GO,NO> Hierarchy;
292 
293  const RCP<const Teuchos::Comm<int> > comm = velCoords->getMap()->getComm();
294 
295  // Pull out Tpetra matrices
296  RCP<Thyra::LinearOpBase<SC> > ThNonConstA11 = rcp_const_cast<Thyra::LinearOpBase<double> >(thA11);
297  RCP<Thyra::LinearOpBase<SC> > ThNonConstA21 = rcp_const_cast<Thyra::LinearOpBase<double> >(thA21);
298  RCP<Thyra::LinearOpBase<SC> > ThNonConstA12 = rcp_const_cast<Thyra::LinearOpBase<double> >(thA12);
299  RCP<Thyra::LinearOpBase<SC> > ThNonConstA11_9Pt = rcp_const_cast<Thyra::LinearOpBase<double> >(thA11_9Pt);
300 
301  RCP<TP_Op> TpetA11 = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA11);
302  RCP<TP_Op> TpetA21 = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA21);
303  RCP<TP_Op> TpetA12 = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA12);
304  RCP<TP_Op> TpetA11_9Pt = Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO>::getTpetraOperator(ThNonConstA11_9Pt);
305 
306  RCP<TP_Crs> TpetCrsA11 = rcp_dynamic_cast<TP_Crs>(TpetA11);
307  RCP<TP_Crs> TpetCrsA21 = rcp_dynamic_cast<TP_Crs>(TpetA21);
308  RCP<TP_Crs> TpetCrsA12 = rcp_dynamic_cast<TP_Crs>(TpetA12);
309  RCP<TP_Crs> TpetCrsA11_9Pt = rcp_dynamic_cast<TP_Crs>(TpetA11_9Pt);
310 
311  RCP<Matrix> A_11 = MueLu::TpetraCrs_To_XpetraMatrix(TpetCrsA11);
312  RCP<Matrix> tmp_A_21 = MueLu::TpetraCrs_To_XpetraMatrix(TpetCrsA21); // needs map modification
313  RCP<Matrix> tmp_A_12 = MueLu::TpetraCrs_To_XpetraMatrix(TpetCrsA12); // needs map modification
314  RCP<Matrix> A_11_9Pt = MueLu::TpetraCrs_To_XpetraMatrix(TpetCrsA11_9Pt);
315 
316  Xpetra::global_size_t numVel = A_11->getRowMap()->getNodeNumElements();
317  Xpetra::global_size_t numPres = tmp_A_21->getRowMap()->getNodeNumElements();
318 
319  // Create new A21 with map so that the global indices of the row map starts
320  // from numVel+1 (where numVel is the number of rows in the A11 block)
321  RCP<const Map> domainMap2 = tmp_A_12->getDomainMap();
322  RCP<const Map> rangeMap2 = tmp_A_21->getRangeMap();
323  Xpetra::global_size_t numRows2 = rangeMap2->getNodeNumElements();
324  Xpetra::global_size_t numCols2 = domainMap2->getNodeNumElements();
325  ArrayView<const GO> rangeElem2 = rangeMap2->getNodeElementList();
326  ArrayView<const GO> domainElem2 = domainMap2->getNodeElementList();
327  ArrayView<const GO> rowElem1 = tmp_A_12->getRowMap()->getNodeElementList();
328  ArrayView<const GO> colElem1 = tmp_A_21->getColMap()->getNodeElementList();
329 
330  Xpetra::UnderlyingLib lib = domainMap2->lib();
331  GO indexBase = domainMap2->getIndexBase();
332 
333  Array<GO> newRowElem2(numRows2, 0);
334  for (Xpetra::global_size_t i = 0; i < numRows2; i++)
335  newRowElem2[i] = numVel + rangeElem2[i];
336 
337  RCP<const Map> newRangeMap2 = MapFactory::Build(lib, numRows2, newRowElem2, indexBase, comm);
338 
339  // maybe should be column map???
340  Array<GO> newColElem2(numCols2, 0);
341  for (Xpetra::global_size_t i = 0; i < numCols2; i++)
342  newColElem2[i] = numVel + domainElem2[i];
343 
344  RCP<const Map> newDomainMap2 = MapFactory::Build(lib, numCols2, newColElem2, indexBase, comm);
345 
346  RCP<Matrix> A_12 = MatrixFactory::Build(tmp_A_12->getRangeMap(), newDomainMap2, tmp_A_12->getNodeMaxNumRowEntries());
347  RCP<Matrix> A_21 = MatrixFactory::Build(newRangeMap2, tmp_A_21->getDomainMap(), tmp_A_21->getNodeMaxNumRowEntries());
348 
349  RCP<CrsMatrix> A_11_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_11) ->getCrsMatrix();
350  RCP<CrsMatrix> A_12_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_12) ->getCrsMatrix();
351  RCP<CrsMatrix> A_21_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_21) ->getCrsMatrix();
352  RCP<CrsMatrix> A_11_crs_9Pt = rcp_dynamic_cast<CrsMatrixWrap>(A_11_9Pt)->getCrsMatrix();
353 
354 #if 0
355  RCP<Matrix> A_22 = MatrixFactory::Build(newRangeMap2, newDomainMap2, 1);
356  RCP<CrsMatrix> A_22_crs = rcp_dynamic_cast<CrsMatrixWrap>(A_22) ->getCrsMatrix();
357 
358  // FIXME: why do we need to perturb A_22?
359  Array<SC> smallVal(1, 1.0e-10);
360 
361  // FIXME: could this be sped up using expertStaticFillComplete?
362  // There was an attempt on doing it, but it did not do the proper thing
363  // with empty columns. See git history
364  ArrayView<const LO> inds;
365  ArrayView<const SC> vals;
366  for (LO row = 0; row < as<LO>(numRows2); ++row) {
367  tmp_A_21->getLocalRowView(row, inds, vals);
368 
369  size_t nnz = inds.size();
370  Array<GO> newInds(nnz, 0);
371  for (LO colID = 0; colID < as<LO>(nnz); colID++)
372  newInds[colID] = colElem1[inds[colID]];
373 
374  A_21_crs->insertGlobalValues(newRowElem2[row], newInds, vals);
375  A_22_crs->insertGlobalValues(newRowElem2[row], Array<LO>(1, newRowElem2[row]), smallVal);
376  }
377  A_21_crs->fillComplete(tmp_A_21->getDomainMap(), newRangeMap2);
378  A_22_crs->fillComplete(newDomainMap2, newRangeMap2);
379 #else
380  RCP<Matrix> A_22 = Teuchos::null;
381  RCP<CrsMatrix> A_22_crs = Teuchos::null;
382 
383  ArrayView<const LO> inds;
384  ArrayView<const SC> vals;
385  for (LO row = 0; row < as<LO>(numRows2); ++row) {
386  tmp_A_21->getLocalRowView(row, inds, vals);
387 
388  size_t nnz = inds.size();
389  Array<GO> newInds(nnz, 0);
390  for (LO colID = 0; colID < as<LO>(nnz); colID++)
391  newInds[colID] = colElem1[inds[colID]];
392 
393  A_21_crs->insertGlobalValues(newRowElem2[row], newInds, vals);
394  }
395  A_21_crs->fillComplete(tmp_A_21->getDomainMap(), newRangeMap2);
396 #endif
397 
398  // Create new A12 with map so that the global indices of the ColMap starts
399  // from numVel+1 (where numVel is the number of rows in the A11 block)
400  for (LO row = 0; row < as<LO>(tmp_A_12->getRowMap()->getNodeNumElements()); ++row) {
401  tmp_A_12->getLocalRowView(row, inds, vals);
402 
403  size_t nnz = inds.size();
404  Array<GO> newInds(nnz, 0);
405  for (LO colID = 0; colID < as<LO>(nnz); colID++)
406  newInds[colID] = newColElem2[inds[colID]];
407 
408  A_12_crs->insertGlobalValues(rowElem1[row], newInds, vals);
409  }
410  A_12_crs->fillComplete(newDomainMap2, tmp_A_12->getRangeMap());
411 
412  RCP<Matrix> A_12_abs = Absolute(*A_12);
413  RCP<Matrix> A_21_abs = Absolute(*A_21);
414 
415  // =========================================================================
416  // Preconditioner construction - I (block)
417  // =========================================================================
418  RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
419  Teuchos::FancyOStream& out = *fancy;
420  out.setOutputToRootOnly(0);
421  RCP<Matrix> BBt = Xpetra::MatrixMatrix<SC,LO,GO,NO>::Multiply(*A_21, false, *A_12, false, out);
422  RCP<Matrix> BBt_abs = Xpetra::MatrixMatrix<SC,LO,GO,NO>::Multiply(*A_21_abs, false, *A_12_abs, false, out);
423 
424  SC dropTol = (paramList.get<int>("useFilters") ? 0.06 : 0.00);
425  RCP<Matrix> filteredA = FilterMatrix(*A_11, *A_11, dropTol);
426  RCP<Matrix> filteredB = FilterMatrix(*BBt, *BBt_abs, dropTol);
427 
428  RCP<Matrix> fA_11_crs = rcp_dynamic_cast<CrsMatrixWrap>(filteredA);
429  RCP<Matrix> fA_12_crs = Teuchos::null;
430  RCP<Matrix> fA_21_crs = Teuchos::null;
431  RCP<Matrix> fA_22_crs = rcp_dynamic_cast<CrsMatrixWrap>(filteredB);
432 
433  // Build the large filtered matrix which requires strided maps
434  std::vector<size_t> stridingInfo(1, 1);
435  int stridedBlockId = -1;
436 
437  Array<GO> elementList(numVel+numPres); // Not RCP ... does this get cleared ?
438  Array<GO> velElem = A_12_crs->getRangeMap()->getNodeElementList();
439  Array<GO> presElem = A_21_crs->getRangeMap()->getNodeElementList();
440 
441  for (Xpetra::global_size_t i = 0 ; i < numVel; i++) elementList[i] = velElem[i];
442  for (Xpetra::global_size_t i = numVel; i < numVel+numPres; i++) elementList[i] = presElem[i-numVel];
443  RCP<const Map> fullMap = StridedMapFactory::Build(Xpetra::UseTpetra, numVel+numPres, elementList(), indexBase, stridingInfo, comm);
444 
445  std::vector<RCP<const Map> > partMaps(2);
446  partMaps[0] = StridedMapFactory::Build(Xpetra::UseTpetra, numVel, velElem, indexBase, stridingInfo, comm);
447  partMaps[1] = StridedMapFactory::Build(Xpetra::UseTpetra, numPres, presElem, indexBase, stridingInfo, comm, stridedBlockId, numVel);
448 
449  // Map extractors are necessary for Xpetra's block operators
450  RCP<const MapExtractor> mapExtractor = MapExtractorFactory::Build(fullMap, partMaps);
451  RCP<BlockedCrsMatrix> fA = rcp(new BlockedCrsMatrix(mapExtractor, mapExtractor, 10));
452  fA->setMatrix(0, 0, fA_11_crs);
453  fA->setMatrix(0, 1, fA_12_crs);
454  fA->setMatrix(1, 0, fA_21_crs);
455  fA->setMatrix(1, 1, fA_22_crs);
456  fA->fillComplete();
457 
458  // -------------------------------------------------------------------------
459  // Preconditioner construction - I.a (filtered hierarchy)
460  // -------------------------------------------------------------------------
462  SetDependencyTree(M, paramList);
463 
464  RCP<Hierarchy> H = rcp(new Hierarchy);
465  RCP<MueLu::Level> finestLevel = H->GetLevel(0);
466  finestLevel->Set("A", rcp_dynamic_cast<Matrix>(fA));
467  finestLevel->Set("p2vMap", p2vMap);
468  finestLevel->Set("CoordinatesVelocity", Xpetra::toXpetra(velCoords));
469  finestLevel->Set("CoordinatesPressure", Xpetra::toXpetra(presCoords));
470  finestLevel->Set("AForPat", A_11_9Pt);
471  H->SetMaxCoarseSize(MUELU_GPD("coarse: max size", int, 1));
472 
473  // The first invocation of Setup() builds the hierarchy using the filtered
474  // matrix. This build includes the grid transfers but not the creation of the
475  // smoothers.
476  // NOTE: we need to indicate what should be kept from the first invocation
477  // for the second invocation, which then focuses on building the smoothers
478  // for the unfiltered matrix.
479  H->Keep("P", M.GetFactory("P") .get());
480  H->Keep("R", M.GetFactory("R") .get());
481  H->Keep("Ptent", M.GetFactory("Ptent").get());
482  H->Setup(M, 0, MUELU_GPD("max levels", int, 3));
483 
484 #if 0
485  for (int i = 1; i < H->GetNumLevels(); i++) {
486  RCP<Matrix> P = H->GetLevel(i)->template Get<RCP<Matrix> >("P");
487  RCP<BlockedCrsMatrix> Pcrs = rcp_dynamic_cast<BlockedCrsMatrix>(P);
488  RCP<Matrix> Pp = Pcrs->getMatrix(1,1);
489  RCP<Matrix> Pv = Pcrs->getMatrix(0,0);
490 
491  Xpetra::IO<SC,LO,GO,NO>::Write("Pp_l" + MueLu::toString(i) + ".mm", *Pp);
492  Xpetra::IO<SC,LO,GO,NO>::Write("Pv_l" + MueLu::toString(i) + ".mm", *Pv);
493  }
494 #endif
495 
496  // -------------------------------------------------------------------------
497  // Preconditioner construction - I.b (smoothers for unfiltered matrix)
498  // -------------------------------------------------------------------------
499  std::string smootherType = MUELU_GPD("smoother: type", std::string, "vanka");
500  ParameterList smootherParams;
501  if (paramList.isSublist("smoother: params"))
502  smootherParams = paramList.sublist("smoother: params");
503  M.SetFactory("Smoother", GetSmoother(smootherType, smootherParams, false/*coarseSolver?*/));
504 
505  std::string coarseType = MUELU_GPD("coarse: type", std::string, "direct");
506  ParameterList coarseParams;
507  if (paramList.isSublist("coarse: params"))
508  coarseParams = paramList.sublist("coarse: params");
509  M.SetFactory("CoarseSolver", GetSmoother(coarseType, coarseParams, true/*coarseSolver?*/));
510 
511 #ifdef HAVE_MUELU_DEBUG
512  M.ResetDebugData();
513 #endif
514 
515  RCP<BlockedCrsMatrix> A = rcp(new BlockedCrsMatrix(mapExtractor, mapExtractor, 10));
516  A->setMatrix(0, 0, A_11);
517  A->setMatrix(0, 1, A_12);
518  A->setMatrix(1, 0, A_21);
519  A->setMatrix(1, 1, A_22);
520  A->fillComplete();
521 
522  H->GetLevel(0)->Set("A", rcp_dynamic_cast<Matrix>(A));
523 
524  H->Setup(M, 0, H->GetNumLevels());
525 
526  return rcp(new MueLu::TpetraOperator<SC,LO,GO,NO>(H));
527  }
528 
529  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
530  RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
532  FilterMatrix(Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Pattern, Scalar dropTol) const {
533  typedef Xpetra::Matrix<SC,LO,GO,NO> Matrix;
534  typedef MueLu::AmalgamationFactory<SC,LO,GO,NO> AmalgamationFactory;
535  typedef MueLu::CoalesceDropFactory<SC,LO,GO,NO> CoalesceDropFactory;
536  typedef MueLu::FactoryManager<SC,LO,GO,NO> FactoryManager;
537  typedef MueLu::FilteredAFactory<SC,LO,GO,NO> FilteredAFactory;
538  typedef MueLu::GraphBase<LO,GO,NO> GraphBase;
539 
540  RCP<GraphBase> filteredGraph;
541  {
542  // Get graph pattern for the pattern matrix
543  MueLu::Level level;
544  level.SetLevelID(1);
545 
546  level.Set<RCP<Matrix> >("A", rcpFromRef(Pattern));
547 
548  RCP<AmalgamationFactory> amalgFactory = rcp(new AmalgamationFactory());
549 
550  RCP<CoalesceDropFactory> dropFactory = rcp(new CoalesceDropFactory());
551  ParameterList dropParams = *(dropFactory->GetValidParameterList());
552  dropParams.set("lightweight wrap", true);
553  dropParams.set("aggregation: drop scheme", "classical");
554  dropParams.set("aggregation: drop tol", dropTol);
555  // dropParams.set("Dirichlet detection threshold", <>);
556  dropFactory->SetParameterList(dropParams);
557  dropFactory->SetFactory("UnAmalgamationInfo", amalgFactory);
558 
559  // Build
560  level.Request("Graph", dropFactory.get());
561  dropFactory->Build(level);
562 
563  level.Get("Graph", filteredGraph, dropFactory.get());
564  }
565 
566  RCP<Matrix> filteredA;
567  {
568  // Filter the original matrix, not the pattern one
569  MueLu::Level level;
570  level.SetLevelID(1);
571 
572  level.Set("A", rcpFromRef(A));
573  level.Set("Graph", filteredGraph);
574  level.Set("Filtering", true);
575 
576  RCP<FilteredAFactory> filterFactory = rcp(new FilteredAFactory());
577  ParameterList filterParams = *(filterFactory->GetValidParameterList());
578  // We need a graph that has proper structure in it. Therefore, we need to
579  // drop older pattern, i.e. not to reuse it
580  filterParams.set("filtered matrix: reuse graph", false);
581  filterParams.set("filtered matrix: use lumping", false);
582  filterFactory->SetParameterList(filterParams);
583 
584  // Build
585  level.Request("A", filterFactory.get());
586  filterFactory->Build(level);
587 
588  level.Get("A", filteredA, filterFactory.get());
589  }
590 
591  // Zero out row sums by fixing the diagonal
592  filteredA->resumeFill();
593  size_t numRows = filteredA->getRowMap()->getNodeNumElements();
594  for (size_t i = 0; i < numRows; i++) {
595  ArrayView<const LO> inds;
596  ArrayView<const SC> vals;
597  filteredA->getLocalRowView(i, inds, vals);
598 
599  size_t nnz = inds.size();
600 
601  Array<SC> valsNew = vals;
602 
603  LO diagIndex = -1;
604  SC diag = Teuchos::ScalarTraits<SC>::zero();
605  for (size_t j = 0; j < inds.size(); j++) {
606  diag += vals[j];
607  if (inds[j] == i)
608  diagIndex = j;
609  }
610  TEUCHOS_TEST_FOR_EXCEPTION(diagIndex == -1, MueLu::Exceptions::RuntimeError,
611  "No diagonal found");
612  if (nnz <= 1)
613  continue;
614 
615  valsNew[diagIndex] -= diag;
616 
617  filteredA->replaceLocalValues(i, inds, valsNew);
618  }
619  filteredA->fillComplete();
620 
621  return filteredA;
622  }
623 
624  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
625  void
628  typedef MueLu::BlockedPFactory <SC,LO,GO,NO> BlockedPFactory;
629  typedef MueLu::GenericRFactory <SC,LO,GO,NO> GenericRFactory;
630  typedef MueLu::BlockedRAPFactory <SC,LO,GO,NO> BlockedRAPFactory;
631  typedef MueLu::SmootherFactory <SC,LO,GO,NO> SmootherFactory;
632  typedef MueLu::BlockedDirectSolver <SC,LO,GO,NO> BlockedDirectSolver;
633  typedef MueLu::FactoryManager <SC,LO,GO,NO> FactoryManager;
634 
635  // Pressure and velocity dependency trees are identical. The only
636  // difference is that pressure has to go first, so that velocity can use
637  // some of pressure data
638  RCP<FactoryManager> M11 = rcp(new FactoryManager()), M22 = rcp(new FactoryManager());
639  SetBlockDependencyTree(*M11, 0, 0, "velocity", paramList);
640  SetBlockDependencyTree(*M22, 1, 1, "pressure", paramList);
641 
642  RCP<BlockedPFactory> PFact = rcp(new BlockedPFactory());
643  ParameterList pParamList = *(PFact->GetValidParameterList());
644  pParamList.set("backwards", true); // do pressure first
645  PFact->SetParameterList(pParamList);
646  PFact->AddFactoryManager(M11);
647  PFact->AddFactoryManager(M22);
648  M.SetFactory("P", PFact);
649 
650  RCP<GenericRFactory> RFact = rcp(new GenericRFactory());
651  RFact->SetFactory("P", PFact);
652  M.SetFactory("R", RFact);
653 
654  RCP<MueLu::Factory > AcFact = rcp(new BlockedRAPFactory());
655  AcFact->SetFactory("R", RFact);
656  AcFact->SetFactory("P", PFact);
657  M.SetFactory("A", AcFact);
658 
659  // Smoothers will be set later
660  M.SetFactory("Smoother", Teuchos::null);
661 
662  RCP<MueLu::Factory> coarseFact = rcp(new SmootherFactory(rcp(new BlockedDirectSolver()), Teuchos::null));
663  // M.SetFactory("CoarseSolver", coarseFact);
664  M.SetFactory("CoarseSolver", Teuchos::null);
665  }
666 
667  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
668  void
670  SetBlockDependencyTree(MueLu::FactoryManager<Scalar,LocalOrdinal,GlobalOrdinal,Node>& M, LocalOrdinal row, LocalOrdinal col, const std::string& mode, const ParameterList& paramList) const {
671  typedef MueLu::ConstraintFactory <SC,LO,GO,NO> ConstraintFactory;
672  typedef MueLu::EminPFactory <SC,LO,GO,NO> EminPFactory;
673  typedef MueLu::GenericRFactory <SC,LO,GO,NO> GenericRFactory;
674  typedef MueLu::PatternFactory <SC,LO,GO,NO> PatternFactory;
675  typedef MueLu::Q2Q1PFactory <SC,LO,GO,NO> Q2Q1PFactory;
676  typedef MueLu::Q2Q1uPFactory <SC,LO,GO,NO> Q2Q1uPFactory;
677  typedef MueLu::SubBlockAFactory <SC,LO,GO,NO> SubBlockAFactory;
678 
679  RCP<SubBlockAFactory> AFact = rcp(new SubBlockAFactory());
680  AFact->SetFactory ("A", MueLu::NoFactory::getRCP());
681  AFact->SetParameter("block row", Teuchos::ParameterEntry(row));
682  AFact->SetParameter("block col", Teuchos::ParameterEntry(col));
683  M.SetFactory("A", AFact);
684 
685  RCP<MueLu::Factory> Q2Q1Fact;
686 
687  const bool isStructured = false;
688 
689  if (isStructured) {
690  Q2Q1Fact = rcp(new Q2Q1PFactory);
691 
692  } else {
693  Q2Q1Fact = rcp(new Q2Q1uPFactory);
694  ParameterList q2q1ParamList = *(Q2Q1Fact->GetValidParameterList());
695  q2q1ParamList.set("mode", mode);
696  if (paramList.isParameter("dump status"))
697  q2q1ParamList.set("dump status", paramList.get<bool>("dump status"));
698  if (paramList.isParameter("phase2"))
699  q2q1ParamList.set("phase2", paramList.get<bool>("phase2"));
700  Q2Q1Fact->SetParameterList(q2q1ParamList);
701  }
702  Q2Q1Fact->SetFactory("A", AFact);
703  M.SetFactory("Ptent", Q2Q1Fact);
704 
705  RCP<PatternFactory> patternFact = rcp(new PatternFactory);
706  ParameterList patternParams = *(patternFact->GetValidParameterList());
707  // Our prolongator constructs the exact pattern we are going to use,
708  // therefore we do not expand it
709  patternParams.set("emin: pattern order", 0);
710  patternFact->SetParameterList(patternParams);
711  patternFact->SetFactory("A", AFact);
712  patternFact->SetFactory("P", Q2Q1Fact);
713  M.SetFactory("Ppattern", patternFact);
714 
715  RCP<ConstraintFactory> CFact = rcp(new ConstraintFactory);
716  CFact->SetFactory("Ppattern", patternFact);
717  M.SetFactory("Constraint", CFact);
718 
719  RCP<EminPFactory> EminPFact = rcp(new EminPFactory());
720  ParameterList eminParams = *(EminPFact->GetValidParameterList());
721  if (paramList.isParameter("emin: num iterations"))
722  eminParams.set("emin: num iterations", paramList.get<int>("emin: num iterations"));
723  if (mode == "pressure") {
724  eminParams.set("emin: iterative method", "cg");
725  } else {
726  eminParams.set("emin: iterative method", "gmres");
727  if (paramList.isParameter("emin: iterative method"))
728  eminParams.set("emin: iterative method", paramList.get<std::string>("emin: iterative method"));
729  }
730  EminPFact->SetParameterList(eminParams);
731  EminPFact->SetFactory("A", AFact);
732  EminPFact->SetFactory("Constraint", CFact);
733  EminPFact->SetFactory("P", Q2Q1Fact);
734  M.SetFactory("P", EminPFact);
735 
736  if (mode == "velocity" && (!paramList.isParameter("velocity: use transpose") || paramList.get<bool>("velocity: use transpose") == false)) {
737  // Pressure system is symmetric, so it does not matter
738  // Velocity system may benefit from running emin in restriction mode (with A^T)
739  RCP<GenericRFactory> RFact = rcp(new GenericRFactory());
740  RFact->SetFactory("P", EminPFact);
741  M.SetFactory("R", RFact);
742  }
743  }
744 
745  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
746  RCP<MueLu::FactoryBase>
748  GetSmoother(const std::string& type, const ParameterList& paramList, bool coarseSolver) const {
749  typedef Teuchos::ParameterEntry ParameterEntry;
750 
751  typedef MueLu::BlockedDirectSolver <SC,LO,GO,NO> BlockedDirectSolver;
752  typedef MueLu::BraessSarazinSmoother <SC,LO,GO,NO> BraessSarazinSmoother;
753  typedef MueLu::DirectSolver <SC,LO,GO,NO> DirectSolver;
754  typedef MueLu::FactoryManager <SC,LO,GO,NO> FactoryManager;
755  typedef MueLu::SchurComplementFactory<SC,LO,GO,NO> SchurComplementFactory;
756  typedef MueLu::SmootherFactory <SC,LO,GO,NO> SmootherFactory;
757  typedef MueLu::SmootherPrototype <SC,LO,GO,NO> SmootherPrototype;
758  typedef MueLu::TrilinosSmoother <SC,LO,GO,NO> TrilinosSmoother;
759 
760  RCP<SmootherPrototype> smootherPrototype;
761  if (type == "none") {
762  return Teuchos::null;
763 
764  } else if (type == "vanka") {
765  // Set up Vanka smoothing via a combination of Schwarz and block relaxation.
766  ParameterList schwarzList;
767  schwarzList.set("schwarz: overlap level", as<int>(0));
768  schwarzList.set("schwarz: zero starting solution", false);
769  schwarzList.set("subdomain solver name", "Block_Relaxation");
770 
771  ParameterList& innerSolverList = schwarzList.sublist("subdomain solver parameters");
772  innerSolverList.set("partitioner: type", "user");
773  innerSolverList.set("partitioner: overlap", MUELU_GPD("partitioner: overlap", int, 1));
774  innerSolverList.set("relaxation: type", MUELU_GPD("relaxation: type", std::string, "Gauss-Seidel"));
775  innerSolverList.set("relaxation: sweeps", MUELU_GPD("relaxation: sweeps", int, 1));
776  innerSolverList.set("relaxation: damping factor", MUELU_GPD("relaxation: damping factor", double, 0.5));
777  innerSolverList.set("relaxation: zero starting solution", false);
778  // innerSolverList.set("relaxation: backward mode", MUELU_GPD("relaxation: backward mode", bool, true); NOT SUPPORTED YET
779 
780  std::string ifpackType = "SCHWARZ";
781 
782  smootherPrototype = rcp(new TrilinosSmoother(ifpackType, schwarzList));
783 
784  } else if (type == "schwarz") {
785 
786  std::string ifpackType = "SCHWARZ";
787 
788  smootherPrototype = rcp(new TrilinosSmoother(ifpackType, paramList));
789 
790  } else if (type == "braess-sarazin") {
791  // Define smoother/solver for BraessSarazin
792  SC omega = MUELU_GPD("bs: omega", double, 1.0);
793  bool lumping = MUELU_GPD("bs: lumping", bool, false);
794 
795  RCP<SchurComplementFactory> schurFact = rcp(new SchurComplementFactory());
796  schurFact->SetParameter("omega", ParameterEntry(omega));
797  schurFact->SetParameter("lumping", ParameterEntry(lumping));
798  schurFact->SetFactory ("A", MueLu::NoFactory::getRCP());
799 
800  // Schur complement solver
801  RCP<SmootherPrototype> schurSmootherPrototype;
802  std::string schurSmootherType = (paramList.isParameter("schur smoother: type") ? paramList.get<std::string>("schur smoother: type") : "RELAXATION");
803  if (schurSmootherType == "RELAXATION") {
804  ParameterList schurSmootherParams = paramList.sublist("schur smoother: params");
805  // schurSmootherParams.set("relaxation: damping factor", omega);
806  schurSmootherPrototype = rcp(new TrilinosSmoother(schurSmootherType, schurSmootherParams));
807  } else {
808  schurSmootherPrototype = rcp(new DirectSolver());
809  }
810  schurSmootherPrototype->SetFactory("A", schurFact);
811 
812  RCP<SmootherFactory> schurSmootherFact = rcp(new SmootherFactory(schurSmootherPrototype));
813 
814  // Define temporary FactoryManager that is used as input for BraessSarazin smoother
815  RCP<FactoryManager> braessManager = rcp(new FactoryManager());
816  braessManager->SetFactory("A", schurFact); // SchurComplement operator for correction step (defined as "A")
817  braessManager->SetFactory("Smoother", schurSmootherFact); // solver/smoother for correction step
818  braessManager->SetFactory("PreSmoother", schurSmootherFact);
819  braessManager->SetFactory("PostSmoother", schurSmootherFact);
820  braessManager->SetIgnoreUserData(true); // always use data from factories defined in factory manager
821 
822  smootherPrototype = rcp(new BraessSarazinSmoother());
823  smootherPrototype->SetParameter("Sweeps", ParameterEntry(MUELU_GPD("bs: sweeps", int, 1)));
824  smootherPrototype->SetParameter("lumping", ParameterEntry(lumping));
825  smootherPrototype->SetParameter("Damping factor", ParameterEntry(omega));
826  smootherPrototype->SetParameter("q2q1 mode", ParameterEntry(true));
827  rcp_dynamic_cast<BraessSarazinSmoother>(smootherPrototype)->AddFactoryManager(braessManager, 0); // set temporary factory manager in BraessSarazin smoother
828 
829  } else if (type == "ilu") {
830  std::string ifpackType = "RILUK";
831 
832  smootherPrototype = rcp(new TrilinosSmoother(ifpackType, paramList));
833 
834  } else if (type == "direct") {
835  smootherPrototype = rcp(new BlockedDirectSolver());
836 
837  } else {
838  throw MueLu::Exceptions::RuntimeError("Unknown smoother type: \"" + type + "\"");
839  }
840 
841  return coarseSolver ? rcp(new SmootherFactory(smootherPrototype, Teuchos::null)) : rcp(new SmootherFactory(smootherPrototype));
842  }
843 
844  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
845  RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
846  MueLuTpetraQ2Q1PreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Absolute(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A) const {
847  typedef Xpetra::CrsMatrix <SC,LO,GO,NO> CrsMatrix;
848  typedef Xpetra::CrsMatrixWrap<SC,LO,GO,NO> CrsMatrixWrap;
849  typedef Xpetra::Matrix <SC,LO,GO,NO> Matrix;
850 
851  const CrsMatrixWrap& Awrap = dynamic_cast<const CrsMatrixWrap&>(A);
852 
853  ArrayRCP<const size_t> iaA;
854  ArrayRCP<const LO> jaA;
855  ArrayRCP<const SC> valA;
856  Awrap.getCrsMatrix()->getAllValues(iaA, jaA, valA);
857 
858  ArrayRCP<size_t> iaB (iaA .size());
859  ArrayRCP<LO> jaB (jaA .size());
860  ArrayRCP<SC> valB(valA.size());
861  for (int i = 0; i < iaA .size(); i++) iaB [i] = iaA[i];
862  for (int i = 0; i < jaA .size(); i++) jaB [i] = jaA[i];
863  for (int i = 0; i < valA.size(); i++) valB[i] = Teuchos::ScalarTraits<SC>::magnitude(valA[i]);
864 
865  RCP<Matrix> B = rcp(new CrsMatrixWrap(A.getRowMap(), A.getColMap(), 0, Xpetra::StaticProfile));
866  RCP<CrsMatrix> Bcrs = rcp_dynamic_cast<CrsMatrixWrap>(B)->getCrsMatrix();
867  Bcrs->setAllValues(iaB, jaB, valB);
868  Bcrs->expertStaticFillComplete(A.getDomainMap(), A.getRangeMap());
869 
870  return B;
871  }
872 
873  // Public functions overridden from Teuchos::Describable
874  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
876  return "Thyra::MueLuTpetraQ2Q1PreconditionerFactory";
877  }
878 
879 } // namespace Thyra
880 
881 #endif
882 #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
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.
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.