Belos Version of the Day
Loading...
Searching...
No Matches
BelosPseudoBlockCGSolMgr.hpp
Go to the documentation of this file.
1//@HEADER
2// ************************************************************************
3//
4// Belos: Block Linear Solvers Package
5// Copyright 2004 Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
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#ifndef BELOS_PSEUDO_BLOCK_CG_SOLMGR_HPP
43#define BELOS_PSEUDO_BLOCK_CG_SOLMGR_HPP
44
49#include "BelosConfigDefs.hpp"
50#include "BelosTypes.hpp"
51
54
57#include "BelosCGIter.hpp"
63#include "Teuchos_LAPACK.hpp"
64#ifdef BELOS_TEUCHOS_TIME_MONITOR
65#include "Teuchos_TimeMonitor.hpp"
66#endif
67
84namespace Belos {
85
87
88
96 PseudoBlockCGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
97 {}};
98
99
100 // Partial specialization for unsupported ScalarType types.
101 // This contains a stub implementation.
102 template<class ScalarType, class MV, class OP,
103 const bool supportsScalarType =
106 public Details::SolverManagerRequiresLapack<ScalarType, MV, OP,
107 Belos::Details::LapackSupportsScalar<ScalarType>::value>
108 {
109 static const bool scalarTypeIsSupported =
111 typedef Details::SolverManagerRequiresLapack<ScalarType, MV, OP,
112 scalarTypeIsSupported> base_type;
113
114 public:
116 base_type ()
117 {}
119 const Teuchos::RCP<Teuchos::ParameterList> &pl) :
120 base_type ()
121 {}
123
124 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> >
125 getResidualStatusTest() const { return Teuchos::null; }
126 };
127
128
129 template<class ScalarType, class MV, class OP>
130 class PseudoBlockCGSolMgr<ScalarType, MV, OP, true> :
131 public Details::SolverManagerRequiresLapack<ScalarType, MV, OP, true>
132 {
133 private:
136 typedef Teuchos::ScalarTraits<ScalarType> SCT;
137 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
138 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
139
140 public:
141
143
144
151
167 PseudoBlockCGSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
168 const Teuchos::RCP<Teuchos::ParameterList> &pl );
169
172
174 Teuchos::RCP<SolverManager<ScalarType, MV, OP> > clone () const override {
175 return Teuchos::rcp(new PseudoBlockCGSolMgr<ScalarType,MV,OP>);
176 }
178
180
181
183 return *problem_;
184 }
185
188 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const override;
189
192 Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const override { return params_; }
193
199 Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
200 return Teuchos::tuple(timerSolve_);
201 }
202
203
214 MagnitudeType achievedTol() const override {
215 return achievedTol_;
216 }
217
219 int getNumIters() const override {
220 return numIters_;
221 }
222
226 bool isLOADetected() const override { return false; }
227
231 ScalarType getConditionEstimate() const {return condEstimate_;}
232 Teuchos::ArrayRCP<MagnitudeType> getEigenEstimates() const {return eigenEstimates_;}
233
235 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> >
236 getResidualStatusTest() const { return convTest_; }
237
239
241
242
244 void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) override { problem_ = problem; }
245
247 void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params ) override;
248
250
252
253
257 void reset( const ResetType type ) override { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
259
261
262
280 ReturnType solve() override;
281
283
286
288 std::string description() const override;
289
291 private:
292 // Compute the condition number estimate
293 void compute_condnum_tridiag_sym(Teuchos::ArrayView<MagnitudeType> diag,
294 Teuchos::ArrayView<MagnitudeType> offdiag,
295 Teuchos::ArrayRCP<MagnitudeType>& lambdas,
296 ScalarType & lambda_min,
297 ScalarType & lambda_max,
298 ScalarType & ConditionNumber );
299
300 // Linear problem.
301 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
302
303 // Output manager.
304 Teuchos::RCP<OutputManager<ScalarType> > printer_;
305 Teuchos::RCP<std::ostream> outputStream_;
306
307 // Status test.
308 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
309 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
310 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
311 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
312
313 // Current parameter list.
314 Teuchos::RCP<Teuchos::ParameterList> params_;
315
321 mutable Teuchos::RCP<const Teuchos::ParameterList> validParams_;
322
323 // Default solver values.
324 static constexpr int maxIters_default_ = 1000;
325 static constexpr bool assertPositiveDefiniteness_default_ = true;
326 static constexpr bool showMaxResNormOnly_default_ = false;
327 static constexpr int verbosity_default_ = Belos::Errors;
328 static constexpr int outputStyle_default_ = Belos::General;
329 static constexpr int outputFreq_default_ = -1;
330 static constexpr int defQuorum_default_ = 1;
331 static constexpr bool foldConvergenceDetectionIntoAllreduce_default_ = false;
332 static constexpr const char * resScale_default_ = "Norm of Initial Residual";
333 static constexpr const char * label_default_ = "Belos";
334 static constexpr bool genCondEst_default_ = false;
335
336 // Current solver values.
337 MagnitudeType convtol_,achievedTol_;
338 int maxIters_, numIters_;
339 int verbosity_, outputStyle_, outputFreq_, defQuorum_;
340 bool assertPositiveDefiniteness_, showMaxResNormOnly_;
341 bool foldConvergenceDetectionIntoAllreduce_;
342 std::string resScale_;
343 bool genCondEst_;
344 ScalarType condEstimate_;
345 Teuchos::ArrayRCP<MagnitudeType> eigenEstimates_;
346
347 // Timers.
348 std::string label_;
349 Teuchos::RCP<Teuchos::Time> timerSolve_;
350
351 // Internal state variables.
352 bool isSet_;
353 };
354
355
356// Empty Constructor
357template<class ScalarType, class MV, class OP>
359 outputStream_(Teuchos::rcpFromRef(std::cout)),
360 convtol_(DefaultSolverParameters::convTol),
361 maxIters_(maxIters_default_),
362 numIters_(0),
363 verbosity_(verbosity_default_),
364 outputStyle_(outputStyle_default_),
365 outputFreq_(outputFreq_default_),
366 defQuorum_(defQuorum_default_),
367 assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
368 showMaxResNormOnly_(showMaxResNormOnly_default_),
369 foldConvergenceDetectionIntoAllreduce_(foldConvergenceDetectionIntoAllreduce_default_),
370 resScale_(resScale_default_),
371 genCondEst_(genCondEst_default_),
372 condEstimate_(-Teuchos::ScalarTraits<ScalarType>::one()),
373 label_(label_default_),
374 isSet_(false)
375{}
376
377// Basic Constructor
378template<class ScalarType, class MV, class OP>
380PseudoBlockCGSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
381 const Teuchos::RCP<Teuchos::ParameterList> &pl ) :
382 problem_(problem),
383 outputStream_(Teuchos::rcpFromRef(std::cout)),
384 convtol_(DefaultSolverParameters::convTol),
385 maxIters_(maxIters_default_),
386 numIters_(0),
387 verbosity_(verbosity_default_),
388 outputStyle_(outputStyle_default_),
389 outputFreq_(outputFreq_default_),
390 defQuorum_(defQuorum_default_),
391 assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
392 showMaxResNormOnly_(showMaxResNormOnly_default_),
393 foldConvergenceDetectionIntoAllreduce_(foldConvergenceDetectionIntoAllreduce_default_),
394 resScale_(resScale_default_),
395 genCondEst_(genCondEst_default_),
396 condEstimate_(-Teuchos::ScalarTraits<ScalarType>::one()),
397 label_(label_default_),
398 isSet_(false)
399{
400 TEUCHOS_TEST_FOR_EXCEPTION(
401 problem_.is_null (), std::invalid_argument,
402 "Belos::PseudoBlockCGSolMgr two-argument constructor: "
403 "'problem' is null. You must supply a non-null Belos::LinearProblem "
404 "instance when calling this constructor.");
405
406 if (! pl.is_null ()) {
407 // Set the parameters using the list that was passed in.
408 setParameters (pl);
409 }
410}
411
412template<class ScalarType, class MV, class OP>
413void
415setParameters (const Teuchos::RCP<Teuchos::ParameterList>& params)
416{
417 using Teuchos::ParameterList;
418 using Teuchos::parameterList;
419 using Teuchos::RCP;
420 using Teuchos::rcp;
421
422 RCP<const ParameterList> defaultParams = this->getValidParameters ();
423
424 // Create the internal parameter list if one doesn't already exist.
425 // Belos' solvers treat the input ParameterList to setParameters as
426 // a "delta" -- that is, a change from the current state -- so the
427 // default parameter list (if the input is null) should be empty.
428 // This explains also why Belos' solvers copy parameters one by one
429 // from the input list to the current list.
430 //
431 // Belos obfuscates the latter, because it takes the input parameter
432 // list by RCP, rather than by (nonconst) reference. The latter
433 // would make more sense, given that it doesn't actually keep the
434 // input parameter list.
435 //
436 // Note, however, that Belos still correctly triggers the "used"
437 // field of each parameter in the input list. While isParameter()
438 // doesn't (apparently) trigger the "used" flag, get() certainly
439 // does.
440
441 if (params_.is_null ()) {
442 // Create an empty list with the same name as the default list.
443 params_ = parameterList (defaultParams->name ());
444 } else {
445 params->validateParameters (*defaultParams);
446 }
447
448 // Check for maximum number of iterations
449 if (params->isParameter ("Maximum Iterations")) {
450 maxIters_ = params->get ("Maximum Iterations", maxIters_default_);
451
452 // Update parameter in our list and in status test.
453 params_->set ("Maximum Iterations", maxIters_);
454 if (! maxIterTest_.is_null ()) {
455 maxIterTest_->setMaxIters (maxIters_);
456 }
457 }
458
459 // Check if positive definiteness assertions are to be performed
460 if (params->isParameter ("Assert Positive Definiteness")) {
461 assertPositiveDefiniteness_ =
462 params->get ("Assert Positive Definiteness",
463 assertPositiveDefiniteness_default_);
464
465 // Update parameter in our list.
466 params_->set ("Assert Positive Definiteness", assertPositiveDefiniteness_);
467 }
468
469 if (params->isParameter("Fold Convergence Detection Into Allreduce")) {
470 foldConvergenceDetectionIntoAllreduce_ = params->get("Fold Convergence Detection Into Allreduce",
471 foldConvergenceDetectionIntoAllreduce_default_);
472 }
473
474 // Check to see if the timer label changed.
475 if (params->isParameter ("Timer Label")) {
476 const std::string tempLabel = params->get ("Timer Label", label_default_);
477
478 // Update parameter in our list and solver timer
479 if (tempLabel != label_) {
480 label_ = tempLabel;
481 params_->set ("Timer Label", label_);
482 const std::string solveLabel =
483 label_ + ": PseudoBlockCGSolMgr total solve time";
484#ifdef BELOS_TEUCHOS_TIME_MONITOR
485 timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
486#endif
487 }
488 }
489
490 // Check for a change in verbosity level
491 if (params->isParameter ("Verbosity")) {
492 if (Teuchos::isParameterType<int> (*params, "Verbosity")) {
493 verbosity_ = params->get ("Verbosity", verbosity_default_);
494 } else {
495 verbosity_ = (int) Teuchos::getParameter<Belos::MsgType> (*params, "Verbosity");
496 }
497
498 // Update parameter in our list.
499 params_->set ("Verbosity", verbosity_);
500 if (! printer_.is_null ()) {
501 printer_->setVerbosity (verbosity_);
502 }
503 }
504
505 // Check for a change in output style
506 if (params->isParameter ("Output Style")) {
507 if (Teuchos::isParameterType<int> (*params, "Output Style")) {
508 outputStyle_ = params->get ("Output Style", outputStyle_default_);
509 } else {
510 // FIXME (mfh 29 Jul 2015) What if the type is wrong?
511 outputStyle_ = (int) Teuchos::getParameter<Belos::OutputType> (*params, "Output Style");
512 }
513
514 // Reconstruct the convergence test if the explicit residual test
515 // is not being used.
516 params_->set ("Output Style", outputStyle_);
517 outputTest_ = Teuchos::null;
518 }
519
520 // output stream
521 if (params->isParameter ("Output Stream")) {
522 outputStream_ = params->get<RCP<std::ostream> > ("Output Stream");
523
524 // Update parameter in our list.
525 params_->set ("Output Stream", outputStream_);
526 if (! printer_.is_null ()) {
527 printer_->setOStream (outputStream_);
528 }
529 }
530
531 // frequency level
532 if (verbosity_ & Belos::StatusTestDetails) {
533 if (params->isParameter ("Output Frequency")) {
534 outputFreq_ = params->get ("Output Frequency", outputFreq_default_);
535 }
536
537 // Update parameter in out list and output status test.
538 params_->set ("Output Frequency", outputFreq_);
539 if (! outputTest_.is_null ()) {
540 outputTest_->setOutputFrequency (outputFreq_);
541 }
542 }
543
544 // Condition estimate
545 if (params->isParameter ("Estimate Condition Number")) {
546 genCondEst_ = params->get ("Estimate Condition Number", genCondEst_default_);
547 }
548
549 // Create output manager if we need to.
550 if (printer_.is_null ()) {
551 printer_ = rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
552 }
553
554 // Convergence
555 typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
556 typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t;
557
558 // Check for convergence tolerance
559 if (params->isParameter ("Convergence Tolerance")) {
560 if (params->isType<MagnitudeType> ("Convergence Tolerance")) {
561 convtol_ = params->get ("Convergence Tolerance",
562 static_cast<MagnitudeType> (DefaultSolverParameters::convTol));
563 }
564 else {
565 convtol_ = params->get ("Convergence Tolerance", DefaultSolverParameters::convTol);
566 }
567
568 // Update parameter in our list and residual tests.
569 params_->set ("Convergence Tolerance", convtol_);
570 if (! convTest_.is_null ()) {
571 convTest_->setTolerance (convtol_);
572 }
573 }
574
575 if (params->isParameter ("Show Maximum Residual Norm Only")) {
576 showMaxResNormOnly_ = params->get<bool> ("Show Maximum Residual Norm Only");
577
578 // Update parameter in our list and residual tests
579 params_->set ("Show Maximum Residual Norm Only", showMaxResNormOnly_);
580 if (! convTest_.is_null ()) {
581 convTest_->setShowMaxResNormOnly (showMaxResNormOnly_);
582 }
583 }
584
585 // Check for a change in scaling, if so we need to build new residual tests.
586 bool newResTest = false;
587 {
588 // "Residual Scaling" is the old parameter name; "Implicit
589 // Residual Scaling" is the new name. We support both options for
590 // backwards compatibility.
591 std::string tempResScale = resScale_;
592 bool implicitResidualScalingName = false;
593 if (params->isParameter ("Residual Scaling")) {
594 tempResScale = params->get<std::string> ("Residual Scaling");
595 }
596 else if (params->isParameter ("Implicit Residual Scaling")) {
597 tempResScale = params->get<std::string> ("Implicit Residual Scaling");
598 implicitResidualScalingName = true;
599 }
600
601 // Only update the scaling if it's different.
602 if (resScale_ != tempResScale) {
603 const Belos::ScaleType resScaleType =
604 convertStringToScaleType (tempResScale);
605 resScale_ = tempResScale;
606
607 // Update parameter in our list and residual tests, using the
608 // given parameter name.
609 if (implicitResidualScalingName) {
610 params_->set ("Implicit Residual Scaling", resScale_);
611 }
612 else {
613 params_->set ("Residual Scaling", resScale_);
614 }
615
616 if (! convTest_.is_null ()) {
617 try {
618 convTest_->defineScaleForm (resScaleType, Belos::TwoNorm);
619 }
620 catch (std::exception& e) {
621 // Make sure the convergence test gets constructed again.
622 newResTest = true;
623 }
624 }
625 }
626 }
627
628 // Get the deflation quorum, or number of converged systems before deflation is allowed
629 if (params->isParameter ("Deflation Quorum")) {
630 defQuorum_ = params->get ("Deflation Quorum", defQuorum_);
631 params_->set ("Deflation Quorum", defQuorum_);
632 if (! convTest_.is_null ()) {
633 convTest_->setQuorum( defQuorum_ );
634 }
635 }
636
637 // Create status tests if we need to.
638
639 // Basic test checks maximum iterations and native residual.
640 if (maxIterTest_.is_null ()) {
641 maxIterTest_ = rcp (new StatusTestMaxIters<ScalarType,MV,OP> (maxIters_));
642 }
643
644 // Implicit residual test, using the native residual to determine if convergence was achieved.
645 if (convTest_.is_null () || newResTest) {
646 convTest_ = rcp (new StatusTestResNorm_t (convtol_, defQuorum_, showMaxResNormOnly_));
647 convTest_->defineScaleForm (convertStringToScaleType (resScale_), Belos::TwoNorm);
648 }
649
650 if (sTest_.is_null () || newResTest) {
651 sTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::OR, maxIterTest_, convTest_));
652 }
653
654 if (outputTest_.is_null () || newResTest) {
655 // Create the status test output class.
656 // This class manages and formats the output from the status test.
657 StatusTestOutputFactory<ScalarType,MV,OP> stoFactory (outputStyle_);
658 outputTest_ = stoFactory.create (printer_, sTest_, outputFreq_,
660
661 // Set the solver string for the output test
662 const std::string solverDesc = " Pseudo Block CG ";
663 outputTest_->setSolverDesc (solverDesc);
664 }
665
666 // Create the timer if we need to.
667 if (timerSolve_.is_null ()) {
668 const std::string solveLabel =
669 label_ + ": PseudoBlockCGSolMgr total solve time";
670#ifdef BELOS_TEUCHOS_TIME_MONITOR
671 timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
672#endif
673 }
674
675 // Inform the solver manager that the current parameters were set.
676 isSet_ = true;
677}
678
679
680template<class ScalarType, class MV, class OP>
681Teuchos::RCP<const Teuchos::ParameterList>
683{
684 using Teuchos::ParameterList;
685 using Teuchos::parameterList;
686 using Teuchos::RCP;
687
688 if (validParams_.is_null()) {
689 // Set all the valid parameters and their default values.
690 RCP<ParameterList> pl = parameterList ();
691 pl->set("Convergence Tolerance", static_cast<MagnitudeType>(DefaultSolverParameters::convTol),
692 "The relative residual tolerance that needs to be achieved by the\n"
693 "iterative solver in order for the linear system to be declared converged.");
694 pl->set("Maximum Iterations", static_cast<int>(maxIters_default_),
695 "The maximum number of block iterations allowed for each\n"
696 "set of RHS solved.");
697 pl->set("Assert Positive Definiteness", static_cast<bool>(assertPositiveDefiniteness_default_),
698 "Whether or not to assert that the linear operator\n"
699 "and the preconditioner are indeed positive definite.");
700 pl->set("Verbosity", static_cast<int>(verbosity_default_),
701 "What type(s) of solver information should be outputted\n"
702 "to the output stream.");
703 pl->set("Output Style", static_cast<int>(outputStyle_default_),
704 "What style is used for the solver information outputted\n"
705 "to the output stream.");
706 pl->set("Output Frequency", static_cast<int>(outputFreq_default_),
707 "How often convergence information should be outputted\n"
708 "to the output stream.");
709 pl->set("Deflation Quorum", static_cast<int>(defQuorum_default_),
710 "The number of linear systems that need to converge before\n"
711 "they are deflated. This number should be <= block size.");
712 pl->set("Output Stream", Teuchos::rcpFromRef(std::cout),
713 "A reference-counted pointer to the output stream where all\n"
714 "solver output is sent.");
715 pl->set("Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
716 "When convergence information is printed, only show the maximum\n"
717 "relative residual norm when the block size is greater than one.");
718 pl->set("Implicit Residual Scaling", resScale_default_,
719 "The type of scaling used in the residual convergence test.");
720 pl->set("Estimate Condition Number", static_cast<bool>(genCondEst_default_),
721 "Whether or not to estimate the condition number of the preconditioned system.");
722 // We leave the old name as a valid parameter for backwards
723 // compatibility (so that validateParametersAndSetDefaults()
724 // doesn't raise an exception if it encounters "Residual
725 // Scaling"). The new name was added for compatibility with other
726 // solvers, none of which use "Residual Scaling".
727 pl->set("Residual Scaling", resScale_default_,
728 "The type of scaling used in the residual convergence test. This "
729 "name is deprecated; the new name is \"Implicit Residual Scaling\".");
730 pl->set("Timer Label", static_cast<const char *>(label_default_),
731 "The string to use as a prefix for the timer labels.");
732 pl->set("Fold Convergence Detection Into Allreduce",static_cast<bool>(foldConvergenceDetectionIntoAllreduce_default_),
733 "Merge the allreduce for convergence detection with the one for CG.\n"
734 "This saves one all-reduce, but incurs more computation.");
735 validParams_ = pl;
736 }
737 return validParams_;
738}
739
740
741// solve()
742template<class ScalarType, class MV, class OP>
744{
745 const char prefix[] = "Belos::PseudoBlockCGSolMgr::solve: ";
746
747 // Set the current parameters if they were not set before.
748 // NOTE: This may occur if the user generated the solver manager with the default constructor and
749 // then didn't set any parameters using setParameters().
750 if (!isSet_) { setParameters( params_ ); }
751
752 TEUCHOS_TEST_FOR_EXCEPTION
753 (! problem_->isProblemSet (), PseudoBlockCGSolMgrLinearProblemFailure,
754 prefix << "The linear problem to solve is not ready. You must call "
755 "setProblem() on the Belos::LinearProblem instance before telling the "
756 "Belos solver to solve it.");
757
758 // Create indices for the linear systems to be solved.
759 int startPtr = 0;
760 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
761 int numCurrRHS = numRHS2Solve;
762
763 std::vector<int> currIdx( numRHS2Solve ), currIdx2( numRHS2Solve );
764 for (int i=0; i<numRHS2Solve; ++i) {
765 currIdx[i] = startPtr+i;
766 currIdx2[i]=i;
767 }
768
769 // Inform the linear problem of the current linear system to solve.
770 problem_->setLSIndex( currIdx );
771
773 // Parameter list (iteration)
774 Teuchos::ParameterList plist;
775
776 plist.set("Assert Positive Definiteness",assertPositiveDefiniteness_);
777 if(genCondEst_) plist.set("Max Size For Condest",maxIters_);
778
779 // Reset the status test.
780 outputTest_->reset();
781
782 // Assume convergence is achieved, then let any failed convergence set this to false.
783 bool isConverged = true;
784
786 // Pseudo-Block CG solver
787 Teuchos::RCP<CGIteration<ScalarType,MV,OP> > block_cg_iter;
788 if (numRHS2Solve == 1) {
789 plist.set("Fold Convergence Detection Into Allreduce",
790 foldConvergenceDetectionIntoAllreduce_);
791 block_cg_iter =
792 Teuchos::rcp (new CGIter<ScalarType,MV,OP> (problem_, printer_, outputTest_, convTest_, plist));
793 } else {
794 block_cg_iter =
795 Teuchos::rcp (new PseudoBlockCGIter<ScalarType,MV,OP> (problem_, printer_, outputTest_, plist));
796 }
797
798 // Setup condition estimate
799 block_cg_iter->setDoCondEst(genCondEst_);
800 bool condEstPerf = false;
801
802 // Enter solve() iterations
803 {
804#ifdef BELOS_TEUCHOS_TIME_MONITOR
805 Teuchos::TimeMonitor slvtimer(*timerSolve_);
806#endif
807
808 while ( numRHS2Solve > 0 ) {
809
810 // Reset the active / converged vectors from this block
811 std::vector<int> convRHSIdx;
812 std::vector<int> currRHSIdx( currIdx );
813 currRHSIdx.resize(numCurrRHS);
814
815 // Reset the number of iterations.
816 block_cg_iter->resetNumIters();
817
818 // Reset the number of calls that the status test output knows about.
819 outputTest_->resetNumCalls();
820
821 // Get the current residual for this block of linear systems.
822 Teuchos::RCP<MV> R_0 = MVT::CloneViewNonConst( *(Teuchos::rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
823
824 // Get a new state struct and initialize the solver.
826 newState.R = R_0;
827 block_cg_iter->initializeCG(newState);
828
829 while(1) {
830
831 // tell block_gmres_iter to iterate
832 try {
833
834 block_cg_iter->iterate();
835
837 //
838 // check convergence first
839 //
841 if ( convTest_->getStatus() == Passed ) {
842
843 // Figure out which linear systems converged.
844 std::vector<int> convIdx = Teuchos::rcp_dynamic_cast<StatusTestGenResNorm<ScalarType,MV,OP> >(convTest_)->convIndices();
845
846 // If the number of converged linear systems is equal to the
847 // number of current linear systems, then we are done with this block.
848 if (convIdx.size() == currRHSIdx.size())
849 break; // break from while(1){block_cg_iter->iterate()}
850
851 // Inform the linear problem that we are finished with this current linear system.
852 problem_->setCurrLS();
853
854 // Reset currRHSIdx to have the right-hand sides that are left to converge for this block.
855 int have = 0;
856 for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
857 bool found = false;
858 for (unsigned int j=0; j<convIdx.size(); ++j) {
859 if (currRHSIdx[i] == convIdx[j]) {
860 found = true;
861 break;
862 }
863 }
864 if (!found) {
865 currIdx2[have] = currIdx2[i];
866 currRHSIdx[have++] = currRHSIdx[i];
867 }
868 }
869 currRHSIdx.resize(have);
870 currIdx2.resize(have);
871
872 // Compute condition estimate if the very first linear system in the block has converged.
873 if (currRHSIdx[0] != 0 && genCondEst_ && !condEstPerf)
874 {
875 // Compute the estimate.
876 ScalarType l_min, l_max;
877 Teuchos::ArrayView<MagnitudeType> diag = block_cg_iter->getDiag();
878 Teuchos::ArrayView<MagnitudeType> offdiag = block_cg_iter->getOffDiag();
879 compute_condnum_tridiag_sym(diag,offdiag,eigenEstimates_,l_min,l_max,condEstimate_);
880
881 // Make sure not to do more condition estimate computations for this solve.
882 block_cg_iter->setDoCondEst(false);
883 condEstPerf = true;
884 }
885
886 // Set the remaining indices after deflation.
887 problem_->setLSIndex( currRHSIdx );
888
889 // Get the current residual vector.
890 std::vector<MagnitudeType> norms;
891 R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
892 for (int i=0; i<have; ++i) { currIdx2[i] = i; }
893
894 // Set the new state and initialize the solver.
896 defstate.R = R_0;
897 block_cg_iter->initializeCG(defstate);
898 }
899
901 //
902 // check for maximum iterations
903 //
905 else if ( maxIterTest_->getStatus() == Passed ) {
906 // we don't have convergence
907 isConverged = false;
908 break; // break from while(1){block_cg_iter->iterate()}
909 }
910
912 //
913 // we returned from iterate(), but none of our status tests Passed.
914 // something is wrong, and it is probably our fault.
915 //
917
918 else {
919 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
920 "Belos::PseudoBlockCGSolMgr::solve(): Invalid return from PseudoBlockCGIter::iterate().");
921 }
922 }
923 catch (const std::exception &e) {
924 printer_->stream(Errors) << "Error! Caught std::exception in PseudoBlockCGIter::iterate() at iteration "
925 << block_cg_iter->getNumIters() << std::endl
926 << e.what() << std::endl;
927 throw;
928 }
929 }
930
931 // Inform the linear problem that we are finished with this block linear system.
932 problem_->setCurrLS();
933
934 // Update indices for the linear systems to be solved.
935 startPtr += numCurrRHS;
936 numRHS2Solve -= numCurrRHS;
937
938 if ( numRHS2Solve > 0 ) {
939
940 numCurrRHS = numRHS2Solve;
941 currIdx.resize( numCurrRHS );
942 currIdx2.resize( numCurrRHS );
943 for (int i=0; i<numCurrRHS; ++i)
944 { currIdx[i] = startPtr+i; currIdx2[i] = i; }
945
946 // Set the next indices.
947 problem_->setLSIndex( currIdx );
948 }
949 else {
950 currIdx.resize( numRHS2Solve );
951 }
952
953 }// while ( numRHS2Solve > 0 )
954
955 }
956
957 // print final summary
958 sTest_->print( printer_->stream(FinalSummary) );
959
960 // print timing information
961#ifdef BELOS_TEUCHOS_TIME_MONITOR
962 // Calling summarize() can be expensive, so don't call unless the
963 // user wants to print out timing details. summarize() will do all
964 // the work even if it's passed a "black hole" output stream.
965 if (verbosity_ & TimingDetails)
966 Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
967#endif
968
969 // get iteration information for this solve
970 numIters_ = maxIterTest_->getNumIters();
971
972 // Save the convergence test value ("achieved tolerance") for this
973 // solve.
974 const std::vector<MagnitudeType>* pTestValues = convTest_->getTestValue();
975 if (pTestValues != NULL && pTestValues->size () > 0) {
976 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
977 }
978
979 // Do condition estimate, if needed
980 if (genCondEst_ && !condEstPerf) {
981 ScalarType l_min, l_max;
982 Teuchos::ArrayView<MagnitudeType> diag = block_cg_iter->getDiag();
983 Teuchos::ArrayView<MagnitudeType> offdiag = block_cg_iter->getOffDiag();
984 compute_condnum_tridiag_sym(diag,offdiag,eigenEstimates_,l_min,l_max,condEstimate_);
985 condEstPerf = true;
986 }
987
988 if (! isConverged) {
989 return Unconverged; // return from PseudoBlockCGSolMgr::solve()
990 }
991 return Converged; // return from PseudoBlockCGSolMgr::solve()
992}
993
994// This method requires the solver manager to return a std::string that describes itself.
995template<class ScalarType, class MV, class OP>
997{
998 std::ostringstream oss;
999 oss << "Belos::PseudoBlockCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
1000 oss << "{";
1001 oss << "}";
1002 return oss.str();
1003}
1004
1005
1006template<class ScalarType, class MV, class OP>
1007void
1009compute_condnum_tridiag_sym (Teuchos::ArrayView<MagnitudeType> diag,
1010 Teuchos::ArrayView<MagnitudeType> offdiag,
1011 Teuchos::ArrayRCP<MagnitudeType>& lambdas,
1012 ScalarType & lambda_min,
1013 ScalarType & lambda_max,
1014 ScalarType & ConditionNumber )
1015{
1016 typedef Teuchos::ScalarTraits<ScalarType> STS;
1017
1018 /* Copied from az_cg.c: compute_condnum_tridiag_sym */
1019 /* diag == ScalarType vector of size N, containing the diagonal
1020 elements of A
1021 offdiag == ScalarType vector of size N-1, containing the offdiagonal
1022 elements of A. Note that A is supposed to be symmatric
1023 */
1024 int info = 0;
1025 const int N = diag.size ();
1026 ScalarType scalar_dummy;
1027 std::vector<MagnitudeType> mag_dummy(4*N);
1028 char char_N = 'N';
1029 Teuchos::LAPACK<int,ScalarType> lapack;
1030
1031 lambdas.resize(N, 0.0);
1032 lambda_min = STS::one ();
1033 lambda_max = STS::one ();
1034 if( N > 2 ) {
1035 lapack.PTEQR (char_N, N, diag.getRawPtr (), offdiag.getRawPtr (),
1036 &scalar_dummy, 1, &mag_dummy[0], &info);
1037 TEUCHOS_TEST_FOR_EXCEPTION
1038 (info < 0, std::logic_error, "Belos::PseudoBlockCGSolMgr::"
1039 "compute_condnum_tridiag_sym: LAPACK's _PTEQR failed with info = "
1040 << info << " < 0. This suggests there might be a bug in the way Belos "
1041 "is calling LAPACK. Please report this to the Belos developers.");
1042 for (int k = 0; k < N; k++) {
1043 lambdas[k] = diag[N - 1 - k];
1044 }
1045 lambda_min = Teuchos::as<ScalarType> (diag[N-1]);
1046 lambda_max = Teuchos::as<ScalarType> (diag[0]);
1047 }
1048
1049 // info > 0 means that LAPACK's eigensolver didn't converge. This
1050 // is unlikely but may be possible. In that case, the best we can
1051 // do is use the eigenvalues that it computes, as long as lambda_max
1052 // >= lambda_min.
1053 if (STS::real (lambda_max) < STS::real (lambda_min)) {
1054 ConditionNumber = STS::one ();
1055 }
1056 else {
1057 // It's OK for the condition number to be Inf.
1058 ConditionNumber = lambda_max / lambda_min;
1059 }
1060
1061} /* compute_condnum_tridiag_sym */
1062
1063
1064
1065
1066
1067} // end Belos namespace
1068
1069#endif /* BELOS_PSEUDO_BLOCK_CG_SOLMGR_HPP */
Belos concrete class for performing the conjugate-gradient (CG) iteration.
Belos concrete class for performing a single-reduction conjugate-gradient (CG) iteration.
Belos header file which uses auto-configuration information to include necessary C++ headers.
Class which describes the linear problem to be solved by the iterative solver.
Class which manages the output and verbosity of the Belos solvers.
Belos concrete class for performing the pseudo-block CG iteration.
Pure virtual base class which describes the basic interface for a solver manager.
Belos::StatusTest for logically combining several status tests.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Belos::StatusTest class for specifying a maximum number of iterations.
A factory class for generating StatusTestOutput objects.
Collection of types and exceptions used within the Belos solvers.
Parent class to all Belos exceptions.
This class implements the preconditioned Conjugate Gradient (CG) iteration.
Type traits class that says whether Teuchos::LAPACK has a valid implementation for the given ScalarTy...
Base class for Belos::SolverManager subclasses which normally can only compile with ScalarType types ...
A linear system to solve, and its associated information.
Traits class which defines basic operations on multivectors.
Class which defines basic traits for the operator type.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
This class implements the pseudo-block CG iteration, where the basic CG algorithm is performed on all...
PseudoBlockCGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i....
PseudoBlockCGSolMgrLinearProblemFailure(const std::string &what_arg)
Teuchos::ArrayRCP< MagnitudeType > getEigenEstimates() const
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
void reset(const ResetType type) override
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Return a reference to the linear problem being solved by this solver manager.
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > getResidualStatusTest() const
Return the residual status test.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
ScalarType getConditionEstimate() const
Gets the estimated condition number.
int getNumIters() const override
Get the iteration count for the most recent call to solve().
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
bool isLOADetected() const override
Return whether a loss of accuracy was detected by this solver during the most current solve.
The Belos::PseudoBlockCGSolMgr provides a powerful and fully-featured solver manager over the pseudo-...
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > getResidualStatusTest() const
PseudoBlockCGSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
A class for extending the status testing capabilities of Belos via logical combinations.
An implementation of StatusTestResNorm using a family of residual norms.
int setTolerance(MagnitudeType tolerance)
Set the value of the tolerance.
A Belos::StatusTest class for specifying a maximum number of iterations.
A factory class for generating StatusTestOutput objects.
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > create(const Teuchos::RCP< OutputManager< ScalarType > > &printer, Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test, int mod, int printStates)
Create the StatusTestOutput object specified by the outputStyle.
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
@ StatusTestDetails
@ FinalSummary
@ TimingDetails
ReturnType
Whether the Belos solve converged for all linear systems.
@ Unconverged
ScaleType
The type of scaling to use on the residual norm value.
ResetType
How to reset the solver.
Structure to contain pointers to CGIteration state variables.
Teuchos::RCP< const MV > R
The current residual.
Default parameters common to most Belos solvers.
static const double convTol
Default convergence tolerance.

Generated for Belos by doxygen 1.10.0