Belos Version of the Day
Loading...
Searching...
No Matches
BelosPCPGSolMgr.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_PCPG_SOLMGR_HPP
43#define BELOS_PCPG_SOLMGR_HPP
44
48
49#include "BelosConfigDefs.hpp"
50#include "BelosTypes.hpp"
53
54#include "BelosPCPGIter.hpp"
55
62#include "Teuchos_LAPACK.hpp"
63#ifdef BELOS_TEUCHOS_TIME_MONITOR
64# include "Teuchos_TimeMonitor.hpp"
65#endif
66#if defined(HAVE_TEUCHOSCORE_CXX11)
67# include <type_traits>
68#endif // defined(HAVE_TEUCHOSCORE_CXX11)
69#include "Teuchos_TypeTraits.hpp"
70
71namespace Belos {
72
74
75
83 PCPGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
84 {}};
85
91 class PCPGSolMgrOrthoFailure : public BelosError {public:
92 PCPGSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
93 {}};
94
101 class PCPGSolMgrLAPACKFailure : public BelosError {public:
102 PCPGSolMgrLAPACKFailure(const std::string& what_arg) : BelosError(what_arg)
103 {}};
104
106
107
131
132 // Partial specialization for complex ScalarType.
133 // This contains a trivial implementation.
134 // See discussion in the class documentation above.
135 //
136 // FIXME (mfh 09 Sep 2015) This also is a stub for types other than
137 // float or double.
138 template<class ScalarType, class MV, class OP,
139 const bool supportsScalarType =
141 ! Teuchos::ScalarTraits<ScalarType>::isComplex>
143 public Details::SolverManagerRequiresRealLapack<ScalarType, MV, OP,
144 Belos::Details::LapackSupportsScalar<ScalarType>::value &&
145 ! Teuchos::ScalarTraits<ScalarType>::isComplex>
146 {
147 static const bool scalarTypeIsSupported =
149 ! Teuchos::ScalarTraits<ScalarType>::isComplex;
150 typedef Details::SolverManagerRequiresRealLapack<ScalarType, MV, OP,
151 scalarTypeIsSupported> base_type;
152
153 public:
155 base_type ()
156 {}
157 PCPGSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
158 const Teuchos::RCP<Teuchos::ParameterList> &pl) :
159 base_type ()
160 {}
161 virtual ~PCPGSolMgr () {}
162
164 Teuchos::RCP<SolverManager<ScalarType, MV, OP> > clone () const override {
166 }
167 };
168
169 template<class ScalarType, class MV, class OP>
170 class PCPGSolMgr<ScalarType, MV, OP, true> :
171 public Details::SolverManagerRequiresRealLapack<ScalarType, MV, OP, true> {
172 private:
175 typedef Teuchos::ScalarTraits<ScalarType> SCT;
176 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
177 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
178
179 public:
181
182
189 PCPGSolMgr();
190
226 PCPGSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
227 const Teuchos::RCP<Teuchos::ParameterList> &pl );
228
230 virtual ~PCPGSolMgr() {};
231
233 virtual Teuchos::RCP<SolverManager<ScalarType, MV, OP> > clone () const {
234 return Teuchos::rcp(new PCPGSolMgr<ScalarType,MV,OP>);
235 }
237
239
240
244 return *problem_;
245 }
246
249 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
250
253 Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
254
260 Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
261 return Teuchos::tuple(timerSolve_);
262 }
263
269 MagnitudeType achievedTol() const {
270 return achievedTol_;
271 }
272
274 int getNumIters() const {
275 return numIters_;
276 }
277
280 bool isLOADetected() const { return false; }
281
283
285
286
288 void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { problem_ = problem; }
289
291 void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
292
294
296
297
301 void reset( const ResetType type ) { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
303
305
306
324 ReturnType solve();
325
327
330
332 std::string description() const;
333
335
336 private:
337
338 // In the A-inner product, perform an RRQR decomposition without using A unless absolutely necessary. Given
339 // the seed space U and C = A U, find U1 and C1 with span(U1)=span(U) such that C1'U1 = I maintaining C=AU.
340 int ARRQR(int numVecs, int numOrthVecs, const Teuchos::SerialDenseMatrix<int,ScalarType>& D);
341
342 // Linear problem.
343 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
344
345 // Output manager.
346 Teuchos::RCP<OutputManager<ScalarType> > printer_;
347 Teuchos::RCP<std::ostream> outputStream_;
348
349 // Status test.
350 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
351 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
352 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
353 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
354
355 // Orthogonalization manager.
356 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
357
358 // Current parameter list.
359 Teuchos::RCP<Teuchos::ParameterList> params_;
360
361 // Default solver values.
362 static constexpr int maxIters_default_ = 1000;
363 static constexpr int deflatedBlocks_default_ = 2;
364 static constexpr int savedBlocks_default_ = 16;
365 static constexpr int verbosity_default_ = Belos::Errors;
366 static constexpr int outputStyle_default_ = Belos::General;
367 static constexpr int outputFreq_default_ = -1;
368 static constexpr const char * label_default_ = "Belos";
369 static constexpr const char * orthoType_default_ = "ICGS";
370
371 //
372 // Current solver values.
373 //
374
376 MagnitudeType convtol_;
377
379 MagnitudeType orthoKappa_;
380
382 MagnitudeType achievedTol_;
383
385 int numIters_;
386
388 int maxIters_;
389
390 int deflatedBlocks_, savedBlocks_, verbosity_, outputStyle_, outputFreq_;
391 std::string orthoType_;
392
393 // Recycled subspace, its image and the residual
394 Teuchos::RCP<MV> U_, C_, R_;
395
396 // Actual dimension of current recycling subspace (<= savedBlocks_ )
397 int dimU_;
398
399 // Timers.
400 std::string label_;
401 Teuchos::RCP<Teuchos::Time> timerSolve_;
402
403 // Internal state variables.
404 bool isSet_;
405 };
406
407
408// Empty Constructor
409template<class ScalarType, class MV, class OP>
411 outputStream_(Teuchos::rcpFromRef(std::cout)),
412 convtol_(DefaultSolverParameters::convTol),
413 orthoKappa_(DefaultSolverParameters::orthoKappa),
414 achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
415 numIters_(0),
416 maxIters_(maxIters_default_),
417 deflatedBlocks_(deflatedBlocks_default_),
418 savedBlocks_(savedBlocks_default_),
419 verbosity_(verbosity_default_),
420 outputStyle_(outputStyle_default_),
421 outputFreq_(outputFreq_default_),
422 orthoType_(orthoType_default_),
423 dimU_(0),
424 label_(label_default_),
425 isSet_(false)
426{}
427
428
429// Basic Constructor
430template<class ScalarType, class MV, class OP>
432 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
433 const Teuchos::RCP<Teuchos::ParameterList> &pl ) :
434 problem_(problem),
435 outputStream_(Teuchos::rcpFromRef(std::cout)),
436
437 convtol_(DefaultSolverParameters::convTol),
438 orthoKappa_(DefaultSolverParameters::orthoKappa),
439 achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
440 numIters_(0),
441 maxIters_(maxIters_default_),
442 deflatedBlocks_(deflatedBlocks_default_),
443 savedBlocks_(savedBlocks_default_),
444 verbosity_(verbosity_default_),
445 outputStyle_(outputStyle_default_),
446 outputFreq_(outputFreq_default_),
447 orthoType_(orthoType_default_),
448 dimU_(0),
449 label_(label_default_),
450 isSet_(false)
451{
452 TEUCHOS_TEST_FOR_EXCEPTION(
453 problem_.is_null (), std::invalid_argument,
454 "Belos::PCPGSolMgr two-argument constructor: "
455 "'problem' is null. You must supply a non-null Belos::LinearProblem "
456 "instance when calling this constructor.");
457
458 if (! pl.is_null ()) {
459 // Set the parameters using the list that was passed in.
460 setParameters (pl);
461 }
462}
463
464
465template<class ScalarType, class MV, class OP>
466void PCPGSolMgr<ScalarType,MV,OP,true>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
467{
468 // Create the internal parameter list if ones doesn't already exist.
469 if (params_ == Teuchos::null) {
470 params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
471 }
472 else {
473 params->validateParameters(*getValidParameters());
474 }
475
476 // Check for maximum number of iterations
477 if (params->isParameter("Maximum Iterations")) {
478 maxIters_ = params->get("Maximum Iterations",maxIters_default_);
479
480 // Update parameter in our list and in status test.
481 params_->set("Maximum Iterations", maxIters_);
482 if (maxIterTest_!=Teuchos::null)
483 maxIterTest_->setMaxIters( maxIters_ );
484 }
485
486 // Check for the maximum numbers of saved and deflated blocks.
487 if (params->isParameter("Num Saved Blocks")) {
488 savedBlocks_ = params->get("Num Saved Blocks",savedBlocks_default_);
489 TEUCHOS_TEST_FOR_EXCEPTION(savedBlocks_ <= 0, std::invalid_argument,
490 "Belos::PCPGSolMgr: \"Num Saved Blocks\" must be strictly positive.");
491
492 // savedBlocks > number of matrix rows and columns, not known in parameters.
493 //TEUCHOS_TEST_FOR_EXCEPTION(savedBlocks_ >= maxIters_, std::invalid_argument,
494 //"Belos::PCPGSolMgr: \"Num Saved Blocks\" must be less than \"Maximum Iterations\".");
495
496 // Update parameter in our list.
497 params_->set("Num Saved Blocks", savedBlocks_);
498 }
499 if (params->isParameter("Num Deflated Blocks")) {
500 deflatedBlocks_ = params->get("Num Deflated Blocks",deflatedBlocks_default_);
501 TEUCHOS_TEST_FOR_EXCEPTION(deflatedBlocks_ < 0, std::invalid_argument,
502 "Belos::PCPGSolMgr: \"Num Deflated Blocks\" must be positive.");
503
504 TEUCHOS_TEST_FOR_EXCEPTION(deflatedBlocks_ > savedBlocks_, std::invalid_argument,
505 "Belos::PCPGSolMgr: \"Num Deflated Blocks\" must be <= \"Num Saved Blocks\".");
506
507 // Update parameter in our list.
508 // The static_cast is for clang link issues with the constexpr before c++17
509 params_->set("Num Deflated Blocks", static_cast<int>(deflatedBlocks_));
510 }
511
512 // Check to see if the timer label changed.
513 if (params->isParameter("Timer Label")) {
514 std::string tempLabel = params->get("Timer Label", label_default_);
515
516 // Update parameter in our list and solver timer
517 if (tempLabel != label_) {
518 label_ = tempLabel;
519 params_->set("Timer Label", label_);
520 std::string solveLabel = label_ + ": PCPGSolMgr total solve time";
521#ifdef BELOS_TEUCHOS_TIME_MONITOR
522 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
523#endif
524 if (ortho_ != Teuchos::null) {
525 ortho_->setLabel( label_ );
526 }
527 }
528 }
529
530 // Check for a change in verbosity level
531 if (params->isParameter("Verbosity")) {
532 if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
533 verbosity_ = params->get("Verbosity", verbosity_default_);
534 } else {
535 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
536 }
537
538 // Update parameter in our list.
539 params_->set("Verbosity", verbosity_);
540 if (printer_ != Teuchos::null)
541 printer_->setVerbosity(verbosity_);
542 }
543
544 // Check for a change in output style
545 if (params->isParameter("Output Style")) {
546 if (Teuchos::isParameterType<int>(*params,"Output Style")) {
547 outputStyle_ = params->get("Output Style", outputStyle_default_);
548 } else {
549 outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
550 }
551
552 // Reconstruct the convergence test if the explicit residual test is not being used.
553 params_->set("Output Style", outputStyle_);
554 outputTest_ = Teuchos::null;
555 }
556
557 // output stream
558 if (params->isParameter("Output Stream")) {
559 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
560
561 // Update parameter in our list.
562 params_->set("Output Stream", outputStream_);
563 if (printer_ != Teuchos::null)
564 printer_->setOStream( outputStream_ );
565 }
566
567 // frequency level
568 if (verbosity_ & Belos::StatusTestDetails) {
569 if (params->isParameter("Output Frequency")) {
570 outputFreq_ = params->get("Output Frequency", outputFreq_default_);
571 }
572
573 // Update parameter in out list and output status test.
574 params_->set("Output Frequency", outputFreq_);
575 if (outputTest_ != Teuchos::null)
576 outputTest_->setOutputFrequency( outputFreq_ );
577 }
578
579 // Create output manager if we need to.
580 if (printer_ == Teuchos::null) {
581 printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
582 }
583
584 // Check if the orthogonalization changed.
585 bool changedOrthoType = false;
586 if (params->isParameter("Orthogonalization")) {
587 std::string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
588 if (tempOrthoType != orthoType_) {
589 orthoType_ = tempOrthoType;
590 changedOrthoType = true;
591 }
592 }
593 params_->set("Orthogonalization", orthoType_);
594
595 // Check which orthogonalization constant to use.
596 if (params->isParameter("Orthogonalization Constant")) {
597 if (params->isType<MagnitudeType> ("Orthogonalization Constant")) {
598 orthoKappa_ = params->get ("Orthogonalization Constant",
599 static_cast<MagnitudeType> (DefaultSolverParameters::orthoKappa));
600 }
601 else {
602 orthoKappa_ = params->get ("Orthogonalization Constant",
604 }
605
606 // Update parameter in our list.
607 params_->set("Orthogonalization Constant",orthoKappa_);
608 if (orthoType_=="DGKS") {
609 if (orthoKappa_ > 0 && ortho_ != Teuchos::null && !changedOrthoType) {
610 Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
611 }
612 }
613 }
614
615 // Create orthogonalization manager if we need to.
616 if (ortho_ == Teuchos::null || changedOrthoType) {
618 Teuchos::RCP<Teuchos::ParameterList> paramsOrtho; // can be null
619 if (orthoType_=="DGKS" && orthoKappa_ > 0) {
620 paramsOrtho->set ("depTol", orthoKappa_ );
621 }
622
623 ortho_ = factory.makeMatOrthoManager (orthoType_, Teuchos::null, printer_, label_, paramsOrtho);
624 }
625
626 // Convergence
627 typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
628 typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t;
629
630 // Check for convergence tolerance
631 if (params->isParameter("Convergence Tolerance")) {
632 if (params->isType<MagnitudeType> ("Convergence Tolerance")) {
633 convtol_ = params->get ("Convergence Tolerance",
634 static_cast<MagnitudeType> (DefaultSolverParameters::convTol));
635 }
636 else {
637 convtol_ = params->get ("Convergence Tolerance", DefaultSolverParameters::convTol);
638 }
639
640 // Update parameter in our list and residual tests.
641 params_->set("Convergence Tolerance", convtol_);
642 if (convTest_ != Teuchos::null)
643 convTest_->setTolerance( convtol_ );
644 }
645
646 // Create status tests if we need to.
647
648 // Basic test checks maximum iterations and native residual.
649 if (maxIterTest_ == Teuchos::null)
650 maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
651
652 if (convTest_ == Teuchos::null)
653 convTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_, 1 ) );
654
655 sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
656
657 // Create the status test output class.
658 // This class manages and formats the output from the status test.
659 StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
660 outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
661
662 // Set the solver string for the output test
663 std::string solverDesc = " PCPG ";
664 outputTest_->setSolverDesc( solverDesc );
665
666 // Create the timer if we need to.
667 if (timerSolve_ == Teuchos::null) {
668 std::string solveLabel = label_ + ": PCPGSolMgr total solve time";
669#ifdef BELOS_TEUCHOS_TIME_MONITOR
670 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
671#endif
672 }
673
674 // Inform the solver manager that the current parameters were set.
675 isSet_ = true;
676}
677
678
679template<class ScalarType, class MV, class OP>
680Teuchos::RCP<const Teuchos::ParameterList>
682{
683 static Teuchos::RCP<const Teuchos::ParameterList> validPL;
684 if (is_null(validPL)) {
685 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
686 // Set all the valid parameters and their default values.
687 pl->set("Convergence Tolerance", static_cast<MagnitudeType>(DefaultSolverParameters::convTol),
688 "The relative residual tolerance that needs to be achieved by the\n"
689 "iterative solver in order for the linear system to be declared converged.");
690 pl->set("Maximum Iterations", static_cast<int>(maxIters_default_),
691 "The maximum number of iterations allowed for each\n"
692 "set of RHS solved.");
693 pl->set("Num Deflated Blocks", static_cast<int>(deflatedBlocks_default_),
694 "The maximum number of vectors in the seed subspace." );
695 pl->set("Num Saved Blocks", static_cast<int>(savedBlocks_default_),
696 "The maximum number of vectors saved from old Krylov subspaces." );
697 pl->set("Verbosity", static_cast<int>(verbosity_default_),
698 "What type(s) of solver information should be outputted\n"
699 "to the output stream.");
700 pl->set("Output Style", static_cast<int>(outputStyle_default_),
701 "What style is used for the solver information outputted\n"
702 "to the output stream.");
703 pl->set("Output Frequency", static_cast<int>(outputFreq_default_),
704 "How often convergence information should be outputted\n"
705 "to the output stream.");
706 pl->set("Output Stream", Teuchos::rcpFromRef(std::cout),
707 "A reference-counted pointer to the output stream where all\n"
708 "solver output is sent.");
709 pl->set("Timer Label", static_cast<const char *>(label_default_),
710 "The string to use as a prefix for the timer labels.");
711 pl->set("Orthogonalization", static_cast<const char *>(orthoType_default_),
712 "The type of orthogonalization to use: DGKS, ICGS, IMGS");
713 pl->set("Orthogonalization Constant",static_cast<MagnitudeType>(DefaultSolverParameters::orthoKappa),
714 "The constant used by DGKS orthogonalization to determine\n"
715 "whether another step of classical Gram-Schmidt is necessary.");
716 validPL = pl;
717 }
718 return validPL;
719}
720
721
722// solve()
723template<class ScalarType, class MV, class OP>
725
726 // Set the current parameters if are not set already.
727 if (!isSet_) { setParameters( params_ ); }
728
729 Teuchos::LAPACK<int,ScalarType> lapack;
730 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
731 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
732
733 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,PCPGSolMgrLinearProblemFailure,
734 "Belos::PCPGSolMgr::solve(): Linear problem is not a valid object.");
735
736 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),PCPGSolMgrLinearProblemFailure,
737 "Belos::PCPGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
738
739 // Create indices for the linear systems to be solved.
740 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
741 std::vector<int> currIdx(1);
742 currIdx[0] = 0;
743
744 bool debug = false;
745
746 // Inform the linear problem of the current linear system to solve.
747 problem_->setLSIndex( currIdx ); // block size == 1
748
749 // Assume convergence is achieved, then let any failed convergence set this to false.
750 bool isConverged = true;
751
753 // PCPG iteration parameter list
754 Teuchos::ParameterList plist;
755 plist.set("Saved Blocks", savedBlocks_);
756 plist.set("Block Size", 1);
757 plist.set("Keep Diagonal", true);
758 plist.set("Initialize Diagonal", true);
759
761 // PCPG solver
762
763 Teuchos::RCP<PCPGIter<ScalarType,MV,OP> > pcpg_iter;
764 pcpg_iter = Teuchos::rcp( new PCPGIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,plist) );
765 // Number of iterations required to generate initial recycle space (if needed)
766
767 // Enter solve() iterations
768 {
769#ifdef BELOS_TEUCHOS_TIME_MONITOR
770 Teuchos::TimeMonitor slvtimer(*timerSolve_);
771#endif
772 while ( numRHS2Solve > 0 ) { // test for quick return
773
774 // Reset the status test.
775 outputTest_->reset();
776
777 // Create the first block in the current Krylov basis (residual).
778 if (R_ == Teuchos::null)
779 R_ = MVT::Clone( *(problem_->getRHS()), 1 );
780
781 problem_->computeCurrResVec( &*R_ );
782
783
784 // Hypothesis: if U_ is not null, then neither is C_ and furthermore U'C= I.
785 // TODO: ensure hypothesis right here ... I have to think about use cases.
786
787 if( U_ != Teuchos::null ){
788 // Hypothesis: if U_ is not null, then neither is C_ and furthermore U'C= I.
789
790 // possibly over solved equation ... I want residual norms
791 // relative to the initial residual, not what I am about to compute.
792 Teuchos::RCP<MV> cur_soln_vec = problem_->getCurrLHSVec();
793 std::vector<MagnitudeType> rnorm0(1);
794 MVT::MvNorm( *R_, rnorm0 ); // rnorm0 = norm(R_);
795
796 // Z := U_'*R_; xo += U_*Z ;R_ -= C_*Z
797 std::cout << "Solver Manager: dimU_ = " << dimU_ << std::endl;
798 Teuchos::SerialDenseMatrix<int,ScalarType> Z( dimU_, 1 );
799
800 Teuchos::RCP<const MV> Uactive, Cactive;
801 std::vector<int> active_columns( dimU_ );
802 for (int i=0; i < dimU_; ++i) active_columns[i] = i;
803 Uactive = MVT::CloneView(*U_, active_columns);
804 Cactive = MVT::CloneView(*C_, active_columns);
805
806 if( debug ){
807 std::cout << " Solver Manager : check duality of seed basis " << std::endl;
808 Teuchos::SerialDenseMatrix<int,ScalarType> H( dimU_, dimU_ );
809 MVT::MvTransMv( one, *Uactive, *Cactive, H );
810 H.print( std::cout );
811 }
812
813 MVT::MvTransMv( one, *Uactive, *R_, Z );
814 Teuchos::RCP<MV> tempU = MVT::Clone( *R_, 1 );
815 MVT::MvTimesMatAddMv( one, *Uactive, Z, zero, *tempU ); // UZ
816 MVT::MvAddMv( one, *tempU, one, *cur_soln_vec, *cur_soln_vec ); // xo += tmp;
817 MVT::MvTimesMatAddMv( one, *Cactive, Z, zero, *tempU ); // CZ
818 MVT::MvAddMv( -one, *tempU, one, *R_, *R_ ); // R_ -= tmp;
819 std::vector<MagnitudeType> rnorm(1);
820 MVT::MvNorm( *R_, rnorm );
821 if( rnorm[0] < rnorm0[0] * .001 ){ //reorthogonalize
822 MVT::MvTransMv( one, *Uactive, *R_, Z );
823 MVT::MvTimesMatAddMv( one, *Uactive, Z, zero, *tempU );
824 MVT::MvAddMv( one, *tempU, one, *cur_soln_vec, *cur_soln_vec ); // xo += UZ;
825 MVT::MvTimesMatAddMv( one, *Cactive, Z, zero, *tempU );
826 MVT::MvAddMv( -one, *tempU, one, *R_, *R_ ); // R_ -= CZ;
827 }
828 Uactive = Teuchos::null;
829 Cactive = Teuchos::null;
830 tempU = Teuchos::null;
831 }
832 else {
833 dimU_ = 0;
834 }
835
836
837 // Set the new state and initialize the solver.
838 PCPGIterState<ScalarType,MV> pcpgState; // fails if R == null.
839
840 pcpgState.R = R_;
841 if( U_ != Teuchos::null ) pcpgState.U = U_;
842 if( C_ != Teuchos::null ) pcpgState.C = C_;
843 if( dimU_ > 0 ) pcpgState.curDim = dimU_;
844 pcpg_iter->initialize(pcpgState);
845
846 // treat initialize() exceptions here? how to use try-catch-throw? DMD
847
848 // Get the current number of deflated blocks with the PCPG iteration
849 dimU_ = pcpgState.curDim;
850 if( !dimU_ ) printer_->stream(Debug) << " No recycled subspace available for RHS index " << currIdx[0] << std::endl << std::endl;
851 pcpg_iter->resetNumIters();
852
853 if( dimU_ > savedBlocks_ )
854 std::cout << "Error: dimU_ = " << dimU_ << " > savedBlocks_ = " << savedBlocks_ << std::endl;
855
856 while(1) { // dummy loop for break
857
858 // tell pcpg_iter to iterate
859 try {
860 if( debug ) printf("********** Calling iterate...\n");
861 pcpg_iter->iterate();
862
864 //
865 // check convergence first
866 //
868 if ( convTest_->getStatus() == Passed ) {
869 // we have convergence
870 break; // break from while(1){pcpg_iter->iterate()}
871 }
873 //
874 // check for maximum iterations
875 //
877 else if ( maxIterTest_->getStatus() == Passed ) {
878 // we don't have convergence
879 isConverged = false;
880 break; // break from while(1){pcpg_iter->iterate()}
881 }
882 else {
883
885 //
886 // we returned from iterate(), but none of our status tests Passed.
887 // Something is wrong, and it is probably the developers fault.
888 //
890
891 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
892 "Belos::PCPGSolMgr::solve(): Invalid return from PCPGIter::iterate().");
893 } // end if
894 } // end try
895 catch (const std::exception &e) {
896 printer_->stream(Errors) << "Error! Caught exception in PCPGIter::iterate() at iteration "
897 << pcpg_iter->getNumIters() << std::endl
898 << e.what() << std::endl;
899 throw;
900 }
901 } // end of while(1)
902
903 // Update the linear problem.
904 Teuchos::RCP<MV> update = pcpg_iter->getCurrentUpdate();
905 problem_->updateSolution( update, true );
906
907 // Inform the linear problem that we are finished with this block linear system.
908 problem_->setCurrLS();
909
910 // Get the state. How did pcpgState die?
911 PCPGIterState<ScalarType,MV> oldState = pcpg_iter->getState();
912
913 dimU_ = oldState.curDim;
914 int q = oldState.prevUdim;
915
916 std::cout << "SolverManager: dimU_ " << dimU_ << " prevUdim= " << q << std::endl;
917
918 if( q > deflatedBlocks_ )
919 std::cout << "SolverManager: Error deflatedBlocks = " << deflatedBlocks_ << std::endl;
920
921 int rank;
922 if( dimU_ > q ){ // Orthogonalize [U;C](:,prevUdim:dimU_)
923 //Given the seed space U and C = A U for some symmetric positive definite A,
924 //find U1 and C1 with span(U1)=span(U) such that C1'U1 = I maintaining C=AU
925
926 //oldState.D->print( std::cout ); D = diag( C'*U )
927
928 U_ = oldState.U; //MVT::MvPrint( *U_, std::cout );
929 C_ = oldState.C; //MVT::MvPrint( *C_, std::cout );
930 rank = ARRQR(dimU_,q, *oldState.D );
931 if( rank < dimU_ ) {
932 std::cout << " rank decreased in ARRQR, something to do? " << std::endl;
933 }
934 dimU_ = rank;
935
936 } // Now U_ and C_ = AU are dual bases.
937
938 if( dimU_ > deflatedBlocks_ ){
939
940 if( !deflatedBlocks_ ){
941 U_ = Teuchos::null;
942 C_ = Teuchos::null;
943 dimU_ = deflatedBlocks_;
944 break;
945 }
946
947 bool Harmonic = false; // (Harmonic) Ritz vectors
948
949 Teuchos::RCP<MV> Uorth;
950
951 std::vector<int> active_cols( dimU_ );
952 for (int i=0; i < dimU_; ++i) active_cols[i] = i;
953
954 if( Harmonic ){
955 Uorth = MVT::CloneCopy(*C_, active_cols);
956 }
957 else{
958 Uorth = MVT::CloneCopy(*U_, active_cols);
959 }
960
961 // Explicitly construct Q and R factors
962 Teuchos::SerialDenseMatrix<int,ScalarType> R(dimU_,dimU_);
963 rank = ortho_->normalize(*Uorth, Teuchos::rcp(&R,false));
964 Uorth = Teuchos::null;
965 // TODO: During the previous solve, the matrix that normalizes U(1:q) was computed and discarded.
966 // One might save it, reuse it here, and just normalize columns U(q+1:dimU_) here.
967
968 // throw an error if U is both A-orthonormal and rank deficient
969 TEUCHOS_TEST_FOR_EXCEPTION(rank < dimU_,PCPGSolMgrOrthoFailure,
970 "Belos::PCPGSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
971
972
973 // R VT' = Ur S,
974 Teuchos::SerialDenseMatrix<int,ScalarType> VT; // Not referenced
975 Teuchos::SerialDenseMatrix<int,ScalarType> Ur; // Not referenced
976 int lwork = 5*dimU_; // minimal, extra computation < 67*dimU_
977 int info = 0; // Hermite
978 int lrwork = 1;
979 if( problem_->isHermitian() ) lrwork = dimU_;
980 std::vector<ScalarType> work(lwork); //
981 std::vector<ScalarType> Svec(dimU_); //
982 std::vector<ScalarType> rwork(lrwork);
983 lapack.GESVD('N', 'O',
984 R.numRows(),R.numCols(),R.values(), R.numRows(),
985 &Svec[0],
986 Ur.values(),1,
987 VT.values(),1, // Output: VT stored in R
988 &work[0], lwork,
989 &rwork[0], &info);
990
991 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, PCPGSolMgrLAPACKFailure,
992 "Belos::PCPGSolMgr::solve(): LAPACK _GESVD failed to compute singular values.");
993
994 if( work[0] != 67. * dimU_ )
995 std::cout << " SVD " << dimU_ << " lwork " << work[0] << std::endl;
996 for( int i=0; i< dimU_; i++)
997 std::cout << i << " " << Svec[i] << std::endl;
998
999 Teuchos::SerialDenseMatrix<int,ScalarType> wholeV( R, Teuchos::TRANS);
1000
1001 int startRow = 0, startCol = 0;
1002 if( Harmonic )
1003 startCol = dimU_ - deflatedBlocks_;
1004
1005 Teuchos::SerialDenseMatrix<int,ScalarType> V(Teuchos::Copy,
1006 wholeV,
1007 wholeV.numRows(),
1008 deflatedBlocks_,
1009 startRow,
1010 startCol);
1011 std::vector<int> active_columns( dimU_ );
1012 std::vector<int> def_cols( deflatedBlocks_ );
1013 for (int i=0; i < dimU_; ++i) active_columns[i] = i;
1014 for (int i=0; i < deflatedBlocks_; ++i) def_cols[i] = i;
1015
1016 Teuchos::RCP<MV> Uactive = MVT::CloneViewNonConst(*U_, def_cols);
1017 Teuchos::RCP<MV> Ucopy = MVT::CloneCopy( *U_, active_columns );
1018 MVT::MvTimesMatAddMv( one, *Ucopy, V, zero, *Uactive ); // U:= U*V
1019 Ucopy = Teuchos::null;
1020 Uactive = Teuchos::null;
1021 Teuchos::RCP<MV> Cactive = MVT::CloneViewNonConst(*C_, def_cols);
1022 Teuchos::RCP<MV> Ccopy = MVT::CloneCopy( *C_, active_columns );
1023 MVT::MvTimesMatAddMv( one, *Ccopy, V, zero, *Cactive ); // C:= C*V
1024 Ccopy = Teuchos::null;
1025 Cactive = Teuchos::null;
1026 dimU_ = deflatedBlocks_;
1027 }
1028 printer_->stream(Debug) << " Generated recycled subspace using RHS index " << currIdx[0] << " of dimension " << dimU_ << std::endl << std::endl;
1029
1030 // Inform the linear problem that we are finished with this block linear system.
1031 problem_->setCurrLS();
1032
1033 // Update indices for the linear systems to be solved.
1034 numRHS2Solve -= 1;
1035 if ( numRHS2Solve > 0 ) {
1036 currIdx[0]++;
1037
1038 // Set the next indices.
1039 problem_->setLSIndex( currIdx );
1040 }
1041 else {
1042 currIdx.resize( numRHS2Solve );
1043 }
1044 }// while ( numRHS2Solve > 0 )
1045 }
1046
1047 // print final summary
1048 sTest_->print( printer_->stream(FinalSummary) );
1049
1050 // print timing information
1051#ifdef BELOS_TEUCHOS_TIME_MONITOR
1052 // Calling summarize() can be expensive, so don't call unless the
1053 // user wants to print out timing details. summarize() will do all
1054 // the work even if it's passed a "black hole" output stream.
1055 if (verbosity_ & TimingDetails)
1056 Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
1057#endif
1058
1059 // Save the convergence test value ("achieved tolerance") for this solve.
1060 {
1061 using Teuchos::rcp_dynamic_cast;
1062 typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type;
1063 // testValues is nonnull and not persistent.
1064 const std::vector<MagnitudeType>* pTestValues =
1065 rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
1066
1067 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1068 "Belos::PCPGSolMgr::solve(): The convergence test's getTestValue() "
1069 "method returned NULL. Please report this bug to the Belos developers.");
1070
1071 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1072 "Belos::PCPGSolMgr::solve(): The convergence test's getTestValue() "
1073 "method returned a vector of length zero. Please report this bug to the "
1074 "Belos developers.");
1075
1076 // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
1077 // achieved tolerances for all vectors in the current solve(), or
1078 // just for the vectors from the last deflation?
1079 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1080 }
1081
1082 // get iteration information for this solve
1083 numIters_ = maxIterTest_->getNumIters();
1084
1085 if (!isConverged) {
1086 return Unconverged; // return from PCPGSolMgr::solve()
1087 }
1088 return Converged; // return from PCPGSolMgr::solve()
1089}
1090
1091// A-orthogonalize the Seed Space
1092// Note that Anasazi::GenOrthoManager provides simplified versions of the algorithm,
1093// that are not rank revealing, and are not designed for PCPG in other ways too.
1094template<class ScalarType, class MV, class OP>
1095int PCPGSolMgr<ScalarType,MV,OP,true>::ARRQR(int p, int q, const Teuchos::SerialDenseMatrix<int,ScalarType>& D)
1096{
1097 using Teuchos::RCP;
1098 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1099 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1100
1101 // Allocate memory for scalars.
1102 Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
1103 Teuchos::SerialDenseMatrix<int,ScalarType> gamma( 1, 1 );
1104 Teuchos::SerialDenseMatrix<int,ScalarType> anorm( 1, 1 );
1105 std::vector<int> curind(1);
1106 std::vector<int> ipiv(p - q); // RRQR Pivot indices
1107 std::vector<ScalarType> Pivots(p); // RRQR Pivots
1108 int i, imax, j, l;
1109 ScalarType rteps = 1.5e-8;
1110
1111 // Scale such that diag( U'C) = I
1112 for( i = q ; i < p ; i++ ){
1113 ipiv[i-q] = i;
1114 curind[0] = i;
1115 RCP<MV> P = MVT::CloneViewNonConst(*U_,curind);
1116 RCP<MV> AP = MVT::CloneViewNonConst(*C_,curind);
1117 anorm(0,0) = one / Teuchos::ScalarTraits<ScalarType>::squareroot( D(i-q,i-q) ) ;
1118 MVT::MvAddMv( anorm(0,0), *P, zero, *AP, *P );
1119 MVT::MvAddMv( zero, *P, anorm(0,0), *AP, *AP );
1120 Pivots[i] = one;
1121 }
1122
1123 for( i = q ; i < p ; i++ ){
1124 if( q < i && i < p-1 ){ // Find the largest pivot
1125 imax = i;
1126 l = ipiv[imax-q];
1127 for( j = i+1 ; j < p ; j++ ){
1128 const int k = ipiv[j-q];
1129 if( Pivots[k] > Pivots[l] ){
1130 imax = j;
1131 l = k;
1132 }
1133 } // end for
1134 if( imax > i ){
1135 l = ipiv[imax-q]; // swap ipiv( imax ) and ipiv(i+1)
1136 ipiv[imax-q] = ipiv[i-q];
1137 ipiv[i-q] = l;
1138 }
1139 } // largest pivot found
1140 int k = ipiv[i-q];
1141
1142 if( Pivots[k] > 1.5625e-2 ){
1143 anorm(0,0) = Pivots[k]; // A-norm of u
1144 }
1145 else{ // anorm(0,0) = sqrt( U(:,k)'*C(:,k) );
1146 curind[0] = k;
1147 RCP<const MV> P = MVT::CloneView(*U_,curind);
1148 RCP<const MV> AP = MVT::CloneView(*C_,curind);
1149 MVT::MvTransMv( one, *P, *AP, anorm );
1150 anorm(0,0) = Teuchos::ScalarTraits<ScalarType>::squareroot( anorm(0,0) ) ;
1151 }
1152 if( rteps <= anorm(0,0) && anorm(0,0) < 9.765625e-4){
1153 /*
1154 C(:,k) = A*U(:,k); % Change C
1155 fixC = U(:, ipiv(1:i-1) )'*C(:,k);
1156 U(:,k) = U(:,k) - U(:, ipiv(1:i-1) )*fixC;
1157 C(:,k) = C(:,k) - C(:, ipiv(1:i-1) )*fixC;
1158 anorm = sqrt( U(:,k)'*C(:,k) );
1159 */
1160 std::cout << "ARRQR: Bad case not implemented" << std::endl;
1161 }
1162 if( anorm(0,0) < rteps ){ // rank [U;C] = i-1
1163 std::cout << "ARRQR : deficient case not implemented " << std::endl;
1164 //U = U(:, ipiv(1:i-1) );
1165 //C = C(:, ipiv(1:i-1) );
1166 p = q + i;
1167 // update curDim_ in State
1168 break;
1169 }
1170 curind[0] = k;
1171 RCP<MV> P = MVT::CloneViewNonConst(*U_,curind);
1172 RCP<MV> AP = MVT::CloneViewNonConst(*C_,curind);
1173 MVT::MvAddMv( anorm(0,0), *P, zero, *AP, *P ); // U(:,k) = U(:,k)/anorm;
1174 MVT::MvAddMv( zero, *P, anorm(0,0), *AP, *AP ); // C(:,k) = C(:,k)/anorm;
1175 P = Teuchos::null;
1176 AP = Teuchos::null;
1177 Pivots[k] = one; // delete, for diagonostic purposes
1178 P = MVT::CloneViewNonConst(*U_,curind); // U(:,k)
1179 AP = MVT::CloneViewNonConst(*C_,curind); // C(:,k)
1180 for( j = i+1 ; j < p ; j++ ){
1181 l = ipiv[j-q]; // ahhh
1182 curind[0] = l;
1183 RCP<MV> Q = MVT::CloneViewNonConst(*U_,curind); // segmentation fault, j=i+1=5
1184 MVT::MvTransMv( one, *Q, *AP, alpha); // alpha(0,0) = U(:,l)'*C(:,k);
1185 MVT::MvAddMv( -alpha(0,0), *P, one, *Q, *Q ); // U(:,l) -= U(:,k) * alpha(0,0);
1186 Q = Teuchos::null;
1187 RCP<MV> AQ = MVT::CloneViewNonConst(*C_,curind);
1188 MVT::MvAddMv( -alpha(0,0), *AP, one, *AQ, *AQ ); // C(:,l) -= C(:,l) - C(:,k) * alpha(0,0);
1189 AQ = Teuchos::null;
1190 gamma(0,0) = ( Pivots[l] - alpha(0,0))*( Pivots[l] + alpha(0,0));
1191 if( gamma(0,0) > 0){
1192 Pivots[l] = Teuchos::ScalarTraits<ScalarType>::squareroot( gamma(0,0) );
1193 }
1194 else {
1195 Pivots[l] = zero; //rank deficiency revealed
1196 }
1197 }
1198 }
1199 return p;
1200}
1201
1202// The method returns a string describing the solver manager.
1203template<class ScalarType, class MV, class OP>
1205{
1206 std::ostringstream oss;
1207 oss << "Belos::PCPGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
1208 oss << "{";
1209 oss << "Ortho Type='"<<orthoType_;
1210 oss << "}";
1211 return oss.str();
1212}
1213
1214} // end Belos namespace
1215
1216#endif /* BELOS_PCPG_SOLMGR_HPP */
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 to iterate Preconditioned Conjugate Projected Gradients.
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
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 real ScalarType t...
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...
This class implements the PCPG iteration, where a single-std::vector Krylov subspace is constructed....
PCPGSolMgrLAPACKFailure is thrown when a nonzero value is retuned from an LAPACK call.
PCPGSolMgrLAPACKFailure(const std::string &what_arg)
PCPGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i.e.
PCPGSolMgrLinearProblemFailure(const std::string &what_arg)
PCPGSolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate orthonormal...
PCPGSolMgrOrthoFailure(const std::string &what_arg)
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get current linear problem being solved for in this object.
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const
Get a parameter list containing the current parameters for this object.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem)
Set the linear problem that needs to be solved.
bool isLOADetected() const
Return whether a loss of accuracy was detected by this solver during the most current solve.
int getNumIters() const
Get the iteration count for the most recent call to solve().
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
virtual Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const
clone for Inverted Injection (DII)
MagnitudeType achievedTol() const
Tolerance achieved by the last solve() invocation.
void reset(const ResetType type)
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
PCPG iterative linear solver.
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
PCPGSolMgr(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.
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.
@ StatusTestDetails
Definition: BelosTypes.hpp:261
@ FinalSummary
Definition: BelosTypes.hpp:259
@ TimingDetails
Definition: BelosTypes.hpp:260
@ Undefined
Definition: BelosTypes.hpp:191
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:155
@ Unconverged
Definition: BelosTypes.hpp:157
@ Converged
Definition: BelosTypes.hpp:156
ResetType
How to reset the solver.
Definition: BelosTypes.hpp:206
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
Structure to contain pointers to PCPGIter state variables.
int prevUdim
Number of block columns in matrices C and U before current iteration.
Teuchos::RCP< MV > U
The recycled subspace.
Teuchos::RCP< MV > C
C = AU, U spans recycled subspace.
int curDim
The current dimension of the reduction.
Teuchos::RCP< MV > R
The current residual.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > D
The current Hessenberg matrix.

Generated for Belos by doxygen 1.9.6