Belos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
BelosBlockCGSolMgr.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_BLOCK_CG_SOLMGR_HPP
43#define BELOS_BLOCK_CG_SOLMGR_HPP
44
49#include "BelosConfigDefs.hpp"
50#include "BelosTypes.hpp"
51
54
55#include "BelosCGIter.hpp"
57#include "BelosBlockCGIter.hpp"
64#include "Teuchos_LAPACK.hpp"
65#ifdef BELOS_TEUCHOS_TIME_MONITOR
66# include "Teuchos_TimeMonitor.hpp"
67#endif
68#if defined(HAVE_TEUCHOSCORE_CXX11)
69# include <type_traits>
70#endif // defined(HAVE_TEUCHOSCORE_CXX11)
71#include <algorithm>
72
89namespace Belos {
90
92
93
101 BlockCGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
102 {}};
103
104 template<class ScalarType, class MV, class OP,
105 const bool lapackSupportsScalarType =
108 public Details::SolverManagerRequiresLapack<ScalarType,MV,OP>
109 {
110 static const bool requiresLapack =
112 typedef Details::SolverManagerRequiresLapack<ScalarType, MV, OP,
114
115 public:
117 base_type ()
118 {}
119 BlockCGSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >& problem,
120 const Teuchos::RCP<Teuchos::ParameterList>& pl) :
121 base_type ()
122 {}
123 virtual ~BlockCGSolMgr () {}
124 };
125
126
127 // Partial specialization for ScalarType types for which
128 // Teuchos::LAPACK has a valid implementation. This contains the
129 // actual working implementation of BlockCGSolMgr.
130 template<class ScalarType, class MV, class OP>
131 class BlockCGSolMgr<ScalarType, MV, OP, true> :
132 public Details::SolverManagerRequiresLapack<ScalarType, MV, OP, true>
133 {
134 // This partial specialization is already chosen for those scalar types
135 // that support lapack, so we don't need to have an additional compile-time
136 // check that the scalar type is float/double/complex.
137// #if defined(HAVE_TEUCHOSCORE_CXX11)
138// # if defined(HAVE_TEUCHOS_COMPLEX)
139// static_assert (std::is_same<ScalarType, std::complex<float> >::value ||
140// std::is_same<ScalarType, std::complex<double> >::value ||
141// std::is_same<ScalarType, float>::value ||
142// std::is_same<ScalarType, double>::value,
143// "Belos::GCRODRSolMgr: ScalarType must be one of the four "
144// "types (S,D,C,Z) supported by LAPACK.");
145// # else
146// static_assert (std::is_same<ScalarType, float>::value ||
147// std::is_same<ScalarType, double>::value,
148// "Belos::GCRODRSolMgr: ScalarType must be float or double. "
149// "Complex arithmetic support is currently disabled. To "
150// "enable it, set Teuchos_ENABLE_COMPLEX=ON.");
151// # endif // defined(HAVE_TEUCHOS_COMPLEX)
152// #endif // defined(HAVE_TEUCHOSCORE_CXX11)
153
154 private:
157 typedef Teuchos::ScalarTraits<ScalarType> SCT;
158 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
159 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
160
161 public:
162
164
165
172
209 BlockCGSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
210 const Teuchos::RCP<Teuchos::ParameterList> &pl );
211
213 virtual ~BlockCGSolMgr() {};
214
216 Teuchos::RCP<SolverManager<ScalarType, MV, OP> > clone () const override {
217 return Teuchos::rcp(new BlockCGSolMgr<ScalarType,MV,OP>);
218 }
220
222
223
225 return *problem_;
226 }
227
230 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const override;
231
234 Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const override { return params_; }
235
241 Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
242 return Teuchos::tuple(timerSolve_);
243 }
244
250 MagnitudeType achievedTol() const override {
251 return achievedTol_;
252 }
253
255 int getNumIters() const override {
256 return numIters_;
257 }
258
261 bool isLOADetected() const override { return false; }
262
264
266
267
269 void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) override { problem_ = problem; }
270
272 void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params ) override;
273
275
277
278
282 void reset( const ResetType type ) override { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
284
286
287
305 ReturnType solve() override;
306
308
311
313 std::string description() const override;
314
316
317 private:
318
320 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
321
323 Teuchos::RCP<OutputManager<ScalarType> > printer_;
325 Teuchos::RCP<std::ostream> outputStream_;
326
331 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
332
334 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
335
337 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
338
340 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
341
343 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
344
346 Teuchos::RCP<Teuchos::ParameterList> params_;
347
348 //
349 // Default solver parameters.
350 //
351 static constexpr int maxIters_default_ = 1000;
352 static constexpr bool adaptiveBlockSize_default_ = true;
353 static constexpr bool showMaxResNormOnly_default_ = false;
354 static constexpr bool useSingleReduction_default_ = false;
355 static constexpr int blockSize_default_ = 1;
356 static constexpr int verbosity_default_ = Belos::Errors;
357 static constexpr int outputStyle_default_ = Belos::General;
358 static constexpr int outputFreq_default_ = -1;
359 static constexpr const char * resNorm_default_ = "TwoNorm";
360 static constexpr bool foldConvergenceDetectionIntoAllreduce_default_ = false;
361 static constexpr const char * resScale_default_ = "Norm of Initial Residual";
362 static constexpr const char * label_default_ = "Belos";
363 static constexpr const char * orthoType_default_ = "ICGS";
364 static constexpr bool assertPositiveDefiniteness_default_ = true;
365
366 //
367 // Current solver parameters and other values.
368 //
369
372
375
382
385
388
390 int blockSize_, verbosity_, outputStyle_, outputFreq_;
391 bool adaptiveBlockSize_, showMaxResNormOnly_, useSingleReduction_;
392 std::string orthoType_, resScale_;
395
397 std::string label_;
398
400 Teuchos::RCP<Teuchos::Time> timerSolve_;
401
403 bool isSet_;
404 };
405
406
407// Empty Constructor
408template<class ScalarType, class MV, class OP>
410 outputStream_(Teuchos::rcpFromRef(std::cout)),
411 convtol_(DefaultSolverParameters::convTol),
412 orthoKappa_(DefaultSolverParameters::orthoKappa),
413 achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
414 maxIters_(maxIters_default_),
415 numIters_(0),
416 blockSize_(blockSize_default_),
417 verbosity_(verbosity_default_),
418 outputStyle_(outputStyle_default_),
419 outputFreq_(outputFreq_default_),
420 adaptiveBlockSize_(adaptiveBlockSize_default_),
421 showMaxResNormOnly_(showMaxResNormOnly_default_),
422 useSingleReduction_(useSingleReduction_default_),
423 orthoType_(orthoType_default_),
424 resScale_(resScale_default_),
425 assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
426 foldConvergenceDetectionIntoAllreduce_(foldConvergenceDetectionIntoAllreduce_default_),
427 label_(label_default_),
428 isSet_(false)
429{}
430
431
432// Basic Constructor
433template<class ScalarType, class MV, class OP>
435BlockCGSolMgr(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
436 const Teuchos::RCP<Teuchos::ParameterList> &pl) :
437 problem_(problem),
438 outputStream_(Teuchos::rcpFromRef(std::cout)),
439 convtol_(DefaultSolverParameters::convTol),
440 orthoKappa_(DefaultSolverParameters::orthoKappa),
441 achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
442 maxIters_(maxIters_default_),
443 numIters_(0),
444 blockSize_(blockSize_default_),
445 verbosity_(verbosity_default_),
446 outputStyle_(outputStyle_default_),
447 outputFreq_(outputFreq_default_),
448 adaptiveBlockSize_(adaptiveBlockSize_default_),
449 showMaxResNormOnly_(showMaxResNormOnly_default_),
450 useSingleReduction_(useSingleReduction_default_),
451 orthoType_(orthoType_default_),
452 resScale_(resScale_default_),
453 assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
454 foldConvergenceDetectionIntoAllreduce_(foldConvergenceDetectionIntoAllreduce_default_),
455 label_(label_default_),
456 isSet_(false)
457{
458 TEUCHOS_TEST_FOR_EXCEPTION(problem_.is_null(), std::invalid_argument,
459 "BlockCGSolMgr's constructor requires a nonnull LinearProblem instance.");
460
461 // If the user passed in a nonnull parameter list, set parameters.
462 // Otherwise, the next solve() call will use default parameters,
463 // unless the user calls setParameters() first.
464 if (! pl.is_null()) {
465 setParameters (pl);
466 }
467}
468
469template<class ScalarType, class MV, class OP>
470void
472setParameters (const Teuchos::RCP<Teuchos::ParameterList> &params)
473{
474 // Create the internal parameter list if one doesn't already exist.
475 if (params_ == Teuchos::null) {
476 params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
477 }
478 else {
479 params->validateParameters(*getValidParameters());
480 }
481
482 // Check for maximum number of iterations
483 if (params->isParameter("Maximum Iterations")) {
484 maxIters_ = params->get("Maximum Iterations",maxIters_default_);
485
486 // Update parameter in our list and in status test.
487 params_->set("Maximum Iterations", maxIters_);
488 if (maxIterTest_!=Teuchos::null)
489 maxIterTest_->setMaxIters( maxIters_ );
490 }
491
492 // Check for blocksize
493 if (params->isParameter("Block Size")) {
494 blockSize_ = params->get("Block Size",blockSize_default_);
495 TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
496 "Belos::BlockCGSolMgr: \"Block Size\" must be strictly positive.");
497
498 // Update parameter in our list.
499 params_->set("Block Size", blockSize_);
500 }
501
502 // Check if the blocksize should be adaptive
503 if (params->isParameter("Adaptive Block Size")) {
504 adaptiveBlockSize_ = params->get("Adaptive Block Size",adaptiveBlockSize_default_);
505
506 // Update parameter in our list.
507 params_->set("Adaptive Block Size", adaptiveBlockSize_);
508 }
509
510 // Check if the user is requesting the single-reduction version of CG (only for blocksize == 1)
511 if (params->isParameter("Use Single Reduction")) {
512 useSingleReduction_ = params->get("Use Single Reduction", useSingleReduction_default_);
513 }
514
515 if (params->isParameter("Fold Convergence Detection Into Allreduce")) {
516 foldConvergenceDetectionIntoAllreduce_ = params->get("Fold Convergence Detection Into Allreduce",
517 foldConvergenceDetectionIntoAllreduce_default_);
518 }
519
520 // Check to see if the timer label changed.
521 if (params->isParameter("Timer Label")) {
522 std::string tempLabel = params->get("Timer Label", label_default_);
523
524 // Update parameter in our list and solver timer
525 if (tempLabel != label_) {
526 label_ = tempLabel;
527 params_->set("Timer Label", label_);
528 std::string solveLabel = label_ + ": BlockCGSolMgr total solve time";
529#ifdef BELOS_TEUCHOS_TIME_MONITOR
530 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
531#endif
532 if (ortho_ != Teuchos::null) {
533 ortho_->setLabel( label_ );
534 }
535 }
536 }
537
538 // Check for a change in verbosity level
539 if (params->isParameter("Verbosity")) {
540 if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
541 verbosity_ = params->get("Verbosity", verbosity_default_);
542 } else {
543 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
544 }
545
546 // Update parameter in our list.
547 params_->set("Verbosity", verbosity_);
548 if (printer_ != Teuchos::null)
549 printer_->setVerbosity(verbosity_);
550 }
551
552 // Check for a change in output style
553 if (params->isParameter("Output Style")) {
554 if (Teuchos::isParameterType<int>(*params,"Output Style")) {
555 outputStyle_ = params->get("Output Style", outputStyle_default_);
556 } else {
557 outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
558 }
559
560 // Update parameter in our list.
561 params_->set("Output Style", outputStyle_);
562 outputTest_ = Teuchos::null;
563 }
564
565 // output stream
566 if (params->isParameter("Output Stream")) {
567 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
568
569 // Update parameter in our list.
570 params_->set("Output Stream", outputStream_);
571 if (printer_ != Teuchos::null)
572 printer_->setOStream( outputStream_ );
573 }
574
575 // frequency level
576 if (verbosity_ & Belos::StatusTestDetails) {
577 if (params->isParameter("Output Frequency")) {
578 outputFreq_ = params->get("Output Frequency", outputFreq_default_);
579 }
580
581 // Update parameter in out list and output status test.
582 params_->set("Output Frequency", outputFreq_);
583 if (outputTest_ != Teuchos::null)
584 outputTest_->setOutputFrequency( outputFreq_ );
585 }
586
587 // Create output manager if we need to.
588 if (printer_ == Teuchos::null) {
589 printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
590 }
591
592 // Check if the orthogonalization changed.
593 bool changedOrthoType = false;
594 if (params->isParameter("Orthogonalization")) {
595 std::string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
596 if (tempOrthoType != orthoType_) {
597 orthoType_ = tempOrthoType;
598 changedOrthoType = true;
599 }
600 }
601 params_->set("Orthogonalization", orthoType_);
602
603 // Check which orthogonalization constant to use.
604 if (params->isParameter("Orthogonalization Constant")) {
605 if (params->isType<MagnitudeType> ("Orthogonalization Constant")) {
606 orthoKappa_ = params->get ("Orthogonalization Constant",
608 }
609 else {
610 orthoKappa_ = params->get ("Orthogonalization Constant",
612 }
613
614 // Update parameter in our list.
615 params_->set("Orthogonalization Constant",orthoKappa_);
616 if (orthoType_=="DGKS") {
617 if (orthoKappa_ > 0 && ortho_ != Teuchos::null && !changedOrthoType) {
618 Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
619 }
620 }
621 }
622
623 // Create orthogonalization manager if we need to.
624 if (ortho_ == Teuchos::null || changedOrthoType) {
626 Teuchos::RCP<Teuchos::ParameterList> paramsOrtho; // can be null
627 if (orthoType_=="DGKS" && orthoKappa_ > 0) {
628 paramsOrtho->set ("depTol", orthoKappa_ );
629 }
630
631 ortho_ = factory.makeMatOrthoManager (orthoType_, Teuchos::null, printer_, label_, paramsOrtho);
632 }
633
634 // Convergence
635 typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
636 typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t;
637
638 // Check for convergence tolerance
639 if (params->isParameter("Convergence Tolerance")) {
640 if (params->isType<MagnitudeType> ("Convergence Tolerance")) {
641 convtol_ = params->get ("Convergence Tolerance",
643 }
644 else {
645 convtol_ = params->get ("Convergence Tolerance", DefaultSolverParameters::convTol);
646 }
647
648 // Update parameter in our list and residual tests.
649 params_->set("Convergence Tolerance", convtol_);
650 if (convTest_ != Teuchos::null)
651 convTest_->setTolerance( convtol_ );
652 }
653
654 if (params->isParameter("Show Maximum Residual Norm Only")) {
655 showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
656
657 // Update parameter in our list and residual tests
658 params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
659 if (convTest_ != Teuchos::null)
660 convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
661 }
662
663 // Check for a change in scaling, if so we need to build new residual tests.
664 bool newResTest = false;
665 {
666 std::string tempResScale = resScale_;
667 if (params->isParameter ("Implicit Residual Scaling")) {
668 tempResScale = params->get<std::string> ("Implicit Residual Scaling");
669 }
670
671 // Only update the scaling if it's different.
672 if (resScale_ != tempResScale) {
673 const Belos::ScaleType resScaleType =
674 convertStringToScaleType (tempResScale);
675 resScale_ = tempResScale;
676
677 // Update parameter in our list and residual tests
678 params_->set ("Implicit Residual Scaling", resScale_);
679
680 if (! convTest_.is_null ()) {
681 try {
682 NormType normType = Belos::TwoNorm;
683 if (params->isParameter("Residual Norm")) {
684 if (params->isType<std::string> ("Residual Norm")) {
685 normType = convertStringToNormType(params->get<std::string> ("Residual Norm"));
686 }
687 }
688 convTest_->defineResForm(StatusTestResNorm_t::Implicit, normType);
689 convTest_->defineScaleForm (resScaleType, Belos::TwoNorm);
690 }
691 catch (std::exception& e) {
692 // Make sure the convergence test gets constructed again.
693 newResTest = true;
694 }
695 }
696 }
697 }
698
699 // Create status tests if we need to.
700
701 // Basic test checks maximum iterations and native residual.
702 if (maxIterTest_ == Teuchos::null)
703 maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
704
705 // Implicit residual test, using the native residual to determine if convergence was achieved.
706 if (convTest_.is_null () || newResTest) {
707
708 NormType normType = Belos::TwoNorm;
709 if (params->isParameter("Residual Norm")) {
710 if (params->isType<std::string> ("Residual Norm")) {
711 normType = convertStringToNormType(params->get<std::string> ("Residual Norm"));
712 }
713 }
714
715 convTest_ = rcp (new StatusTestResNorm_t (convtol_, 1, showMaxResNormOnly_));
716 convTest_->defineResForm(StatusTestResNorm_t::Implicit, normType);
717 convTest_->defineScaleForm (convertStringToScaleType (resScale_), Belos::TwoNorm);
718 }
719
720 if (sTest_.is_null () || newResTest)
721 sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
722
723 if (outputTest_.is_null () || newResTest) {
724
725 // Create the status test output class.
726 // This class manages and formats the output from the status test.
727 StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
728 outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
729
730 // Set the solver string for the output test
731 std::string solverDesc = " Block CG ";
732 outputTest_->setSolverDesc( solverDesc );
733
734 }
735
736 // BelosCgIter accepts a parameter specifying whether to assert for the positivity of p^H*A*p in the CG iteration
737 if (params->isParameter("Assert Positive Definiteness")) {
738 assertPositiveDefiniteness_ = Teuchos::getParameter<bool>(*params,"Assert Positive Definiteness");
739 params_->set("Assert Positive Definiteness", assertPositiveDefiniteness_);
740 }
741
742 // Create the timer if we need to.
743 if (timerSolve_ == Teuchos::null) {
744 std::string solveLabel = label_ + ": BlockCGSolMgr total solve time";
745#ifdef BELOS_TEUCHOS_TIME_MONITOR
746 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
747#endif
748 }
749
750 // Inform the solver manager that the current parameters were set.
751 isSet_ = true;
752}
753
754
755template<class ScalarType, class MV, class OP>
756Teuchos::RCP<const Teuchos::ParameterList>
758{
759 static Teuchos::RCP<const Teuchos::ParameterList> validPL;
760
761 // Set all the valid parameters and their default values.
762 if(is_null(validPL)) {
763 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
764 pl->set("Convergence Tolerance", static_cast<MagnitudeType>(DefaultSolverParameters::convTol),
765 "The relative residual tolerance that needs to be achieved by the\n"
766 "iterative solver in order for the linear system to be declared converged.");
767 pl->set("Maximum Iterations", static_cast<int>(maxIters_default_),
768 "The maximum number of block iterations allowed for each\n"
769 "set of RHS solved.");
770 pl->set("Block Size", static_cast<int>(blockSize_default_),
771 "The number of vectors in each block.");
772 pl->set("Adaptive Block Size", static_cast<bool>(adaptiveBlockSize_default_),
773 "Whether the solver manager should adapt to the block size\n"
774 "based on the number of RHS to solve.");
775 pl->set("Verbosity", static_cast<int>(verbosity_default_),
776 "What type(s) of solver information should be outputted\n"
777 "to the output stream.");
778 pl->set("Output Style", static_cast<int>(outputStyle_default_),
779 "What style is used for the solver information outputted\n"
780 "to the output stream.");
781 pl->set("Output Frequency", static_cast<int>(outputFreq_default_),
782 "How often convergence information should be outputted\n"
783 "to the output stream.");
784 pl->set("Output Stream", Teuchos::rcpFromRef(std::cout),
785 "A reference-counted pointer to the output stream where all\n"
786 "solver output is sent.");
787 pl->set("Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
788 "When convergence information is printed, only show the maximum\n"
789 "relative residual norm when the block size is greater than one.");
790 pl->set("Use Single Reduction", static_cast<bool>(useSingleReduction_default_),
791 "Use single reduction iteration when the block size is one.");
792 pl->set("Implicit Residual Scaling", resScale_default_,
793 "The type of scaling used in the residual convergence test.");
794 pl->set("Timer Label", static_cast<const char *>(label_default_),
795 "The string to use as a prefix for the timer labels.");
796 pl->set("Orthogonalization", static_cast<const char *>(orthoType_default_),
797 "The type of orthogonalization to use: DGKS, ICGS, or IMGS.");
798 pl->set("Assert Positive Definiteness",static_cast<bool>(assertPositiveDefiniteness_default_),
799 "Assert for positivity of p^H*A*p in CG iteration.");
800 pl->set("Orthogonalization Constant",static_cast<MagnitudeType>(DefaultSolverParameters::orthoKappa),
801 "The constant used by DGKS orthogonalization to determine\n"
802 "whether another step of classical Gram-Schmidt is necessary.");
803 pl->set("Residual Norm",static_cast<const char *>(resNorm_default_),
804 "Norm used for the convergence check on the residual.");
805 pl->set("Fold Convergence Detection Into Allreduce",static_cast<bool>(foldConvergenceDetectionIntoAllreduce_default_),
806 "Merge the allreduce for convergence detection with the one for CG.\n"
807 "This saves one all-reduce, but incurs more computation.");
808 validPL = pl;
809 }
810 return validPL;
811}
812
813
814// solve()
815template<class ScalarType, class MV, class OP>
817 using Teuchos::RCP;
818 using Teuchos::rcp;
819 using Teuchos::rcp_const_cast;
820 using Teuchos::rcp_dynamic_cast;
821
822 // Set the current parameters if they were not set before. NOTE:
823 // This may occur if the user generated the solver manager with the
824 // default constructor and then didn't set any parameters using
825 // setParameters().
826 if (!isSet_) {
827 setParameters(Teuchos::parameterList(*getValidParameters()));
828 }
829
830 Teuchos::LAPACK<int,ScalarType> lapack;
831
832 TEUCHOS_TEST_FOR_EXCEPTION( !problem_->isProblemSet(),
834 "Belos::BlockCGSolMgr::solve(): Linear problem is not ready, setProblem() "
835 "has not been called.");
836
837 // Create indices for the linear systems to be solved.
838 int startPtr = 0;
839 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
840 int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
841
842 std::vector<int> currIdx, currIdx2;
843 // If an adaptive block size is allowed then only the linear
844 // systems that need to be solved are solved. Otherwise, the index
845 // set is generated that informs the linear problem that some
846 // linear systems are augmented.
847 if ( adaptiveBlockSize_ ) {
848 blockSize_ = numCurrRHS;
849 currIdx.resize( numCurrRHS );
850 currIdx2.resize( numCurrRHS );
851 for (int i=0; i<numCurrRHS; ++i)
852 { currIdx[i] = startPtr+i; currIdx2[i]=i; }
853
854 }
855 else {
856 currIdx.resize( blockSize_ );
857 currIdx2.resize( blockSize_ );
858 for (int i=0; i<numCurrRHS; ++i)
859 { currIdx[i] = startPtr+i; currIdx2[i]=i; }
860 for (int i=numCurrRHS; i<blockSize_; ++i)
861 { currIdx[i] = -1; currIdx2[i] = i; }
862 }
863
864 // Inform the linear problem of the current linear system to solve.
865 problem_->setLSIndex( currIdx );
866
868 // Set up the parameter list for the Iteration subclass.
869 Teuchos::ParameterList plist;
870 plist.set("Block Size",blockSize_);
871
872 // Reset the output status test (controls all the other status tests).
873 outputTest_->reset();
874
875 // Assume convergence is achieved, then let any failed convergence
876 // set this to false. "Innocent until proven guilty."
877 bool isConverged = true;
878
880 // Set up the BlockCG Iteration subclass.
881
882 plist.set("Assert Positive Definiteness", assertPositiveDefiniteness_);
883
884 RCP<CGIteration<ScalarType,MV,OP> > block_cg_iter;
885 if (blockSize_ == 1) {
886 // Standard (nonblock) CG is faster for the special case of a
887 // block size of 1. A single reduction iteration can also be used
888 // if collectives are more expensive than vector updates.
889 plist.set("Fold Convergence Detection Into Allreduce",
890 foldConvergenceDetectionIntoAllreduce_);
891 if (useSingleReduction_) {
892 block_cg_iter =
893 rcp (new CGSingleRedIter<ScalarType,MV,OP> (problem_, printer_,
894 outputTest_, convTest_, plist));
895 }
896 else {
897 block_cg_iter =
898 rcp (new CGIter<ScalarType,MV,OP> (problem_, printer_,
899 outputTest_, convTest_, plist));
900 }
901 } else {
902 block_cg_iter =
903 rcp (new BlockCGIter<ScalarType,MV,OP> (problem_, printer_, outputTest_,
904 ortho_, plist));
905 }
906
907
908 // Enter solve() iterations
909 {
910#ifdef BELOS_TEUCHOS_TIME_MONITOR
911 Teuchos::TimeMonitor slvtimer(*timerSolve_);
912#endif
913
914 while ( numRHS2Solve > 0 ) {
915 //
916 // Reset the active / converged vectors from this block
917 std::vector<int> convRHSIdx;
918 std::vector<int> currRHSIdx( currIdx );
919 currRHSIdx.resize(numCurrRHS);
920
921 // Reset the number of iterations.
922 block_cg_iter->resetNumIters();
923
924 // Reset the number of calls that the status test output knows about.
925 outputTest_->resetNumCalls();
926
927 // Get the current residual for this block of linear systems.
928 RCP<MV> R_0 = MVT::CloneViewNonConst( *(rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
929
930 // Set the new state and initialize the solver.
932 newstate.R = R_0;
933 block_cg_iter->initializeCG(newstate);
934
935 while(1) {
936
937 // tell block_cg_iter to iterate
938 try {
939 block_cg_iter->iterate();
940 //
941 // Check whether any of the linear systems converged.
942 //
943 if (convTest_->getStatus() == Passed) {
944 // At least one of the linear system(s) converged.
945 //
946 // Get the column indices of the linear systems that converged.
947 typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type;
948 std::vector<int> convIdx =
949 rcp_dynamic_cast<conv_test_type>(convTest_)->convIndices();
950
951 // If the number of converged linear systems equals the
952 // number of linear systems currently being solved, then
953 // we are done with this block.
954 if (convIdx.size() == currRHSIdx.size())
955 break; // break from while(1){block_cg_iter->iterate()}
956
957 // Inform the linear problem that we are finished with
958 // this current linear system.
959 problem_->setCurrLS();
960
961 // Reset currRHSIdx to contain the right-hand sides that
962 // are left to converge for this block.
963 int have = 0;
964 std::vector<int> unconvIdx(currRHSIdx.size());
965 for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
966 bool found = false;
967 for (unsigned int j=0; j<convIdx.size(); ++j) {
968 if (currRHSIdx[i] == convIdx[j]) {
969 found = true;
970 break;
971 }
972 }
973 if (!found) {
974 currIdx2[have] = currIdx2[i];
975 currRHSIdx[have++] = currRHSIdx[i];
976 }
977 else {
978 }
979 }
980 currRHSIdx.resize(have);
981 currIdx2.resize(have);
982
983 // Set the remaining indices after deflation.
984 problem_->setLSIndex( currRHSIdx );
985
986 // Get the current residual vector.
987 std::vector<MagnitudeType> norms;
988 R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
989 for (int i=0; i<have; ++i) { currIdx2[i] = i; }
990
991 // Set the new blocksize for the solver.
992 block_cg_iter->setBlockSize( have );
993
994 // Set the new state and initialize the solver.
996 defstate.R = R_0;
997 block_cg_iter->initializeCG(defstate);
998 }
999 //
1000 // None of the linear systems converged. Check whether the
1001 // maximum iteration count was reached.
1002 //
1003 else if (maxIterTest_->getStatus() == Passed) {
1004 isConverged = false; // None of the linear systems converged.
1005 break; // break from while(1){block_cg_iter->iterate()}
1006 }
1007 //
1008 // iterate() returned, but none of our status tests Passed.
1009 // This indicates a bug.
1010 //
1011 else {
1012 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
1013 "Belos::BlockCGSolMgr::solve(): Neither the convergence test nor "
1014 "the maximum iteration count test passed. Please report this bug "
1015 "to the Belos developers.");
1016 }
1017 }
1018 catch (const std::exception &e) {
1019 std::ostream& err = printer_->stream (Errors);
1020 err << "Error! Caught std::exception in CGIteration::iterate() at "
1021 << "iteration " << block_cg_iter->getNumIters() << std::endl
1022 << e.what() << std::endl;
1023 throw;
1024 }
1025 }
1026
1027 // Inform the linear problem that we are finished with this
1028 // block linear system.
1029 problem_->setCurrLS();
1030
1031 // Update indices for the linear systems to be solved.
1032 startPtr += numCurrRHS;
1033 numRHS2Solve -= numCurrRHS;
1034 if ( numRHS2Solve > 0 ) {
1035 numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1036
1037 if ( adaptiveBlockSize_ ) {
1038 blockSize_ = numCurrRHS;
1039 currIdx.resize( numCurrRHS );
1040 currIdx2.resize( numCurrRHS );
1041 for (int i=0; i<numCurrRHS; ++i)
1042 { currIdx[i] = startPtr+i; currIdx2[i] = i; }
1043 }
1044 else {
1045 currIdx.resize( blockSize_ );
1046 currIdx2.resize( blockSize_ );
1047 for (int i=0; i<numCurrRHS; ++i)
1048 { currIdx[i] = startPtr+i; currIdx2[i] = i; }
1049 for (int i=numCurrRHS; i<blockSize_; ++i)
1050 { currIdx[i] = -1; currIdx2[i] = i; }
1051 }
1052 // Set the next indices.
1053 problem_->setLSIndex( currIdx );
1054
1055 // Set the new blocksize for the solver.
1056 block_cg_iter->setBlockSize( blockSize_ );
1057 }
1058 else {
1059 currIdx.resize( numRHS2Solve );
1060 }
1061
1062 }// while ( numRHS2Solve > 0 )
1063
1064 }
1065
1066 // print final summary
1067 sTest_->print( printer_->stream(FinalSummary) );
1068
1069 // print timing information
1070#ifdef BELOS_TEUCHOS_TIME_MONITOR
1071 // Calling summarize() requires communication in general, so don't
1072 // call it unless the user wants to print out timing details.
1073 // summarize() will do all the work even if it's passed a "black
1074 // hole" output stream.
1075 if (verbosity_ & TimingDetails) {
1076 Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
1077 }
1078#endif
1079
1080 // Save the iteration count for this solve.
1081 numIters_ = maxIterTest_->getNumIters();
1082
1083 // Save the convergence test value ("achieved tolerance") for this solve.
1084 {
1085 typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type;
1086 // testValues is nonnull and not persistent.
1087 const std::vector<MagnitudeType>* pTestValues =
1088 rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
1089
1090 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1091 "Belos::BlockCGSolMgr::solve(): The convergence test's getTestValue() "
1092 "method returned NULL. Please report this bug to the Belos developers.");
1093
1094 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1095 "Belos::BlockCGSolMgr::solve(): The convergence test's getTestValue() "
1096 "method returned a vector of length zero. Please report this bug to the "
1097 "Belos developers.");
1098
1099 // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
1100 // achieved tolerances for all vectors in the current solve(), or
1101 // just for the vectors from the last deflation?
1102 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1103 }
1104
1105 if (!isConverged) {
1106 return Unconverged; // return from BlockCGSolMgr::solve()
1107 }
1108 return Converged; // return from BlockCGSolMgr::solve()
1109}
1110
1111// This method requires the solver manager to return a std::string that describes itself.
1112template<class ScalarType, class MV, class OP>
1114{
1115 std::ostringstream oss;
1116 oss << "Belos::BlockCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
1117 oss << "{";
1118 oss << "Ortho Type='"<<orthoType_<<"\', Block Size=" << blockSize_;
1119 oss << "}";
1120 return oss.str();
1121}
1122
1123} // end Belos namespace
1124
1125#endif /* BELOS_BLOCK_CG_SOLMGR_HPP */
Belos concrete class for performing the block conjugate-gradient (CG) iteration.
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.
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.
Definition: BelosTypes.hpp:60
Stub implementation of BlockCGIter, for ScalarType types for which Teuchos::LAPACK does NOT have a va...
BlockCGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i....
BlockCGSolMgrLinearProblemFailure(const std::string &what_arg)
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Return a reference to the linear problem being solved by this solver manager.
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
The linear problem to solve.
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
Teuchos::RCP< Teuchos::ParameterList > params_
Current parameter list.
Teuchos::RCP< Teuchos::Time > timerSolve_
Solve timer.
MagnitudeType orthoKappa_
Orthogonalization parameter (read from parameter list).
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
Teuchos::ScalarTraits< MagnitudeType > MT
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
int maxIters_
Maximum iteration count (read from parameter list).
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
Maximum iteration count stopping criterion.
Teuchos::RCP< std::ostream > outputStream_
Output stream to which the output manager prints.
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
Aggregate stopping criterion.
void reset(const ResetType type) override
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
bool isSet_
Whether or not the parameters have been set (via setParameters()).
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > ortho_
Orthogonalization manager.
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
Output "status test" that controls all the other status tests.
int getNumIters() const override
Get the iteration count for the most recent call to solve().
std::string label_
Prefix label for all the timers.
MagnitudeType achievedTol_
Tolerance achieved by the last solve() invocation.
MagnitudeType convtol_
Convergence tolerance (read from parameter list).
bool isLOADetected() const override
Return whether a loss of accuracy was detected by this solver during the most current solve.
Teuchos::RCP< OutputManager< ScalarType > > printer_
Output manager, that handles printing of different kinds of messages.
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > convTest_
Convergence stopping criterion.
int numIters_
Number of iterations taken by the last solve() invocation.
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.
The Belos::BlockCGSolMgr provides a powerful and fully-featured solver manager over the CG and BlockC...
static const bool requiresLapack
BlockCGSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
Details::SolverManagerRequiresLapack< ScalarType, MV, OP, requiresLapack > base_type
This class implements the preconditioned Conjugate Gradient (CG) iteration.
Definition: BelosCGIter.hpp:79
This class implements the preconditioned single-reduction 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.
Enumeration of all valid Belos (Mat)OrthoManager classes.
Teuchos::RCP< Belos::MatOrthoManager< Scalar, MV, OP > > makeMatOrthoManager(const std::string &ortho, const Teuchos::RCP< const OP > &M, const Teuchos::RCP< OutputManager< Scalar > > &, const std::string &label, const Teuchos::RCP< Teuchos::ParameterList > &params)
Return an instance of the specified MatOrthoManager subclass.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
A class for extending the status testing capabilities of Belos via logical combinations.
An implementation of StatusTestResNorm using a family of residual norms.
std::vector< int > convIndices()
Returns the std::vector containing the indices of the residuals that passed the test.
const std::vector< MagnitudeType > * getTestValue() const
Returns the test value, , computed in most recent call to CheckStatus.
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.
Definition: BelosTypes.cpp:106
NormType
The type of vector norm to compute.
Definition: BelosTypes.hpp:97
@ TwoNorm
Definition: BelosTypes.hpp:98
@ StatusTestDetails
Definition: BelosTypes.hpp:261
@ FinalSummary
Definition: BelosTypes.hpp:259
@ TimingDetails
Definition: BelosTypes.hpp:260
@ Undefined
Definition: BelosTypes.hpp:191
NormType convertStringToNormType(const std::string &normType)
Convert the given string to its NormType enum value.
Definition: BelosTypes.cpp:88
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:155
@ Unconverged
Definition: BelosTypes.hpp:157
@ Converged
Definition: BelosTypes.hpp:156
ScaleType
The type of scaling to use on the residual norm value.
Definition: BelosTypes.hpp:120
ResetType
How to reset the solver.
Definition: BelosTypes.hpp:206
Structure to contain pointers to CGIteration state variables.
Teuchos::RCP< const MV > R
The current residual.
Default parameters common to most Belos solvers.
Definition: BelosTypes.hpp:283
static const double convTol
Default convergence tolerance.
Definition: BelosTypes.hpp:293
static const double orthoKappa
DGKS orthogonalization constant.
Definition: BelosTypes.hpp:299