42#ifndef BELOS_MINRES_ITER_HPP
43#define BELOS_MINRES_ITER_HPP
72#include "Teuchos_SerialDenseMatrix.hpp"
73#include "Teuchos_SerialDenseVector.hpp"
74#include "Teuchos_ScalarTraits.hpp"
75#include "Teuchos_ParameterList.hpp"
76#include "Teuchos_TimeMonitor.hpp"
93template<
class ScalarType,
class MV,
class OP>
103 typedef Teuchos::ScalarTraits< ScalarType >
SCT;
105 typedef Teuchos::ScalarTraits< MagnitudeType >
SMT;
121 const Teuchos::ParameterList& params);
182 throw std::logic_error(
"getState() cannot be called unless "
183 "the state has been initialized");
208 Teuchos::RCP<const MV>
213 std::vector<MagnitudeType>& theNorms = *norms;
214 if (theNorms.size() < 1)
218 return Teuchos::null;
227 void symOrtho( ScalarType a, ScalarType b, ScalarType *c, ScalarType *s, ScalarType *r );
242 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
243 "Belos::MinresIter::setBlockSize(): Cannot use a block size that is not one.");
263 const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > >
lp_;
264 const Teuchos::RCP< OutputManager< ScalarType > >
om_;
265 const Teuchos::RCP< StatusTest< ScalarType, MV, OP > >
stest_;
299 Teuchos::RCP< MV >
Y_;
305 Teuchos::RCP< MV >
W_;
318 Teuchos::SerialDenseMatrix<int,ScalarType>
beta1_;
324 template<
class ScalarType,
class MV,
class OP>
328 const Teuchos::ParameterList & ):
333 stateStorageInitialized_(false),
341 template <
class ScalarType,
class MV,
class OP>
344 if (!stateStorageInitialized_) {
347 Teuchos::RCP< const MV > lhsMV = lp_->getLHS();
348 Teuchos::RCP< const MV > rhsMV = lp_->getRHS();
349 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
350 stateStorageInitialized_ =
false;
357 if (Y_ == Teuchos::null) {
359 Teuchos::RCP< const MV > tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
360 TEUCHOS_TEST_FOR_EXCEPTION( tmp == Teuchos::null,
361 std::invalid_argument,
362 "Belos::MinresIter::setStateSize(): linear problem does not specify multivectors to clone from.");
363 Y_ = MVT::Clone( *tmp, 1 );
364 R1_ = MVT::Clone( *tmp, 1 );
365 R2_ = MVT::Clone( *tmp, 1 );
366 W_ = MVT::Clone( *tmp, 1 );
367 W1_ = MVT::Clone( *tmp, 1 );
368 W2_ = MVT::Clone( *tmp, 1 );
371 stateStorageInitialized_ =
true;
379 template <
class ScalarType,
class MV,
class OP>
383 if (!stateStorageInitialized_)
386 TEUCHOS_TEST_FOR_EXCEPTION( !stateStorageInitialized_,
387 std::invalid_argument,
388 "Belos::MinresIter::initialize(): Cannot initialize state storage!" );
390 TEUCHOS_TEST_FOR_EXCEPTION( newstate.
Y == Teuchos::null,
391 std::invalid_argument,
392 "Belos::MinresIter::initialize(): MinresIterationState does not have initial residual.");
394 std::string errstr(
"Belos::MinresIter::initialize(): Specified multivectors must have a consistent length and width.");
395 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
Y) != MVT::GetGlobalLength(*Y_),
396 std::invalid_argument,
398 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
Y) != 1,
399 std::invalid_argument,
403 const ScalarType one = SCT::one();
409 MVT::Assign( *newstate.
Y, *R2_ );
410 MVT::Assign( *newstate.
Y, *R1_ );
414 MVT::MvInit ( *W2_ );
416 if ( lp_->getLeftPrec() != Teuchos::null ) {
417 lp_->applyLeftPrec( *newstate.
Y, *Y_ );
418 if ( lp_->getRightPrec() != Teuchos::null ) {
419 Teuchos::RCP<MV> tmp = MVT::CloneCopy( *Y_ );
420 lp_->applyRightPrec( *tmp, *Y_ );
423 else if ( lp_->getRightPrec() != Teuchos::null ) {
424 lp_->applyRightPrec( *newstate.
Y, *Y_ );
427 if (newstate.
Y != Y_) {
429 MVT::Assign( *newstate.
Y, *Y_ );
434 beta1_ = Teuchos::SerialDenseMatrix<int,ScalarType>( 1, 1 );
435 MVT::MvTransMv( one, *newstate.
Y, *Y_, beta1_ );
437 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(beta1_(0,0)) < m_zero,
438 std::invalid_argument,
439 "The preconditioner is not positive definite." );
441 if( SCT::magnitude(beta1_(0,0)) == m_zero )
444 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
445 MVT::MvInit( *cur_soln_vec );
448 beta1_(0,0) = SCT::squareroot( beta1_(0,0) );
457 template <
class ScalarType,
class MV,
class OP>
463 if (initialized_ ==
false) {
468 const ScalarType one = SCT::one();
469 const ScalarType zero = SCT::zero();
473 Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
474 Teuchos::SerialDenseMatrix<int,ScalarType> beta( beta1_ );
475 phibar_ = Teuchos::ScalarTraits<ScalarType>::magnitude( beta1_(0,0) );
478 ScalarType oldBeta = zero;
479 ScalarType epsln = zero;
480 ScalarType cs = -one;
481 ScalarType sn = zero;
482 ScalarType dbar = zero;
492 Teuchos::RCP<MV> V = MVT::Clone( *Y_, 1 );
493 Teuchos::RCP<MV> tmpY, tmpW;
496 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
499 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1,
501 "Belos::MinresIter::iterate(): current linear system has more than one vector!" );
506 while (stest_->checkStatus(
this) !=
Passed) {
513 MVT::MvAddMv (one / beta(0,0), *Y_, zero, *Y_, *V);
516 lp_->applyOp (*V, *Y_);
519 MVT::MvAddMv (one, *Y_, -beta(0,0)/oldBeta, *R1_, *Y_);
522 MVT::MvTransMv (one, *V, *Y_, alpha);
525 MVT::MvAddMv (one, *Y_, -alpha(0,0)/beta(0,0), *R2_, *Y_);
535 if ( lp_->getLeftPrec() != Teuchos::null ) {
536 lp_->applyLeftPrec( *R2_, *Y_ );
537 if ( lp_->getRightPrec() != Teuchos::null ) {
538 Teuchos::RCP<MV> tmp = MVT::CloneCopy( *Y_ );
539 lp_->applyRightPrec( *tmp, *Y_ );
542 else if ( lp_->getRightPrec() != Teuchos::null ) {
543 lp_->applyRightPrec( *R2_, *Y_ );
546 MVT::Assign( *R2_, *Y_ );
551 MVT::MvTransMv( one, *R2_, *Y_, beta );
563 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(beta(0,0)) < m_zero,
565 "Belos::MinresIter::iterate(): Encountered negative "
566 "value " << beta(0,0) <<
" for r2^H*M*r2 at itera"
567 "tion " << iter_ <<
": MINRES cannot continue." );
568 beta(0,0) = SCT::squareroot( beta(0,0) );
576 delta = cs*dbar + sn*alpha(0,0);
577 gbar = sn*dbar - cs*alpha(0,0);
578 epsln = sn*beta(0,0);
579 dbar = - cs*beta(0,0);
582 this->symOrtho(gbar, beta(0,0), &cs, &sn, &gamma);
585 phibar_ = Teuchos::ScalarTraits<ScalarType>::magnitude( sn * phibar_ );
589 MVT::Assign( *W_, *W1_ );
596 MVT::MvAddMv( one, *V, -oldeps, *W1_, *W_ );
597 MVT::MvAddMv( one, *W_, -delta, *W2_, *W_ );
598 MVT::MvScale( *W_, one / gamma );
602 MVT::MvAddMv( one, *cur_soln_vec, phi, *W_, *cur_soln_vec );
603 lp_->updateSolution();
613 template <
class ScalarType,
class MV,
class OP>
615 ScalarType *c, ScalarType *s, ScalarType *r
618 const ScalarType one = SCT::one();
619 const ScalarType zero = SCT::zero();
623 if ( absB == m_zero ) {
626 if ( absA == m_zero )
630 }
else if ( absA == m_zero ) {
634 }
else if ( absB >= absA ) {
635 ScalarType tau = a / b;
636 if ( Teuchos::ScalarTraits<ScalarType>::real(b) < m_zero )
637 *s = -one / SCT::squareroot( one+tau*tau );
639 *s = one / SCT::squareroot( one+tau*tau );
643 ScalarType tau = b / a;
644 if ( Teuchos::ScalarTraits<ScalarType>::real(a) < m_zero )
645 *c = -one / SCT::squareroot( one+tau*tau );
647 *c = one / SCT::squareroot( one+tau*tau );
Belos header file which uses auto-configuration information to include necessary C++ headers.
Class which describes the linear problem to be solved by the iterative solver.
Pure virtual base class which augments the basic interface for a minimal residual linear solver itera...
Declaration of basic traits for the multivector type.
Class which defines basic traits for the operator type.
Class which manages the output and verbosity of the Belos solvers.
Pure virtual base class for defining the status testing capabilities of Belos.
Collection of types and exceptions used within the Belos solvers.
A linear system to solve, and its associated information.
void initializeMinres(const MinresIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
Teuchos::SerialDenseMatrix< int, ScalarType > beta1_
Coefficient in the MINRES iteration.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
Teuchos::RCP< MV > W1_
Previous direction vector.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
MinresIter(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< OutputManager< ScalarType > > &printer, const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &tester, const Teuchos::ParameterList ¶ms)
Constructor.
int iter_
Current number of iterations performed.
bool isInitialized()
States whether the solver has been initialized or not.
int getNumIters() const
Get the current iteration count.
bool isInitialized() const
States whether the solver has been initialized or not.
void setBlockSize(int blockSize)
Set the blocksize to be used by the iterative solver in solving this linear problem.
Teuchos::RCP< MV > Y_
Preconditioned residual.
MinresIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Teuchos::RCP< MV > W_
Direction vector.
MultiVecTraits< ScalarType, MV > MVT
Teuchos::ScalarTraits< MagnitudeType > SMT
const Teuchos::RCP< OutputManager< ScalarType > > om_
Teuchos::RCP< MV > R1_
Previous residual.
void initialize()
Initialize the solver.
void symOrtho(ScalarType a, ScalarType b, ScalarType *c, ScalarType *s, ScalarType *r)
Teuchos::RCP< MV > R2_
Previous residual.
void iterate()
Perform MINRES iterations until convergence or error.
SCT::magnitudeType MagnitudeType
Teuchos::ScalarTraits< ScalarType > SCT
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
OperatorTraits< ScalarType, MV, OP > OPT
bool initialized_
Whether the solver has been initialized.
void resetNumIters(int iter=0)
Reset the iteration count.
void setStateSize()
Method for initalizing the state storage needed by MINRES.
Teuchos::RCP< MV > W2_
Previous direction vector.
bool stateStorageInitialized_
Whether the state storage has been initialized.
virtual ~MinresIter()
Destructor.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
MinresIterateFailure is thrown when the MinresIteration object is unable to compute the next iterate ...
Traits class which defines basic operations on multivectors.
Class which defines basic traits for the operator type.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
A pure virtual class for defining the status tests for the Belos iterative solvers.
Structure to contain pointers to MinresIteration state variables.
Teuchos::RCP< const MV > W
The current direction vector.
Teuchos::RCP< const MV > R2
Previous residual.
Teuchos::RCP< const MV > W1
Previous direction vector.
Teuchos::RCP< const MV > Y
The current residual.
Teuchos::RCP< const MV > W2
Previous direction vector.
Teuchos::RCP< const MV > R1
Previous residual.