42 #ifndef BELOS_MINRES_ITER_HPP
43 #define BELOS_MINRES_ITER_HPP
93 template<
class ScalarType,
class MV,
class OP>
182 throw std::logic_error(
"getState() cannot be called unless "
183 "the state has been initialized");
213 std::vector<MagnitudeType>& theNorms = *norms;
214 if (theNorms.size() < 1)
216 theNorms[0] = phibar_;
218 return Teuchos::null;
227 void symOrtho( ScalarType a, ScalarType b, ScalarType *c, ScalarType *s, ScalarType *r );
243 "Belos::MinresIter::setBlockSize(): Cannot use a block size that is not one.");
283 bool stateStorageInitialized_;
324 template<
class ScalarType,
class MV,
class OP>
333 stateStorageInitialized_(false),
341 template <
class ScalarType,
class MV,
class OP>
344 if (!stateStorageInitialized_) {
349 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
350 stateStorageInitialized_ =
false;
357 if (Y_ == 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_)
387 std::invalid_argument,
388 "Belos::MinresIter::initialize(): Cannot initialize state storage!" );
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.");
396 std::invalid_argument,
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 ) {
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_ );
435 MVT::MvTransMv( one, *newstate.
Y, *Y_, beta1_ );
438 std::invalid_argument,
439 "The preconditioner is not positive definite." );
441 if( SCT::magnitude(beta1_(0,0)) == m_zero )
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();
478 ScalarType oldBeta = zero;
479 ScalarType epsln = zero;
480 ScalarType cs = -one;
481 ScalarType sn = zero;
482 ScalarType dbar = zero;
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 ) {
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 );
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);
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;
637 *s = -one / SCT::squareroot( one+tau*tau );
639 *s = one / SCT::squareroot( one+tau*tau );
643 ScalarType tau = b / a;
645 *c = -one / SCT::squareroot( one+tau*tau );
647 *c = one / SCT::squareroot( one+tau*tau );
Teuchos::RCP< const MV > R1
Previous residual.
Collection of types and exceptions used within the Belos solvers.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
Class which manages the output and verbosity of the Belos solvers.
Pure virtual base class for defining the status testing capabilities of Belos.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void resetNumIters(int iter=0)
Reset the iteration count.
Declaration of basic traits for the multivector type.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
A pure virtual class for defining the status tests for the Belos iterative solvers.
SCT::magnitudeType MagnitudeType
Class which defines basic traits for the operator type.
virtual ~MinresIter()
Destructor.
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.
Traits class which defines basic operations on multivectors.
bool isInitialized()
States whether the solver has been initialized or not.
void iterate()
Perform MINRES iterations until convergence or error.
Teuchos::RCP< const MV > R2
Previous residual.
MinresIterateFailure is thrown when the MinresIteration object is unable to compute the next iterate ...
void symOrtho(ScalarType a, ScalarType b, ScalarType *c, ScalarType *s, ScalarType *r)
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
Structure to contain pointers to MinresIteration state variables.
Teuchos::ScalarTraits< MagnitudeType > SMT
MultiVecTraits< ScalarType, MV > MVT
MinresIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Teuchos::RCP< const MV > W2
Previous direction vector.
static magnitudeType magnitude(T a)
Pure virtual base class which augments the basic interface for a minimal residual linear solver itera...
OperatorTraits< ScalarType, MV, OP > OPT
void initialize()
Initialize the solver.
Teuchos::RCP< const MV > Y
The current residual.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
void initializeMinres(const MinresIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
int getNumIters() const
Get the current iteration count.
void setBlockSize(int blockSize)
Set the blocksize to be used by the iterative solver in solving this linear problem.
Class which defines basic traits for the operator type.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Teuchos::ScalarTraits< ScalarType > SCT
bool isInitialized() const
States whether the solver has been initialized or not.
Belos header file which uses auto-configuration information to include necessary C++ headers...
Teuchos::RCP< const MV > W1
Previous direction vector.
Teuchos::RCP< const MV > W
The current direction vector.