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