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)
217 theNorms[0] = phibar_;
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_;
284 bool stateStorageInitialized_;
293 MagnitudeType phibar_;
300 Teuchos::RCP< MV > Y_;
302 Teuchos::RCP< MV > R1_;
304 Teuchos::RCP< MV > R2_;
306 Teuchos::RCP< MV > W_;
308 Teuchos::RCP< MV > W1_;
310 Teuchos::RCP< MV > W2_;
319 Teuchos::SerialDenseMatrix<int,ScalarType> beta1_;
325 template<
class ScalarType,
class MV,
class OP>
329 const Teuchos::ParameterList ¶ms ):
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.");
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.");
397 std::invalid_argument,
400 std::invalid_argument,
404 const ScalarType one = SCT::one();
417 if ( lp_->getLeftPrec() != Teuchos::null ) {
418 lp_->applyLeftPrec( *newstate.
Y, *Y_ );
421 if (newstate.
Y != Y_) {
428 beta1_ = Teuchos::SerialDenseMatrix<int,ScalarType>( 1, 1 );
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();
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;
489 Teuchos::RCP<MV> tmpY, tmpW;
492 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
497 "Belos::MinresIter::iterate(): current linear system has more than one vector!" );
502 while (stest_->checkStatus(
this) !=
Passed) {
512 lp_->applyOp (*V, *Y_);
525 MVT::MvAddMv (one, *Y_, -alpha(0,0)/beta(0,0), *R2_, *Y_);
535 if ( lp_->getLeftPrec() != Teuchos::null ) {
536 lp_->applyLeftPrec( *R2_, *Y_ );
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_ );
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.
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.
static void MvInit(MV &mv, const ScalarType alpha=Teuchos::ScalarTraits< ScalarType >::zero())
Replace each element of the vectors in mv with alpha.
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &mv, Teuchos::SerialDenseMatrix< int, ScalarType > &B)
Compute a dense matrix B through the matrix-matrix multiply .
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.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
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.
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
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.
Pure virtual base class which augments the basic interface for a minimal residual linear solver itera...
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
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.
static void MvScale(MV &mv, const ScalarType alpha)
Scale each element of the vectors in mv with alpha.
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.