MueLu Version of the Day
Loading...
Searching...
No Matches
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
103namespace 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()->getLocalNumElements();
315 Xpetra::global_size_t numPres = tmp_A_21->getRowMap()->getLocalNumElements();
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->getLocalNumElements();
322 Xpetra::global_size_t numCols2 = domainMap2->getLocalNumElements();
323 ArrayView<const GO> rangeElem2 = rangeMap2->getLocalElementList();
324 ArrayView<const GO> domainElem2 = domainMap2->getLocalElementList();
325 ArrayView<const GO> rowElem1 = tmp_A_12->getRowMap()->getLocalElementList();
326 ArrayView<const GO> colElem1 = tmp_A_21->getColMap()->getLocalElementList();
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->getLocalMaxNumRowEntries());
345 RCP<Matrix> A_21 = MatrixFactory::Build(newRangeMap2, tmp_A_21->getDomainMap(), tmp_A_21->getLocalMaxNumRowEntries());
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()->getLocalNumElements()); ++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()->getLocalElementList();
437 Array<GO> presElem = A_21_crs->getRangeMap()->getLocalElementList();
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()->getLocalNumElements();
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
Various adapters that will create a MueLu preconditioner that is a Tpetra::Operator.
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultScalar Scalar
#define MUELU_GPD(name, type, defaultValue)
AmalgamationFactory for subblocks of strided map based amalgamation data.
Factory for creating a graph based on a given matrix.
Exception throws to report errors in the internal logical of the program.
This class specifies the default factory that should generate some data on a Level if the data does n...
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.
Factory for building filtered matrices using filtered graphs.
MueLu representation of a graph.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
void SetLevelID(int levelID)
Set level number.
Definition: MueLu_Level.cpp:78
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
void Request(const FactoryBase &factory)
Increment the storage counter for all the inputs of a factory.
static const RCP< const NoFactory > getRCP()
Static Get() functions.
Factory for building the Schur Complement for a 2x2 block matrix.
Wraps an existing MueLu::Hierarchy as a Tpetra::Operator.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
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 SetDependencyTree(MueLu::FactoryManager< SC, LO, GO, NO > &M, const ParameterList &paramList) const
Teuchos::RCP< Xpetra::Matrix< SC, LO, GO, NO > > Absolute(const Xpetra::Matrix< SC, LO, GO, NO > &A) const
void uninitializePrec(PreconditionerBase< SC > *prec, Teuchos::RCP< const LinearOpSourceBase< SC > > *fwdOp, ESupportSolveUse *supportSolveUse) const
void SetBlockDependencyTree(MueLu::FactoryManager< SC, LO, GO, NO > &M, LO row, LO col, const std::string &mode, const ParameterList &paramList) 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
RCP< MueLu::FactoryBase > GetSmoother(const std::string &type, const ParameterList &paramList, bool coarseSolver) const
void initializePrec(const Teuchos::RCP< const LinearOpSourceBase< SC > > &fwdOp, PreconditionerBase< SC > *prec, const ESupportSolveUse supportSolveUse) const
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &paramList)
Teuchos::RCP< const Teuchos::ParameterList > getParameterList() const
bool isCompatible(const LinearOpSourceBase< SC > &fwdOp) const
RCP< Xpetra::Matrix< SC, LO, GO, NO > > TpetraCrs_To_XpetraMatrix(const Teuchos::RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > &Atpetra)
std::string toString(const T &what)
Little helper function to convert non-string types to strings.