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 
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:
109  typedef Teuchos::ScalarTraits<ScalarType> SCT;
110  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
111  typedef Teuchos::ScalarTraits<MagnitudeType> MT;
112 
113  public:
114 
116 
117 
124 
141  PseudoBlockTFQMRSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
142  const Teuchos::RCP<Teuchos::ParameterList> &pl );
143 
147 
149 
150 
152  return *problem_;
153  }
154 
157  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
158 
161  Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
162 
168  Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
169  return Teuchos::tuple(timerSolve_);
170  }
171 
177  MagnitudeType achievedTol() const {
178  return achievedTol_;
179  }
180 
182  int getNumIters() const {
183  return numIters_;
184  }
185 
193  bool isLOADetected() const { return false; }
194 
196 
198 
199 
201  void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { problem_ = problem; }
202 
204  void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
205 
207 
209 
210 
214  void reset( const ResetType type ) { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
216 
218 
219 
237  ReturnType solve();
238 
240 
243 
245  std::string description() const;
246 
248 
249  private:
250 
251  // Method for checking current status test against defined linear problem.
252  bool checkStatusTest();
253 
254  // Linear problem.
255  Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
256 
257  // Output manager.
258  Teuchos::RCP<OutputManager<ScalarType> > printer_;
259  Teuchos::RCP<std::ostream> outputStream_;
260 
261  // Status test.
262  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
263  Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
264  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
265  Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > expConvTest_, impConvTest_;
266  Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
267 
268  // Current parameter list.
269  Teuchos::RCP<Teuchos::ParameterList> params_;
270 
271  // Default solver values.
272  static const MagnitudeType convtol_default_;
273  static const int maxIters_default_;
274  static const bool expResTest_default_;
275  static const int verbosity_default_;
276  static const int outputStyle_default_;
277  static const int outputFreq_default_;
278  static const int defQuorum_default_;
279  static const std::string impResScale_default_;
280  static const std::string expResScale_default_;
281  static const std::string label_default_;
282  static const Teuchos::RCP<std::ostream> outputStream_default_;
283 
284  // Current solver values.
285  MagnitudeType convtol_, achievedTol_;
286  int maxIters_, numIters_;
287  int verbosity_, outputStyle_, outputFreq_, defQuorum_;
288  bool expResTest_;
289  std::string impResScale_, expResScale_;
290 
291  // Timers.
292  std::string label_;
293  Teuchos::RCP<Teuchos::Time> timerSolve_;
294 
295  // Internal state variables.
296  bool isSet_, isSTSet_;
297  };
298 
299 
300 // Default solver values.
301 template<class ScalarType, class MV, class OP>
302 const typename PseudoBlockTFQMRSolMgr<ScalarType,MV,OP>::MagnitudeType PseudoBlockTFQMRSolMgr<ScalarType,MV,OP>::convtol_default_ = 1e-8;
303 
304 template<class ScalarType, class MV, class OP>
306 
307 template<class ScalarType, class MV, class OP>
309 
310 template<class ScalarType, class MV, class OP>
312 
313 template<class ScalarType, class MV, class OP>
315 
316 template<class ScalarType, class MV, class OP>
318 
319 template<class ScalarType, class MV, class OP>
321 
322 template<class ScalarType, class MV, class OP>
323 const std::string PseudoBlockTFQMRSolMgr<ScalarType,MV,OP>::impResScale_default_ = "Norm of Preconditioned Initial Residual";
324 
325 template<class ScalarType, class MV, class OP>
326 const std::string PseudoBlockTFQMRSolMgr<ScalarType,MV,OP>::expResScale_default_ = "Norm of Initial Residual";
327 
328 template<class ScalarType, class MV, class OP>
330 
331 template<class ScalarType, class MV, class OP>
332 const Teuchos::RCP<std::ostream> PseudoBlockTFQMRSolMgr<ScalarType,MV,OP>::outputStream_default_ = Teuchos::rcp(&std::cout,false);
333 
334 
335 // Empty Constructor
336 template<class ScalarType, class MV, class OP>
338  outputStream_(outputStream_default_),
339  convtol_(convtol_default_),
340  achievedTol_(Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
341  maxIters_(maxIters_default_),
342  numIters_(0),
343  verbosity_(verbosity_default_),
344  outputStyle_(outputStyle_default_),
345  outputFreq_(outputFreq_default_),
346  defQuorum_(defQuorum_default_),
347  expResTest_(expResTest_default_),
348  impResScale_(impResScale_default_),
349  expResScale_(expResScale_default_),
350  label_(label_default_),
351  isSet_(false),
352  isSTSet_(false)
353 {}
354 
355 
356 // Basic Constructor
357 template<class ScalarType, class MV, class OP>
359  const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
360  const Teuchos::RCP<Teuchos::ParameterList> &pl ) :
361  problem_(problem),
362  outputStream_(outputStream_default_),
363  convtol_(convtol_default_),
364  achievedTol_(Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
365  maxIters_(maxIters_default_),
366  numIters_(0),
367  verbosity_(verbosity_default_),
368  outputStyle_(outputStyle_default_),
369  outputFreq_(outputFreq_default_),
370  defQuorum_(defQuorum_default_),
371  expResTest_(expResTest_default_),
372  impResScale_(impResScale_default_),
373  expResScale_(expResScale_default_),
374  label_(label_default_),
375  isSet_(false),
376  isSTSet_(false)
377 {
378  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
379 
380  // If the parameter list pointer is null, then set the current parameters to the default parameter list.
381  if ( !is_null(pl) ) {
382  setParameters( pl );
383  }
384 }
385 
386 template<class ScalarType, class MV, class OP>
387 void PseudoBlockTFQMRSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
388 {
389  // Create the internal parameter list if ones doesn't already exist.
390  if (params_ == Teuchos::null) {
391  params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
392  }
393  else {
394  params->validateParameters(*getValidParameters());
395  }
396 
397  // Check for maximum number of iterations
398  if (params->isParameter("Maximum Iterations")) {
399  maxIters_ = params->get("Maximum Iterations",maxIters_default_);
400 
401  // Update parameter in our list and in status test.
402  params_->set("Maximum Iterations", maxIters_);
403  if (maxIterTest_!=Teuchos::null)
404  maxIterTest_->setMaxIters( maxIters_ );
405  }
406 
407  // Check to see if the timer label changed.
408  if (params->isParameter("Timer Label")) {
409  std::string tempLabel = params->get("Timer Label", label_default_);
410 
411  // Update parameter in our list and solver timer
412  if (tempLabel != label_) {
413  label_ = tempLabel;
414  params_->set("Timer Label", label_);
415  std::string solveLabel = label_ + ": PseudoBlockTFQMRSolMgr total solve time";
416 #ifdef BELOS_TEUCHOS_TIME_MONITOR
417  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
418 #endif
419  }
420  }
421 
422  // Check for a change in verbosity level
423  if (params->isParameter("Verbosity")) {
424  if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
425  verbosity_ = params->get("Verbosity", verbosity_default_);
426  } else {
427  verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
428  }
429 
430  // Update parameter in our list.
431  params_->set("Verbosity", verbosity_);
432  if (printer_ != Teuchos::null)
433  printer_->setVerbosity(verbosity_);
434  }
435 
436  // Check for a change in output style
437  if (params->isParameter("Output Style")) {
438  if (Teuchos::isParameterType<int>(*params,"Output Style")) {
439  outputStyle_ = params->get("Output Style", outputStyle_default_);
440  } else {
441  outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
442  }
443 
444  // Reconstruct the convergence test if the explicit residual test is not being used.
445  params_->set("Output Style", outputStyle_);
446  isSTSet_ = false;
447  }
448 
449  // output stream
450  if (params->isParameter("Output Stream")) {
451  outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
452 
453  // Update parameter in our list.
454  params_->set("Output Stream", outputStream_);
455  if (printer_ != Teuchos::null)
456  printer_->setOStream( outputStream_ );
457  }
458 
459  // frequency level
460  if (verbosity_ & Belos::StatusTestDetails) {
461  if (params->isParameter("Output Frequency")) {
462  outputFreq_ = params->get("Output Frequency", outputFreq_default_);
463  }
464 
465  // Update parameter in out list and output status test.
466  params_->set("Output Frequency", outputFreq_);
467  if (outputTest_ != Teuchos::null)
468  outputTest_->setOutputFrequency( outputFreq_ );
469  }
470 
471  // Create output manager if we need to.
472  if (printer_ == Teuchos::null) {
473  printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
474  }
475 
476  // Check for convergence tolerance
477  if (params->isParameter("Convergence Tolerance")) {
478  convtol_ = params->get("Convergence Tolerance",convtol_default_);
479 
480  // Update parameter in our list and residual tests.
481  params_->set("Convergence Tolerance", convtol_);
482  if (impConvTest_ != Teuchos::null)
483  impConvTest_->setTolerance( convtol_ );
484  if (expConvTest_ != Teuchos::null)
485  expConvTest_->setTolerance( convtol_ );
486  }
487 
488  // Check for a change in scaling, if so we need to build new residual tests.
489  if (params->isParameter("Implicit Residual Scaling")) {
490  std::string tempImpResScale = Teuchos::getParameter<std::string>( *params, "Implicit Residual Scaling" );
491 
492  // Only update the scaling if it's different.
493  if (impResScale_ != tempImpResScale) {
494  Belos::ScaleType impResScaleType = convertStringToScaleType( tempImpResScale );
495  impResScale_ = tempImpResScale;
496 
497  // Update parameter in our list and residual tests
498  params_->set("Implicit Residual Scaling", impResScale_);
499  if (impConvTest_ != Teuchos::null) {
500  try {
501  impConvTest_->defineScaleForm( impResScaleType, Belos::TwoNorm );
502  }
503  catch (std::exception& e) {
504  // Make sure the convergence test gets constructed again.
505  isSTSet_ = false;
506  }
507  }
508  }
509  }
510 
511  if (params->isParameter("Explicit Residual Scaling")) {
512  std::string tempExpResScale = Teuchos::getParameter<std::string>( *params, "Explicit Residual Scaling" );
513 
514  // Only update the scaling if it's different.
515  if (expResScale_ != tempExpResScale) {
516  Belos::ScaleType expResScaleType = convertStringToScaleType( tempExpResScale );
517  expResScale_ = tempExpResScale;
518 
519  // Update parameter in our list and residual tests
520  params_->set("Explicit Residual Scaling", expResScale_);
521  if (expConvTest_ != Teuchos::null) {
522  try {
523  expConvTest_->defineScaleForm( expResScaleType, Belos::TwoNorm );
524  }
525  catch (std::exception& e) {
526  // Make sure the convergence test gets constructed again.
527  isSTSet_ = false;
528  }
529  }
530  }
531  }
532 
533  if (params->isParameter("Explicit Residual Test")) {
534  expResTest_ = Teuchos::getParameter<bool>( *params,"Explicit Residual Test" );
535 
536  // Reconstruct the convergence test if the explicit residual test is not being used.
537  params_->set("Explicit Residual Test", expResTest_);
538  if (expConvTest_ == Teuchos::null) {
539  isSTSet_ = false;
540  }
541  }
542 
543  // Get the deflation quorum, or number of converged systems before deflation is allowed
544  if (params->isParameter("Deflation Quorum")) {
545  defQuorum_ = params->get("Deflation Quorum", defQuorum_);
546  params_->set ("Deflation Quorum", defQuorum_);
547  if (! impConvTest_.is_null ()) {
548  impConvTest_->setQuorum (defQuorum_);
549  }
550  if (! expConvTest_.is_null ()) {
551  expConvTest_->setQuorum (defQuorum_);
552  }
553  }
554 
555  // Create the timer if we need to.
556  if (timerSolve_ == Teuchos::null) {
557  std::string solveLabel = label_ + ": PseudoBlockTFQMRSolMgr total solve time";
558 #ifdef BELOS_TEUCHOS_TIME_MONITOR
559  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
560 #endif
561  }
562 
563  // Inform the solver manager that the current parameters were set.
564  isSet_ = true;
565 }
566 
567 
568 // Check the status test versus the defined linear problem
569 template<class ScalarType, class MV, class OP>
571 
572  typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
573  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestGenResNorm_t;
574 
575  // Basic test checks maximum iterations and native residual.
576  maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
577 
578  if (expResTest_) {
579 
580  // Implicit residual test, using the native residual to determine if convergence was achieved.
581  Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
582  Teuchos::rcp( new StatusTestGenResNorm_t( convtol_, defQuorum_ ) );
583  tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
584  impConvTest_ = tmpImpConvTest;
585 
586  // Explicit residual test once the native residual is below the tolerance
587  Teuchos::RCP<StatusTestGenResNorm_t> tmpExpConvTest =
588  Teuchos::rcp( new StatusTestGenResNorm_t( convtol_, defQuorum_ ) );
589  tmpExpConvTest->defineResForm( StatusTestGenResNorm_t::Explicit, Belos::TwoNorm );
590  tmpExpConvTest->defineScaleForm( convertStringToScaleType(expResScale_), Belos::TwoNorm );
591  expConvTest_ = tmpExpConvTest;
592 
593  // The convergence test is a combination of the "cheap" implicit test and explicit test.
594  convTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::SEQ, impConvTest_, expConvTest_ ) );
595  }
596  else {
597 
598  // Implicit residual test, using the native residual to determine if convergence was achieved.
599  Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
600  Teuchos::rcp( new StatusTestGenResNorm_t( convtol_, defQuorum_ ) );
601  tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
602  impConvTest_ = tmpImpConvTest;
603 
604  // Set the explicit and total convergence test to this implicit test that checks for accuracy loss.
605  expConvTest_ = impConvTest_;
606  convTest_ = impConvTest_;
607  }
608  sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
609 
610  // Create the status test output class.
611  // This class manages and formats the output from the status test.
612  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
613  outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
614 
615  // Set the solver string for the output test
616  std::string solverDesc = " Pseudo Block TFQMR ";
617  outputTest_->setSolverDesc( solverDesc );
618 
619 
620  // The status test is now set.
621  isSTSet_ = true;
622 
623  return false;
624 }
625 
626 
627 template<class ScalarType, class MV, class OP>
628 Teuchos::RCP<const Teuchos::ParameterList>
630 {
631  static Teuchos::RCP<const Teuchos::ParameterList> validPL;
632 
633  // Set all the valid parameters and their default values.
634  if(is_null(validPL)) {
635  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
636  pl->set("Convergence Tolerance", convtol_default_,
637  "The relative residual tolerance that needs to be achieved by the\n"
638  "iterative solver in order for the linear system to be declared converged.");
639  pl->set("Maximum Iterations", maxIters_default_,
640  "The maximum number of block iterations allowed for each\n"
641  "set of RHS solved.");
642  pl->set("Verbosity", verbosity_default_,
643  "What type(s) of solver information should be outputted\n"
644  "to the output stream.");
645  pl->set("Output Style", outputStyle_default_,
646  "What style is used for the solver information outputted\n"
647  "to the output stream.");
648  pl->set("Output Frequency", outputFreq_default_,
649  "How often convergence information should be outputted\n"
650  "to the output stream.");
651  pl->set("Deflation Quorum", defQuorum_default_,
652  "The number of linear systems that need to converge before they are deflated.");
653  pl->set("Output Stream", outputStream_default_,
654  "A reference-counted pointer to the output stream where all\n"
655  "solver output is sent.");
656  pl->set("Explicit Residual Test", expResTest_default_,
657  "Whether the explicitly computed residual should be used in the convergence test.");
658  pl->set("Implicit Residual Scaling", impResScale_default_,
659  "The type of scaling used in the implicit residual convergence test.");
660  pl->set("Explicit Residual Scaling", expResScale_default_,
661  "The type of scaling used in the explicit residual convergence test.");
662  pl->set("Timer Label", label_default_,
663  "The string to use as a prefix for the timer labels.");
664  // pl->set("Restart Timers", restartTimers_);
665  validPL = pl;
666  }
667  return validPL;
668 }
669 
670 
671 // solve()
672 template<class ScalarType, class MV, class OP>
674 
675  // Set the current parameters if they were not set before.
676  // NOTE: This may occur if the user generated the solver manager with the default constructor and
677  // then didn't set any parameters using setParameters().
678  if (!isSet_) {
679  setParameters(Teuchos::parameterList(*getValidParameters()));
680  }
681 
682  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,PseudoBlockTFQMRSolMgrLinearProblemFailure,
683  "Belos::PseudoBlockTFQMRSolMgr::solve(): Linear problem is not a valid object.");
684 
685  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),PseudoBlockTFQMRSolMgrLinearProblemFailure,
686  "Belos::PseudoBlockTFQMRSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
687 
688  if (!isSTSet_) {
689  TEUCHOS_TEST_FOR_EXCEPTION( checkStatusTest(),PseudoBlockTFQMRSolMgrLinearProblemFailure,
690  "Belos::PseudoBlockTFQMRSolMgr::solve(): Linear problem and requested status tests are incompatible.");
691  }
692 
693  // Create indices for the linear systems to be solved.
694  int startPtr = 0;
695  int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
696  int numCurrRHS = numRHS2Solve;
697 
698  std::vector<int> currIdx( numRHS2Solve ), currIdx2( numRHS2Solve );
699  for (int i=0; i<numRHS2Solve; ++i) {
700  currIdx[i] = startPtr+i;
701  currIdx2[i]=i;
702  }
703 
704  // Inform the linear problem of the current linear system to solve.
705  problem_->setLSIndex( currIdx );
706 
708  // Parameter list
709  Teuchos::ParameterList plist;
710 
711  // Reset the status test.
712  outputTest_->reset();
713 
714  // Assume convergence is achieved, then let any failed convergence set this to false.
715  bool isConverged = true;
716 
718  // TFQMR solver
719 
720  Teuchos::RCP<PseudoBlockTFQMRIter<ScalarType,MV,OP> > block_tfqmr_iter =
721  Teuchos::rcp( new PseudoBlockTFQMRIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,plist) );
722 
723  // Enter solve() iterations
724  {
725 #ifdef BELOS_TEUCHOS_TIME_MONITOR
726  Teuchos::TimeMonitor slvtimer(*timerSolve_);
727 #endif
728 
729  while ( numRHS2Solve > 0 ) {
730  //
731  // Reset the active / converged vectors from this block
732  std::vector<int> convRHSIdx;
733  std::vector<int> currRHSIdx( currIdx );
734  currRHSIdx.resize(numCurrRHS);
735 
736  // Reset the number of iterations.
737  block_tfqmr_iter->resetNumIters();
738 
739  // Reset the number of calls that the status test output knows about.
740  outputTest_->resetNumCalls();
741 
742  // Get the current residual for this block of linear systems.
743  Teuchos::RCP<MV> R_0 = MVT::CloneViewNonConst( *(Teuchos::rcp_const_cast<MV>(problem_->getInitPrecResVec())), currIdx );
744 
745  // Set the new state and initialize the solver.
747  newstate.R = R_0;
748  block_tfqmr_iter->initializeTFQMR(newstate);
749 
750  while(1) {
751 
752  // tell block_tfqmr_iter to iterate
753  try {
754  block_tfqmr_iter->iterate();
755 
757  //
758  // check convergence first
759  //
761  if ( convTest_->getStatus() == Passed ) {
762 
763  // Figure out which linear systems converged.
764  std::vector<int> convIdx = expConvTest_->convIndices();
765 
766  // If the number of converged linear systems is equal to the
767  // number of current linear systems, then we are done with this block.
768  if (convIdx.size() == currRHSIdx.size())
769  break; // break from while(1){block_tfqmr_iter->iterate()}
770 
771  // Update the current solution with the update computed by the iteration object.
772  problem_->updateSolution( block_tfqmr_iter->getCurrentUpdate(), true );
773 
774  // Inform the linear problem that we are finished with this current linear system.
775  problem_->setCurrLS();
776 
777  // Reset currRHSIdx to have the right-hand sides that are left to converge for this block.
778  int have = 0;
779  std::vector<int> unconvIdx(currRHSIdx.size());
780  for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
781  bool found = false;
782  for (unsigned int j=0; j<convIdx.size(); ++j) {
783  if (currRHSIdx[i] == convIdx[j]) {
784  found = true;
785  break;
786  }
787  }
788  if (!found) {
789  currIdx2[have] = currIdx2[i];
790  currRHSIdx[have++] = currRHSIdx[i];
791  }
792  }
793  currRHSIdx.resize(have);
794  currIdx2.resize(have);
795 
796  // Set the remaining indices after deflation.
797  problem_->setLSIndex( currRHSIdx );
798 
799  // Get the current residual vector.
800  // Set the new state and initialize the solver.
801  R_0 = MVT::Clone( *(problem_->getInitPrecResVec()), currRHSIdx.size() );
802  problem_->computeCurrPrecResVec( &*R_0 );
803 
804  // Set the new state and initialize the solver.
806  defstate.R = R_0;
807  block_tfqmr_iter->initializeTFQMR(defstate);
808  }
810  //
811  // check for maximum iterations
812  //
814  else if ( maxIterTest_->getStatus() == Passed ) {
815  // we don't have convergence
816  isConverged = false;
817  break; // break from while(1){block_tfqmr_iter->iterate()}
818  }
819 
821  //
822  // we returned from iterate(), but none of our status tests Passed.
823  // something is wrong, and it is probably our fault.
824  //
826 
827  else {
828  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
829  "Belos::PseudoBlockTFQMRSolMgr::solve(): Invalid return from PseudoBlockTFQMRIter::iterate().");
830  }
831  }
832  catch (const std::exception &e) {
833  printer_->stream(Errors) << "Error! Caught std::exception in PseudoBlockTFQMRIter::iterate() at iteration "
834  << block_tfqmr_iter->getNumIters() << std::endl
835  << e.what() << std::endl;
836  throw;
837  }
838  }
839 
840  // Update the current solution with the update computed by the iteration object.
841  problem_->updateSolution( block_tfqmr_iter->getCurrentUpdate(), true );
842 
843  // Inform the linear problem that we are finished with this block linear system.
844  problem_->setCurrLS();
845 
846  // Update indices for the linear systems to be solved.
847  startPtr += numCurrRHS;
848  numRHS2Solve -= numCurrRHS;
849  if ( numRHS2Solve > 0 ) {
850  numCurrRHS = numRHS2Solve;
851  currIdx.resize( numCurrRHS );
852  currIdx2.resize( numCurrRHS );
853  for (int i=0; i<numCurrRHS; ++i)
854  { currIdx[i] = startPtr+i; currIdx2[i] = i; }
855 
856  // Adapt the status test quorum if we need to.
857  if (defQuorum_ > numCurrRHS) {
858  if (impConvTest_ != Teuchos::null)
859  impConvTest_->setQuorum( numCurrRHS );
860  if (expConvTest_ != Teuchos::null)
861  expConvTest_->setQuorum( numCurrRHS );
862  }
863 
864  // Set the next indices.
865  problem_->setLSIndex( currIdx );
866  }
867  else {
868  currIdx.resize( numRHS2Solve );
869  }
870 
871  }// while ( numRHS2Solve > 0 )
872 
873  }
874 
875  // print final summary
876  sTest_->print( printer_->stream(FinalSummary) );
877 
878  // print timing information
879 #ifdef BELOS_TEUCHOS_TIME_MONITOR
880  // Calling summarize() can be expensive, so don't call unless the
881  // user wants to print out timing details. summarize() will do all
882  // the work even if it's passed a "black hole" output stream.
883  if (verbosity_ & TimingDetails)
884  Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
885 #endif
886 
887  // get iteration information for this solve
888  numIters_ = maxIterTest_->getNumIters();
889 
890  // Save the convergence test value ("achieved tolerance") for this
891  // solve. For this solver, convTest_ may either be a single
892  // (implicit) residual norm test, or a combination of two residual
893  // norm tests. In the latter case, the master convergence test
894  // convTest_ is a SEQ combo of the implicit resp. explicit tests.
895  // If the implicit test never passes, then the explicit test won't
896  // ever be executed. This manifests as
897  // expConvTest_->getTestValue()->size() < 1. We deal with this case
898  // by using the values returned by impConvTest_->getTestValue().
899  {
900  // We'll fetch the vector of residual norms one way or the other.
901  const std::vector<MagnitudeType>* pTestValues = NULL;
902  if (expResTest_) {
903  pTestValues = expConvTest_->getTestValue();
904  if (pTestValues == NULL || pTestValues->size() < 1) {
905  pTestValues = impConvTest_->getTestValue();
906  }
907  }
908  else {
909  // Only the implicit residual norm test is being used.
910  pTestValues = impConvTest_->getTestValue();
911  }
912  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
913  "Belos::PseudoBlockTFQMRSolMgr::solve(): The implicit convergence test's "
914  "getTestValue() method returned NULL. Please report this bug to the "
915  "Belos developers.");
916  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
917  "Belos::TMQMRSolMgr::solve(): The implicit convergence test's "
918  "getTestValue() method returned a vector of length zero. Please report "
919  "this bug to the Belos developers.");
920 
921  // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
922  // achieved tolerances for all vectors in the current solve(), or
923  // just for the vectors from the last deflation?
924  achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
925  }
926 
927  if (!isConverged) {
928  return Unconverged; // return from PseudoBlockTFQMRSolMgr::solve()
929  }
930  return Converged; // return from PseudoBlockTFQMRSolMgr::solve()
931 }
932 
933 // This method requires the solver manager to return a std::string that describes itself.
934 template<class ScalarType, class MV, class OP>
936 {
937  std::ostringstream oss;
938  oss << "Belos::PseudoBlockTFQMRSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
939  oss << "{}";
940  return oss.str();
941 }
942 
943 } // end Belos namespace
944 
945 #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:87
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)
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...
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem)
Set the linear problem that needs to be solved.
ScaleType
The type of scaling to use on the residual norm value.
Definition: BelosTypes.hpp:119
Structure to contain pointers to PseudoBlockTFQMRIter state variables.
A factory class for generating StatusTestOutput objects.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Return a reference to the linear problem being solved by this solver manager.
Belos concrete class for generating iterations with the preconditioned tranpose-free QMR (TFQMR) meth...
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const
Get a parameter list containing the current parameters for this object.
Teuchos::RCP< const MV > R
The current residual basis.
bool isLOADetected() const
Whether loss of accuracy was detected during the last solve() invocation.
An implementation of StatusTestResNorm using a family of residual norms.
Belos::StatusTest class for specifying a maximum number of iterations.
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > &params)
Set the parameters the solver manager should use to solve the linear problem.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
A factory class for generating StatusTestOutput objects.
MagnitudeType achievedTol() const
Tolerance achieved by the last solve() invocation.
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...
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
PseudoBlockTFQMRSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i...
ResetType
How to reset the solver.
Definition: BelosTypes.hpp:200
Pure virtual base class which describes the basic interface for a solver manager. ...
std::string description() const
Method to return description of the pseudo-block TFQMR solver manager.
int getNumIters() const
Get the iteration count for the most recent call to solve().
PseudoBlockTFQMRSolMgrOrthoFailure(const std::string &what_arg)
ReturnType solve()
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:149
The Belos::SolverManager is a templated virtual base class that defines the basic interface that any ...
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Get a parameter list containing the valid parameters for this object.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
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.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
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
Belos header file which uses auto-configuration information to include necessary C++ headers...

Generated on Thu Jul 21 2016 14:43:57 for Belos by doxygen 1.8.11