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" 77 #include "Teuchos_BLAS.hpp" 94 template<
class ScalarType,
class MV,
class OP>
104 typedef Teuchos::ScalarTraits< ScalarType >
SCT;
106 typedef Teuchos::ScalarTraits< MagnitudeType >
SMT;
122 const Teuchos::ParameterList& params);
183 throw std::logic_error(
"getState() cannot be called unless " 184 "the state has been initialized");
209 Teuchos::RCP<const MV>
214 std::vector<MagnitudeType>& theNorms = *norms;
215 if (theNorms.size() < 1)
219 return Teuchos::null;
228 void symOrtho( ScalarType a, ScalarType b, ScalarType *c, ScalarType *s, ScalarType *r );
243 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
244 "Belos::MinresIter::setBlockSize(): Cannot use a block size that is not one.");
264 const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > >
lp_;
265 const Teuchos::RCP< OutputManager< ScalarType > >
om_;
266 const Teuchos::RCP< StatusTest< ScalarType, MV, OP > >
stest_;
300 Teuchos::RCP< MV >
Y_;
306 Teuchos::RCP< MV >
W_;
319 Teuchos::SerialDenseMatrix<int,ScalarType>
beta1_;
325 template<
class ScalarType,
class MV,
class OP>
329 const Teuchos::ParameterList & ):
334 stateStorageInitialized_(false),
342 template <
class ScalarType,
class MV,
class OP>
345 if (!stateStorageInitialized_) {
348 Teuchos::RCP< const MV > lhsMV = lp_->getLHS();
349 Teuchos::RCP< const MV > rhsMV = lp_->getRHS();
350 if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
351 stateStorageInitialized_ =
false;
358 if (Y_ == Teuchos::null) {
360 Teuchos::RCP< const MV > tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
361 TEUCHOS_TEST_FOR_EXCEPTION( tmp == Teuchos::null,
362 std::invalid_argument,
363 "Belos::MinresIter::setStateSize(): linear problem does not specify multivectors to clone from.");
364 Y_ = MVT::Clone( *tmp, 1 );
365 R1_ = MVT::Clone( *tmp, 1 );
366 R2_ = MVT::Clone( *tmp, 1 );
367 W_ = MVT::Clone( *tmp, 1 );
368 W1_ = MVT::Clone( *tmp, 1 );
369 W2_ = MVT::Clone( *tmp, 1 );
372 stateStorageInitialized_ =
true;
380 template <
class ScalarType,
class MV,
class OP>
384 if (!stateStorageInitialized_)
387 TEUCHOS_TEST_FOR_EXCEPTION( !stateStorageInitialized_,
388 std::invalid_argument,
389 "Belos::MinresIter::initialize(): Cannot initialize state storage!" );
391 TEUCHOS_TEST_FOR_EXCEPTION( newstate.
Y == Teuchos::null,
392 std::invalid_argument,
393 "Belos::MinresIter::initialize(): MinresIterationState does not have initial residual.");
395 std::string errstr(
"Belos::MinresIter::initialize(): Specified multivectors must have a consistent length and width.");
396 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.
Y) != MVT::GetGlobalLength(*Y_),
397 std::invalid_argument,
399 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
Y) != 1,
400 std::invalid_argument,
404 const ScalarType one = SCT::one();
410 MVT::MvAddMv( one, *newstate.
Y, zero, *newstate.
Y, *R2_ );
411 MVT::MvAddMv( one, *newstate.
Y, zero, *newstate.
Y, *R1_ );
415 MVT::MvInit ( *W2_ );
417 if ( lp_->getLeftPrec() != Teuchos::null ) {
418 lp_->applyLeftPrec( *newstate.
Y, *Y_ );
421 if (newstate.
Y != Y_) {
423 MVT::MvAddMv( one, *newstate.
Y, zero, *newstate.
Y, *Y_ );
428 beta1_ = Teuchos::SerialDenseMatrix<int,ScalarType>( 1, 1 );
429 MVT::MvTransMv( one, *newstate.
Y, *Y_, beta1_ );
431 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(beta1_(0,0)) < zero,
432 std::invalid_argument,
433 "The preconditioner is not positive definite." );
435 if( SCT::magnitude(beta1_(0,0)) == zero )
438 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
439 MVT::MvInit( *cur_soln_vec );
442 beta1_(0,0) = SCT::squareroot( beta1_(0,0) );
451 template <
class ScalarType,
class MV,
class OP>
457 if (initialized_ ==
false) {
461 Teuchos::BLAS<int,ScalarType> blas;
464 const ScalarType one = SCT::one();
468 Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
469 Teuchos::SerialDenseMatrix<int,ScalarType> beta( beta1_ );
470 phibar_ = Teuchos::ScalarTraits<ScalarType>::magnitude( beta1_(0,0) );
471 ScalarType shift = zero;
474 ScalarType oldBeta = zero;
475 ScalarType epsln = zero;
476 ScalarType cs = -one;
477 ScalarType sn = zero;
478 ScalarType dbar = zero;
488 Teuchos::RCP<MV> V = MVT::Clone( *Y_, 1 );
489 Teuchos::RCP<MV> tmpY, tmpW;
492 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
495 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1,
497 "Belos::MinresIter::iterate(): current linear system has more than one vector!" );
502 while (stest_->checkStatus(
this) !=
Passed) {
509 MVT::MvAddMv (one / beta(0,0), *Y_, zero, *Y_, *V);
512 lp_->applyOp (*V, *Y_);
516 MVT::MvAddMv (one, *Y_, -shift, *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_ );
539 MVT::MvAddMv( one, *R2_, zero, *R2_, *Y_ );
544 MVT::MvTransMv( one, *R2_, *Y_, beta );
556 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(beta(0,0)) <= zero,
558 "Belos::MinresIter::iterate(): Encountered nonpositi" 559 "ve value " << beta(0,0) <<
" for r2^H*M*r2 at itera" 560 "tion " << iter_ <<
": MINRES cannot continue." );
561 beta(0,0) = SCT::squareroot( beta(0,0) );
569 delta = cs*dbar + sn*alpha(0,0);
570 gbar = sn*dbar - cs*alpha(0,0);
571 epsln = sn*beta(0,0);
572 dbar = - cs*beta(0,0);
575 this->symOrtho(gbar, beta(0,0), &cs, &sn, &gamma);
578 phibar_ = Teuchos::ScalarTraits<ScalarType>::magnitude( sn * phibar_ );
582 MVT::MvAddMv( one, *W_, zero, *W_, *W1_ );
589 MVT::MvAddMv( one, *V, -oldeps, *W1_, *W_ );
590 MVT::MvAddMv( one, *W_, -delta, *W2_, *W_ );
591 MVT::MvScale( *W_, one / gamma );
595 MVT::MvAddMv( one, *cur_soln_vec, phi, *W_, *cur_soln_vec );
596 lp_->updateSolution();
606 template <
class ScalarType,
class MV,
class OP>
608 ScalarType *c, ScalarType *s, ScalarType *r
611 const ScalarType one = SCT::one();
612 const ScalarType zero = SCT::zero();
616 if ( absB == m_zero ) {
619 if ( absA == m_zero )
623 }
else if ( absA == m_zero ) {
627 }
else if ( absB >= absA ) {
628 ScalarType tau = a / b;
629 if ( Teuchos::ScalarTraits<ScalarType>::real(b) < m_zero )
630 *s = -one / SCT::squareroot( one+tau*tau );
632 *s = one / SCT::squareroot( one+tau*tau );
636 ScalarType tau = b / a;
637 if ( Teuchos::ScalarTraits<ScalarType>::real(a) < m_zero )
638 *c = -one / SCT::squareroot( one+tau*tau );
640 *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.