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)
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.");
324 template<
class ScalarType,
class MV,
class OP>
333 stateStorageInitialized_(false),
341 template <
class ScalarType,
class MV,
class OP>
344 if (!stateStorageInitialized_) {
350 stateStorageInitialized_ =
false;
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_ );
417 lp_->applyLeftPrec( *newstate.
Y, *Y_ );
420 lp_->applyRightPrec( *tmp, *Y_ );
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_);
536 lp_->applyLeftPrec( *R2_, *Y_ );
539 lp_->applyRightPrec( *tmp, *Y_ );
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.
Teuchos::RCP< MV > R1_
Previous residual.
Teuchos::RCP< MV > W2_
Previous direction vector.
const Teuchos::RCP< OutputManager< ScalarType > > om_
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.
bool initialized_
Whether the solver has been initialized.
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.
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
void iterate()
Perform MINRES iterations until convergence or error.
Teuchos::SerialDenseMatrix< int, ScalarType > beta1_
Coefficient in the MINRES iteration.
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
int iter_
Current number of iterations performed.
bool stateStorageInitialized_
Whether the state storage has been initialized.
MinresIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Teuchos::RCP< const MV > W2
Previous direction vector.
Teuchos::RCP< MV > W1_
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 setStateSize()
Method for initalizing the state storage needed by MINRES.
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< MV > R2_
Previous residual.
Teuchos::RCP< const MV > W
The current direction vector.
Teuchos::RCP< MV > W_
Direction vector.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
Teuchos::RCP< MV > Y_
Preconditioned residual.