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

Generated for Belos by doxygen 1.8.14