42 #ifndef BELOS_LSQR_ITER_HPP
43 #define BELOS_LSQR_ITER_HPP
76 template<
class ScalarType,
class MV,
class OP>
206 "LSQRIter::setBlockSize(): Cannot use a block size that is not one.");
288 template<
class ScalarType,
class MV,
class OP>
297 stateStorageInitialized_(false),
303 mat_resid_norm_(0.0),
311 template <
class ScalarType,
class MV,
class OP>
314 if (!stateStorageInitialized_) {
319 stateStorageInitialized_ =
false;
331 U_ = MVT::Clone( *rhsMV, 1 );
332 V_ = MVT::Clone( *lhsMV, 1 );
333 W_ = MVT::Clone( *lhsMV, 1 );
337 stateStorageInitialized_ =
true;
345 template <
class ScalarType,
class MV,
class OP>
351 if (!stateStorageInitialized_)
355 "LSQRIter::initialize(): Cannot initialize state storage!");
357 std::string errstr(
"LSQRIter::initialize(): Specified multivectors must have a consistent length and width.");
362 RCP<const MV> lhsMV = lp_->getLHS();
364 const bool debugSerialLSQR =
false;
366 if( debugSerialLSQR )
368 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> lhsNorm(1);
369 MVT::MvNorm( *lhsMV, lhsNorm );
370 std::cout <<
"initializeLSQR lhsNorm " << lhsNorm[0] << std::endl;
374 RCP<const MV> rhsMV = lp_->getInitPrecResVec();
377 MVT::MvAddMv( one, *rhsMV, zero, *U_, *U_);
379 RCP<const OP> M_left = lp_->getLeftPrec();
380 RCP<const OP>
A = lp_->getOperator();
381 RCP<const OP> M_right = lp_->getRightPrec();
383 if( debugSerialLSQR )
385 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> rhsNorm(1);
386 MVT::MvNorm( *U_, rhsNorm );
387 std::cout <<
"initializeLSQR | U_ | : " << rhsNorm[0] << std::endl;
399 if ( M_left.is_null())
407 RCP<MV> tempInRangeOfA = MVT::CloneCopy (*U_);
408 OPT::Apply (*M_left, *U_, *tempInRangeOfA,
CONJTRANS);
409 OPT::Apply (*A, *tempInRangeOfA, *V_,
CONJTRANS);
412 if (! M_right.is_null())
414 RCP<MV> tempInDomainOfA = MVT::CloneCopy (*V_);
416 OPT::Apply (*M_right, *tempInDomainOfA, *V_,
CONJTRANS);
421 MVT::MvAddMv( one, *V_, zero, *V_, *W_);
423 frob_mat_norm_ = zero;
434 template <
class ScalarType,
class MV,
class OP>
440 if (initialized_ ==
false) {
450 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> alpha(1);
451 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> beta(1);
452 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> xi(1);
455 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> wnorm2(1);
456 ScalarType rhobar, phibar, cs1, phi, rho, cs, sn, theta, xxnorm = MTzero, common;
457 ScalarType zetabar, sn1, psi, res = zero, ddnorm = zero, gamma, tau;
458 ScalarType anorm2 = zero;
459 ScalarType cs2 = -one, sn2 = zero, gammabar, zeta = zero, delta;
463 AV = MVT::Clone( *U_, 1);
465 AtU = MVT::Clone( *V_, 1);
466 const bool debugSerialLSQR =
false;
474 "LSQRIter::iterate(): current linear system has more than one vector!" );
478 MVT::MvNorm( *U_, beta );
479 if( SCT::real(beta[0]) > zero )
481 MVT::MvScale( *U_, one / beta[0] );
486 MVT::MvScale( *V_, one / beta[0] );
491 MVT::MvNorm( *V_, alpha );
492 if( debugSerialLSQR )
496 std::cout << iter_ <<
" First alpha " << alpha[0] <<
" beta " << beta[0] <<
" lambda " << lambda_ << std::endl;
498 if( SCT::real(alpha[0]) > zero )
500 MVT::MvScale( *V_, one / alpha[0] );
502 if( beta[0] * alpha[0] > zero )
504 MVT::MvScale( *W_, one / (beta[0] * alpha[0]) );
508 MVT::MvScale( *W_, zero );
512 RCP<const OP> M_left = lp_->getLeftPrec();
513 RCP<const OP>
A = lp_->getOperator();
514 RCP<const OP> M_right = lp_->getRightPrec();
519 resid_norm_ = beta[0];
520 mat_resid_norm_ = alpha[0] * beta[0];
534 if ( M_right.is_null() )
536 OPT::Apply(*A, *V_, *AV);
540 RCP<MV> tempInDomainOfA = MVT::CloneCopy (*V_);
541 OPT::Apply (*M_right, *V_, *tempInDomainOfA);
542 OPT::Apply(*A, *tempInDomainOfA, *AV);
545 if (! M_left.is_null())
547 RCP<MV> tempInRangeOfA = MVT::CloneCopy (*AV);
548 OPT::Apply (*M_left, *tempInRangeOfA, *AV);
552 if ( !( M_left.is_null() && M_right.is_null() )
553 && debugSerialLSQR && iter_ == 1)
559 if (! M_left.is_null())
561 RCP<MV> tempInRangeOfA = MVT::CloneCopy (*U_);
562 OPT::Apply (*M_left, *U_, *tempInRangeOfA,
CONJTRANS);
563 OPT::Apply (*A, *tempInRangeOfA, *AtU,
CONJTRANS);
569 if ( !( M_right.is_null() ) )
571 RCP<MV> tempInDomainOfA = MVT::CloneCopy (*AtU);
572 OPT::Apply (*M_right, *tempInDomainOfA, *AtU,
CONJTRANS);
575 MVT::MvAddMv( one, *AtU, -alpha[0], *V_, *AtU );
576 MVT::MvNorm( *AtU, xi );
577 std::cout <<
"| V alpha - A' u |= " << xi[0] << std::endl;
580 MVT::MvTransMv( one, *U_, *U_, uotuo );
581 std::cout <<
"<U, U> = " << uotuo << std::endl;
583 std::cout <<
"alpha = " << alpha[0] << std::endl;
586 MVT::MvTransMv( one, *AV, *U_, utav );
587 std::cout <<
"<AV, U> = alpha = " << utav << std::endl;
590 MVT::MvAddMv( one, *AV, -alpha[0], *U_, *U_ );
591 MVT::MvNorm( *U_, beta);
593 anorm2 += alpha[0]*alpha[0] + beta[0]*beta[0] + lambda_*lambda_;
596 if ( SCT::real(beta[0]) > zero )
599 MVT::MvScale( *U_, one / beta[0] );
601 if (M_left.is_null())
607 RCP<MV> tempInRangeOfA = MVT::CloneCopy (*U_);
608 OPT::Apply (*M_left, *U_, *tempInRangeOfA,
CONJTRANS);
609 OPT::Apply(*A, *tempInRangeOfA, *AtU,
CONJTRANS);
611 if (! M_right.is_null())
613 RCP<MV> tempInDomainOfA = MVT::CloneCopy (*AtU);
614 OPT::Apply (*M_right, *tempInDomainOfA, *AtU,
CONJTRANS);
617 MVT::MvAddMv( one, *AtU, -beta[0], *V_, *V_ );
618 MVT::MvNorm( *V_, alpha );
625 if ( SCT::real(alpha[0]) > zero )
627 MVT::MvScale( *V_, one / alpha[0] );
633 cs1 = rhobar / common;
634 sn1 = lambda_ / common;
636 phibar = cs1 * phibar;
643 theta = sn * alpha[0];
644 rhobar = -cs * alpha[0];
646 phibar = sn * phibar;
650 gammabar = -cs2 * rho;
651 zetabar = (phi - delta*zeta) / gammabar;
654 cs2 = gammabar / gamma;
656 zeta = (phi - delta*zeta) / gamma;
660 if ( M_right.is_null())
662 MVT::MvAddMv( phi / rho, *W_, one, *cur_soln_vec, *cur_soln_vec);
666 RCP<MV> tempInDomainOfA = MVT::CloneCopy (*W_);
667 OPT::Apply (*M_right, *W_, *tempInDomainOfA);
668 MVT::MvAddMv( phi / rho, *tempInDomainOfA, one, *cur_soln_vec, *cur_soln_vec);
671 MVT::MvNorm( *W_, wnorm2 );
672 ddnorm += (one / rho)*(one / rho) * wnorm2[0]*wnorm2[0];
673 MVT::MvAddMv( one, *V_, -theta / rho, *W_, *W_ );
Teuchos::RCP< const MV > U
Bidiagonalization vector.
Collection of types and exceptions used within the Belos solvers.
IterationState contains the data that defines the state of the LSQR solver at any given time...
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
const Teuchos::RCP< Belos::StatusTest< ScalarType, MV, OP > > stest_
ScalarType mat_resid_norm_
Class which manages the output and verbosity of the Belos solvers.
const Teuchos::RCP< Belos::OutputManager< ScalarType > > om_
ScalarType resid_norm
The current residual norm.
Teuchos::ScalarTraits< ScalarType >::magnitudeType lambda
The damping value.
void resetNumIters(int iter=0)
Reset the iteration count.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
Pure virtual base class for defining the status testing capabilities of Belos.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Declaration of basic traits for the multivector type.
void initializeLSQR(LSQRIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, completing the initial state.
bool stateStorageInitialized_
ScalarType sol_norm
An estimate of the norm of the solution.
ScalarType mat_cond_num
An approximation to the condition number of A.
A pure virtual class for defining the status tests for the Belos iterative solvers.
Belos::OperatorTraits< ScalarType, MV, OP > OPT
Class which defines basic traits for the operator type.
Traits class which defines basic operations on multivectors.
void iterate()
This method performs LSQR iterations until the status test indicates the need to stop or an error occ...
ScalarType frob_mat_norm
An approximation to the Frobenius norm of A.
ScalarType frob_mat_norm_
ScalarType mat_resid_norm
An estimate of the norm of A^T*resid.
Teuchos::RCP< const MV > W
The search direction vector.
Teuchos::RCP< const MV > V
Bidiagonalization vector.
virtual ~LSQRIter()
Destructor.
A linear system to solve, and its associated information.
void setStateSize()
Method for initalizing the state storage needed by LSQR.
Class which describes the linear problem to be solved by the iterative solver.
int getNumIters() const
Get the current iteration count.
LSQRIter(const Teuchos::RCP< Belos::LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Belos::OutputManager< ScalarType > > &printer, const Teuchos::RCP< Belos::StatusTest< ScalarType, MV, OP > > &tester, Teuchos::ParameterList ¶ms)
LSQRIter constructor with linear problem, solver utilities, and parameter list of solver options...
Teuchos::ScalarTraits< ScalarType > SCT
const Teuchos::RCP< Belos::LinearProblem< ScalarType, MV, OP > > lp_
static magnitudeType magnitude(T a)
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
Belos::MultiVecTraits< ScalarType, MV > MVT
SCT::magnitudeType MagnitudeType
LSQRIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
void initialize()
The solver is initialized using initializeLSQR.
Class which defines basic traits for the operator type.
Implementation of the LSQR iteration.
ScalarType bnorm
The norm of the RHS vector b.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
bool isInitialized()
States whether the solver has been initialized or not.
Structure to contain pointers to LSQRIteration state variables, ...
Belos header file which uses auto-configuration information to include necessary C++ headers...
const Belos::LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
LSQRIterateFailure is thrown when the LSQRIteration object is unable to compute the next iterate in t...
void setBlockSize(int blockSize)
Set the blocksize to be used by the iterative solver to solve this linear problem.