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" 93 template<
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 );
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.
bool isInitialized() const
States whether the solver has been initialized or not.
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.
void resetNumIters(int iter=0)
Reset the iteration count.
Declaration of basic traits for the multivector type.
int getNumIters() const
Get the current iteration count.
bool initialized_
Whether the solver has been initialized.
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.
MinresIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
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)
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
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.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Teuchos::RCP< const MV > W2
Previous direction vector.
Teuchos::RCP< MV > W1_
Previous direction vector.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
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.
void initializeMinres(const MinresIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
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.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Teuchos::ScalarTraits< ScalarType > SCT
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.