MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_Ifpack2Smoother_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 MUELU_IFPACK2SMOOTHER_DEF_HPP
47#define MUELU_IFPACK2SMOOTHER_DEF_HPP
48
49#include "MueLu_ConfigDefs.hpp"
50
51#if defined(HAVE_MUELU_TPETRA) && defined(HAVE_MUELU_IFPACK2)
52
53#include <Teuchos_ParameterList.hpp>
54
55#include <Tpetra_RowMatrix.hpp>
56
57#include <Ifpack2_Chebyshev.hpp>
58#include <Ifpack2_RILUK.hpp>
59#include <Ifpack2_Relaxation.hpp>
60#include <Ifpack2_ILUT.hpp>
61#include <Ifpack2_BlockRelaxation.hpp>
62#include <Ifpack2_Factory.hpp>
63#include <Ifpack2_Parameters.hpp>
64
65#include <Xpetra_BlockedCrsMatrix.hpp>
66#include <Xpetra_CrsMatrix.hpp>
67#include <Xpetra_CrsMatrixWrap.hpp>
68#include <Xpetra_TpetraBlockCrsMatrix.hpp>
69#include <Xpetra_Matrix.hpp>
70#include <Xpetra_MatrixMatrix.hpp>
71#include <Xpetra_MultiVectorFactory.hpp>
72#include <Xpetra_TpetraMultiVector.hpp>
73
74#include <Tpetra_BlockCrsMatrix_Helpers.hpp>
75
77#include "MueLu_Level.hpp"
79#include "MueLu_Utilities.hpp"
80#include "MueLu_Monitor.hpp"
81#include "MueLu_Aggregates.hpp"
82
83
84#ifdef HAVE_MUELU_INTREPID2
87#include "Intrepid2_Basis.hpp"
88#include "Kokkos_DynRankView.hpp"
89#endif
90
91
92namespace MueLu {
93
94 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
95 Ifpack2Smoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Ifpack2Smoother(const std::string& type, const Teuchos::ParameterList& paramList, const LO& overlap)
96 : type_(type), overlap_(overlap)
97 {
98 typedef Tpetra::RowMatrix<SC,LO,GO,NO> tRowMatrix;
99 bool isSupported = Ifpack2::Factory::isSupported<tRowMatrix>(type_) || (type_ == "LINESMOOTHING_TRIDI_RELAXATION" ||
100 type_ == "LINESMOOTHING_TRIDI RELAXATION" ||
101 type_ == "LINESMOOTHING_TRIDIRELAXATION" ||
102 type_ == "LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
103 type_ == "LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
104 type_ == "LINESMOOTHING_TRIDIAGONALRELAXATION" ||
105 type_ == "LINESMOOTHING_BANDED_RELAXATION" ||
106 type_ == "LINESMOOTHING_BANDED RELAXATION" ||
107 type_ == "LINESMOOTHING_BANDEDRELAXATION" ||
108 type_ == "LINESMOOTHING_BLOCK_RELAXATION" ||
109 type_ == "LINESMOOTHING_BLOCK RELAXATION" ||
110 type_ == "LINESMOOTHING_BLOCKRELAXATION" ||
111 type_ == "TOPOLOGICAL" ||
112 type_ == "AGGREGATE");
113 this->declareConstructionOutcome(!isSupported, "Ifpack2 does not provide the smoother '" + type_ + "'.");
114 if (isSupported)
115 SetParameterList(paramList);
116 }
117
118 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
120 Factory::SetParameterList(paramList);
121
123 // It might be invalid to change parameters after the setup, but it depends entirely on Ifpack implementation.
124 // TODO: I don't know if Ifpack returns an error code or exception or ignore parameters modification in this case...
125 SetPrecParameters();
126 }
127 }
128
129 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
131 std::string prefix = this->ShortClassName() + ": SetPrecParameters";
132 RCP<TimeMonitor> tM = rcp(new TimeMonitor(*this, prefix, Timings0));
133 ParameterList& paramList = const_cast<ParameterList&>(this->GetParameterList());
134
135 paramList.setParameters(list);
136
137 RCP<ParameterList> precList = this->RemoveFactoriesFromList(this->GetParameterList());
138
139 if(!precList.is_null() && precList->isParameter("partitioner: type") && precList->get<std::string>("partitioner: type") == "linear" &&
140 !precList->isParameter("partitioner: local parts")) {
141 LO matrixBlockSize = 1;
142 int lclSize = A_->getRangeMap()->getLocalNumElements();
143 RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_);
144 if (!matA.is_null()) {
145 lclSize = matA->getLocalNumRows();
146 matrixBlockSize = matA->GetFixedBlockSize();
147 }
148 precList->set("partitioner: local parts", lclSize / matrixBlockSize);
149 }
150
151 prec_->setParameters(*precList);
152
153 paramList.setParameters(*precList); // what about that??
154 }
155
156 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
158 this->Input(currentLevel, "A");
159
160 if (type_ == "LINESMOOTHING_TRIDI_RELAXATION" ||
161 type_ == "LINESMOOTHING_TRIDI RELAXATION" ||
162 type_ == "LINESMOOTHING_TRIDIRELAXATION" ||
163 type_ == "LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
164 type_ == "LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
165 type_ == "LINESMOOTHING_TRIDIAGONALRELAXATION" ||
166 type_ == "LINESMOOTHING_BANDED_RELAXATION" ||
167 type_ == "LINESMOOTHING_BANDED RELAXATION" ||
168 type_ == "LINESMOOTHING_BANDEDRELAXATION" ||
169 type_ == "LINESMOOTHING_BLOCK_RELAXATION" ||
170 type_ == "LINESMOOTHING_BLOCK RELAXATION" ||
171 type_ == "LINESMOOTHING_BLOCKRELAXATION") {
172 this->Input(currentLevel, "CoarseNumZLayers"); // necessary for fallback criterion
173 this->Input(currentLevel, "LineDetection_VertLineIds"); // necessary to feed block smoother
174 }
175 else if (type_ == "BLOCK RELAXATION" ||
176 type_ == "BLOCK_RELAXATION" ||
177 type_ == "BLOCKRELAXATION" ||
178 // Banded
179 type_ == "BANDED_RELAXATION" ||
180 type_ == "BANDED RELAXATION" ||
181 type_ == "BANDEDRELAXATION" ||
182 // Tridiagonal
183 type_ == "TRIDI_RELAXATION" ||
184 type_ == "TRIDI RELAXATION" ||
185 type_ == "TRIDIRELAXATION" ||
186 type_ == "TRIDIAGONAL_RELAXATION" ||
187 type_ == "TRIDIAGONAL RELAXATION" ||
188 type_ == "TRIDIAGONALRELAXATION")
189 {
190 //We need to check for the "partitioner type" = "line"
191 ParameterList precList = this->GetParameterList();
192 if(precList.isParameter("partitioner: type") &&
193 precList.get<std::string>("partitioner: type") == "line") {
194 this->Input(currentLevel, "Coordinates");
195 }
196 }
197 else if (type_ == "TOPOLOGICAL")
198 {
199 // for the topological smoother, we require an element to node map:
200 this->Input(currentLevel, "pcoarsen: element to node map");
201 }
202 else if (type_ == "AGGREGATE")
203 {
204 // Aggregate smoothing needs aggregates
205 this->Input(currentLevel,"Aggregates");
206 }
207 else if (type_ == "HIPTMAIR") {
208 // Hiptmair needs D0 and NodeMatrix
209 this->Input(currentLevel,"NodeMatrix");
210 this->Input(currentLevel,"D0");
211 }
212
213 }
214
215 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
217 FactoryMonitor m(*this, "Setup Smoother", currentLevel);
218 A_ = Factory::Get< RCP<Operator> >(currentLevel, "A");
219 ParameterList& paramList = const_cast<ParameterList&>(this->GetParameterList());
220
221 // If the user asked us to convert the matrix into BlockCrsMatrix form, we do that now
222
223 if(paramList.isParameter("smoother: use blockcrsmatrix storage") && paramList.get<bool>("smoother: use blockcrsmatrix storage")) {
224 int blocksize = 1;
225 RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_);
226 if (!matA.is_null())
227 blocksize = matA->GetFixedBlockSize();
228 if (blocksize) {
229 // NOTE: Don't think you can move this out of the if block. You can't. The test MueLu_MeshTyingBlocked_SimpleSmoother_2dof_medium_MPI_1 will fail
230
231 RCP<CrsMatrixWrap> AcrsWrap = rcp_dynamic_cast<CrsMatrixWrap>(A_);
232 if(AcrsWrap.is_null())
233 throw std::runtime_error("Ifpack2Smoother: Cannot convert matrix A to CrsMatrixWrap object.");
234 RCP<CrsMatrix> Acrs = AcrsWrap->getCrsMatrix();
235 if(Acrs.is_null())
236 throw std::runtime_error("Ifpack2Smoother: Cannot extract CrsMatrix from matrix A.");
237 RCP<TpetraCrsMatrix> At = rcp_dynamic_cast<TpetraCrsMatrix>(Acrs);
238 if(At.is_null()) {
239 if(!Xpetra::Helpers<Scalar,LO,GO,Node>::isTpetraBlockCrs(matA))
240 throw std::runtime_error("Ifpack2Smoother: Cannot extract CrsMatrix or BlockCrsMatrix from matrix A.");
241 this->GetOStream(Statistics0) << "Ifpack2Smoother: Using (native) BlockCrsMatrix storage with blocksize "<<blocksize<<std::endl;
242 paramList.remove("smoother: use blockcrsmatrix storage");
243 }
244 else {
245 RCP<Tpetra::BlockCrsMatrix<Scalar, LO, GO, Node> > blockCrs = Tpetra::convertToBlockCrsMatrix(*At->getTpetra_CrsMatrix(),blocksize);
246 RCP<CrsMatrix> blockCrs_as_crs = rcp(new TpetraBlockCrsMatrix(blockCrs));
247 RCP<CrsMatrixWrap> blockWrap = rcp(new CrsMatrixWrap(blockCrs_as_crs));
248 A_ = blockWrap;
249 this->GetOStream(Statistics0) << "Ifpack2Smoother: Using BlockCrsMatrix storage with blocksize "<<blocksize<<std::endl;
250
251 paramList.remove("smoother: use blockcrsmatrix storage");
252 }
253 }
254 }
255
256 if (type_ == "SCHWARZ")
257 SetupSchwarz(currentLevel);
258
259 else if (type_ == "LINESMOOTHING_TRIDI_RELAXATION" ||
260 type_ == "LINESMOOTHING_TRIDI RELAXATION" ||
261 type_ == "LINESMOOTHING_TRIDIRELAXATION" ||
262 type_ == "LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
263 type_ == "LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
264 type_ == "LINESMOOTHING_TRIDIAGONALRELAXATION" ||
265 type_ == "LINESMOOTHING_BANDED_RELAXATION" ||
266 type_ == "LINESMOOTHING_BANDED RELAXATION" ||
267 type_ == "LINESMOOTHING_BANDEDRELAXATION" ||
268 type_ == "LINESMOOTHING_BLOCK_RELAXATION" ||
269 type_ == "LINESMOOTHING_BLOCK RELAXATION" ||
270 type_ == "LINESMOOTHING_BLOCKRELAXATION")
271 SetupLineSmoothing(currentLevel);
272
273 else if (type_ == "BLOCK_RELAXATION" ||
274 type_ == "BLOCK RELAXATION" ||
275 type_ == "BLOCKRELAXATION" ||
276 // Banded
277 type_ == "BANDED_RELAXATION" ||
278 type_ == "BANDED RELAXATION" ||
279 type_ == "BANDEDRELAXATION" ||
280 // Tridiagonal
281 type_ == "TRIDI_RELAXATION" ||
282 type_ == "TRIDI RELAXATION" ||
283 type_ == "TRIDIRELAXATION" ||
284 type_ == "TRIDIAGONAL_RELAXATION" ||
285 type_ == "TRIDIAGONAL RELAXATION" ||
286 type_ == "TRIDIAGONALRELAXATION")
287 SetupBlockRelaxation(currentLevel);
288
289 else if (type_ == "CHEBYSHEV")
290 SetupChebyshev(currentLevel);
291
292 else if (type_ == "TOPOLOGICAL")
293 {
294#ifdef HAVE_MUELU_INTREPID2
295 SetupTopological(currentLevel);
296#else
297 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "'TOPOLOGICAL' smoother choice requires Intrepid2");
298#endif
299 }
300 else if (type_ == "AGGREGATE")
301 SetupAggregate(currentLevel);
302
303 else if (type_ == "HIPTMAIR")
304 SetupHiptmair(currentLevel);
305
306 else
307 {
308 SetupGeneric(currentLevel);
309 }
310
312
313 this->GetOStream(Statistics1) << description() << std::endl;
314 }
315
316 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
318 typedef Tpetra::RowMatrix<SC,LO,GO,NO> tRowMatrix;
319
320 bool reusePreconditioner = false;
321 if (this->IsSetup() == true) {
322 // Reuse the constructed preconditioner
323 this->GetOStream(Runtime1) << "MueLu::Ifpack2Smoother::SetupSchwarz(): Setup() has already been called, assuming reuse" << std::endl;
324
325 bool isTRowMatrix = true;
326 RCP<const tRowMatrix> tA;
327 try {
329 } catch (Exceptions::BadCast&) {
330 isTRowMatrix = false;
331 }
332
333 RCP<Ifpack2::Details::CanChangeMatrix<tRowMatrix> > prec = rcp_dynamic_cast<Ifpack2::Details::CanChangeMatrix<tRowMatrix> >(prec_);
334 if (!prec.is_null() && isTRowMatrix) {
335 prec->setMatrix(tA);
336 reusePreconditioner = true;
337 } else {
338 this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupSchwarz(): reuse of this type is not available "
339 "(either failed cast to CanChangeMatrix, or to Tpetra Row Matrix), reverting to full construction" << std::endl;
340 }
341 }
342
343 if (!reusePreconditioner) {
344 ParameterList& paramList = const_cast<ParameterList&>(this->GetParameterList());
345
346 bool isBlockedMatrix = false;
347 RCP<Matrix> merged2Mat;
348
349 std::string sublistName = "subdomain solver parameters";
350 if (paramList.isSublist(sublistName)) {
351 // If we are doing "user" partitioning, we assume that what the user
352 // really wants to do is make tiny little subdomains with one row
353 // assigned to each subdomain. The rows used for these little
354 // subdomains correspond to those in the 2nd block row. Then,
355 // if we overlap these mini-subdomains, we will do something that
356 // looks like Vanka (grabbing all velocities associated with each
357 // each pressure unknown). In addition, we put all Dirichlet points
358 // as a little mini-domain.
359 ParameterList& subList = paramList.sublist(sublistName);
360
361 std::string partName = "partitioner: type";
362 // Pretty sure no one has been using this. Unfortunately, old if
363 // statement (which checked for equality with "user") prevented
364 // anyone from employing other types of Ifpack2 user partition
365 // options. Leaving this and switching if to "vanka user" just
366 // in case some day someone might want to use this.
367 if (subList.isParameter(partName) && subList.get<std::string>(partName) == "vanka user") {
368 isBlockedMatrix = true;
369
370 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
371 TEUCHOS_TEST_FOR_EXCEPTION(bA.is_null(), Exceptions::BadCast,
372 "Matrix A must be of type BlockedCrsMatrix.");
373
374 size_t numVels = bA->getMatrix(0,0)->getLocalNumRows();
375 size_t numPres = bA->getMatrix(1,0)->getLocalNumRows();
376 size_t numRows = rcp_dynamic_cast<Matrix>(A_, true)->getLocalNumRows();
377
378 ArrayRCP<LocalOrdinal> blockSeeds(numRows, Teuchos::OrdinalTraits<LocalOrdinal>::invalid());
379
380 size_t numBlocks = 0;
381 for (size_t rowOfB = numVels; rowOfB < numVels+numPres; ++rowOfB)
382 blockSeeds[rowOfB] = numBlocks++;
383
384 RCP<BlockedCrsMatrix> bA2 = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
385 TEUCHOS_TEST_FOR_EXCEPTION(bA2.is_null(), Exceptions::BadCast,
386 "Matrix A must be of type BlockedCrsMatrix.");
387
388 merged2Mat = bA2->Merge();
389
390 // Add Dirichlet rows to the list of seeds
391 ArrayRCP<const bool> boundaryNodes = Utilities::DetectDirichletRows(*merged2Mat, 0.0);
392 bool haveBoundary = false;
393 for (LO i = 0; i < boundaryNodes.size(); i++)
394 if (boundaryNodes[i]) {
395 // FIXME:
396 // 1. would not this [] overlap with some in the previos blockSeed loop?
397 // 2. do we need to distinguish between pressure and velocity Dirichlet b.c.
398 blockSeeds[i] = numBlocks;
399 haveBoundary = true;
400 }
401 if (haveBoundary)
402 numBlocks++;
403
404 subList.set("partitioner: type","user");
405 subList.set("partitioner: map", blockSeeds);
406 subList.set("partitioner: local parts", as<int>(numBlocks));
407
408 } else {
409 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
410 if (!bA.is_null()) {
411 isBlockedMatrix = true;
412 merged2Mat = bA->Merge();
413 }
414 }
415 }
416
417 RCP<const tRowMatrix> tA;
418 if (isBlockedMatrix == true) tA = Utilities::Op2NonConstTpetraRow(merged2Mat);
420
421 prec_ = Ifpack2::Factory::create(type_, tA, overlap_);
422 SetPrecParameters();
423
424 prec_->initialize();
425 }
426
427 prec_->compute();
428 }
429
430
431
432 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
434
435 ParameterList& paramList = const_cast<ParameterList&>(this->GetParameterList());
436
437 if (this->IsSetup() == true) {
438 this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupAggregate(): Setup() has already been called" << std::endl;
439 this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupAggregate(): reuse of this type is not available, reverting to full construction" << std::endl;
440 }
441
442 this->GetOStream(Statistics0) << "Ifpack2Smoother: Using Aggregate Smoothing"<<std::endl;
443
444 RCP<Aggregates> aggregates = Factory::Get<RCP<Aggregates> >(currentLevel,"Aggregates");
445
446 RCP<const LOMultiVector> vertex2AggId = aggregates->GetVertex2AggId();
447 ArrayRCP<LO> aggregate_ids = rcp_const_cast<LOMultiVector>(vertex2AggId)->getDataNonConst(0);
448 ArrayRCP<LO> dof_ids;
449
450 // We need to unamalgamate, if the FixedBlockSize > 1
451 LO blocksize = 1;
452 RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_);
453 if (!matA.is_null())
454 blocksize = matA->GetFixedBlockSize();
455 if(blocksize > 1) {
456 dof_ids.resize(aggregate_ids.size() * blocksize);
457 for(LO i=0; i<(LO)aggregate_ids.size(); i++) {
458 for(LO j=0; j<(LO)blocksize; j++)
459 dof_ids[i*blocksize+j] = aggregate_ids[i];
460 }
461 }
462 else {
463 dof_ids = aggregate_ids;
464 }
465
466
467 paramList.set("partitioner: map", dof_ids);
468 paramList.set("partitioner: type", "user");
469 paramList.set("partitioner: overlap", 0);
470 paramList.set("partitioner: local parts", (int)aggregates->GetNumAggregates());
471
472 RCP<const Tpetra::RowMatrix<SC, LO, GO, NO> > tA = Utilities::Op2NonConstTpetraRow(A_);
473
474 type_ = "BLOCKRELAXATION";
475 prec_ = Ifpack2::Factory::create(type_, tA, overlap_);
476 SetPrecParameters();
477 prec_->initialize();
478 prec_->compute();
479 }
480
481#ifdef HAVE_MUELU_INTREPID2
482 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
484 /*
485
486 basic notion:
487
488 Look for user input indicating topo dimension, something like "topological domain type: {node|edge|face}"
489 Call something like what you can find in Poisson example line 1180 to set seeds for a smoother
490
491 */
492 if (this->IsSetup() == true) {
493 this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupTopological(): Setup() has already been called" << std::endl;
494 this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupTopological(): reuse of this type is not available, reverting to full construction" << std::endl;
495 }
496
497 ParameterList& paramList = const_cast<ParameterList&>(this->GetParameterList());
498
499 typedef typename Node::device_type::execution_space ES;
500
501 typedef Kokkos::DynRankView<LocalOrdinal,typename Node::device_type> FCO; //
502
503 LocalOrdinal lo_invalid = Teuchos::OrdinalTraits<LO>::invalid();
504
505 using namespace std;
506
507 const Teuchos::RCP<FCO> elemToNode = Factory::Get<Teuchos::RCP<FCO> >(currentLevel,"pcoarsen: element to node map");
508
509 string basisString = paramList.get<string>("pcoarsen: hi basis");
510 int degree;
511 // NOTE: To make sure Stokhos works we only instantiate these guys with double. There's a lot
512 // of stuff in the guts of Intrepid2 that doesn't play well with Stokhos as of yet. Here, we only
513 // care about the assignment of basis ordinals to topological entities, so this code is actually
514 // independent of the Scalar type--hard-coding double here won't hurt us.
515 auto basis = MueLuIntrepid::BasisFactory<double,ES>(basisString, degree);
516
517 string topologyTypeString = paramList.get<string>("smoother: neighborhood type");
518 int dimension;
519 if (topologyTypeString == "node")
520 dimension = 0;
521 else if (topologyTypeString == "edge")
522 dimension = 1;
523 else if (topologyTypeString == "face")
524 dimension = 2;
525 else if (topologyTypeString == "cell")
526 dimension = basis->getBaseCellTopology().getDimension();
527 else
528 TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Unrecognized smoother neighborhood type. Supported types are node, edge, face.");
529 RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_, true);
530 vector<vector<LocalOrdinal>> seeds;
531 MueLuIntrepid::FindGeometricSeedOrdinals(basis, *elemToNode, seeds, *matA->getRowMap(), *matA->getColMap());
532
533 // Ifpack2 wants the seeds in an array of the same length as the number of local elements,
534 // with local partition #s marked for the ones that are seeds, and invalid for the rest
535 int myNodeCount = matA->getRowMap()->getLocalNumElements();
536 ArrayRCP<LocalOrdinal> nodeSeeds(myNodeCount,lo_invalid);
537 int localPartitionNumber = 0;
538 for (LocalOrdinal seed : seeds[dimension])
539 {
540 nodeSeeds[seed] = localPartitionNumber++;
541 }
542
543 paramList.remove("smoother: neighborhood type");
544 paramList.remove("pcoarsen: hi basis");
545
546 paramList.set("partitioner: map", nodeSeeds);
547 paramList.set("partitioner: type", "user");
548 paramList.set("partitioner: overlap", 1);
549 paramList.set("partitioner: local parts", int(seeds[dimension].size()));
550
551 RCP<const Tpetra::RowMatrix<SC, LO, GO, NO> > tA = Utilities::Op2NonConstTpetraRow(A_);
552
553 type_ = "BLOCKRELAXATION";
554 prec_ = Ifpack2::Factory::create(type_, tA, overlap_);
555 SetPrecParameters();
556 prec_->initialize();
557 prec_->compute();
558 }
559#endif
560
561 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
563 if (this->IsSetup() == true) {
564 this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupLineSmoothing(): Setup() has already been called" << std::endl;
565 this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupLineSmoothing(): reuse of this type is not available, reverting to full construction" << std::endl;
566 }
567
568 ParameterList& myparamList = const_cast<ParameterList&>(this->GetParameterList());
569
570 LO CoarseNumZLayers = Factory::Get<LO>(currentLevel,"CoarseNumZLayers");
571 if (CoarseNumZLayers > 0) {
572 Teuchos::ArrayRCP<LO> TVertLineIdSmoo = Factory::Get< Teuchos::ArrayRCP<LO> >(currentLevel, "LineDetection_VertLineIds");
573
574 // determine number of local parts
575 LO maxPart = 0;
576 for(size_t k = 0; k < Teuchos::as<size_t>(TVertLineIdSmoo.size()); k++) {
577 if(maxPart < TVertLineIdSmoo[k]) maxPart = TVertLineIdSmoo[k];
578 }
579 RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_, true);
580 size_t numLocalRows = matA->getLocalNumRows();
581
582 TEUCHOS_TEST_FOR_EXCEPTION(numLocalRows % TVertLineIdSmoo.size() != 0, Exceptions::RuntimeError,
583 "MueLu::Ifpack2Smoother::Setup(): the number of local nodes is incompatible with the TVertLineIdsSmoo.");
584
585 //actualDofsPerNode is the actual number of matrix rows per mesh element.
586 //It is encoded in either the MueLu Level, or in the Xpetra matrix block size.
587 //This value is needed by Ifpack2 to do decoupled block relaxation.
588 int actualDofsPerNode = numLocalRows / TVertLineIdSmoo.size();
589 LO matrixBlockSize = matA->GetFixedBlockSize();
590 if(matrixBlockSize > 1 && actualDofsPerNode > 1)
591 {
592 TEUCHOS_TEST_FOR_EXCEPTION(actualDofsPerNode != matrixBlockSize, Exceptions::RuntimeError,
593 "MueLu::Ifpack2Smoother::Setup(): A is a block matrix but its block size and DOFs/node from partitioner disagree");
594 }
595 else if(matrixBlockSize > 1)
596 {
597 actualDofsPerNode = matrixBlockSize;
598 }
599 myparamList.set("partitioner: PDE equations", actualDofsPerNode);
600
601 if (numLocalRows == Teuchos::as<size_t>(TVertLineIdSmoo.size())) {
602 myparamList.set("partitioner: type","user");
603 myparamList.set("partitioner: map",TVertLineIdSmoo);
604 myparamList.set("partitioner: local parts",maxPart+1);
605 } else {
606 // we assume a constant number of DOFs per node
607 size_t numDofsPerNode = numLocalRows / TVertLineIdSmoo.size();
608
609 // Create a new Teuchos::ArrayRCP<LO> of size numLocalRows and fill it with the corresponding information
610 Teuchos::ArrayRCP<LO> partitionerMap(numLocalRows, Teuchos::OrdinalTraits<LocalOrdinal>::invalid());
611 for (size_t blockRow = 0; blockRow < Teuchos::as<size_t>(TVertLineIdSmoo.size()); ++blockRow)
612 for (size_t dof = 0; dof < numDofsPerNode; dof++)
613 partitionerMap[blockRow * numDofsPerNode + dof] = TVertLineIdSmoo[blockRow];
614 myparamList.set("partitioner: type","user");
615 myparamList.set("partitioner: map",partitionerMap);
616 myparamList.set("partitioner: local parts",maxPart + 1);
617 }
618
619 if (type_ == "LINESMOOTHING_BANDED_RELAXATION" ||
620 type_ == "LINESMOOTHING_BANDED RELAXATION" ||
621 type_ == "LINESMOOTHING_BANDEDRELAXATION")
622 type_ = "BANDEDRELAXATION";
623 else if (type_ == "LINESMOOTHING_TRIDI_RELAXATION" ||
624 type_ == "LINESMOOTHING_TRIDI RELAXATION" ||
625 type_ == "LINESMOOTHING_TRIDIRELAXATION" ||
626 type_ == "LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
627 type_ == "LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
628 type_ == "LINESMOOTHING_TRIDIAGONALRELAXATION")
629 type_ = "TRIDIAGONALRELAXATION";
630 else
631 type_ = "BLOCKRELAXATION";
632 } else {
633 // line detection failed -> fallback to point-wise relaxation
634 this->GetOStream(Runtime0) << "Line detection failed: fall back to point-wise relaxation" << std::endl;
635 myparamList.remove("partitioner: type",false);
636 myparamList.remove("partitioner: map", false);
637 myparamList.remove("partitioner: local parts",false);
638 type_ = "RELAXATION";
639 }
640
641 RCP<const Tpetra::RowMatrix<SC, LO, GO, NO> > tA = Utilities::Op2NonConstTpetraRow(A_);
642
643 prec_ = Ifpack2::Factory::create(type_, tA, overlap_);
644 SetPrecParameters();
645 prec_->initialize();
646 prec_->compute();
647 }
648
649 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
651 typedef Tpetra::RowMatrix<SC,LO,GO,NO> tRowMatrix;
652
653 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
654 if (!bA.is_null())
655 A_ = bA->Merge();
656
657 RCP<const tRowMatrix> tA = Utilities::Op2NonConstTpetraRow(A_);
658
659 bool reusePreconditioner = false;
660 if (this->IsSetup() == true) {
661 // Reuse the constructed preconditioner
662 this->GetOStream(Runtime1) << "MueLu::Ifpack2Smoother::SetupBlockRelaxation(): Setup() has already been called, assuming reuse" << std::endl;
663
664 RCP<Ifpack2::Details::CanChangeMatrix<tRowMatrix> > prec = rcp_dynamic_cast<Ifpack2::Details::CanChangeMatrix<tRowMatrix> >(prec_);
665 if (!prec.is_null()) {
666 prec->setMatrix(tA);
667 reusePreconditioner = true;
668 } else {
669 this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupBlockRelaxation(): reuse of this type is not available (failed cast to CanChangeMatrix), "
670 "reverting to full construction" << std::endl;
671 }
672 }
673
674 if (!reusePreconditioner) {
675 ParameterList& myparamList = const_cast<ParameterList&>(this->GetParameterList());
676 myparamList.print();
677 if(myparamList.isParameter("partitioner: type") &&
678 myparamList.get<std::string>("partitioner: type") == "line") {
679 Teuchos::RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> > xCoordinates =
680 Factory::Get<Teuchos::RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> > >(currentLevel, "Coordinates");
681 Teuchos::RCP<Tpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> > coordinates = Teuchos::rcpFromRef(Xpetra::toTpetra<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>(*xCoordinates));
682
683 RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_);
684 size_t lclSize = A_->getRangeMap()->getLocalNumElements();
685 if (!matA.is_null())
686 lclSize = matA->getLocalNumRows();
687 size_t numDofsPerNode = lclSize / xCoordinates->getMap()->getLocalNumElements();
688 myparamList.set("partitioner: coordinates", coordinates);
689 myparamList.set("partitioner: PDE equations", (int) numDofsPerNode);
690 }
691
692 prec_ = Ifpack2::Factory::create(type_, tA, overlap_);
693 SetPrecParameters();
694 prec_->initialize();
695 }
696
697 prec_->compute();
698 }
699
700 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
702 typedef Tpetra::RowMatrix<SC,LO,GO,NO> tRowMatrix;
703 using STS = Teuchos::ScalarTraits<SC>;
704 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
705 if (!bA.is_null())
706 A_ = bA->Merge();
707
708 RCP<const tRowMatrix> tA = Utilities::Op2NonConstTpetraRow(A_);
709
710 bool reusePreconditioner = false;
711
712 if (this->IsSetup() == true) {
713 // Reuse the constructed preconditioner
714 this->GetOStream(Runtime1) << "MueLu::Ifpack2Smoother::SetupChebyshev(): Setup() has already been called, assuming reuse" << std::endl;
715
716 RCP<Ifpack2::Details::CanChangeMatrix<tRowMatrix> > prec = rcp_dynamic_cast<Ifpack2::Details::CanChangeMatrix<tRowMatrix> >(prec_);
717 if (!prec.is_null()) {
718 prec->setMatrix(tA);
719 reusePreconditioner = true;
720 } else {
721 this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupChebyshev(): reuse of this type is not available (failed cast to CanChangeMatrix), "
722 "reverting to full construction" << std::endl;
723 }
724 }
725
726 // Take care of the eigenvalues
727 ParameterList& paramList = const_cast<ParameterList&>(this->GetParameterList());
728 SC negone = -STS::one();
729 SC lambdaMax = SetupChebyshevEigenvalues(currentLevel,"A","",paramList);
730
731
732 if (!reusePreconditioner) {
733 prec_ = Ifpack2::Factory::create(type_, tA, overlap_);
734 SetPrecParameters();
735 {
736 SubFactoryMonitor(*this, "Preconditioner init", currentLevel);
737 prec_->initialize();
738 }
739 } else
740 SetPrecParameters();
741
742 {
743 SubFactoryMonitor(*this, "Preconditioner compute", currentLevel);
744 prec_->compute();
745 }
746
747 if (lambdaMax == negone) {
748 typedef Tpetra::RowMatrix<SC, LO, GO, NO> MatrixType;
749
750 Teuchos::RCP<Ifpack2::Chebyshev<MatrixType> > chebyPrec = rcp_dynamic_cast<Ifpack2::Chebyshev<MatrixType> >(prec_);
751 if (chebyPrec != Teuchos::null) {
752 lambdaMax = chebyPrec->getLambdaMaxForApply();
753 RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_);
754 if (!matA.is_null())
755 matA->SetMaxEigenvalueEstimate(lambdaMax);
756 this->GetOStream(Statistics1) << "chebyshev: max eigenvalue (calculated by Ifpack2)" << " = " << lambdaMax << std::endl;
757 }
758 TEUCHOS_TEST_FOR_EXCEPTION(lambdaMax == negone, Exceptions::RuntimeError, "MueLu::Ifpack2Smoother::Setup(): no maximum eigenvalue estimate");
759 }
760 }
761
762
763 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
765 typedef Tpetra::RowMatrix<SC,LO,GO,NO> tRowMatrix;
766 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
767 if (!bA.is_null())
768 A_ = bA->Merge();
769
770 RCP<const tRowMatrix> tA = Utilities::Op2NonConstTpetraRow(A_);
771
772 bool reusePreconditioner = false;
773 if (this->IsSetup() == true) {
774 // Reuse the constructed preconditioner
775 this->GetOStream(Runtime1) << "MueLu::Ifpack2Smoother::SetupHiptmair(): Setup() has already been called, assuming reuse" << std::endl;
776
777 RCP<Ifpack2::Details::CanChangeMatrix<tRowMatrix> > prec = rcp_dynamic_cast<Ifpack2::Details::CanChangeMatrix<tRowMatrix> >(prec_);
778 if (!prec.is_null()) {
779 prec->setMatrix(tA);
780 reusePreconditioner = true;
781 } else {
782 this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupHiptmair(): reuse of this type is not available (failed cast to CanChangeMatrix), "
783 "reverting to full construction" << std::endl;
784 }
785 }
786
787 // If we're doing Chebyshev subsmoothing, we'll need to get the eigenvalues
788 ParameterList& paramList = const_cast<ParameterList&>(this->GetParameterList());
789 std::string smoother1 = paramList.get("hiptmair: smoother type 1","CHEBYSHEV");
790 std::string smoother2 = paramList.get("hiptmair: smoother type 2","CHEBYSHEV");
791 // SC lambdaMax11,lambdaMax22;
792
793 if(smoother1 == "CHEBYSHEV") {
794 ParameterList & list1 = paramList.sublist("hiptmair: smoother list 1");
795 //lambdaMax11 =
796 SetupChebyshevEigenvalues(currentLevel,"A","EdgeMatrix ",list1);
797 }
798 if(smoother2 == "CHEBYSHEV") {
799 ParameterList & list2 = paramList.sublist("hiptmair: smoother list 2");
800 //lambdaMax22 =
801 SetupChebyshevEigenvalues(currentLevel,"A","EdgeMatrix ",list2);
802 }
803
804 // FIXME: Should really add some checks to make sure the eigenvalue calcs worked like in
805 // the regular SetupChebyshev
806
807 // Grab the auxillary matrices and stick them on the list
808 RCP<Operator> NodeMatrix = currentLevel.Get<RCP<Operator> >("NodeMatrix");
809 RCP<Operator> D0 = currentLevel.Get<RCP<Operator> >("D0");
810
811 RCP<tRowMatrix> tNodeMatrix = Utilities::Op2NonConstTpetraRow(NodeMatrix);
812 RCP<tRowMatrix> tD0 = Utilities::Op2NonConstTpetraRow(D0);
813
814
815 Teuchos::ParameterList newlist;
816 newlist.set("P",tD0);
817 newlist.set("PtAP",tNodeMatrix);
818 if (!reusePreconditioner) {
819 prec_ = Ifpack2::Factory::create(type_, tA, overlap_);
820 SetPrecParameters(newlist);
821 prec_->initialize();
822 }
823
824 prec_->compute();
825 }
826
827
828
829 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
830 Scalar Ifpack2Smoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::SetupChebyshevEigenvalues(Level & currentLevel, const std::string & matrixName, const std::string & label, ParameterList & paramList) const {
831 // Helper: This gets used for smoothers that want to set up Chebyhev
832 typedef Teuchos::ScalarTraits<SC> STS;
833 SC negone = -STS::one();
834 RCP<const Matrix> currentA = currentLevel.Get<RCP<Matrix> >(matrixName);
835 SC lambdaMax = negone;
836
837 std::string maxEigString = "chebyshev: max eigenvalue";
838 std::string eigRatioString = "chebyshev: ratio eigenvalue";
839
840 // Get/calculate the maximum eigenvalue
841 if (paramList.isParameter(maxEigString)) {
842 if (paramList.isType<double>(maxEigString))
843 lambdaMax = paramList.get<double>(maxEigString);
844 else
845 lambdaMax = paramList.get<SC>(maxEigString);
846 this->GetOStream(Statistics1) << label << maxEigString << " (cached with smoother parameter list) = " << lambdaMax << std::endl;
847 RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_);
848 if (!matA.is_null())
849 matA->SetMaxEigenvalueEstimate(lambdaMax);
850
851 } else {
852 lambdaMax = currentA->GetMaxEigenvalueEstimate();
853 if (lambdaMax != negone) {
854 this->GetOStream(Statistics1) << label << maxEigString << " (cached with matrix) = " << lambdaMax << std::endl;
855 paramList.set(maxEigString, lambdaMax);
856 }
857 }
858
859 // Calculate the eigenvalue ratio
860 const SC defaultEigRatio = 20;
861
862 SC ratio = defaultEigRatio;
863 if (paramList.isParameter(eigRatioString)) {
864 if (paramList.isType<double>(eigRatioString))
865 ratio = paramList.get<double>(eigRatioString);
866 else
867 ratio = paramList.get<SC>(eigRatioString);
868 }
869 if (currentLevel.GetLevelID()) {
870 // Update ratio to be
871 // ratio = max(number of fine DOFs / number of coarse DOFs, defaultValue)
872 //
873 // NOTE: We don't need to request previous level matrix as we know for sure it was constructed
874 RCP<const Matrix> fineA = currentLevel.GetPreviousLevel()->Get<RCP<Matrix> >(matrixName);
875 size_t nRowsFine = fineA->getGlobalNumRows();
876 size_t nRowsCoarse = currentA->getGlobalNumRows();
877
878 SC levelRatio = as<SC>(as<float>(nRowsFine)/nRowsCoarse);
879 if (STS::magnitude(levelRatio) > STS::magnitude(ratio))
880 ratio = levelRatio;
881 }
882
883 this->GetOStream(Statistics1) << label << eigRatioString << " (computed) = " << ratio << std::endl;
884 paramList.set(eigRatioString, ratio);
885
886 if (paramList.isParameter("chebyshev: use rowsumabs diagonal scaling")) {
887 this->GetOStream(Runtime1) << "chebyshev: using rowsumabs diagonal scaling" << std::endl;
888 bool doScale = false;
889 doScale = paramList.get<bool>("chebyshev: use rowsumabs diagonal scaling");
890 paramList.remove("chebyshev: use rowsumabs diagonal scaling");
891 double chebyReplaceTol = Teuchos::ScalarTraits<Scalar>::eps()*100;
892 std::string paramName = "chebyshev: rowsumabs diagonal replacement tolerance";
893 if (paramList.isParameter(paramName)) {
894 chebyReplaceTol = paramList.get<double>(paramName);
895 paramList.remove(paramName);
896 }
897 double chebyReplaceVal = Teuchos::ScalarTraits<double>::zero();
898 paramName = "chebyshev: rowsumabs diagonal replacement value";
899 if (paramList.isParameter(paramName)) {
900 chebyReplaceVal = paramList.get<double>(paramName);
901 paramList.remove(paramName);
902 }
903 bool chebyReplaceSingleEntryRowWithZero = false;
904 paramName = "chebyshev: rowsumabs replace single entry row with zero";
905 if (paramList.isParameter(paramName)) {
906 chebyReplaceSingleEntryRowWithZero = paramList.get<bool>(paramName);
907 paramList.remove(paramName);
908 }
909 bool useAverageAbsDiagVal = false;
910 paramName = "chebyshev: rowsumabs use automatic diagonal tolerance";
911 if (paramList.isParameter(paramName)) {
912 useAverageAbsDiagVal = paramList.get<bool>(paramName);
913 paramList.remove(paramName);
914 }
915 if (doScale) {
916 const bool doReciprocal = true;
917 RCP<Vector> lumpedDiagonal = Utilities::GetLumpedMatrixDiagonal(*currentA, doReciprocal, chebyReplaceTol, chebyReplaceVal, chebyReplaceSingleEntryRowWithZero, useAverageAbsDiagVal);
918 const Xpetra::TpetraVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmpVec = dynamic_cast<const Xpetra::TpetraVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(*lumpedDiagonal);
919 paramList.set("chebyshev: operator inv diagonal",tmpVec.getTpetra_Vector());
920 }
921 }
922
923 return lambdaMax;
924 }
925
926
927
928
929 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
931 typedef Tpetra::RowMatrix<SC,LO,GO,NO> tRowMatrix;
932 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
933 if (!bA.is_null())
934 A_ = bA->Merge();
935
936 RCP<const tRowMatrix> tA = Utilities::Op2NonConstTpetraRow(A_);
937
938 bool reusePreconditioner = false;
939 if (this->IsSetup() == true) {
940 // Reuse the constructed preconditioner
941 this->GetOStream(Runtime1) << "MueLu::Ifpack2Smoother::SetupGeneric(): Setup() has already been called, assuming reuse" << std::endl;
942
943 RCP<Ifpack2::Details::CanChangeMatrix<tRowMatrix> > prec = rcp_dynamic_cast<Ifpack2::Details::CanChangeMatrix<tRowMatrix> >(prec_);
944 if (!prec.is_null()) {
945 prec->setMatrix(tA);
946 reusePreconditioner = true;
947 } else {
948 this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::SetupGeneric(): reuse of this type is not available (failed cast to CanChangeMatrix), "
949 "reverting to full construction" << std::endl;
950 }
951 }
952
953 if (!reusePreconditioner) {
954 prec_ = Ifpack2::Factory::create(type_, tA, overlap_);
955 SetPrecParameters();
956 prec_->initialize();
957 }
958
959 prec_->compute();
960 }
961
962 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
963 void Ifpack2Smoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Apply(MultiVector& X, const MultiVector& B, bool InitialGuessIsZero) const {
964 TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError, "MueLu::Ifpack2Smoother::Apply(): Setup() has not been called");
965
966 // Forward the InitialGuessIsZero option to Ifpack2
967 // TODO: It might be nice to switch back the internal
968 // "zero starting solution" option of the ifpack2 object prec_ to his
969 // initial value at the end but there is no way right now to get
970 // the current value of the "zero starting solution" in ifpack2.
971 // It's not really an issue, as prec_ can only be used by this method.
972 Teuchos::ParameterList paramList;
973 bool supportInitialGuess = false;
974 const Teuchos::ParameterList params = this->GetParameterList();
975
976 if (prec_->supportsZeroStartingSolution()) {
977 prec_->setZeroStartingSolution(InitialGuessIsZero);
978 supportInitialGuess = true;
979 } else if (type_ == "SCHWARZ") {
980 paramList.set("schwarz: zero starting solution", InitialGuessIsZero);
981 //Because additive Schwarz has "delta" semantics, it's sufficient to
982 //toggle only the zero initial guess flag, and not pass in already
983 //set parameters. If we call SetPrecParameters, the subdomain solver
984 //will be destroyed.
985 prec_->setParameters(paramList);
986 supportInitialGuess = true;
987 }
988
989 //TODO JJH 30Apr2014 Calling SetPrecParameters(paramList) when the smoother
990 //is Ifpack2::AdditiveSchwarz::setParameterList() will destroy the subdomain
991 //(aka inner) solver. This behavior is documented but a departure from what
992 //it previously did, and what other Ifpack2 solvers currently do. So I have
993 //moved SetPrecParameters(paramList) into the if-else block above.
994
995 // Apply
996 if (InitialGuessIsZero || supportInitialGuess) {
997 Tpetra::MultiVector<SC,LO,GO,NO>& tpX = Utilities::MV2NonConstTpetraMV(X);
998 const Tpetra::MultiVector<SC,LO,GO,NO>& tpB = Utilities::MV2TpetraMV(B);
999 prec_->apply(tpB, tpX);
1000 } else {
1001 typedef Teuchos::ScalarTraits<Scalar> TST;
1002
1003 RCP<MultiVector> Residual;
1004 {
1005 std::string prefix = this->ShortClassName() + ": Apply: ";
1006 RCP<TimeMonitor> tM = rcp(new TimeMonitor(*this, prefix + "residual calculation", Timings0));
1007 Residual = Utilities::Residual(*A_, X, B);
1008 }
1009
1010 RCP<MultiVector> Correction = MultiVectorFactory::Build(A_->getDomainMap(), X.getNumVectors());
1011
1012 Tpetra::MultiVector<SC,LO,GO,NO>& tpX = Utilities::MV2NonConstTpetraMV(*Correction);
1013 const Tpetra::MultiVector<SC,LO,GO,NO>& tpB = Utilities::MV2TpetraMV(*Residual);
1014
1015 prec_->apply(tpB, tpX);
1016
1017 X.update(TST::one(), *Correction, TST::one());
1018 }
1019 }
1020
1021 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
1022 RCP<MueLu::SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Ifpack2Smoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Copy() const {
1023 RCP<Ifpack2Smoother> smoother = rcp(new Ifpack2Smoother(*this) );
1024 smoother->SetParameterList(this->GetParameterList());
1025 return smoother;
1026 }
1027
1028 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
1030 std::ostringstream out;
1032 out << prec_->description();
1033 } else {
1035 out << "{type = " << type_ << "}";
1036 }
1037 return out.str();
1038 }
1039
1040 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
1041 void Ifpack2Smoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::print(Teuchos::FancyOStream &out, const VerbLevel verbLevel) const {
1043
1044 if (verbLevel & Parameters0)
1045 out0 << "Prec. type: " << type_ << std::endl;
1046
1047 if (verbLevel & Parameters1) {
1048 out0 << "Parameter list: " << std::endl;
1049 Teuchos::OSTab tab2(out);
1050 out << this->GetParameterList();
1051 out0 << "Overlap: " << overlap_ << std::endl;
1052 }
1053
1054 if (verbLevel & External)
1055 if (prec_ != Teuchos::null) {
1056 Teuchos::OSTab tab2(out);
1057 out << *prec_ << std::endl;
1058 }
1059
1060 if (verbLevel & Debug) {
1061 out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl
1062 << "-" << std::endl
1063 << "RCP<prec_>: " << prec_ << std::endl;
1064 }
1065 }
1066
1067 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
1069 typedef Tpetra::RowMatrix<SC,LO,GO,NO> MatrixType;
1070 // NOTE: Only works for a subset of Ifpack2's smoothers
1071 RCP<Ifpack2::Relaxation<MatrixType> > pr = rcp_dynamic_cast<Ifpack2::Relaxation<MatrixType> >(prec_);
1072 if(!pr.is_null()) return pr->getNodeSmootherComplexity();
1073
1074 RCP<Ifpack2::Chebyshev<MatrixType> > pc = rcp_dynamic_cast<Ifpack2::Chebyshev<MatrixType> >(prec_);
1075 if(!pc.is_null()) return pc->getNodeSmootherComplexity();
1076
1077 RCP<Ifpack2::BlockRelaxation<MatrixType> > pb = rcp_dynamic_cast<Ifpack2::BlockRelaxation<MatrixType> >(prec_);
1078 if(!pb.is_null()) return pb->getNodeSmootherComplexity();
1079
1080 RCP<Ifpack2::ILUT<MatrixType> > pi = rcp_dynamic_cast<Ifpack2::ILUT<MatrixType> >(prec_);
1081 if(!pi.is_null()) return pi->getNodeSmootherComplexity();
1082
1083 RCP<Ifpack2::RILUK<MatrixType> > pk = rcp_dynamic_cast<Ifpack2::RILUK<MatrixType> >(prec_);
1084 if(!pk.is_null()) return pk->getNodeSmootherComplexity();
1085
1086
1087 return Teuchos::OrdinalTraits<size_t>::invalid();
1088 }
1089
1090
1091} // namespace MueLu
1092
1093#endif // HAVE_MUELU_TPETRA && HAVE_MUELU_IFPACK2
1094#endif // MUELU_IFPACK2SMOOTHER_DEF_HPP
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultScalar Scalar
static Teuchos::RCP< Preconditioner< typename MatrixType::scalar_type, typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type > > create(const std::string &precType, const Teuchos::RCP< const MatrixType > &matrix)
virtual std::string description() const
Return a simple one-line description of this object.
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that encapsulates Ifpack2 smoothers.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
void SetupSchwarz(Level &currentLevel)
void Setup(Level &currentLevel)
Set up the smoother.
Scalar SetupChebyshevEigenvalues(Level &currentLevel, const std::string &matrixName, const std::string &label, Teuchos::ParameterList &paramList) const
void SetupAggregate(Level &currentLevel)
void SetupBlockRelaxation(Level &currentLevel)
void SetupChebyshev(Level &currentLevel)
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
void SetParameterList(const Teuchos::ParameterList &paramList)
Set parameters from a parameter list and return with default values.
RCP< SmootherPrototype > Copy() const
std::string description() const
Return a simple one-line description of this object.
void SetupHiptmair(Level &currentLevel)
void SetPrecParameters(const Teuchos::ParameterList &list=Teuchos::ParameterList()) const
void SetupTopological(Level &currentLevel)
std::string type_
ifpack2-specific key phrase that denote smoother type
friend class Ifpack2Smoother
Constructor.
void DeclareInput(Level &currentLevel) const
Input.
void SetupLineSmoothing(Level &currentLevel)
void SetupGeneric(Level &currentLevel)
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the preconditioner.
Class that holds all level-specific information.
RCP< Level > & GetPreviousLevel()
Previous level.
int GetLevelID() const
Return level number.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
virtual void SetParameterList(const Teuchos::ParameterList &paramList)=0
Set parameters from a parameter list and return with default values.
void declareConstructionOutcome(bool fail, std::string msg)
bool IsSetup() const
Get the state of a smoother prototype.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Integrates Teuchos::TimeMonitor with MueLu verbosity system.
static RCP< Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraRow(RCP< Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static Teuchos::RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetLumpedMatrixDiagonal(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > const &A, const bool doReciprocal=false, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100, Scalar tolReplacement=Teuchos::ScalarTraits< Scalar >::zero(), const bool replaceSingleEntryRowWithZero=false, const bool useAverageAbsDiagVal=false)
static RCP< const Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2TpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > const vec)
Helper utility to pull out the underlying Tpetra objects from an Xpetra object.
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
static RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &RHS)
static Teuchos::ArrayRCP< const bool > DetectDirichletRows(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::magnitude(0.), const bool count_twos_as_dirichlet=false)
void FindGeometricSeedOrdinals(Teuchos::RCP< Basis > basis, const LOFieldContainer &elementToNodeMap, std::vector< std::vector< LocalOrdinal > > &seeds, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &rowMap, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &columnMap)
Namespace for MueLu classes and methods.
@ Warnings0
Important warning messages (one line)
@ Debug
Print additional debugging information.
@ Statistics1
Print more statistics.
@ External
Print external lib objects.
@ Timings0
High level timing information (use Teuchos::TimeMonitor::summarize() to print)
@ Runtime0
One-liner description of what is happening.
@ Runtime1
Description of what is happening (more verbose)
@ Parameters0
Print class parameters.
@ Statistics0
Print statistics that do not involve significant additional computation.
@ Parameters1
Print class parameters (more parameters, more verbose)