Ifpack2 Templated Preconditioning Package Version 1.0
Loading...
Searching...
No Matches
Ifpack2_AdditiveSchwarz_def.hpp
Go to the documentation of this file.
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41*/
42
54
55#ifndef IFPACK2_ADDITIVESCHWARZ_DEF_HPP
56#define IFPACK2_ADDITIVESCHWARZ_DEF_HPP
57
58#include "Trilinos_Details_LinearSolverFactory.hpp"
59// We need Ifpack2's implementation of LinearSolver, because we use it
60// to wrap the user-provided Ifpack2::Preconditioner in
61// Ifpack2::AdditiveSchwarz::setInnerPreconditioner.
62#include "Ifpack2_Details_LinearSolver.hpp"
63#include "Ifpack2_Details_getParamTryingTypes.hpp"
64
65#if defined(HAVE_IFPACK2_XPETRA) && defined(HAVE_IFPACK2_ZOLTAN2)
66#include "Zoltan2_TpetraRowGraphAdapter.hpp"
67#include "Zoltan2_OrderingProblem.hpp"
68#include "Zoltan2_OrderingSolution.hpp"
69#endif
70
72#include "Ifpack2_Parameters.hpp"
73#include "Ifpack2_LocalFilter.hpp"
74#include "Ifpack2_ReorderFilter.hpp"
75#include "Ifpack2_SingletonFilter.hpp"
76#include "Ifpack2_Details_AdditiveSchwarzFilter.hpp"
77
78#ifdef HAVE_MPI
79#include "Teuchos_DefaultMpiComm.hpp"
80#endif
81
82#include "Teuchos_StandardParameterEntryValidators.hpp"
83#include <locale> // std::toupper
84
85
86// FIXME (mfh 25 Aug 2015) Work-around for Bug 6392. This doesn't
87// need to be a weak symbol because it only refers to a function in
88// the Ifpack2 package.
89namespace Ifpack2 {
90namespace Details {
91 extern void registerLinearSolverFactory ();
92} // namespace Details
93} // namespace Ifpack2
94
95#ifdef HAVE_IFPACK2_DEBUG
96
97namespace { // (anonymous)
98
99 template<class MV>
100 bool
101 anyBad (const MV& X)
102 {
103 using STS = Teuchos::ScalarTraits<typename MV::scalar_type>;
104 using magnitude_type = typename STS::magnitudeType;
105 using STM = Teuchos::ScalarTraits<magnitude_type>;
106
107 Teuchos::Array<magnitude_type> norms (X.getNumVectors ());
108 X.norm2 (norms ());
109 bool good = true;
110 for (size_t j = 0; j < X.getNumVectors (); ++j) {
111 if (STM::isnaninf (norms[j])) {
112 good = false;
113 break;
114 }
115 }
116 return ! good;
117 }
118
119} // namespace (anonymous)
120
121#endif // HAVE_IFPACK2_DEBUG
122
123namespace Ifpack2 {
124
125template<class MatrixType, class LocalInverseType>
126bool
127AdditiveSchwarz<MatrixType, LocalInverseType>::hasInnerPrecName () const
128{
129 const char* options[4] = {
130 "inner preconditioner name",
131 "subdomain solver name",
132 "schwarz: inner preconditioner name",
133 "schwarz: subdomain solver name"
134 };
135 const int numOptions = 4;
136 bool match = false;
137 for (int k = 0; k < numOptions && ! match; ++k) {
138 if (List_.isParameter (options[k])) {
139 return true;
140 }
141 }
142 return false;
143}
144
145
146template<class MatrixType, class LocalInverseType>
147void
148AdditiveSchwarz<MatrixType, LocalInverseType>::removeInnerPrecName ()
149{
150 const char* options[4] = {
151 "inner preconditioner name",
152 "subdomain solver name",
153 "schwarz: inner preconditioner name",
154 "schwarz: subdomain solver name"
155 };
156 const int numOptions = 4;
157 for (int k = 0; k < numOptions; ++k) {
158 List_.remove (options[k], false);
159 }
160}
161
162
163template<class MatrixType, class LocalInverseType>
164std::string
165AdditiveSchwarz<MatrixType, LocalInverseType>::innerPrecName () const
166{
167 const char* options[4] = {
168 "inner preconditioner name",
169 "subdomain solver name",
170 "schwarz: inner preconditioner name",
171 "schwarz: subdomain solver name"
172 };
173 const int numOptions = 4;
174 std::string newName;
175 bool match = false;
176
177 // As soon as one parameter option matches, ignore all others.
178 for (int k = 0; k < numOptions && ! match; ++k) {
179 const Teuchos::ParameterEntry* paramEnt =
180 List_.getEntryPtr (options[k]);
181 if (paramEnt != nullptr && paramEnt->isType<std::string> ()) {
182 newName = Teuchos::getValue<std::string> (*paramEnt);
183 match = true;
184 }
185 }
186 return match ? newName : defaultInnerPrecName ();
187}
188
189
190template<class MatrixType, class LocalInverseType>
191void
192AdditiveSchwarz<MatrixType, LocalInverseType>::removeInnerPrecParams ()
193{
194 const char* options[4] = {
195 "inner preconditioner parameters",
196 "subdomain solver parameters",
197 "schwarz: inner preconditioner parameters",
198 "schwarz: subdomain solver parameters"
199 };
200 const int numOptions = 4;
201
202 // As soon as one parameter option matches, ignore all others.
203 for (int k = 0; k < numOptions; ++k) {
204 List_.remove (options[k], false);
205 }
206}
207
208
209template<class MatrixType, class LocalInverseType>
210std::pair<Teuchos::ParameterList, bool>
211AdditiveSchwarz<MatrixType, LocalInverseType>::innerPrecParams () const
212{
213 const char* options[4] = {
214 "inner preconditioner parameters",
215 "subdomain solver parameters",
216 "schwarz: inner preconditioner parameters",
217 "schwarz: subdomain solver parameters"
218 };
219 const int numOptions = 4;
220 Teuchos::ParameterList params;
221
222 // As soon as one parameter option matches, ignore all others.
223 bool match = false;
224 for (int k = 0; k < numOptions && ! match; ++k) {
225 if (List_.isSublist (options[k])) {
226 params = List_.sublist (options[k]);
227 match = true;
228 }
229 }
230 // Default is an empty list of parameters.
231 return std::make_pair (params, match);
232}
233
234template<class MatrixType, class LocalInverseType>
235std::string
236AdditiveSchwarz<MatrixType, LocalInverseType>::defaultInnerPrecName ()
237{
238 // The default inner preconditioner is "ILUT", for backwards
239 // compatibility with the original AdditiveSchwarz implementation.
240 return "ILUT";
241}
242
243template<class MatrixType, class LocalInverseType>
245AdditiveSchwarz (const Teuchos::RCP<const row_matrix_type>& A) :
246 Matrix_ (A)
247{}
248
249template<class MatrixType, class LocalInverseType>
251AdditiveSchwarz (const Teuchos::RCP<const row_matrix_type>& A,
252 const int overlapLevel) :
253 Matrix_ (A),
254 OverlapLevel_ (overlapLevel)
255{}
256
257template<class MatrixType,class LocalInverseType>
258Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type > >
260getDomainMap () const
261{
262 TEUCHOS_TEST_FOR_EXCEPTION(
263 Matrix_.is_null (), std::runtime_error, "Ifpack2::AdditiveSchwarz::"
264 "getDomainMap: The matrix to precondition is null. You must either pass "
265 "a nonnull matrix to the constructor, or call setMatrix() with a nonnull "
266 "input, before you may call this method.");
267 return Matrix_->getDomainMap ();
268}
269
270
271template<class MatrixType,class LocalInverseType>
272Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
274{
275 TEUCHOS_TEST_FOR_EXCEPTION(
276 Matrix_.is_null (), std::runtime_error, "Ifpack2::AdditiveSchwarz::"
277 "getRangeMap: The matrix to precondition is null. You must either pass "
278 "a nonnull matrix to the constructor, or call setMatrix() with a nonnull "
279 "input, before you may call this method.");
280 return Matrix_->getRangeMap ();
281}
282
283
284template<class MatrixType,class LocalInverseType>
285Teuchos::RCP<const Tpetra::RowMatrix<typename MatrixType::scalar_type, typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> > AdditiveSchwarz<MatrixType,LocalInverseType>::getMatrix() const
286{
287 return Matrix_;
288}
289
290
291template<class MatrixType,class LocalInverseType>
292void
294apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &B,
295 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &Y,
296 Teuchos::ETransp mode,
297 scalar_type alpha,
298 scalar_type beta) const
299{
300 using Teuchos::Time;
301 using Teuchos::TimeMonitor;
302 using Teuchos::RCP;
303 using Teuchos::rcp;
304 using Teuchos::rcp_dynamic_cast;
305 typedef Teuchos::ScalarTraits<scalar_type> STS;
306 const char prefix[] = "Ifpack2::AdditiveSchwarz::apply: ";
307
308 TEUCHOS_TEST_FOR_EXCEPTION
309 (! IsComputed_, std::runtime_error,
310 prefix << "isComputed() must be true before you may call apply().");
311 TEUCHOS_TEST_FOR_EXCEPTION
312 (Matrix_.is_null (), std::logic_error, prefix <<
313 "The input matrix A is null, but the preconditioner says that it has "
314 "been computed (isComputed() is true). This should never happen, since "
315 "setMatrix() should always mark the preconditioner as not computed if "
316 "its argument is null. "
317 "Please report this bug to the Ifpack2 developers.");
318 TEUCHOS_TEST_FOR_EXCEPTION
319 (Inverse_.is_null (), std::runtime_error,
320 prefix << "The subdomain solver is null. "
321 "This can only happen if you called setInnerPreconditioner() with a null "
322 "input, after calling initialize() or compute(). If you choose to call "
323 "setInnerPreconditioner() with a null input, you must then call it with "
324 "a nonnull input before you may call initialize() or compute().");
325 TEUCHOS_TEST_FOR_EXCEPTION
326 (B.getNumVectors() != Y.getNumVectors(), std::invalid_argument,
327 prefix << "B and Y must have the same number of columns. B has " <<
328 B.getNumVectors () << " columns, but Y has " << Y.getNumVectors() << ".");
329 TEUCHOS_TEST_FOR_EXCEPTION
330 (IsOverlapping_ && OverlappingMatrix_.is_null (), std::logic_error,
331 prefix << "The overlapping matrix is null. "
332 "This should never happen if IsOverlapping_ is true. "
333 "Please report this bug to the Ifpack2 developers.");
334 TEUCHOS_TEST_FOR_EXCEPTION
335 (! IsOverlapping_ && localMap_.is_null (), std::logic_error,
336 prefix << "localMap_ is null. "
337 "This should never happen if IsOverlapping_ is false. "
338 "Please report this bug to the Ifpack2 developers.");
339 TEUCHOS_TEST_FOR_EXCEPTION
340 (alpha != STS::one (), std::logic_error,
341 prefix << "Not implemented for alpha != 1.");
342 TEUCHOS_TEST_FOR_EXCEPTION
343 (beta != STS::zero (), std::logic_error,
344 prefix << "Not implemented for beta != 0.");
345
346#ifdef HAVE_IFPACK2_DEBUG
347 {
348 const bool bad = anyBad (B);
349 TEUCHOS_TEST_FOR_EXCEPTION
350 (bad, std::runtime_error, "Ifpack2::AdditiveSchwarz::apply: "
351 "The 2-norm of the input B is NaN or Inf.");
352 }
353#endif // HAVE_IFPACK2_DEBUG
354
355#ifdef HAVE_IFPACK2_DEBUG
356 if (! ZeroStartingSolution_) {
357 const bool bad = anyBad (Y);
358 TEUCHOS_TEST_FOR_EXCEPTION
359 (bad, std::runtime_error, "Ifpack2::AdditiveSchwarz::apply: "
360 "On input, the initial guess Y has 2-norm NaN or Inf "
361 "(ZeroStartingSolution_ is false).");
362 }
363#endif // HAVE_IFPACK2_DEBUG
364
365 const std::string timerName ("Ifpack2::AdditiveSchwarz::apply");
366 RCP<Time> timer = TimeMonitor::lookupCounter (timerName);
367 if (timer.is_null ()) {
368 timer = TimeMonitor::getNewCounter (timerName);
369 }
370 double startTime = timer->wallTime();
371
372 { // Start timing here.
373 TimeMonitor timeMon (*timer);
374
375 const scalar_type ZERO = Teuchos::ScalarTraits<scalar_type>::zero ();
376 const size_t numVectors = B.getNumVectors ();
377
378 // mfh 25 Apr 2015: Fix for currently failing
379 // Ifpack2_AdditiveSchwarz_RILUK test.
380 if (ZeroStartingSolution_) {
381 Y.putScalar (ZERO);
382 }
383
384 // set up for overlap communication
385 MV* OverlappingB = nullptr;
386 MV* OverlappingY = nullptr;
387 {
388 RCP<const map_type> B_and_Y_map = IsOverlapping_ ?
389 OverlappingMatrix_->getRowMap () : localMap_;
390 if (overlapping_B_.get () == nullptr ||
391 overlapping_B_->getNumVectors () != numVectors) {
392 overlapping_B_.reset (new MV (B_and_Y_map, numVectors, false));
393 }
394 if (overlapping_Y_.get () == nullptr ||
395 overlapping_Y_->getNumVectors () != numVectors) {
396 overlapping_Y_.reset (new MV (B_and_Y_map, numVectors, false));
397 }
398 OverlappingB = overlapping_B_.get ();
399 OverlappingY = overlapping_Y_.get ();
400 // FIXME (mfh 25 Jun 2019) It's not clear whether we really need
401 // to fill with zeros here, but that's what was happening before.
402 OverlappingB->putScalar (ZERO);
403 OverlappingY->putScalar (ZERO);
404 }
405
406 RCP<MV> globalOverlappingB;
407 if (! IsOverlapping_) {
408 globalOverlappingB =
409 OverlappingB->offsetViewNonConst (Matrix_->getRowMap (), 0);
410
411 // Create Import object on demand, if necessary.
412 if (DistributedImporter_.is_null ()) {
413 // FIXME (mfh 15 Apr 2014) Why can't we just ask the Matrix
414 // for its Import object? Of course a general RowMatrix might
415 // not necessarily have one.
416 DistributedImporter_ =
417 rcp (new import_type (Matrix_->getRowMap (),
418 Matrix_->getDomainMap ()));
419 }
420 }
421
422 if (R_.get () == nullptr || R_->getNumVectors () != numVectors) {
423 R_.reset (new MV (B.getMap (), numVectors, false));
424 }
425 if (C_.get () == nullptr || C_->getNumVectors () != numVectors) {
426 C_.reset (new MV (Y.getMap (), numVectors, false));
427 }
428 // If taking averages in overlap region, we need to compute
429 // the number of procs who have a copy of each overlap dof
430 Teuchos::ArrayRCP<scalar_type> dataNumOverlapCopies;
431 if (IsOverlapping_ && AvgOverlap_) {
432 if (num_overlap_copies_.get() == nullptr) {
433 num_overlap_copies_.reset (new MV (Y.getMap (), 1, false));
434 RCP<MV> onesVec( new MV(OverlappingMatrix_->getRowMap(), 1, false) );
435 onesVec->putScalar(Teuchos::ScalarTraits<scalar_type>::one());
436 rcp_dynamic_cast<OverlappingRowMatrix<row_matrix_type>> (OverlappingMatrix_)->exportMultiVector (*onesVec, *(num_overlap_copies_.get ()), CombineMode_);
437 }
438 dataNumOverlapCopies = num_overlap_copies_.get ()->getDataNonConst(0);
439 }
440
441 MV* R = R_.get ();
442 MV* C = C_.get ();
443
444 // FIXME (mfh 25 Jun 2019) It was never clear whether C had to be
445 // initialized to zero. R definitely should not need this.
446 C->putScalar (ZERO);
447
448 for (int ni=0; ni<NumIterations_; ++ni)
449 {
450#ifdef HAVE_IFPACK2_DEBUG
451 {
452 const bool bad = anyBad (Y);
453 TEUCHOS_TEST_FOR_EXCEPTION
454 (bad, std::runtime_error, "Ifpack2::AdditiveSchwarz::apply: "
455 "At top of iteration " << ni << ", the 2-norm of Y is NaN or Inf.");
456 }
457#endif // HAVE_IFPACK2_DEBUG
458
459 Tpetra::deep_copy(*R, B);
460
461 // if (ZeroStartingSolution_ && ni == 0) {
462 // Y.putScalar (STS::zero ());
463 // }
464 if (!ZeroStartingSolution_ || ni > 0) {
465 //calculate residual
466 Matrix_->apply (Y, *R, mode, -STS::one(), STS::one());
467
468#ifdef HAVE_IFPACK2_DEBUG
469 {
470 const bool bad = anyBad (*R);
471 TEUCHOS_TEST_FOR_EXCEPTION
472 (bad, std::runtime_error, "Ifpack2::AdditiveSchwarz::apply: "
473 "At iteration " << ni << ", the 2-norm of R (result of computing "
474 "residual with Y) is NaN or Inf.");
475 }
476#endif // HAVE_IFPACK2_DEBUG
477 }
478
479 // do communication if necessary
480 if (IsOverlapping_) {
481 TEUCHOS_TEST_FOR_EXCEPTION
482 (OverlappingMatrix_.is_null (), std::logic_error, prefix <<
483 "IsOverlapping_ is true, but OverlappingMatrix_, while nonnull, is "
484 "not an OverlappingRowMatrix<row_matrix_type>. Please report this "
485 "bug to the Ifpack2 developers.");
486 OverlappingMatrix_->importMultiVector (*R, *OverlappingB, Tpetra::INSERT);
487
488 //JJH We don't need to import the solution Y we are always solving AY=R with initial guess zero
489 //if (ZeroStartingSolution_ == false)
490 // overlapMatrix->importMultiVector (Y, *OverlappingY, Tpetra::INSERT);
491 /*
492 FIXME from Ifpack1: Will not work with non-zero starting solutions.
493 TODO JJH 3/20/15 I don't know whether this comment is still valid.
494
495 Here is the log for the associated commit 720b2fa4 to Ifpack1:
496
497 "Added a note to recall that the nonzero starting solution will not
498 work properly if reordering, filtering or wider overlaps are used. This only
499 applied to methods like Jacobi, Gauss-Seidel, and SGS (in both point and block
500 version), and not to ILU-type preconditioners."
501 */
502
503#ifdef HAVE_IFPACK2_DEBUG
504 {
505 const bool bad = anyBad (*OverlappingB);
506 TEUCHOS_TEST_FOR_EXCEPTION
507 (bad, std::runtime_error, "Ifpack2::AdditiveSchwarz::apply: "
508 "At iteration " << ni << ", result of importMultiVector from R "
509 "to OverlappingB, has 2-norm NaN or Inf.");
510 }
511#endif // HAVE_IFPACK2_DEBUG
512 } else {
513 globalOverlappingB->doImport (*R, *DistributedImporter_, Tpetra::INSERT);
514
515#ifdef HAVE_IFPACK2_DEBUG
516 {
517 const bool bad = anyBad (*globalOverlappingB);
518 TEUCHOS_TEST_FOR_EXCEPTION
519 (bad, std::runtime_error, "Ifpack2::AdditiveSchwarz::apply: "
520 "At iteration " << ni << ", result of doImport from R, has 2-norm "
521 "NaN or Inf.");
522 }
523#endif // HAVE_IFPACK2_DEBUG
524 }
525
526#ifdef HAVE_IFPACK2_DEBUG
527 {
528 const bool bad = anyBad (*OverlappingB);
529 TEUCHOS_TEST_FOR_EXCEPTION
530 (bad, std::runtime_error, "Ifpack2::AdditiveSchwarz::apply: "
531 "At iteration " << ni << ", right before localApply, the 2-norm of "
532 "OverlappingB is NaN or Inf.");
533 }
534#endif // HAVE_IFPACK2_DEBUG
535
536 // local solve
537 localApply(*OverlappingB, *OverlappingY);
538
539#ifdef HAVE_IFPACK2_DEBUG
540 {
541 const bool bad = anyBad (*OverlappingY);
542 TEUCHOS_TEST_FOR_EXCEPTION
543 (bad, std::runtime_error, "Ifpack2::AdditiveSchwarz::apply: "
544 "At iteration " << ni << ", after localApply and before export / "
545 "copy, the 2-norm of OverlappingY is NaN or Inf.");
546 }
547#endif // HAVE_IFPACK2_DEBUG
548
549#ifdef HAVE_IFPACK2_DEBUG
550 {
551 const bool bad = anyBad (*C);
552 TEUCHOS_TEST_FOR_EXCEPTION
553 (bad, std::runtime_error, "Ifpack2::AdditiveSchwarz::apply: "
554 "At iteration " << ni << ", before export / copy, the 2-norm of C "
555 "is NaN or Inf.");
556 }
557#endif // HAVE_IFPACK2_DEBUG
558
559 // do communication if necessary
560 if (IsOverlapping_) {
561 TEUCHOS_TEST_FOR_EXCEPTION
562 (OverlappingMatrix_.is_null (), std::logic_error, prefix
563 << "OverlappingMatrix_ is null when it shouldn't be. "
564 "Please report this bug to the Ifpack2 developers.");
565 OverlappingMatrix_->exportMultiVector (*OverlappingY, *C, CombineMode_);
566
567 // average solution in overlap regions if requested via "schwarz: combine mode" "AVG"
568 if (AvgOverlap_) {
569 Teuchos::ArrayRCP<scalar_type> dataC = C->getDataNonConst(0);
570 for (int i = 0; i < (int) C->getMap()->getLocalNumElements(); i++) {
571 dataC[i] = dataC[i]/dataNumOverlapCopies[i];
572 }
573 }
574 }
575 else {
576 // mfh 16 Apr 2014: Make a view of Y with the same Map as
577 // OverlappingY, so that we can copy OverlappingY into Y. This
578 // replaces code that iterates over all entries of OverlappingY,
579 // copying them one at a time into Y. That code assumed that
580 // the rows of Y and the rows of OverlappingY have the same
581 // global indices in the same order; see Bug 5992.
582 RCP<MV> C_view = C->offsetViewNonConst (OverlappingY->getMap (), 0);
583 Tpetra::deep_copy (*C_view, *OverlappingY);
584 }
585
586#ifdef HAVE_IFPACK2_DEBUG
587 {
588 const bool bad = anyBad (*C);
589 TEUCHOS_TEST_FOR_EXCEPTION
590 (bad, std::runtime_error, "Ifpack2::AdditiveSchwarz::apply: "
591 "At iteration " << ni << ", before Y := C + Y, the 2-norm of C "
592 "is NaN or Inf.");
593 }
594#endif // HAVE_IFPACK2_DEBUG
595
596#ifdef HAVE_IFPACK2_DEBUG
597 {
598 const bool bad = anyBad (Y);
599 TEUCHOS_TEST_FOR_EXCEPTION
600 (bad, std::runtime_error, "Ifpack2::AdditiveSchwarz::apply: "
601 "Before Y := C + Y, at iteration " << ni << ", the 2-norm of Y "
602 "is NaN or Inf.");
603 }
604#endif // HAVE_IFPACK2_DEBUG
605
606 Y.update(UpdateDamping_, *C, STS::one());
607
608#ifdef HAVE_IFPACK2_DEBUG
609 {
610 const bool bad = anyBad (Y);
611 TEUCHOS_TEST_FOR_EXCEPTION
612 (bad, std::runtime_error, "Ifpack2::AdditiveSchwarz::apply: "
613 "At iteration " << ni << ", after Y := C + Y, the 2-norm of Y "
614 "is NaN or Inf.");
615 }
616#endif // HAVE_IFPACK2_DEBUG
617 } // for each iteration
618
619 } // Stop timing here
620
621#ifdef HAVE_IFPACK2_DEBUG
622 {
623 const bool bad = anyBad (Y);
624 TEUCHOS_TEST_FOR_EXCEPTION
625 (bad, std::runtime_error, "Ifpack2::AdditiveSchwarz::apply: "
626 "The 2-norm of the output Y is NaN or Inf.");
627 }
628#endif // HAVE_IFPACK2_DEBUG
629
630 ++NumApply_;
631
632 ApplyTime_ += (timer->wallTime() - startTime);
633}
634
635template<class MatrixType,class LocalInverseType>
636void
638localApply (MV& OverlappingB, MV& OverlappingY) const
639{
640 using Teuchos::RCP;
641 using Teuchos::rcp_dynamic_cast;
642
643 const size_t numVectors = OverlappingB.getNumVectors ();
644
645 auto additiveSchwarzFilter = rcp_dynamic_cast<Details::AdditiveSchwarzFilter<MatrixType>>(innerMatrix_);
646 if(additiveSchwarzFilter)
647 {
648 //Create the reduced system innerMatrix_ * ReducedY = ReducedB.
649 //This effectively fuses 3 tasks:
650 // -SingletonFilter::SolveSingletons (solve entries of OverlappingY corresponding to singletons)
651 // -SingletonFilter::CreateReducedRHS (fill ReducedReorderedB from OverlappingB, with entries in singleton columns eliminated)
652 // -ReorderFilter::permuteOriginalToReordered (apply permutation to ReducedReorderedB)
653 MV ReducedReorderedB (additiveSchwarzFilter->getRowMap(), numVectors);
654 MV ReducedReorderedY (additiveSchwarzFilter->getRowMap(), numVectors);
655 additiveSchwarzFilter->CreateReducedProblem(OverlappingB, OverlappingY, ReducedReorderedB);
656 //Apply inner solver
657 Inverse_->solve (ReducedReorderedY, ReducedReorderedB);
658 //Scatter ReducedY back to non-singleton rows of OverlappingY, according to the reordering.
659 additiveSchwarzFilter->UpdateLHS(ReducedReorderedY, OverlappingY);
660 }
661 else
662 {
663 if (FilterSingletons_) {
664 // process singleton filter
665 MV ReducedB (SingletonMatrix_->getRowMap (), numVectors);
666 MV ReducedY (SingletonMatrix_->getRowMap (), numVectors);
667
668 RCP<SingletonFilter<row_matrix_type> > singletonFilter =
669 rcp_dynamic_cast<SingletonFilter<row_matrix_type> > (SingletonMatrix_);
670 TEUCHOS_TEST_FOR_EXCEPTION
671 (! SingletonMatrix_.is_null () && singletonFilter.is_null (),
672 std::logic_error, "Ifpack2::AdditiveSchwarz::localApply: "
673 "SingletonFilter_ is nonnull but is not a SingletonFilter"
674 "<row_matrix_type>. This should never happen. Please report this bug "
675 "to the Ifpack2 developers.");
676 singletonFilter->SolveSingletons (OverlappingB, OverlappingY);
677 singletonFilter->CreateReducedRHS (OverlappingY, OverlappingB, ReducedB);
678
679 // process reordering
680 if (! UseReordering_) {
681 Inverse_->solve (ReducedY, ReducedB);
682 }
683 else {
684 RCP<ReorderFilter<row_matrix_type> > rf =
685 rcp_dynamic_cast<ReorderFilter<row_matrix_type> > (ReorderedLocalizedMatrix_);
686 TEUCHOS_TEST_FOR_EXCEPTION
687 (! ReorderedLocalizedMatrix_.is_null () && rf.is_null (), std::logic_error,
688 "Ifpack2::AdditiveSchwarz::localApply: ReorderedLocalizedMatrix_ is "
689 "nonnull but is not a ReorderFilter<row_matrix_type>. This should "
690 "never happen. Please report this bug to the Ifpack2 developers.");
691 MV ReorderedB (ReducedB, Teuchos::Copy);
692 MV ReorderedY (ReducedY, Teuchos::Copy);
693 rf->permuteOriginalToReordered (ReducedB, ReorderedB);
694 Inverse_->solve (ReorderedY, ReorderedB);
695 rf->permuteReorderedToOriginal (ReorderedY, ReducedY);
696 }
697
698 // finish up with singletons
699 singletonFilter->UpdateLHS (ReducedY, OverlappingY);
700 }
701 else {
702 // process reordering
703 if (! UseReordering_) {
704 Inverse_->solve (OverlappingY, OverlappingB);
705 }
706 else {
707 MV ReorderedB (OverlappingB, Teuchos::Copy);
708 MV ReorderedY (OverlappingY, Teuchos::Copy);
709
710 RCP<ReorderFilter<row_matrix_type> > rf =
711 rcp_dynamic_cast<ReorderFilter<row_matrix_type> > (ReorderedLocalizedMatrix_);
712 TEUCHOS_TEST_FOR_EXCEPTION
713 (! ReorderedLocalizedMatrix_.is_null () && rf.is_null (), std::logic_error,
714 "Ifpack2::AdditiveSchwarz::localApply: ReorderedLocalizedMatrix_ is "
715 "nonnull but is not a ReorderFilter<row_matrix_type>. This should "
716 "never happen. Please report this bug to the Ifpack2 developers.");
717 rf->permuteOriginalToReordered (OverlappingB, ReorderedB);
718 Inverse_->solve (ReorderedY, ReorderedB);
719 rf->permuteReorderedToOriginal (ReorderedY, OverlappingY);
720 }
721 }
722 }
723}
724
725
726template<class MatrixType,class LocalInverseType>
728setParameters (const Teuchos::ParameterList& plist)
729{
730 // mfh 18 Nov 2013: Ifpack2's setParameters() method passes in the
731 // input list as const. This means that we have to copy it before
732 // validation or passing into setParameterList().
733 List_ = plist;
734 this->setParameterList (Teuchos::rcpFromRef (List_));
735}
736
737
738
739template<class MatrixType,class LocalInverseType>
741setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& plist)
742{
743 using Tpetra::CombineMode;
744 using Teuchos::ParameterEntry;
745 using Teuchos::ParameterEntryValidator;
746 using Teuchos::ParameterList;
747 using Teuchos::RCP;
748 using Teuchos::rcp;
749 using Teuchos::rcp_dynamic_cast;
750 using Teuchos::StringToIntegralParameterEntryValidator;
751 using Details::getParamTryingTypes;
752 const char prefix[] = "Ifpack2::AdditiveSchwarz: ";
753
754 if (plist.is_null ()) {
755 // Assume that the user meant to set default parameters by passing
756 // in an empty list.
757 this->setParameterList (rcp (new ParameterList ()));
758 }
759 // FIXME (mfh 26 Aug 2015) It's not necessarily true that plist is
760 // nonnull at this point.
761
762 // At this point, plist should be nonnull.
763 TEUCHOS_TEST_FOR_EXCEPTION(
764 plist.is_null (), std::logic_error, "Ifpack2::AdditiveSchwarz::"
765 "setParameterList: plist is null. This should never happen, since the "
766 "method should have replaced a null input list with a nonnull empty list "
767 "by this point. Please report this bug to the Ifpack2 developers.");
768
769 // TODO JJH 24March2015 The list needs to be validated. Not sure why this is commented out.
770 // try {
771 // List_.validateParameters (* getValidParameters ());
772 // }
773 // catch (std::exception& e) {
774 // std::cerr << "Ifpack2::AdditiveSchwarz::setParameterList: Validation failed with the following error message: " << e.what () << std::endl;
775 // throw e;
776 // }
777
778 // mfh 18 Nov 2013: Supplying the current value as the default value
779 // when calling ParameterList::get() ensures "delta" behavior when
780 // users pass in new parameters: any unspecified parameters in the
781 // new list retain their values in the old list. This preserves
782 // backwards compatiblity with this class' previous behavior. Note
783 // that validateParametersAndSetDefaults() would have different
784 // behavior: any parameters not in the new list would get default
785 // values, which could be different than their values in the
786 // original list.
787
788 const std::string cmParamName ("schwarz: combine mode");
789 const ParameterEntry* cmEnt = plist->getEntryPtr (cmParamName);
790 if (cmEnt != nullptr) {
791 if (cmEnt->isType<CombineMode> ()) {
792 CombineMode_ = Teuchos::getValue<CombineMode> (*cmEnt);
793 }
794 else if (cmEnt->isType<int> ()) {
795 const int cm = Teuchos::getValue<int> (*cmEnt);
796 CombineMode_ = static_cast<CombineMode> (cm);
797 }
798 else if (cmEnt->isType<std::string> ()) {
799 // Try to get the combine mode as a string. If this works, use
800 // the validator to convert to int. This is painful, but
801 // necessary in order to do validation, since the input list may
802 // not necessarily come with a validator.
803 const ParameterEntry& validEntry =
804 getValidParameters ()->getEntry (cmParamName);
805 RCP<const ParameterEntryValidator> v = validEntry.validator ();
806 using vs2e_type = StringToIntegralParameterEntryValidator<CombineMode>;
807 RCP<const vs2e_type> vs2e = rcp_dynamic_cast<const vs2e_type> (v, true);
808
809 ParameterEntry& inputEntry = plist->getEntry (cmParamName);
810 // As AVG is only a Schwarz option and does not exist in Tpetra's
811 // version of CombineMode, we use a separate boolean local to
812 // Schwarz in conjunction with CombineMode_ == ADD to handle
813 // averaging. Here, we change input entry to ADD and set the boolean.
814 if (strncmp(Teuchos::getValue<std::string>(inputEntry).c_str(),"AVG",3) == 0) {
815 inputEntry.template setValue<std::string>("ADD");
816 AvgOverlap_ = true;
817 }
818 CombineMode_ = vs2e->getIntegralValue (inputEntry, cmParamName);
819 }
820 }
821 // If doing user partitioning with Block Jacobi relaxation and overlapping blocks, we might
822 // later need to know whether or not the overlapping Schwarz scheme is "ADD" or "ZERO" (which
823 // is really RAS Schwarz. If it is "ADD", communication will be necessary when computing the
824 // proper weights needed to combine solution values in overlap regions
825 if (plist->isParameter("subdomain solver name")) {
826 if (plist->get<std::string>("subdomain solver name") == "BLOCK_RELAXATION") {
827 if (plist->isSublist("subdomain solver parameters")) {
828 if (plist->sublist("subdomain solver parameters").isParameter("relaxation: type")) {
829 if (plist->sublist("subdomain solver parameters").get<std::string>("relaxation: type") == "Jacobi" ) {
830 if (plist->sublist("subdomain solver parameters").isParameter("partitioner: type")) {
831 if (plist->sublist("subdomain solver parameters").get<std::string>("partitioner: type") == "user") {
832 if (CombineMode_ == Tpetra::ADD) plist->sublist("subdomain solver parameters").set("partitioner: combine mode","ADD");
833 if (CombineMode_ == Tpetra::ZERO) plist->sublist("subdomain solver parameters").set("partitioner: combine mode","ZERO");
834 AvgOverlap_ = false; // averaging already taken care of by the partitioner: nonsymmetric overlap combine option
835 }
836 }
837 }
838 }
839 }
840 }
841 }
842
843 OverlapLevel_ = plist->get ("schwarz: overlap level", OverlapLevel_);
844
845 // We set IsOverlapping_ in initialize(), once we know that Matrix_ is nonnull.
846
847 // Will we be doing reordering? Unlike Ifpack, we'll use a
848 // "schwarz: reordering list" to give to Zoltan2.
849 UseReordering_ = plist->get ("schwarz: use reordering", UseReordering_);
850
851#if !defined(HAVE_IFPACK2_XPETRA) || !defined(HAVE_IFPACK2_ZOLTAN2)
852 TEUCHOS_TEST_FOR_EXCEPTION(
853 UseReordering_, std::invalid_argument, "Ifpack2::AdditiveSchwarz::"
854 "setParameters: You specified \"schwarz: use reordering\" = true. "
855 "This is only valid when Trilinos was built with Ifpack2, Xpetra, and "
856 "Zoltan2 enabled. Either Xpetra or Zoltan2 was not enabled in your build "
857 "of Trilinos.");
858#endif
859
860 // FIXME (mfh 18 Nov 2013) Now would be a good time to validate the
861 // "schwarz: reordering list" parameter list. Currently, that list
862 // gets extracted in setup().
863
864 // if true, filter singletons. NOTE: the filtered matrix can still have
865 // singletons! A simple example: upper triangular matrix, if I remove
866 // the lower node, I still get a matrix with a singleton! However, filter
867 // singletons should help for PDE problems with Dirichlet BCs.
868 FilterSingletons_ = plist->get ("schwarz: filter singletons", FilterSingletons_);
869
870 // Allow for damped Schwarz updates
871 getParamTryingTypes<scalar_type, scalar_type, double>
872 (UpdateDamping_, *plist, "schwarz: update damping", prefix);
873
874 // If the inner solver doesn't exist yet, don't create it.
875 // initialize() creates it.
876 //
877 // If the inner solver _does_ exist, there are three cases,
878 // depending on what the user put in the input ParameterList.
879 //
880 // 1. The user did /not/ provide a parameter specifying the inner
881 // solver's type, nor did the user specify a sublist of
882 // parameters for the inner solver
883 // 2. The user did /not/ provide a parameter specifying the inner
884 // solver's type, but /did/ specify a sublist of parameters for
885 // the inner solver
886 // 3. The user provided a parameter specifying the inner solver's
887 // type (it does not matter in this case whether the user gave
888 // a sublist of parameters for the inner solver)
889 //
890 // AdditiveSchwarz has "delta" (relative) semantics for setting
891 // parameters. This means that if the user did not specify the
892 // inner solver's type, we presume that the type has not changed.
893 // Thus, if the inner solver exists, we don't need to recreate it.
894 //
895 // In Case 3, if the user bothered to specify the inner solver's
896 // type, then we must assume it may differ than the current inner
897 // solver's type. Thus, we have to recreate the inner solver. We
898 // achieve this here by assigning null to Inverse_; initialize()
899 // will recreate the solver when it is needed. Our assumption here
900 // is necessary because Ifpack2::Preconditioner does not have a
901 // method for querying a preconditioner's "type" (i.e., name) as a
902 // string. Remember that the user may have previously set an
903 // arbitrary inner solver by calling setInnerPreconditioner().
904 //
905 // See note at the end of setInnerPreconditioner().
906
907 if (! Inverse_.is_null ()) {
908 // "CUSTOM" explicitly indicates that the user called or plans to
909 // call setInnerPreconditioner.
910 if (hasInnerPrecName () && innerPrecName () != "CUSTOM") {
911 // Wipe out the current inner solver. initialize() will
912 // recreate it with the correct type.
913 Inverse_ = Teuchos::null;
914 }
915 else {
916 // Extract and apply the sublist of parameters to give to the
917 // inner solver, if there is such a sublist of parameters.
918 std::pair<ParameterList, bool> result = innerPrecParams ();
919 if (result.second) {
920 // FIXME (mfh 26 Aug 2015) Rewrite innerPrecParams() so this
921 // isn't another deep copy.
922 Inverse_->setParameters (rcp (new ParameterList (result.first)));
923 }
924 }
925 }
926
927 NumIterations_ = plist->get ("schwarz: num iterations", NumIterations_);
928 ZeroStartingSolution_ =
929 plist->get ("schwarz: zero starting solution", ZeroStartingSolution_);
930}
931
932
933
934template<class MatrixType,class LocalInverseType>
935Teuchos::RCP<const Teuchos::ParameterList>
937getValidParameters () const
938{
939 using Teuchos::ParameterList;
940 using Teuchos::parameterList;
941 using Teuchos::RCP;
942 using Teuchos::rcp_const_cast;
943
944 if (validParams_.is_null ()) {
945 const int overlapLevel = 0;
946 const bool useReordering = false;
947 const bool filterSingletons = false;
948 const int numIterations = 1;
949 const bool zeroStartingSolution = true;
950 const scalar_type updateDamping = Teuchos::ScalarTraits<scalar_type>::one ();
951 ParameterList reorderingSublist;
952 reorderingSublist.set ("order_method", std::string ("rcm"));
953
954 RCP<ParameterList> plist = parameterList ("Ifpack2::AdditiveSchwarz");
955
956 Tpetra::setCombineModeParameter (*plist, "schwarz: combine mode");
957 plist->set ("schwarz: overlap level", overlapLevel);
958 plist->set ("schwarz: use reordering", useReordering);
959 plist->set ("schwarz: reordering list", reorderingSublist);
960 // mfh 24 Mar 2015: We accept this for backwards compatibility
961 // ONLY. It is IGNORED.
962 plist->set ("schwarz: compute condest", false);
963 plist->set ("schwarz: filter singletons", filterSingletons);
964 plist->set ("schwarz: num iterations", numIterations);
965 plist->set ("schwarz: zero starting solution", zeroStartingSolution);
966 plist->set ("schwarz: update damping", updateDamping);
967
968 // FIXME (mfh 18 Nov 2013) Get valid parameters from inner solver.
969 // JJH The inner solver should handle its own validation.
970 //
971 // FIXME (mfh 18 Nov 2013) Get valid parameters from Zoltan2, if
972 // Zoltan2 was enabled in the build.
973 // JJH Zoltan2 should handle its own validation.
974 //
975
976 validParams_ = rcp_const_cast<const ParameterList> (plist);
977 }
978 return validParams_;
979}
980
981
982template<class MatrixType,class LocalInverseType>
984{
985 using Tpetra::global_size_t;
986 using Teuchos::RCP;
987 using Teuchos::rcp;
988 using Teuchos::SerialComm;
989 using Teuchos::Time;
990 using Teuchos::TimeMonitor;
991
992 const std::string timerName ("Ifpack2::AdditiveSchwarz::initialize");
993 RCP<Time> timer = TimeMonitor::lookupCounter (timerName);
994 if (timer.is_null ()) {
995 timer = TimeMonitor::getNewCounter (timerName);
996 }
997 double startTime = timer->wallTime();
998
999 { // Start timing here.
1000 TimeMonitor timeMon (*timer);
1001
1002 TEUCHOS_TEST_FOR_EXCEPTION(
1003 Matrix_.is_null (), std::runtime_error, "Ifpack2::AdditiveSchwarz::"
1004 "initialize: The matrix to precondition is null. You must either pass "
1005 "a nonnull matrix to the constructor, or call setMatrix() with a nonnull "
1006 "input, before you may call this method.");
1007
1008 IsInitialized_ = false;
1009 IsComputed_ = false;
1010 overlapping_B_.reset (nullptr);
1011 overlapping_Y_.reset (nullptr);
1012 R_.reset (nullptr);
1013 C_.reset (nullptr);
1014
1015 RCP<const Teuchos::Comm<int> > comm = Matrix_->getComm ();
1016 RCP<const map_type> rowMap = Matrix_->getRowMap ();
1017 const global_size_t INVALID =
1018 Teuchos::OrdinalTraits<global_size_t>::invalid ();
1019
1020 // If there's only one process in the matrix's communicator,
1021 // then there's no need to compute overlap.
1022 if (comm->getSize () == 1) {
1023 OverlapLevel_ = 0;
1024 IsOverlapping_ = false;
1025 } else if (OverlapLevel_ != 0) {
1026 IsOverlapping_ = true;
1027 }
1028
1029 if (OverlapLevel_ == 0) {
1030 const global_ordinal_type indexBase = rowMap->getIndexBase ();
1031 RCP<const SerialComm<int> > localComm (new SerialComm<int> ());
1032 // FIXME (mfh 15 Apr 2014) What if indexBase isn't the least
1033 // global index in the list of GIDs on this process?
1034 localMap_ =
1035 rcp (new map_type (INVALID, rowMap->getLocalNumElements (),
1036 indexBase, localComm));
1037 }
1038
1039 // compute the overlapping matrix if necessary
1040 if (IsOverlapping_) {
1041 Teuchos::TimeMonitor t(*Teuchos::TimeMonitor::getNewTimer("OverlappingRowMatrix construction"));
1042 OverlappingMatrix_ = rcp (new OverlappingRowMatrix<row_matrix_type> (Matrix_, OverlapLevel_));
1043 }
1044
1045 setup (); // This does a lot of the initialization work.
1046
1047 if (! Inverse_.is_null ()) {
1048 Inverse_->symbolic (); // Initialize subdomain solver.
1049 }
1050
1051 } // Stop timing here.
1052
1053 IsInitialized_ = true;
1054 ++NumInitialize_;
1055
1056 InitializeTime_ += (timer->wallTime() - startTime);
1057}
1058
1059template<class MatrixType,class LocalInverseType>
1061{
1062 return IsInitialized_;
1063}
1064
1065
1066template<class MatrixType,class LocalInverseType>
1068{
1069 using Teuchos::RCP;
1070 using Teuchos::Time;
1071 using Teuchos::TimeMonitor;
1072
1073 if (! IsInitialized_) {
1074 initialize ();
1075 }
1076
1077 TEUCHOS_TEST_FOR_EXCEPTION(
1078 ! isInitialized (), std::logic_error, "Ifpack2::AdditiveSchwarz::compute: "
1079 "The preconditioner is not yet initialized, "
1080 "even though initialize() supposedly has been called. "
1081 "This should never happen. "
1082 "Please report this bug to the Ifpack2 developers.");
1083
1084 TEUCHOS_TEST_FOR_EXCEPTION(
1085 Inverse_.is_null (), std::runtime_error,
1086 "Ifpack2::AdditiveSchwarz::compute: The subdomain solver is null. "
1087 "This can only happen if you called setInnerPreconditioner() with a null "
1088 "input, after calling initialize() or compute(). If you choose to call "
1089 "setInnerPreconditioner() with a null input, you must then call it with a "
1090 "nonnull input before you may call initialize() or compute().");
1091
1092 const std::string timerName ("Ifpack2::AdditiveSchwarz::compute");
1093 RCP<Time> timer = TimeMonitor::lookupCounter (timerName);
1094 if (timer.is_null ()) {
1095 timer = TimeMonitor::getNewCounter (timerName);
1096 }
1097 TimeMonitor timeMon (*timer);
1098 double startTime = timer->wallTime();
1099
1100 // compute () assumes that the values of Matrix_ (aka A) have changed.
1101 // If this has overlap, do an import from the input matrix to the halo.
1102 if (IsOverlapping_) {
1103 Teuchos::TimeMonitor t(*Teuchos::TimeMonitor::getNewTimer("Halo Import"));
1104 OverlappingMatrix_->doExtImport();
1105 }
1106 // At this point, either Matrix_ or OverlappingMatrix_ (depending on whether this is overlapping)
1107 // has new values and unchanged structure. If we are using AdditiveSchwarzFilter, update the local matrix.
1108 //
1109 if(auto asf = Teuchos::rcp_dynamic_cast<Details::AdditiveSchwarzFilter<MatrixType>>(innerMatrix_))
1110 {
1111 Teuchos::TimeMonitor t(*Teuchos::TimeMonitor::getNewTimer("Fill Local Matrix"));
1112 //NOTE: if this compute() call comes right after the initialize() with no intervening matrix changes, this call is redundant.
1113 //initialize() already filled the local matrix. However, we have no way to tell if this is the case.
1114 asf->updateMatrixValues();
1115 }
1116 //Now, whether the Inverse_'s matrix is the AdditiveSchwarzFilter's local matrix or simply Matrix_/OverlappingMatrix_,
1117 //it will be able to see the new values and update itself accordingly.
1118
1119 { // Start timing here.
1120
1121 IsComputed_ = false;
1122 Inverse_->numeric ();
1123 } // Stop timing here.
1124
1125 IsComputed_ = true;
1126 ++NumCompute_;
1127
1128 ComputeTime_ += (timer->wallTime() - startTime);
1129}
1130
1131//==============================================================================
1132// Returns true if the preconditioner has been successfully computed, false otherwise.
1133template<class MatrixType,class LocalInverseType>
1135{
1136 return IsComputed_;
1137}
1138
1139
1140template<class MatrixType,class LocalInverseType>
1142{
1143 return NumInitialize_;
1144}
1145
1146
1147template<class MatrixType,class LocalInverseType>
1149{
1150 return NumCompute_;
1151}
1152
1153
1154template<class MatrixType,class LocalInverseType>
1156{
1157 return NumApply_;
1158}
1159
1160
1161template<class MatrixType,class LocalInverseType>
1163{
1164 return InitializeTime_;
1165}
1166
1167
1168template<class MatrixType,class LocalInverseType>
1170{
1171 return ComputeTime_;
1172}
1173
1174
1175template<class MatrixType,class LocalInverseType>
1177{
1178 return ApplyTime_;
1179}
1180
1181
1182template<class MatrixType,class LocalInverseType>
1184{
1185 std::ostringstream out;
1186
1187 out << "\"Ifpack2::AdditiveSchwarz\": {";
1188 if (this->getObjectLabel () != "") {
1189 out << "Label: \"" << this->getObjectLabel () << "\", ";
1190 }
1191 out << "Initialized: " << (isInitialized () ? "true" : "false")
1192 << ", Computed: " << (isComputed () ? "true" : "false")
1193 << ", Iterations: " << NumIterations_
1194 << ", Overlap level: " << OverlapLevel_
1195 << ", Subdomain reordering: \"" << ReorderingAlgorithm_ << "\"";
1196 out << ", Combine mode: \"";
1197 if (CombineMode_ == Tpetra::INSERT) {
1198 out << "INSERT";
1199 } else if (CombineMode_ == Tpetra::ADD) {
1200 out << "ADD";
1201 } else if (CombineMode_ == Tpetra::REPLACE) {
1202 out << "REPLACE";
1203 } else if (CombineMode_ == Tpetra::ABSMAX) {
1204 out << "ABSMAX";
1205 } else if (CombineMode_ == Tpetra::ZERO) {
1206 out << "ZERO";
1207 }
1208 out << "\"";
1209 if (Matrix_.is_null ()) {
1210 out << ", Matrix: null";
1211 }
1212 else {
1213 out << ", Global matrix dimensions: ["
1214 << Matrix_->getGlobalNumRows () << ", "
1215 << Matrix_->getGlobalNumCols () << "]";
1216 }
1217 out << ", Inner solver: ";
1218 if (! Inverse_.is_null ()) {
1219 Teuchos::RCP<Teuchos::Describable> inv =
1220 Teuchos::rcp_dynamic_cast<Teuchos::Describable> (Inverse_);
1221 if (! inv.is_null ()) {
1222 out << "{" << inv->description () << "}";
1223 } else {
1224 out << "{" << "Some inner solver" << "}";
1225 }
1226 } else {
1227 out << "null";
1228 }
1229
1230 out << "}";
1231 return out.str ();
1232}
1233
1234
1235template<class MatrixType,class LocalInverseType>
1236void
1238describe (Teuchos::FancyOStream& out,
1239 const Teuchos::EVerbosityLevel verbLevel) const
1240{
1241 using Teuchos::OSTab;
1242 using Teuchos::TypeNameTraits;
1243 using std::endl;
1244
1245 const int myRank = Matrix_->getComm ()->getRank ();
1246 const int numProcs = Matrix_->getComm ()->getSize ();
1247 const Teuchos::EVerbosityLevel vl =
1248 (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
1249
1250 if (vl > Teuchos::VERB_NONE) {
1251 // describe() starts with a tab, by convention.
1252 OSTab tab0 (out);
1253 if (myRank == 0) {
1254 out << "\"Ifpack2::AdditiveSchwarz\":";
1255 }
1256 OSTab tab1 (out);
1257 if (myRank == 0) {
1258 out << "MatrixType: " << TypeNameTraits<MatrixType>::name () << endl;
1259 out << "LocalInverseType: " << TypeNameTraits<LocalInverseType>::name () << endl;
1260 if (this->getObjectLabel () != "") {
1261 out << "Label: \"" << this->getObjectLabel () << "\"" << endl;
1262 }
1263
1264 out << "Overlap level: " << OverlapLevel_ << endl
1265 << "Combine mode: \"";
1266 if (CombineMode_ == Tpetra::INSERT) {
1267 out << "INSERT";
1268 } else if (CombineMode_ == Tpetra::ADD) {
1269 out << "ADD";
1270 } else if (CombineMode_ == Tpetra::REPLACE) {
1271 out << "REPLACE";
1272 } else if (CombineMode_ == Tpetra::ABSMAX) {
1273 out << "ABSMAX";
1274 } else if (CombineMode_ == Tpetra::ZERO) {
1275 out << "ZERO";
1276 }
1277 out << "\"" << endl
1278 << "Subdomain reordering: \"" << ReorderingAlgorithm_ << "\"" << endl;
1279 }
1280
1281 if (Matrix_.is_null ()) {
1282 if (myRank == 0) {
1283 out << "Matrix: null" << endl;
1284 }
1285 }
1286 else {
1287 if (myRank == 0) {
1288 out << "Matrix:" << endl;
1289 std::flush (out);
1290 }
1291 Matrix_->getComm ()->barrier (); // wait for output to finish
1292 Matrix_->describe (out, Teuchos::VERB_LOW);
1293 }
1294
1295 if (myRank == 0) {
1296 out << "Number of initialize calls: " << getNumInitialize () << endl
1297 << "Number of compute calls: " << getNumCompute () << endl
1298 << "Number of apply calls: " << getNumApply () << endl
1299 << "Total time in seconds for initialize: " << getInitializeTime () << endl
1300 << "Total time in seconds for compute: " << getComputeTime () << endl
1301 << "Total time in seconds for apply: " << getApplyTime () << endl;
1302 }
1303
1304 if (Inverse_.is_null ()) {
1305 if (myRank == 0) {
1306 out << "Subdomain solver: null" << endl;
1307 }
1308 }
1309 else {
1310 if (vl < Teuchos::VERB_EXTREME) {
1311 if (myRank == 0) {
1312 out << "Subdomain solver: not null" << endl;
1313 }
1314 }
1315 else { // vl >= Teuchos::VERB_EXTREME
1316 for (int p = 0; p < numProcs; ++p) {
1317 if (p == myRank) {
1318 out << "Subdomain solver on Process " << myRank << ":";
1319 if (Inverse_.is_null ()) {
1320 out << "null" << endl;
1321 } else {
1322 Teuchos::RCP<Teuchos::Describable> inv =
1323 Teuchos::rcp_dynamic_cast<Teuchos::Describable> (Inverse_);
1324 if (! inv.is_null ()) {
1325 out << endl;
1326 inv->describe (out, vl);
1327 } else {
1328 out << "null" << endl;
1329 }
1330 }
1331 }
1332 Matrix_->getComm ()->barrier ();
1333 Matrix_->getComm ()->barrier ();
1334 Matrix_->getComm ()->barrier (); // wait for output to finish
1335 }
1336 }
1337 }
1338
1339 Matrix_->getComm ()->barrier (); // wait for output to finish
1340 }
1341}
1342
1343
1344template<class MatrixType,class LocalInverseType>
1345std::ostream& AdditiveSchwarz<MatrixType,LocalInverseType>::print(std::ostream& os) const
1346{
1347 Teuchos::FancyOStream fos(Teuchos::rcp(&os,false));
1348 fos.setOutputToRootOnly(0);
1349 describe(fos);
1350 return(os);
1351}
1352
1353
1354template<class MatrixType,class LocalInverseType>
1356{
1357 return OverlapLevel_;
1358}
1359
1360
1361template<class MatrixType,class LocalInverseType>
1363{
1364#ifdef HAVE_MPI
1365 using Teuchos::MpiComm;
1366#endif // HAVE_MPI
1367 using Teuchos::ArrayRCP;
1368 using Teuchos::ParameterList;
1369 using Teuchos::RCP;
1370 using Teuchos::rcp;
1371 using Teuchos::rcp_dynamic_cast;
1372 using Teuchos::rcpFromRef;
1373
1374 TEUCHOS_TEST_FOR_EXCEPTION(
1375 Matrix_.is_null (), std::runtime_error, "Ifpack2::AdditiveSchwarz::"
1376 "initialize: The matrix to precondition is null. You must either pass "
1377 "a nonnull matrix to the constructor, or call setMatrix() with a nonnull "
1378 "input, before you may call this method.");
1379
1380 //If the matrix is a CrsMatrix or OverlappingRowMatrix, use the high-performance
1381 //AdditiveSchwarzFilter. Otherwise, use composition of Reordered/Singleton/LocalFilter.
1382 auto matrixCrs = rcp_dynamic_cast<const crs_matrix_type>(Matrix_);
1383 if(!OverlappingMatrix_.is_null() || !matrixCrs.is_null())
1384 {
1385 ArrayRCP<local_ordinal_type> perm;
1386 ArrayRCP<local_ordinal_type> revperm;
1387 if (UseReordering_) {
1388 Teuchos::TimeMonitor t(*Teuchos::TimeMonitor::getNewTimer("Reordering"));
1389#if defined(HAVE_IFPACK2_XPETRA) && defined(HAVE_IFPACK2_ZOLTAN2)
1390 // Unlike Ifpack, Zoltan2 does all the dirty work here.
1391 Teuchos::ParameterList zlist = List_.sublist ("schwarz: reordering list");
1392 ReorderingAlgorithm_ = zlist.get<std::string> ("order_method", "rcm");
1393
1394 if(ReorderingAlgorithm_ == "user") {
1395 // User-provided reordering
1396 perm = zlist.get<Teuchos::ArrayRCP<local_ordinal_type> >("user ordering");
1397 revperm = zlist.get<Teuchos::ArrayRCP<local_ordinal_type> >("user reverse ordering");
1398 }
1399 else {
1400 // Zoltan2 reordering
1401 typedef Tpetra::RowGraph
1402 <local_ordinal_type, global_ordinal_type, node_type> row_graph_type;
1403 typedef Zoltan2::TpetraRowGraphAdapter<row_graph_type> z2_adapter_type;
1404 auto constActiveGraph = Teuchos::rcp_const_cast<const row_graph_type>(
1405 IsOverlapping_ ? OverlappingMatrix_->getGraph() : Matrix_->getGraph());
1406 z2_adapter_type Zoltan2Graph (constActiveGraph);
1407
1408 typedef Zoltan2::OrderingProblem<z2_adapter_type> ordering_problem_type;
1409#ifdef HAVE_MPI
1410 // Grab the MPI Communicator and build the ordering problem with that
1411 MPI_Comm myRawComm;
1412
1413 RCP<const MpiComm<int> > mpicomm =
1414 rcp_dynamic_cast<const MpiComm<int> > (Matrix_->getComm ());
1415 if (mpicomm == Teuchos::null) {
1416 myRawComm = MPI_COMM_SELF;
1417 } else {
1418 myRawComm = * (mpicomm->getRawMpiComm ());
1419 }
1420 ordering_problem_type MyOrderingProblem (&Zoltan2Graph, &zlist, myRawComm);
1421#else
1422 ordering_problem_type MyOrderingProblem (&Zoltan2Graph, &zlist);
1423#endif
1424 MyOrderingProblem.solve ();
1425
1426 {
1427 typedef Zoltan2::LocalOrderingSolution<local_ordinal_type>
1428 ordering_solution_type;
1429
1430 ordering_solution_type sol (*MyOrderingProblem.getLocalOrderingSolution());
1431
1432 // perm[i] gives the where OLD index i shows up in the NEW
1433 // ordering. revperm[i] gives the where NEW index i shows
1434 // up in the OLD ordering. Note that perm is actually the
1435 // "inverse permutation," in Zoltan2 terms.
1436 perm = sol.getPermutationRCPConst (true);
1437 revperm = sol.getPermutationRCPConst ();
1438 }
1439 }
1440#else
1441 // This is a logic_error, not a runtime_error, because
1442 // setParameters() should have excluded this case already.
1443 TEUCHOS_TEST_FOR_EXCEPTION(
1444 true, std::logic_error, "Ifpack2::AdditiveSchwarz::setup: "
1445 "The Zoltan2 and Xpetra packages must be enabled in order "
1446 "to support reordering.");
1447#endif
1448 }
1449 else
1450 {
1451 local_ordinal_type numLocalRows = OverlappingMatrix_.is_null() ? matrixCrs->getLocalNumRows() : OverlappingMatrix_->getLocalNumRows();
1452 //Use an identity ordering.
1453 //TODO: create a non-permuted code path in AdditiveSchwarzFilter, in the case that neither
1454 //reordering nor singleton filtering are enabled. In this situation it's like LocalFilter.
1455 perm = ArrayRCP<local_ordinal_type>(numLocalRows);
1456 revperm = ArrayRCP<local_ordinal_type>(numLocalRows);
1457 for(local_ordinal_type i = 0; i < numLocalRows; i++)
1458 {
1459 perm[i] = i;
1460 revperm[i] = i;
1461 }
1462 }
1463 //Now, construct the filter
1464 {
1465 Teuchos::TimeMonitor t(*Teuchos::TimeMonitor::getNewTimer("Filter construction"));
1466 RCP<Details::AdditiveSchwarzFilter<MatrixType>> asf;
1467 if(OverlappingMatrix_.is_null())
1468 asf = rcp(new Details::AdditiveSchwarzFilter<MatrixType>(matrixCrs, perm, revperm, FilterSingletons_));
1469 else
1470 asf = rcp(new Details::AdditiveSchwarzFilter<MatrixType>(OverlappingMatrix_, perm, revperm, FilterSingletons_));
1471 innerMatrix_ = asf;
1472 }
1473 }
1474 else
1475 {
1476 // Localized version of Matrix_ or OverlappingMatrix_.
1477 RCP<row_matrix_type> LocalizedMatrix;
1478
1479 // The "most current local matrix." At the end of this method, this
1480 // will be handed off to the inner solver.
1481 RCP<row_matrix_type> ActiveMatrix;
1482
1483 // Create localized matrix.
1484 if (! OverlappingMatrix_.is_null ()) {
1485 LocalizedMatrix = rcp (new LocalFilter<row_matrix_type> (OverlappingMatrix_));
1486 }
1487 else {
1488 LocalizedMatrix = rcp (new LocalFilter<row_matrix_type> (Matrix_));
1489 }
1490
1491 // Sanity check; I don't trust the logic above to have created LocalizedMatrix.
1492 TEUCHOS_TEST_FOR_EXCEPTION(
1493 LocalizedMatrix.is_null (), std::logic_error,
1494 "Ifpack2::AdditiveSchwarz::setup: LocalizedMatrix is null, after the code "
1495 "that claimed to have created it. This should never be the case. Please "
1496 "report this bug to the Ifpack2 developers.");
1497
1498 // Mark localized matrix as active
1499 ActiveMatrix = LocalizedMatrix;
1500
1501 // Singleton Filtering
1502 if (FilterSingletons_) {
1503 SingletonMatrix_ = rcp (new SingletonFilter<row_matrix_type> (LocalizedMatrix));
1504 ActiveMatrix = SingletonMatrix_;
1505 }
1506
1507 // Do reordering
1508 if (UseReordering_) {
1509#if defined(HAVE_IFPACK2_XPETRA) && defined(HAVE_IFPACK2_ZOLTAN2)
1510 // Unlike Ifpack, Zoltan2 does all the dirty work here.
1511 typedef ReorderFilter<row_matrix_type> reorder_filter_type;
1512 Teuchos::ParameterList zlist = List_.sublist ("schwarz: reordering list");
1513 ReorderingAlgorithm_ = zlist.get<std::string> ("order_method", "rcm");
1514
1515 ArrayRCP<local_ordinal_type> perm;
1516 ArrayRCP<local_ordinal_type> revperm;
1517
1518 if(ReorderingAlgorithm_ == "user") {
1519 // User-provided reordering
1520 perm = zlist.get<Teuchos::ArrayRCP<local_ordinal_type> >("user ordering");
1521 revperm = zlist.get<Teuchos::ArrayRCP<local_ordinal_type> >("user reverse ordering");
1522 }
1523 else {
1524 // Zoltan2 reordering
1525 typedef Tpetra::RowGraph
1526 <local_ordinal_type, global_ordinal_type, node_type> row_graph_type;
1527 typedef Zoltan2::TpetraRowGraphAdapter<row_graph_type> z2_adapter_type;
1528 RCP<const row_graph_type> constActiveGraph =
1529 Teuchos::rcp_const_cast<const row_graph_type>(ActiveMatrix->getGraph());
1530 z2_adapter_type Zoltan2Graph (constActiveGraph);
1531
1532 typedef Zoltan2::OrderingProblem<z2_adapter_type> ordering_problem_type;
1533#ifdef HAVE_MPI
1534 // Grab the MPI Communicator and build the ordering problem with that
1535 MPI_Comm myRawComm;
1536
1537 RCP<const MpiComm<int> > mpicomm =
1538 rcp_dynamic_cast<const MpiComm<int> > (ActiveMatrix->getComm ());
1539 if (mpicomm == Teuchos::null) {
1540 myRawComm = MPI_COMM_SELF;
1541 } else {
1542 myRawComm = * (mpicomm->getRawMpiComm ());
1543 }
1544 ordering_problem_type MyOrderingProblem (&Zoltan2Graph, &zlist, myRawComm);
1545#else
1546 ordering_problem_type MyOrderingProblem (&Zoltan2Graph, &zlist);
1547#endif
1548 MyOrderingProblem.solve ();
1549
1550 {
1551 typedef Zoltan2::LocalOrderingSolution<local_ordinal_type>
1552 ordering_solution_type;
1553
1554 ordering_solution_type sol (*MyOrderingProblem.getLocalOrderingSolution());
1555
1556 // perm[i] gives the where OLD index i shows up in the NEW
1557 // ordering. revperm[i] gives the where NEW index i shows
1558 // up in the OLD ordering. Note that perm is actually the
1559 // "inverse permutation," in Zoltan2 terms.
1560 perm = sol.getPermutationRCPConst (true);
1561 revperm = sol.getPermutationRCPConst ();
1562 }
1563 }
1564 // All reorderings here...
1565 ReorderedLocalizedMatrix_ = rcp (new reorder_filter_type (ActiveMatrix, perm, revperm));
1566
1567
1568 ActiveMatrix = ReorderedLocalizedMatrix_;
1569#else
1570 // This is a logic_error, not a runtime_error, because
1571 // setParameters() should have excluded this case already.
1572 TEUCHOS_TEST_FOR_EXCEPTION(
1573 true, std::logic_error, "Ifpack2::AdditiveSchwarz::setup: "
1574 "The Zoltan2 and Xpetra packages must be enabled in order "
1575 "to support reordering.");
1576#endif
1577 }
1578 innerMatrix_ = ActiveMatrix;
1579 }
1580
1581 TEUCHOS_TEST_FOR_EXCEPTION(
1582 innerMatrix_.is_null (), std::logic_error, "Ifpack2::AdditiveSchwarz::"
1583 "setup: Inner matrix is null right before constructing inner solver. "
1584 "Please report this bug to the Ifpack2 developers.");
1585
1586 // Construct the inner solver if necessary.
1587 if (Inverse_.is_null ()) {
1588 const std::string innerName = innerPrecName ();
1589 TEUCHOS_TEST_FOR_EXCEPTION(
1590 innerName == "INVALID", std::logic_error,
1591 "Ifpack2::AdditiveSchwarz::initialize: AdditiveSchwarz doesn't "
1592 "know how to create an instance of your LocalInverseType \""
1593 << Teuchos::TypeNameTraits<LocalInverseType>::name () << "\". "
1594 "Please talk to the Ifpack2 developers for details.");
1595
1596 TEUCHOS_TEST_FOR_EXCEPTION(
1597 innerName == "CUSTOM", std::runtime_error, "Ifpack2::AdditiveSchwarz::"
1598 "initialize: If the \"inner preconditioner name\" parameter (or any "
1599 "alias thereof) has the value \"CUSTOM\", then you must first call "
1600 "setInnerPreconditioner with a nonnull inner preconditioner input before "
1601 "you may call initialize().");
1602
1603 // FIXME (mfh 26 Aug 2015) Once we fix Bug 6392, the following
1604 // three lines of code can and SHOULD go away.
1605 if (! Trilinos::Details::Impl::registeredSomeLinearSolverFactory ("Ifpack2")) {
1607 }
1608
1609 // FIXME (mfh 26 Aug 2015) Provide the capability to get inner
1610 // solvers from packages other than Ifpack2.
1611 typedef typename MV::mag_type MT;
1612 RCP<inner_solver_type> innerPrec =
1613 Trilinos::Details::getLinearSolver<MV, OP, MT> ("Ifpack2", innerName);
1614 TEUCHOS_TEST_FOR_EXCEPTION(
1615 innerPrec.is_null (), std::logic_error,
1616 "Ifpack2::AdditiveSchwarz::setup: Failed to create inner preconditioner "
1617 "with name \"" << innerName << "\".");
1618 innerPrec->setMatrix (innerMatrix_);
1619
1620 // Extract and apply the sublist of parameters to give to the
1621 // inner solver, if there is such a sublist of parameters.
1622 std::pair<Teuchos::ParameterList, bool> result = innerPrecParams ();
1623 if (result.second) {
1624 // FIXME (mfh 26 Aug 2015) We don't really want to use yet
1625 // another deep copy of the ParameterList here.
1626 innerPrec->setParameters (rcp (new ParameterList (result.first)));
1627 }
1628 Inverse_ = innerPrec; // "Commit" the inner solver.
1629 }
1630 else if (Inverse_->getMatrix ().getRawPtr () != innerMatrix_.getRawPtr ()) {
1631 // The new inner matrix is different from the inner
1632 // preconditioner's current matrix, so give the inner
1633 // preconditioner the new inner matrix.
1634 Inverse_->setMatrix (innerMatrix_);
1635 }
1636 TEUCHOS_TEST_FOR_EXCEPTION(
1637 Inverse_.is_null (), std::logic_error, "Ifpack2::AdditiveSchwarz::"
1638 "setup: Inverse_ is null right after we were supposed to have created it."
1639 " Please report this bug to the Ifpack2 developers.");
1640
1641 // We don't have to call setInnerPreconditioner() here, because we
1642 // had the inner matrix (innerMatrix_) before creation of the inner
1643 // preconditioner. Calling setInnerPreconditioner here would be
1644 // legal, but it would require an unnecessary reset of the inner
1645 // preconditioner (i.e., calling initialize() and compute() again).
1646}
1647
1648
1649template<class MatrixType, class LocalInverseType>
1654 node_type> >& innerPrec)
1655{
1656 if (! innerPrec.is_null ()) {
1657 // Make sure that the new inner solver knows how to have its matrix changed.
1658 typedef Details::CanChangeMatrix<row_matrix_type> can_change_type;
1659 can_change_type* innerSolver = dynamic_cast<can_change_type*> (&*innerPrec);
1660 TEUCHOS_TEST_FOR_EXCEPTION(
1661 innerSolver == NULL, std::invalid_argument, "Ifpack2::AdditiveSchwarz::"
1662 "setInnerPreconditioner: The input preconditioner does not implement the "
1663 "setMatrix() feature. Only input preconditioners that inherit from "
1664 "Ifpack2::Details::CanChangeMatrix implement this feature.");
1665
1666 // If users provide an inner solver, we assume that
1667 // AdditiveSchwarz's current inner solver parameters no longer
1668 // apply. (In fact, we will remove those parameters from
1669 // AdditiveSchwarz's current list below.) Thus, we do /not/ apply
1670 // the current sublist of inner solver parameters to the input
1671 // inner solver.
1672
1673 // mfh 03 Jan 2014: Thanks to Paul Tsuji for pointing out that
1674 // it's perfectly legal for innerMatrix_ to be null here. This
1675 // can happen if initialize() has not been called yet. For
1676 // example, when Ifpack2::Factory creates an AdditiveSchwarz
1677 // instance, it calls setInnerPreconditioner() without first
1678 // calling initialize().
1679
1680 // Give the local matrix to the new inner solver.
1681 if(auto asf = Teuchos::rcp_dynamic_cast<Details::AdditiveSchwarzFilter<MatrixType>>(innerMatrix_))
1682 innerSolver->setMatrix (asf->getFilteredMatrix());
1683 else
1684 innerSolver->setMatrix (innerMatrix_);
1685
1686 // If the user previously specified a parameter for the inner
1687 // preconditioner's type, then clear out that parameter and its
1688 // associated sublist. Replace the inner preconditioner's type with
1689 // "CUSTOM", to make it obvious that AdditiveSchwarz's ParameterList
1690 // does not necessarily describe the current inner preconditioner.
1691 // We have to remove all allowed aliases of "inner preconditioner
1692 // name" before we may set it to "CUSTOM". Users may also set this
1693 // parameter to "CUSTOM" themselves, but this is not required.
1694 removeInnerPrecName ();
1695 removeInnerPrecParams ();
1696 List_.set ("inner preconditioner name", "CUSTOM");
1697
1698 // Bring the new inner solver's current status (initialized or
1699 // computed) in line with AdditiveSchwarz's current status.
1700 if (isInitialized ()) {
1701 innerPrec->initialize ();
1702 }
1703 if (isComputed ()) {
1704 innerPrec->compute ();
1705 }
1706 }
1707
1708 // If the new inner solver is null, we don't change the initialized
1709 // or computed status of AdditiveSchwarz. That way, AdditiveSchwarz
1710 // won't have to recompute innerMatrix_ if the inner solver changes.
1711 // This does introduce a new error condition in compute() and
1712 // apply(), but that's OK.
1713
1714 // Set the new inner solver.
1716 global_ordinal_type, node_type> inner_solver_impl_type;
1717 Inverse_ = Teuchos::rcp (new inner_solver_impl_type (innerPrec, "CUSTOM"));
1718}
1719
1720template<class MatrixType, class LocalInverseType>
1722setMatrix (const Teuchos::RCP<const row_matrix_type>& A)
1723{
1724 // Don't set the matrix unless it is different from the current one.
1725 if (A.getRawPtr () != Matrix_.getRawPtr ()) {
1726 IsInitialized_ = false;
1727 IsComputed_ = false;
1728
1729 // Reset all the state computed in initialize() and compute().
1730 OverlappingMatrix_ = Teuchos::null;
1731 ReorderedLocalizedMatrix_ = Teuchos::null;
1732 innerMatrix_ = Teuchos::null;
1733 SingletonMatrix_ = Teuchos::null;
1734 localMap_ = Teuchos::null;
1735 overlapping_B_.reset (nullptr);
1736 overlapping_Y_.reset (nullptr);
1737 R_.reset (nullptr);
1738 C_.reset (nullptr);
1739 DistributedImporter_ = Teuchos::null;
1740
1741 Matrix_ = A;
1742 }
1743}
1744
1745} // namespace Ifpack2
1746
1747// NOTE (mfh 26 Aug 2015) There's no need to instantiate for CrsMatrix
1748// too. All Ifpack2 preconditioners can and should do dynamic casts
1749// internally, if they need a type more specific than RowMatrix.
1750#define IFPACK2_ADDITIVESCHWARZ_INSTANT(S,LO,GO,N) \
1751 template class Ifpack2::AdditiveSchwarz< Tpetra::RowMatrix<S, LO, GO, N> >;
1752
1753#endif // IFPACK2_ADDITIVESCHWARZ_DECL_HPP
1754
void registerLinearSolverFactory()
Register Ifpack2's LinearSolverFactory with the central repository, for all enabled combinations of t...
Definition: Ifpack2_Details_registerLinearSolverFactory.cpp:68
Declaration of interface for preconditioners that can change their matrix after construction.
Additive Schwarz domain decomposition for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:296
typename MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:317
std::string description() const
Return a simple one-line description of this object.
Definition: Ifpack2_AdditiveSchwarz_def.hpp:1183
virtual bool isInitialized() const
Returns true if the preconditioner has been successfully initialized, false otherwise.
Definition: Ifpack2_AdditiveSchwarz_def.hpp:1060
virtual int getOverlapLevel() const
Returns the level of overlap.
Definition: Ifpack2_AdditiveSchwarz_def.hpp:1355
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Get a list of the preconditioner's default parameters.
Definition: Ifpack2_AdditiveSchwarz_def.hpp:937
typename MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input MatrixType.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:320
virtual int getNumCompute() const
Returns the number of calls to compute().
Definition: Ifpack2_AdditiveSchwarz_def.hpp:1148
virtual int getNumApply() const
Returns the number of calls to apply().
Definition: Ifpack2_AdditiveSchwarz_def.hpp:1155
virtual double getComputeTime() const
Returns the time spent in compute().
Definition: Ifpack2_AdditiveSchwarz_def.hpp:1169
virtual double getApplyTime() const
Returns the time spent in apply().
Definition: Ifpack2_AdditiveSchwarz_def.hpp:1176
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Change the matrix to be preconditioned.
Definition: Ifpack2_AdditiveSchwarz_def.hpp:1722
typename MatrixType::node_type node_type
The Node type used by the input MatrixType.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:323
virtual void initialize()
Computes all (graph-related) data necessary to initialize the preconditioner.
Definition: Ifpack2_AdditiveSchwarz_def.hpp:983
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &plist)
Set the preconditioner's parameters.
Definition: Ifpack2_AdditiveSchwarz_def.hpp:741
virtual void setInnerPreconditioner(const Teuchos::RCP< Preconditioner< scalar_type, local_ordinal_type, global_ordinal_type, node_type > > &innerPrec)
Set the inner preconditioner.
Definition: Ifpack2_AdditiveSchwarz_def.hpp:1651
virtual double getInitializeTime() const
Returns the time spent in initialize().
Definition: Ifpack2_AdditiveSchwarz_def.hpp:1162
virtual Teuchos::RCP< const Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > > getRangeMap() const
The range Map of this operator.
Definition: Ifpack2_AdditiveSchwarz_def.hpp:273
virtual Teuchos::RCP< const Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > > getDomainMap() const
The domain Map of this operator.
Definition: Ifpack2_AdditiveSchwarz_def.hpp:260
virtual void compute()
Computes all (coefficient) data necessary to apply the preconditioner.
Definition: Ifpack2_AdditiveSchwarz_def.hpp:1067
virtual std::ostream & print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
Definition: Ifpack2_AdditiveSchwarz_def.hpp:1345
typename MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:314
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object with some verbosity level to an FancyOStream object.
Definition: Ifpack2_AdditiveSchwarz_def.hpp:1238
virtual void setParameters(const Teuchos::ParameterList &plist)
Set the preconditioner's parameters.
Definition: Ifpack2_AdditiveSchwarz_def.hpp:728
virtual Teuchos::RCP< const row_matrix_type > getMatrix() const
The input matrix.
Definition: Ifpack2_AdditiveSchwarz_def.hpp:285
virtual void apply(const Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &X, Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, scalar_type alpha=Teuchos::ScalarTraits< scalar_type >::one(), scalar_type beta=Teuchos::ScalarTraits< scalar_type >::zero()) const
Apply the preconditioner to X, putting the result in Y.
Definition: Ifpack2_AdditiveSchwarz_def.hpp:294
virtual bool isComputed() const
Returns true if the preconditioner has been successfully computed, false otherwise.
Definition: Ifpack2_AdditiveSchwarz_def.hpp:1134
virtual int getNumInitialize() const
Returns the number of calls to initialize().
Definition: Ifpack2_AdditiveSchwarz_def.hpp:1141
AdditiveSchwarz(const Teuchos::RCP< const row_matrix_type > &A)
Constructor that takes a matrix.
Definition: Ifpack2_AdditiveSchwarz_def.hpp:245
Mix-in interface for preconditioners that can change their matrix after construction.
Definition: Ifpack2_Details_CanChangeMatrix.hpp:93
Ifpack2's implementation of Trilinos::Details::LinearSolver interface.
Definition: Ifpack2_Details_LinearSolver_decl.hpp:110
Sparse matrix (Tpetra::RowMatrix subclass) with ghost rows.
Definition: Ifpack2_OverlappingRowMatrix_decl.hpp:59
Interface for all Ifpack2 preconditioners.
Definition: Ifpack2_Preconditioner.hpp:108
void registerLinearSolverFactory()
Ifpack2 implementation details.
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:74
void getValidParameters(Teuchos::ParameterList &params)
Fills a list which contains all the parameters possibly used by Ifpack2.
Definition: Ifpack2_Parameters.cpp:51