Belos  Version of the Day
BelosBiCGStabSolMgr.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_BICGSTAB_SOLMGR_HPP
43 #define BELOS_BICGSTAB_SOLMGR_HPP
44 
49 #include "BelosConfigDefs.hpp"
50 #include "BelosTypes.hpp"
51 
52 #include "BelosLinearProblem.hpp"
53 #include "BelosSolverManager.hpp"
54 
55 #include "BelosBiCGStabIter.hpp"
58 #include "BelosStatusTestCombo.hpp"
60 #include "BelosOutputManager.hpp"
61 #include "Teuchos_BLAS.hpp"
62 #include "Teuchos_LAPACK.hpp"
63 #ifdef BELOS_TEUCHOS_TIME_MONITOR
64 #include "Teuchos_TimeMonitor.hpp"
65 #endif
66 
83 namespace Belos {
84 
86 
87 
95  BiCGStabSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
96  {}};
97 
104  class BiCGStabSolMgrOrthoFailure : public BelosError {public:
105  BiCGStabSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
106  {}};
107 
108  template<class ScalarType, class MV, class OP>
109  class BiCGStabSolMgr : public SolverManager<ScalarType,MV,OP> {
110 
111  private:
114  typedef Teuchos::ScalarTraits<ScalarType> SCT;
115  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
116  typedef Teuchos::ScalarTraits<MagnitudeType> MT;
117 
118  public:
119 
121 
122 
128  BiCGStabSolMgr();
129 
139  BiCGStabSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
140  const Teuchos::RCP<Teuchos::ParameterList> &pl );
141 
143  virtual ~BiCGStabSolMgr() {};
145 
147 
148 
150  return *problem_;
151  }
152 
155  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
156 
159  Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
160 
166  Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
167  return Teuchos::tuple(timerSolve_);
168  }
169 
170 
181  MagnitudeType achievedTol() const {
182  return achievedTol_;
183  }
184 
186  int getNumIters() const {
187  return numIters_;
188  }
189 
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  private:
249 
250  // Linear problem.
251  Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
252 
253  // Output manager.
254  Teuchos::RCP<OutputManager<ScalarType> > printer_;
255  Teuchos::RCP<std::ostream> outputStream_;
256 
257  // Status test.
258  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
259  Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
260  Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
261  Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
262 
263  // Orthogonalization manager.
264  Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
265 
266  // Current parameter list.
267  Teuchos::RCP<Teuchos::ParameterList> params_;
268 
274  mutable Teuchos::RCP<const Teuchos::ParameterList> validParams_;
275 
276  // Default solver values.
277  static const MagnitudeType convtol_default_;
278  static const int maxIters_default_;
279  static const bool showMaxResNormOnly_default_;
280  static const int verbosity_default_;
281  static const int outputStyle_default_;
282  static const int outputFreq_default_;
283  static const int defQuorum_default_;
284  static const std::string resScale_default_;
285  static const std::string label_default_;
286  static const Teuchos::RCP<std::ostream> outputStream_default_;
287 
288  // Current solver values.
289  MagnitudeType convtol_,achievedTol_;
290  int maxIters_, numIters_;
291  int verbosity_, outputStyle_, outputFreq_, defQuorum_;
292  bool showMaxResNormOnly_;
293  std::string resScale_;
294 
295  // Timers.
296  std::string label_;
297  Teuchos::RCP<Teuchos::Time> timerSolve_;
298 
299  // Internal state variables.
300  bool isSet_;
301  };
302 
303 
304 // Default solver values.
305 template<class ScalarType, class MV, class OP>
306 const typename BiCGStabSolMgr<ScalarType,MV,OP>::MagnitudeType BiCGStabSolMgr<ScalarType,MV,OP>::convtol_default_ = 1e-8;
307 
308 template<class ScalarType, class MV, class OP>
310 
311 template<class ScalarType, class MV, class OP>
313 
314 template<class ScalarType, class MV, class OP>
316 
317 template<class ScalarType, class MV, class OP>
319 
320 template<class ScalarType, class MV, class OP>
322 
323 template<class ScalarType, class MV, class OP>
325 
326 template<class ScalarType, class MV, class OP>
327 const std::string BiCGStabSolMgr<ScalarType,MV,OP>::resScale_default_ = "Norm of Initial Residual";
328 
329 template<class ScalarType, class MV, class OP>
330 const std::string BiCGStabSolMgr<ScalarType,MV,OP>::label_default_ = "Belos";
331 
332 template<class ScalarType, class MV, class OP>
333 const Teuchos::RCP<std::ostream> BiCGStabSolMgr<ScalarType,MV,OP>::outputStream_default_ = Teuchos::rcp(&std::cout,false);
334 
335 // Empty Constructor
336 template<class ScalarType, class MV, class OP>
338  outputStream_(outputStream_default_),
339  convtol_(convtol_default_),
340  maxIters_(maxIters_default_),
341  numIters_(0),
342  verbosity_(verbosity_default_),
343  outputStyle_(outputStyle_default_),
344  outputFreq_(outputFreq_default_),
345  defQuorum_(defQuorum_default_),
346  showMaxResNormOnly_(showMaxResNormOnly_default_),
347  resScale_(resScale_default_),
348  label_(label_default_),
349  isSet_(false)
350 {}
351 
352 // Basic Constructor
353 template<class ScalarType, class MV, class OP>
355 BiCGStabSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
356  const Teuchos::RCP<Teuchos::ParameterList> &pl ) :
357  problem_(problem),
358  outputStream_(outputStream_default_),
359  convtol_(convtol_default_),
360  maxIters_(maxIters_default_),
361  numIters_(0),
362  verbosity_(verbosity_default_),
363  outputStyle_(outputStyle_default_),
364  outputFreq_(outputFreq_default_),
365  defQuorum_(defQuorum_default_),
366  showMaxResNormOnly_(showMaxResNormOnly_default_),
367  resScale_(resScale_default_),
368  label_(label_default_),
369  isSet_(false)
370 {
371  TEUCHOS_TEST_FOR_EXCEPTION(
372  problem_.is_null (), std::invalid_argument,
373  "Belos::BiCGStabSolMgr two-argument constructor: "
374  "'problem' is null. You must supply a non-null Belos::LinearProblem "
375  "instance when calling this constructor.");
376 
377  if (! pl.is_null ()) {
378  // Set the parameters using the list that was passed in.
379  setParameters (pl);
380  }
381 }
382 
383 template<class ScalarType, class MV, class OP>
384 void BiCGStabSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
385 {
386  using Teuchos::ParameterList;
387  using Teuchos::parameterList;
388  using Teuchos::RCP;
389 
390  RCP<const ParameterList> defaultParams = getValidParameters();
391 
392  // Create the internal parameter list if one doesn't already exist.
393  if (params_.is_null()) {
394  params_ = parameterList (*defaultParams);
395  } else {
396  params->validateParameters (*defaultParams);
397  }
398 
399  // Check for maximum number of iterations
400  if (params->isParameter("Maximum Iterations")) {
401  maxIters_ = params->get("Maximum Iterations",maxIters_default_);
402 
403  // Update parameter in our list and in status test.
404  params_->set("Maximum Iterations", maxIters_);
405  if (maxIterTest_!=Teuchos::null)
406  maxIterTest_->setMaxIters( maxIters_ );
407  }
408 
409  // Check to see if the timer label changed.
410  if (params->isParameter("Timer Label")) {
411  std::string tempLabel = params->get("Timer Label", label_default_);
412 
413  // Update parameter in our list and solver timer
414  if (tempLabel != label_) {
415  label_ = tempLabel;
416  params_->set("Timer Label", label_);
417  std::string solveLabel = label_ + ": BiCGStabSolMgr total solve time";
418 #ifdef BELOS_TEUCHOS_TIME_MONITOR
419  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
420 #endif
421  if (ortho_ != Teuchos::null) {
422  ortho_->setLabel( label_ );
423  }
424  }
425  }
426 
427  // Check for a change in verbosity level
428  if (params->isParameter("Verbosity")) {
429  if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
430  verbosity_ = params->get("Verbosity", verbosity_default_);
431  } else {
432  verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
433  }
434 
435  // Update parameter in our list.
436  params_->set("Verbosity", verbosity_);
437  if (printer_ != Teuchos::null)
438  printer_->setVerbosity(verbosity_);
439  }
440 
441  // Check for a change in output style
442  if (params->isParameter("Output Style")) {
443  if (Teuchos::isParameterType<int>(*params,"Output Style")) {
444  outputStyle_ = params->get("Output Style", outputStyle_default_);
445  } else {
446  outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
447  }
448 
449  // Reconstruct the convergence test if the explicit residual test is not being used.
450  params_->set("Output Style", outputStyle_);
451  outputTest_ = Teuchos::null;
452  }
453 
454  // output stream
455  if (params->isParameter("Output Stream")) {
456  outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
457 
458  // Update parameter in our list.
459  params_->set("Output Stream", outputStream_);
460  if (printer_ != Teuchos::null)
461  printer_->setOStream( outputStream_ );
462  }
463 
464  // frequency level
465  if (verbosity_ & Belos::StatusTestDetails) {
466  if (params->isParameter("Output Frequency")) {
467  outputFreq_ = params->get("Output Frequency", outputFreq_default_);
468  }
469 
470  // Update parameter in out list and output status test.
471  params_->set("Output Frequency", outputFreq_);
472  if (outputTest_ != Teuchos::null)
473  outputTest_->setOutputFrequency( outputFreq_ );
474  }
475 
476  // Create output manager if we need to.
477  if (printer_ == Teuchos::null) {
478  printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
479  }
480 
481  // Convergence
482  typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
483  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t;
484 
485  // Check for convergence tolerance
486  if (params->isParameter("Convergence Tolerance")) {
487  convtol_ = params->get("Convergence Tolerance",convtol_default_);
488 
489  // Update parameter in our list and residual tests.
490  params_->set("Convergence Tolerance", convtol_);
491  if (convTest_ != Teuchos::null)
492  convTest_->setTolerance( convtol_ );
493  }
494 
495  if (params->isParameter("Show Maximum Residual Norm Only")) {
496  showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
497 
498  // Update parameter in our list and residual tests
499  params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
500  if (convTest_ != Teuchos::null)
501  convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
502  }
503 
504  // Check for a change in scaling, if so we need to build new residual tests.
505  bool newResTest = false;
506  {
507  // "Residual Scaling" is the old parameter name; "Implicit
508  // Residual Scaling" is the new name. We support both options for
509  // backwards compatibility.
510  std::string tempResScale = resScale_;
511  bool implicitResidualScalingName = false;
512  if (params->isParameter ("Residual Scaling")) {
513  tempResScale = params->get<std::string> ("Residual Scaling");
514  }
515  else if (params->isParameter ("Implicit Residual Scaling")) {
516  tempResScale = params->get<std::string> ("Implicit Residual Scaling");
517  implicitResidualScalingName = true;
518  }
519 
520  // Only update the scaling if it's different.
521  if (resScale_ != tempResScale) {
522  Belos::ScaleType resScaleType = convertStringToScaleType( tempResScale );
523  resScale_ = tempResScale;
524 
525  // Update parameter in our list and residual tests, using the
526  // given parameter name.
527  if (implicitResidualScalingName) {
528  params_->set ("Implicit Residual Scaling", resScale_);
529  }
530  else {
531  params_->set ("Residual Scaling", resScale_);
532  }
533 
534  if (! convTest_.is_null()) {
535  try {
536  convTest_->defineScaleForm( resScaleType, Belos::TwoNorm );
537  }
538  catch (std::exception& e) {
539  // Make sure the convergence test gets constructed again.
540  newResTest = true;
541  }
542  }
543  }
544  }
545 
546  // Get the deflation quorum, or number of converged systems before deflation is allowed
547  if (params->isParameter("Deflation Quorum")) {
548  defQuorum_ = params->get("Deflation Quorum", defQuorum_);
549  params_->set("Deflation Quorum", defQuorum_);
550  if (convTest_ != Teuchos::null)
551  convTest_->setQuorum( defQuorum_ );
552  }
553 
554  // Create status tests if we need to.
555 
556  // Basic test checks maximum iterations and native residual.
557  if (maxIterTest_ == Teuchos::null)
558  maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
559 
560  // Implicit residual test, using the native residual to determine if convergence was achieved.
561  if (convTest_ == Teuchos::null || newResTest) {
562  convTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_, defQuorum_, showMaxResNormOnly_ ) );
563  convTest_->defineScaleForm( convertStringToScaleType( resScale_ ), Belos::TwoNorm );
564  }
565 
566  if (sTest_ == Teuchos::null || newResTest)
567  sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
568 
569  if (outputTest_ == Teuchos::null || newResTest) {
570 
571  // Create the status test output class.
572  // This class manages and formats the output from the status test.
573  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
574  outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
575 
576  // Set the solver string for the output test
577  std::string solverDesc = " Pseudo Block BiCGStab ";
578  outputTest_->setSolverDesc( solverDesc );
579 
580  }
581 
582  // Create the timer if we need to.
583  if (timerSolve_ == Teuchos::null) {
584  std::string solveLabel = label_ + ": BiCGStabSolMgr total solve time";
585 #ifdef BELOS_TEUCHOS_TIME_MONITOR
586  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
587 #endif
588  }
589 
590  // Inform the solver manager that the current parameters were set.
591  isSet_ = true;
592 }
593 
594 
595 template<class ScalarType, class MV, class OP>
596 Teuchos::RCP<const Teuchos::ParameterList>
598 {
599  using Teuchos::ParameterList;
600  using Teuchos::parameterList;
601  using Teuchos::RCP;
602 
603  if (validParams_.is_null()) {
604  // Set all the valid parameters and their default values.
605  RCP<ParameterList> pl = parameterList ();
606  pl->set("Convergence Tolerance", convtol_default_,
607  "The relative residual tolerance that needs to be achieved by the\n"
608  "iterative solver in order for the linera system to be declared converged.");
609  pl->set("Maximum Iterations", maxIters_default_,
610  "The maximum number of block iterations allowed for each\n"
611  "set of RHS solved.");
612  pl->set("Verbosity", verbosity_default_,
613  "What type(s) of solver information should be outputted\n"
614  "to the output stream.");
615  pl->set("Output Style", outputStyle_default_,
616  "What style is used for the solver information outputted\n"
617  "to the output stream.");
618  pl->set("Output Frequency", outputFreq_default_,
619  "How often convergence information should be outputted\n"
620  "to the output stream.");
621  pl->set("Deflation Quorum", defQuorum_default_,
622  "The number of linear systems that need to converge before\n"
623  "they are deflated. This number should be <= block size.");
624  pl->set("Output Stream", outputStream_default_,
625  "A reference-counted pointer to the output stream where all\n"
626  "solver output is sent.");
627  pl->set("Show Maximum Residual Norm Only", showMaxResNormOnly_default_,
628  "When convergence information is printed, only show the maximum\n"
629  "relative residual norm when the block size is greater than one.");
630  pl->set("Implicit Residual Scaling", resScale_default_,
631  "The type of scaling used in the residual convergence test.");
632  // We leave the old name as a valid parameter for backwards
633  // compatibility (so that validateParametersAndSetDefaults()
634  // doesn't raise an exception if it encounters "Residual
635  // Scaling"). The new name was added for compatibility with other
636  // solvers, none of which use "Residual Scaling".
637  pl->set("Residual Scaling", resScale_default_,
638  "The type of scaling used in the residual convergence test. This "
639  "name is deprecated; the new name is \"Implicit Residual Scaling\".");
640  pl->set("Timer Label", label_default_,
641  "The string to use as a prefix for the timer labels.");
642  // defaultParams_->set("Restart Timers", restartTimers_);
643  validParams_ = pl;
644  }
645  return validParams_;
646 }
647 
648 
649 template<class ScalarType, class MV, class OP>
651 {
652  // Set the current parameters if they were not set before.
653  // NOTE: This may occur if the user generated the solver manager with the default constructor and
654  // then didn't set any parameters using setParameters().
655  if (! isSet_) {
656  setParameters (params_);
657  }
658 
659  Teuchos::BLAS<int,ScalarType> blas;
660 
661  TEUCHOS_TEST_FOR_EXCEPTION
662  (! problem_->isProblemSet (), BiCGStabSolMgrLinearProblemFailure,
663  "Belos::BiCGStabSolMgr::solve: Linear problem is not ready. "
664  "You must call setProblem() on the LinearProblem before you may solve it.");
665  TEUCHOS_TEST_FOR_EXCEPTION
666  (problem_->isLeftPrec (), std::logic_error, "Belos::BiCGStabSolMgr::solve: "
667  "The left-preconditioned case has not yet been implemented. Please use "
668  "right preconditioning for now. If you need to use left preconditioning, "
669  "please contact the Belos developers. Left preconditioning is more "
670  "interesting in BiCGStab because whether it works depends on the initial "
671  "guess (e.g., an initial guess of all zeros might NOT work).");
672 
673  // Create indices for the linear systems to be solved.
674  int startPtr = 0;
675  int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
676  int numCurrRHS = numRHS2Solve;
677 
678  std::vector<int> currIdx( numRHS2Solve ), currIdx2( numRHS2Solve );
679  for (int i=0; i<numRHS2Solve; ++i) {
680  currIdx[i] = startPtr+i;
681  currIdx2[i]=i;
682  }
683 
684  // Inform the linear problem of the current linear system to solve.
685  problem_->setLSIndex( currIdx );
686 
688  // Parameter list (iteration)
689  Teuchos::ParameterList plist;
690 
691  // Reset the status test.
692  outputTest_->reset();
693 
694  // Assume convergence is achieved, then let any failed convergence set this to false.
695  bool isConverged = true;
696 
698  // Pseudo-Block BiCGStab solver
699 
700  Teuchos::RCP<BiCGStabIter<ScalarType,MV,OP> > block_cg_iter
701  = Teuchos::rcp( new BiCGStabIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,plist) );
702 
703  // Enter solve() iterations
704  {
705 #ifdef BELOS_TEUCHOS_TIME_MONITOR
706  Teuchos::TimeMonitor slvtimer(*timerSolve_);
707 #endif
708 
709  //bool first_time=true;
710  while ( numRHS2Solve > 0 ) {
711  // Reset the active / converged vectors from this block
712  std::vector<int> convRHSIdx;
713  std::vector<int> currRHSIdx( currIdx );
714  currRHSIdx.resize(numCurrRHS);
715 
716  // Reset the number of iterations.
717  block_cg_iter->resetNumIters();
718 
719  // Reset the number of calls that the status test output knows about.
720  outputTest_->resetNumCalls();
721 
722  // Get the current residual for this block of linear systems.
723  Teuchos::RCP<MV> R_0 = MVT::CloneViewNonConst( *(Teuchos::rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
724 
725  // Get a new state struct and initialize the solver.
727  newState.R = R_0;
728  block_cg_iter->initializeBiCGStab(newState);
729 
730  while(1) {
731 
732  // tell block_gmres_iter to iterate
733  try {
734 
735  block_cg_iter->iterate();
736 
738  //
739  // check convergence first
740  //
742  if ( convTest_->getStatus() == Passed ) {
743 
744  // Figure out which linear systems converged.
745  std::vector<int> convIdx = Teuchos::rcp_dynamic_cast<StatusTestGenResNorm<ScalarType,MV,OP> >(convTest_)->convIndices();
746 
747  // If the number of converged linear systems is equal to the
748  // number of current linear systems, then we are done with this block.
749  if (convIdx.size() == currRHSIdx.size())
750  break; // break from while(1){block_cg_iter->iterate()}
751 
752  // Inform the linear problem that we are finished with this current linear system.
753  problem_->setCurrLS();
754 
755  // Reset currRHSIdx to have the right-hand sides that are left to converge for this block.
756  int have = 0;
757  std::vector<int> unconvIdx(currRHSIdx.size());
758  for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
759  bool found = false;
760  for (unsigned int j=0; j<convIdx.size(); ++j) {
761  if (currRHSIdx[i] == convIdx[j]) {
762  found = true;
763  break;
764  }
765  }
766  if (!found) {
767  currIdx2[have] = currIdx2[i];
768  currRHSIdx[have++] = currRHSIdx[i];
769  }
770  }
771  currRHSIdx.resize(have);
772  currIdx2.resize(have);
773 
774  // Set the remaining indices after deflation.
775  problem_->setLSIndex( currRHSIdx );
776 
777  // Get the current residual vector.
778  std::vector<MagnitudeType> norms;
779  R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
780  for (int i=0; i<have; ++i) { currIdx2[i] = i; }
781 
782  // Set the new state and initialize the solver.
784  defstate.R = R_0;
785  block_cg_iter->initializeBiCGStab(defstate);
786  }
787 
789  //
790  // check for maximum iterations
791  //
793  else if ( maxIterTest_->getStatus() == Passed ) {
794  // we don't have convergence
795  isConverged = false;
796  break; // break from while(1){block_cg_iter->iterate()}
797  }
798 
800  //
801  // we returned from iterate(), but none of our status tests Passed.
802  // something is wrong, and it is probably our fault.
803  //
805 
806  else {
807  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
808  "Belos::BiCGStabSolMgr::solve(): Invalid return from BiCGStabIter::iterate().");
809  }
810  }
811  catch (const std::exception &e) {
812  printer_->stream(Errors) << "Error! Caught std::exception in BiCGStabIter::iterate() at iteration "
813  << block_cg_iter->getNumIters() << std::endl
814  << e.what() << std::endl;
815  throw;
816  }
817  }
818 
819  // Inform the linear problem that we are finished with this block linear system.
820  problem_->setCurrLS();
821 
822  // Update indices for the linear systems to be solved.
823  startPtr += numCurrRHS;
824  numRHS2Solve -= numCurrRHS;
825 
826  if ( numRHS2Solve > 0 ) {
827 
828  numCurrRHS = numRHS2Solve;
829  currIdx.resize( numCurrRHS );
830  currIdx2.resize( numCurrRHS );
831  for (int i=0; i<numCurrRHS; ++i)
832  { currIdx[i] = startPtr+i; currIdx2[i] = i; }
833 
834  // Set the next indices.
835  problem_->setLSIndex( currIdx );
836  }
837  else {
838  currIdx.resize( numRHS2Solve );
839  }
840 
841  //first_time=false;
842  }// while ( numRHS2Solve > 0 )
843 
844  }
845 
846  // print final summary
847  sTest_->print( printer_->stream(FinalSummary) );
848 
849  // print timing information
850 #ifdef BELOS_TEUCHOS_TIME_MONITOR
851  // Calling summarize() can be expensive, so don't call unless the
852  // user wants to print out timing details. summarize() will do all
853  // the work even if it's passed a "black hole" output stream.
854  if (verbosity_ & TimingDetails)
855  Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
856 #endif
857 
858  // get iteration information for this solve
859  numIters_ = maxIterTest_->getNumIters();
860 
861 
862  // Save the convergence test value ("achieved tolerance") for this
863  // solve.
864  const std::vector<MagnitudeType>* pTestValues = convTest_->getTestValue();
865  achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
866 
867 
868  if (!isConverged ) {
869  return Unconverged; // return from BiCGStabSolMgr::solve()
870  }
871  return Converged; // return from BiCGStabSolMgr::solve()
872 }
873 
874 // This method requires the solver manager to return a std::string that describes itself.
875 template<class ScalarType, class MV, class OP>
877 {
878  std::ostringstream oss;
879  oss << "Belos::BiCGStabSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
880  oss << "{";
881  oss << "}";
882  return oss.str();
883 }
884 
885 
886 
887 } // end Belos namespace
888 
889 #endif /* BELOS_BICGSTAB_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.
BiCGStabSolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate orthono...
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
This class implements the pseudo-block BiCGStab iteration, where the basic BiCGStab algorithm is perf...
Class which manages the output and verbosity of the Belos solvers.
MagnitudeType achievedTol() const
Tolerance achieved by the last solve() invocation.
Teuchos::RCP< const MV > R
The current residual.
ScaleType
The type of scaling to use on the residual norm value.
Definition: BelosTypes.hpp:119
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const
Get a parameter list containing the current parameters for this object.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Return a reference to the linear problem being solved by this solver manager.
A factory class for generating StatusTestOutput objects.
Structure to contain pointers to BiCGStabIteration state variables.
bool isLOADetected() const
Return whether a loss of accuracy was detected by this solver during the most current solve...
An implementation of StatusTestResNorm using a family of residual norms.
std::string description() const
Method to return description of the block BiCGStab solver manager.
Belos::StatusTest class for specifying a maximum number of iterations.
int getNumIters() const
Get the iteration count for the most recent call to solve().
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Get a parameter list containing the valid parameters for this object.
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > &params)
Set the parameters the solver manager should use to solve the linear problem.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem)
Set the linear problem that needs to be solved.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
A factory class for generating StatusTestOutput objects.
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
BiCGStabSolMgrLinearProblemFailure(const std::string &what_arg)
A Belos::StatusTest class for specifying a maximum number of iterations.
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).
ResetType
How to reset the solver.
Definition: BelosTypes.hpp:200
Pure virtual base class which describes the basic interface for a solver manager. ...
BiCGStabSolMgr()
Empty constructor for BiCGStabSolMgr. This constructor takes no arguments and sets the default values...
virtual ~BiCGStabSolMgr()
Destructor.
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
BiCGStabSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i...
ReturnType solve()
This method performs possibly repeated calls to the underlying linear solver&#39;s iterate() routine unti...
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
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 ...
void reset(const ResetType type)
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
The Belos::BiCGStabSolMgr provides a powerful and fully-featured solver manager over the pseudo-block...
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...
BiCGStabSolMgrOrthoFailure(const std::string &what_arg)
Belos concrete class for performing the pseudo-block BiCGStab iteration.

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