Belos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
BelosPseudoBlockTFQMRSolMgr.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_TFQMR_SOLMGR_HPP
43#define BELOS_PSEUDO_BLOCK_TFQMR_SOLMGR_HPP
44
49#include "BelosConfigDefs.hpp"
50#include "BelosTypes.hpp"
51
54
61#ifdef BELOS_TEUCHOS_TIME_MONITOR
62#include "Teuchos_TimeMonitor.hpp"
63#endif
64
78namespace Belos {
79
81
82
90 PseudoBlockTFQMRSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
91 {}};
92
93 template<class ScalarType, class MV, class OP>
94 class PseudoBlockTFQMRSolMgr : public SolverManager<ScalarType,MV,OP> {
95
96 private:
99 typedef Teuchos::ScalarTraits<ScalarType> SCT;
100 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
101 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
102
103 public:
104
106
107
114
131 PseudoBlockTFQMRSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
132 const Teuchos::RCP<Teuchos::ParameterList> &pl );
133
136
138 Teuchos::RCP<SolverManager<ScalarType, MV, OP> > clone () const override {
139 return Teuchos::rcp(new PseudoBlockTFQMRSolMgr<ScalarType,MV,OP>);
140 }
142
144
145
147 return *problem_;
148 }
149
152 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const override;
153
156 Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const override { return params_; }
157
163 Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
164 return Teuchos::tuple(timerSolve_);
165 }
166
172 MagnitudeType achievedTol() const override {
173 return achievedTol_;
174 }
175
177 int getNumIters() const override {
178 return numIters_;
179 }
180
188 bool isLOADetected() const override { return false; }
190
192
193
195 void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) override { problem_ = problem; }
196
198 void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params ) override;
199
201
203
208 void reset( const ResetType type ) override { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
210
212
213
231 ReturnType solve() override;
232
234
236
238 std::string description() const override;
239
241 private:
242
243 // Method for checking current status test against defined linear problem.
244 bool checkStatusTest();
245
246 // Linear problem.
247 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
248
249 // Output manager.
250 Teuchos::RCP<OutputManager<ScalarType> > printer_;
251 Teuchos::RCP<std::ostream> outputStream_;
252
253 // Status test.
254 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
255 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
256 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
257 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > expConvTest_, impConvTest_;
258 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
259
260 // Current parameter list.
261 Teuchos::RCP<Teuchos::ParameterList> params_;
262
263 // Default solver values.
264 static constexpr int maxIters_default_ = 1000;
265 static constexpr bool expResTest_default_ = false;
266 static constexpr int verbosity_default_ = Belos::Errors;
267 static constexpr int outputStyle_default_ = Belos::General;
268 static constexpr int outputFreq_default_ = -1;
269 static constexpr int defQuorum_default_ = 1;
270 static constexpr const char * impResScale_default_ = "Norm of Preconditioned Initial Residual";
271 static constexpr const char * expResScale_default_ = "Norm of Initial Residual";
272 static constexpr const char * label_default_ = "Belos";
273
274 // Current solver values.
280
281 // Timers.
282 std::string label_;
283 Teuchos::RCP<Teuchos::Time> timerSolve_;
284
285 // Internal state variables.
287 };
288
289
290// Empty Constructor
291template<class ScalarType, class MV, class OP>
293 outputStream_(Teuchos::rcpFromRef(std::cout)),
294 convtol_(DefaultSolverParameters::convTol),
295 impTolScale_(DefaultSolverParameters::impTolScale),
296 achievedTol_(Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
297 maxIters_(maxIters_default_),
298 numIters_(0),
299 verbosity_(verbosity_default_),
300 outputStyle_(outputStyle_default_),
301 outputFreq_(outputFreq_default_),
302 defQuorum_(defQuorum_default_),
303 expResTest_(expResTest_default_),
304 impResScale_(impResScale_default_),
305 expResScale_(expResScale_default_),
306 label_(label_default_),
307 isSet_(false),
308 isSTSet_(false)
309{}
310
311
312// Basic Constructor
313template<class ScalarType, class MV, class OP>
315 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
316 const Teuchos::RCP<Teuchos::ParameterList> &pl ) :
317 problem_(problem),
318 outputStream_(Teuchos::rcpFromRef(std::cout)),
319 convtol_(DefaultSolverParameters::convTol),
320 impTolScale_(DefaultSolverParameters::impTolScale),
321 achievedTol_(Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
322 maxIters_(maxIters_default_),
323 numIters_(0),
324 verbosity_(verbosity_default_),
325 outputStyle_(outputStyle_default_),
326 outputFreq_(outputFreq_default_),
327 defQuorum_(defQuorum_default_),
328 expResTest_(expResTest_default_),
329 impResScale_(impResScale_default_),
330 expResScale_(expResScale_default_),
331 label_(label_default_),
332 isSet_(false),
333 isSTSet_(false)
334{
335 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
336
337 // If the parameter list pointer is null, then set the current parameters to the default parameter list.
338 if ( !is_null(pl) ) {
339 setParameters( pl );
340 }
341}
342
343template<class ScalarType, class MV, class OP>
344void PseudoBlockTFQMRSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
345{
346 // Create the internal parameter list if ones doesn't already exist.
347 if (params_ == Teuchos::null) {
348 params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
349 }
350 else {
351 params->validateParameters(*getValidParameters());
352 }
353
354 // Check for maximum number of iterations
355 if (params->isParameter("Maximum Iterations")) {
356 maxIters_ = params->get("Maximum Iterations",maxIters_default_);
357
358 // Update parameter in our list and in status test.
359 params_->set("Maximum Iterations", maxIters_);
360 if (maxIterTest_!=Teuchos::null)
361 maxIterTest_->setMaxIters( maxIters_ );
362 }
363
364 // Check to see if the timer label changed.
365 if (params->isParameter("Timer Label")) {
366 std::string tempLabel = params->get("Timer Label", label_default_);
367
368 // Update parameter in our list and solver timer
369 if (tempLabel != label_) {
370 label_ = tempLabel;
371 params_->set("Timer Label", label_);
372 std::string solveLabel = label_ + ": PseudoBlockTFQMRSolMgr total solve time";
373#ifdef BELOS_TEUCHOS_TIME_MONITOR
374 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
375#endif
376 }
377 }
378
379 // Check for a change in verbosity level
380 if (params->isParameter("Verbosity")) {
381 if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
382 verbosity_ = params->get("Verbosity", verbosity_default_);
383 } else {
384 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
385 }
386
387 // Update parameter in our list.
388 params_->set("Verbosity", verbosity_);
389 if (printer_ != Teuchos::null)
390 printer_->setVerbosity(verbosity_);
391 }
392
393 // Check for a change in output style
394 if (params->isParameter("Output Style")) {
395 if (Teuchos::isParameterType<int>(*params,"Output Style")) {
396 outputStyle_ = params->get("Output Style", outputStyle_default_);
397 } else {
398 outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
399 }
400
401 // Reconstruct the convergence test if the explicit residual test is not being used.
402 params_->set("Output Style", outputStyle_);
403 isSTSet_ = false;
404 }
405
406 // output stream
407 if (params->isParameter("Output Stream")) {
408 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
409
410 // Update parameter in our list.
411 params_->set("Output Stream", outputStream_);
412 if (printer_ != Teuchos::null)
413 printer_->setOStream( outputStream_ );
414 }
415
416 // frequency level
417 if (verbosity_ & Belos::StatusTestDetails) {
418 if (params->isParameter("Output Frequency")) {
419 outputFreq_ = params->get("Output Frequency", outputFreq_default_);
420 }
421
422 // Update parameter in out list and output status test.
423 params_->set("Output Frequency", outputFreq_);
424 if (outputTest_ != Teuchos::null)
425 outputTest_->setOutputFrequency( outputFreq_ );
426 }
427
428 // Create output manager if we need to.
429 if (printer_ == Teuchos::null) {
430 printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
431 }
432
433 // Check for convergence tolerance
434 if (params->isParameter("Convergence Tolerance")) {
435 if (params->isType<MagnitudeType> ("Convergence Tolerance")) {
436 convtol_ = params->get ("Convergence Tolerance",
438 }
439 else {
440 convtol_ = params->get ("Convergence Tolerance", DefaultSolverParameters::convTol);
441 }
442
443 // Update parameter in our list and residual tests.
444 params_->set("Convergence Tolerance", convtol_);
445 isSTSet_ = false;
446 }
447
448 if (params->isParameter("Implicit Tolerance Scale Factor")) {
449 if (params->isType<MagnitudeType> ("Implicit Tolerance Scale Factor")) {
450 impTolScale_ = params->get ("Implicit Tolerance Scale Factor",
452
453 }
454 else {
455 impTolScale_ = params->get ("Implicit Tolerance Scale Factor",
457 }
458
459 // Update parameter in our list.
460 params_->set("Implicit Tolerance Scale Factor", impTolScale_);
461 isSTSet_ = false;
462 }
463
464 if (params->isParameter("Implicit Residual Scaling")) {
465 std::string tempImpResScale = Teuchos::getParameter<std::string>( *params, "Implicit Residual Scaling" );
466
467 // Only update the scaling if it's different.
468 if (impResScale_ != tempImpResScale) {
469 impResScale_ = tempImpResScale;
470
471 // Update parameter in our list.
472 params_->set("Implicit Residual Scaling", impResScale_);
473 isSTSet_ = false;
474 }
475 }
476
477 if (params->isParameter("Explicit Residual Scaling")) {
478 std::string tempExpResScale = Teuchos::getParameter<std::string>( *params, "Explicit Residual Scaling" );
479
480 // Only update the scaling if it's different.
481 if (expResScale_ != tempExpResScale) {
482 expResScale_ = tempExpResScale;
483
484 // Update parameter in our list.
485 params_->set("Explicit Residual Scaling", expResScale_);
486 isSTSet_ = false;
487 }
488 }
489
490 if (params->isParameter("Explicit Residual Test")) {
491 expResTest_ = Teuchos::getParameter<bool>( *params,"Explicit Residual Test" );
492
493 // Reconstruct the convergence test if the explicit residual test is not being used.
494 params_->set("Explicit Residual Test", expResTest_);
495 if (expConvTest_ == Teuchos::null) {
496 isSTSet_ = false;
497 }
498 }
499
500 // Get the deflation quorum, or number of converged systems before deflation is allowed
501 if (params->isParameter("Deflation Quorum")) {
502 defQuorum_ = params->get("Deflation Quorum", defQuorum_);
503 params_->set ("Deflation Quorum", defQuorum_);
504 if (! impConvTest_.is_null ()) {
505 impConvTest_->setQuorum (defQuorum_);
506 }
507 if (! expConvTest_.is_null ()) {
508 expConvTest_->setQuorum (defQuorum_);
509 }
510 }
511
512 // Create the timer if we need to.
513 if (timerSolve_ == Teuchos::null) {
514 std::string solveLabel = label_ + ": PseudoBlockTFQMRSolMgr total solve time";
515#ifdef BELOS_TEUCHOS_TIME_MONITOR
516 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
517#endif
518 }
519
520 // Inform the solver manager that the current parameters were set.
521 isSet_ = true;
522}
523
524
525// Check the status test versus the defined linear problem
526template<class ScalarType, class MV, class OP>
528
529 typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
530 typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestGenResNorm_t;
531
532 // Basic test checks maximum iterations and native residual.
533 maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
534
535 if (expResTest_) {
536
537 // Implicit residual test, using the native residual to determine if convergence was achieved.
538 Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
539 Teuchos::rcp( new StatusTestGenResNorm_t( impTolScale_*convtol_, defQuorum_ ) );
540 tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
541 impConvTest_ = tmpImpConvTest;
542
543 // Explicit residual test once the native residual is below the tolerance
544 Teuchos::RCP<StatusTestGenResNorm_t> tmpExpConvTest =
545 Teuchos::rcp( new StatusTestGenResNorm_t( convtol_, defQuorum_ ) );
546 tmpExpConvTest->defineResForm( StatusTestGenResNorm_t::Explicit, Belos::TwoNorm );
547 tmpExpConvTest->defineScaleForm( convertStringToScaleType(expResScale_), Belos::TwoNorm );
548 expConvTest_ = tmpExpConvTest;
549
550 // The convergence test is a combination of the "cheap" implicit test and explicit test.
551 convTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::SEQ, impConvTest_, expConvTest_ ) );
552 }
553 else {
554
555 // Implicit residual test, using the native residual to determine if convergence was achieved.
556 Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
557 Teuchos::rcp( new StatusTestGenResNorm_t( convtol_, defQuorum_ ) );
558 tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
559 impConvTest_ = tmpImpConvTest;
560
561 // Set the explicit and total convergence test to this implicit test that checks for accuracy loss.
562 expConvTest_ = impConvTest_;
563 convTest_ = impConvTest_;
564 }
565 sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
566
567 // Create the status test output class.
568 // This class manages and formats the output from the status test.
569 StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
570 outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
571
572 // Set the solver string for the output test
573 std::string solverDesc = " Pseudo Block TFQMR ";
574 outputTest_->setSolverDesc( solverDesc );
575
576
577 // The status test is now set.
578 isSTSet_ = true;
579
580 return false;
581}
582
583
584template<class ScalarType, class MV, class OP>
585Teuchos::RCP<const Teuchos::ParameterList>
587{
588 static Teuchos::RCP<const Teuchos::ParameterList> validPL;
589
590 // Set all the valid parameters and their default values.
591 if(is_null(validPL)) {
592 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
593
594 // The static_cast is to resolve an issue with older clang versions which
595 // would cause the constexpr to link fail. With c++17 the problem is resolved.
596 pl->set("Convergence Tolerance", static_cast<MagnitudeType>(DefaultSolverParameters::convTol),
597 "The relative residual tolerance that needs to be achieved by the\n"
598 "iterative solver in order for the linear system to be declared converged.");
599 pl->set("Implicit Tolerance Scale Factor", static_cast<MagnitudeType>(DefaultSolverParameters::impTolScale),
600 "The scale factor used by the implicit residual test when explicit residual\n"
601 "testing is used. May enable faster convergence when TFQMR bound is too loose.");
602 pl->set("Maximum Iterations", static_cast<int>(maxIters_default_),
603 "The maximum number of block iterations allowed for each\n"
604 "set of RHS solved.");
605 pl->set("Verbosity", static_cast<int>(verbosity_default_),
606 "What type(s) of solver information should be outputted\n"
607 "to the output stream.");
608 pl->set("Output Style", static_cast<int>(outputStyle_default_),
609 "What style is used for the solver information outputted\n"
610 "to the output stream.");
611 pl->set("Output Frequency", static_cast<int>(outputFreq_default_),
612 "How often convergence information should be outputted\n"
613 "to the output stream.");
614 pl->set("Deflation Quorum", static_cast<int>(defQuorum_default_),
615 "The number of linear systems that need to converge before they are deflated.");
616 pl->set("Output Stream", Teuchos::rcpFromRef(std::cout),
617 "A reference-counted pointer to the output stream where all\n"
618 "solver output is sent.");
619 pl->set("Explicit Residual Test", static_cast<bool>(expResTest_default_),
620 "Whether the explicitly computed residual should be used in the convergence test.");
621 pl->set("Implicit Residual Scaling", static_cast<const char *>(impResScale_default_),
622 "The type of scaling used in the implicit residual convergence test.");
623 pl->set("Explicit Residual Scaling", static_cast<const char *>(expResScale_default_),
624 "The type of scaling used in the explicit residual convergence test.");
625 pl->set("Timer Label", static_cast<const char *>(label_default_),
626 "The string to use as a prefix for the timer labels.");
627 validPL = pl;
628 }
629 return validPL;
630}
631
632
633// solve()
634template<class ScalarType, class MV, class OP>
636
637 // Set the current parameters if they were not set before.
638 // NOTE: This may occur if the user generated the solver manager with the default constructor and
639 // then didn't set any parameters using setParameters().
640 if (!isSet_) {
641 setParameters(Teuchos::parameterList(*getValidParameters()));
642 }
643
644 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,PseudoBlockTFQMRSolMgrLinearProblemFailure,
645 "Belos::PseudoBlockTFQMRSolMgr::solve(): Linear problem is not a valid object.");
646
647 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),PseudoBlockTFQMRSolMgrLinearProblemFailure,
648 "Belos::PseudoBlockTFQMRSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
649
650 if (!isSTSet_) {
651 TEUCHOS_TEST_FOR_EXCEPTION( checkStatusTest(),PseudoBlockTFQMRSolMgrLinearProblemFailure,
652 "Belos::PseudoBlockTFQMRSolMgr::solve(): Linear problem and requested status tests are incompatible.");
653 }
654
655 // Create indices for the linear systems to be solved.
656 int startPtr = 0;
657 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
658 int numCurrRHS = numRHS2Solve;
659
660 std::vector<int> currIdx( numRHS2Solve );
661 for (int i=0; i<numRHS2Solve; ++i) {
662 currIdx[i] = startPtr+i;
663 }
664
665 // Inform the linear problem of the current linear system to solve.
666 problem_->setLSIndex( currIdx );
667
669 // Parameter list
670 Teuchos::ParameterList plist;
671
672 // Reset the status test.
673 outputTest_->reset();
674
675 // Assume convergence is achieved, then let any failed convergence set this to false.
676 bool isConverged = true;
677
679 // TFQMR solver
680
681 Teuchos::RCP<PseudoBlockTFQMRIter<ScalarType,MV,OP> > block_tfqmr_iter =
682 Teuchos::rcp( new PseudoBlockTFQMRIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,plist) );
683
684 // Enter solve() iterations
685 {
686#ifdef BELOS_TEUCHOS_TIME_MONITOR
687 Teuchos::TimeMonitor slvtimer(*timerSolve_);
688#endif
689
690 while ( numRHS2Solve > 0 ) {
691 //
692 // Reset the active / converged vectors from this block
693 std::vector<int> convRHSIdx;
694 std::vector<int> currRHSIdx( currIdx );
695 currRHSIdx.resize(numCurrRHS);
696
697 // Reset the number of iterations.
698 block_tfqmr_iter->resetNumIters();
699
700 // Reset the number of calls that the status test output knows about.
701 outputTest_->resetNumCalls();
702
703 // Get the current residual for this block of linear systems.
704 Teuchos::RCP<MV> R_0 = MVT::CloneViewNonConst( *(Teuchos::rcp_const_cast<MV>(problem_->getInitPrecResVec())), currIdx );
705
706 // Set the new state and initialize the solver.
708 newstate.Rtilde = R_0;
709 block_tfqmr_iter->initializeTFQMR(newstate);
710
711 while(1) {
712
713 // tell block_tfqmr_iter to iterate
714 try {
715 block_tfqmr_iter->iterate();
716
718 //
719 // check convergence first
720 //
722 if ( convTest_->getStatus() == Passed ) {
723
724 // Figure out which linear systems converged.
725 std::vector<int> convIdx = expConvTest_->convIndices();
726
727 // If the number of converged linear systems is equal to the
728 // number of current linear systems, then we are done with this block.
729 if (convIdx.size() == currRHSIdx.size())
730 break; // break from while(1){block_tfqmr_iter->iterate()}
731
732 // Update the current solution with the update computed by the iteration object.
733 problem_->updateSolution( block_tfqmr_iter->getCurrentUpdate(), true );
734
735 // Inform the linear problem that we are finished with this current linear system.
736 problem_->setCurrLS();
737
738 // Reset currRHSIdx to have the right-hand sides that are left to converge for this block.
739 int have = 0;
740 std::vector<int> unconvIdx( currRHSIdx.size() );
741 for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
742 bool found = false;
743 for (unsigned int j=0; j<convIdx.size(); ++j) {
744 if (currRHSIdx[i] == convIdx[j]) {
745 found = true;
746 break;
747 }
748 }
749 if (!found) {
750 unconvIdx[have] = i;
751 currRHSIdx[have++] = currRHSIdx[i];
752 }
753 }
754 unconvIdx.resize(have);
755 currRHSIdx.resize(have);
756
757 // Set the remaining indices after deflation.
758 problem_->setLSIndex( currRHSIdx );
759
760 // Get the current residual vector.
761 // Set the new state and initialize the solver.
762 PseudoBlockTFQMRIterState<ScalarType,MV> currentState = block_tfqmr_iter->getState();
763
764 // Set the new state and initialize the solver.
766
767 // Copy over the vectors.
768 defstate.Rtilde = MVT::CloneView( *currentState.Rtilde, unconvIdx);
769 defstate.U = MVT::CloneView( *currentState.U, unconvIdx );
770 defstate.AU = MVT::CloneView( *currentState.AU, unconvIdx );
771 defstate.V = MVT::CloneView( *currentState.V, unconvIdx );
772 defstate.W = MVT::CloneView( *currentState.W, unconvIdx );
773 defstate.D = MVT::CloneView( *currentState.D, unconvIdx );
774
775 // Copy over the scalars.
776 for (std::vector<int>::iterator uIter = unconvIdx.begin(); uIter != unconvIdx.end(); uIter++)
777 {
778 defstate.alpha.push_back( currentState.alpha[ *uIter ] );
779 defstate.eta.push_back( currentState.eta[ *uIter ] );
780 defstate.rho.push_back( currentState.rho[ *uIter ] );
781 defstate.tau.push_back( currentState.tau[ *uIter ] );
782 defstate.theta.push_back( currentState.theta[ *uIter ] );
783 }
784
785 block_tfqmr_iter->initializeTFQMR(defstate);
786 }
788 //
789 // check for maximum iterations
790 //
792 else if ( maxIterTest_->getStatus() == Passed ) {
793 // we don't have convergence
794 isConverged = false;
795 break; // break from while(1){block_tfqmr_iter->iterate()}
796 }
797
799 //
800 // we returned from iterate(), but none of our status tests Passed.
801 // something is wrong, and it is probably our fault.
802 //
804
805 else {
806 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
807 "Belos::PseudoBlockTFQMRSolMgr::solve(): Invalid return from PseudoBlockTFQMRIter::iterate().");
808 }
809 }
810 catch (const std::exception &e) {
811 printer_->stream(Errors) << "Error! Caught std::exception in PseudoBlockTFQMRIter::iterate() at iteration "
812 << block_tfqmr_iter->getNumIters() << std::endl
813 << e.what() << std::endl;
814 throw;
815 }
816 }
817
818 // Update the current solution with the update computed by the iteration object.
819 problem_->updateSolution( block_tfqmr_iter->getCurrentUpdate(), true );
820
821 // Inform the linear problem that we are finished with this block linear system.
822 problem_->setCurrLS();
823
824 // Update indices for the linear systems to be solved.
825 startPtr += numCurrRHS;
826 numRHS2Solve -= numCurrRHS;
827 if ( numRHS2Solve > 0 ) {
828 numCurrRHS = numRHS2Solve;
829 currIdx.resize( numCurrRHS );
830 for (int i=0; i<numCurrRHS; ++i)
831 { currIdx[i] = startPtr+i; }
832
833 // Adapt the status test quorum if we need to.
834 if (defQuorum_ > numCurrRHS) {
835 if (impConvTest_ != Teuchos::null)
836 impConvTest_->setQuorum( numCurrRHS );
837 if (expConvTest_ != Teuchos::null)
838 expConvTest_->setQuorum( numCurrRHS );
839 }
840
841 // Set the next indices.
842 problem_->setLSIndex( currIdx );
843 }
844 else {
845 currIdx.resize( numRHS2Solve );
846 }
847
848 }// while ( numRHS2Solve > 0 )
849
850 }
851
852 // print final summary
853 sTest_->print( printer_->stream(FinalSummary) );
854
855 // print timing information
856#ifdef BELOS_TEUCHOS_TIME_MONITOR
857 // Calling summarize() can be expensive, so don't call unless the
858 // user wants to print out timing details. summarize() will do all
859 // the work even if it's passed a "black hole" output stream.
860 if (verbosity_ & TimingDetails)
861 Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
862#endif
863
864 // get iteration information for this solve
865 numIters_ = maxIterTest_->getNumIters();
866
867 // Save the convergence test value ("achieved tolerance") for this
868 // solve. For this solver, convTest_ may either be a single
869 // (implicit) residual norm test, or a combination of two residual
870 // norm tests. In the latter case, the master convergence test
871 // convTest_ is a SEQ combo of the implicit resp. explicit tests.
872 // If the implicit test never passes, then the explicit test won't
873 // ever be executed. This manifests as
874 // expConvTest_->getTestValue()->size() < 1. We deal with this case
875 // by using the values returned by impConvTest_->getTestValue().
876 {
877 // We'll fetch the vector of residual norms one way or the other.
878 const std::vector<MagnitudeType>* pTestValues = NULL;
879 if (expResTest_) {
880 pTestValues = expConvTest_->getTestValue();
881 if (pTestValues == NULL || pTestValues->size() < 1) {
882 pTestValues = impConvTest_->getTestValue();
883 }
884 }
885 else {
886 // Only the implicit residual norm test is being used.
887 pTestValues = impConvTest_->getTestValue();
888 }
889 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
890 "Belos::PseudoBlockTFQMRSolMgr::solve(): The implicit convergence test's "
891 "getTestValue() method returned NULL. Please report this bug to the "
892 "Belos developers.");
893 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
894 "Belos::TMQMRSolMgr::solve(): The implicit convergence test's "
895 "getTestValue() method returned a vector of length zero. Please report "
896 "this bug to the Belos developers.");
897
898 // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
899 // achieved tolerances for all vectors in the current solve(), or
900 // just for the vectors from the last deflation?
901 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
902 }
903
904 if (!isConverged) {
905 return Unconverged; // return from PseudoBlockTFQMRSolMgr::solve()
906 }
907 return Converged; // return from PseudoBlockTFQMRSolMgr::solve()
908}
909
910// This method requires the solver manager to return a std::string that describes itself.
911template<class ScalarType, class MV, class OP>
913{
914 std::ostringstream oss;
915 oss << "Belos::PseudoBlockTFQMRSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
916 oss << "{}";
917 return oss.str();
918}
919
920} // end Belos namespace
921
922#endif /* BELOS_PSEUDO_BLOCK_TFQMR_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 for generating iterations with the preconditioned tranpose-free QMR (TFQMR) meth...
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
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 preconditioned transpose-free QMR algorithm for solving non-Hermitian linea...
PseudoBlockTFQMRSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i....
The Belos::PseudoBlockTFQMRSolMgr provides a powerful and fully-featured solver manager over the pseu...
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
static constexpr const char * label_default_
Teuchos::RCP< Teuchos::Time > timerSolve_
static constexpr const char * impResScale_default_
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > expConvTest_
Teuchos::RCP< std::ostream > outputStream_
bool isLOADetected() const override
Whether loss of accuracy was detected during the last solve() invocation.
void reset(const ResetType type) override
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
static constexpr const char * expResScale_default_
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
Teuchos::ScalarTraits< MagnitudeType > MT
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override
Get a parameter list containing the valid parameters for this object.
Teuchos::RCP< OutputManager< ScalarType > > printer_
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > impConvTest_
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Return a reference to the linear problem being solved by this solver manager.
OperatorTraits< ScalarType, MV, OP > OPT
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > convTest_
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
std::string description() const override
Method to return description of the pseudo-block TFQMR solver manager.
Teuchos::ScalarTraits< ScalarType > SCT
PseudoBlockTFQMRSolMgr()
Empty constructor for PseudoBlockTFQMRSolMgr. This constructor takes no arguments and sets the defaul...
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > &params) override
Set the parameters the solver manager should use to solve the linear problem.
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
MultiVecTraits< ScalarType, MV > MVT
int getNumIters() const override
Get the iteration count for the most recent call to solve().
Teuchos::RCP< Teuchos::ParameterList > params_
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
ReturnType solve() override
This method performs possibly repeated calls to the underlying linear solver's iterate() routine unti...
The Belos::SolverManager is a templated virtual base class that defines the basic interface that any ...
A class for extending the status testing capabilities of Belos via logical combinations.
An implementation of StatusTestResNorm using a family of residual norms.
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
@ 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
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 impTolScale
"Implicit Tolerance Scale Factor"
Definition: BelosTypes.hpp:305
static const double convTol
Default convergence tolerance.
Definition: BelosTypes.hpp:293
Structure to contain pointers to PseudoBlockTFQMRIter state variables.
Teuchos::RCP< const MV > W
The current residual basis.