42 #ifndef BELOS_LSQR_ITER_HPP 43 #define BELOS_LSQR_ITER_HPP 61 #include "Teuchos_SerialDenseMatrix.hpp" 62 #include "Teuchos_SerialDenseVector.hpp" 63 #include "Teuchos_ScalarTraits.hpp" 64 #include "Teuchos_ParameterList.hpp" 65 #include "Teuchos_TimeMonitor.hpp" 76 template<
class ScalarType,
class MV,
class OP>
86 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
101 Teuchos::ParameterList ¶ms );
165 state.
bnorm = bnorm_;
183 Teuchos::RCP<const MV>
getNativeResiduals( std::vector<MagnitudeType> *norms )
const {
return Teuchos::null; }
205 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
206 "LSQRIter::setBlockSize(): Cannot use a block size that is not one.");
225 const Teuchos::RCP<Belos::LinearProblem<ScalarType,MV,OP> > lp_;
226 const Teuchos::RCP<Belos::OutputManager<ScalarType> > om_;
227 const Teuchos::RCP<Belos::StatusTest<ScalarType,MV,OP> > stest_;
242 bool stateStorageInitialized_;
261 MagnitudeType lambda_;
264 ScalarType resid_norm_;
267 ScalarType frob_mat_norm_;
270 ScalarType mat_cond_num_;
273 ScalarType mat_resid_norm_;
276 ScalarType sol_norm_;
288 template<
class ScalarType,
class MV,
class OP>
292 Teuchos::ParameterList ¶ms ):
297 stateStorageInitialized_(false),
303 mat_resid_norm_(0.0),
311 template <
class ScalarType,
class MV,
class OP>
314 if (!stateStorageInitialized_) {
316 Teuchos::RCP<const MV> rhsMV = lp_->getInitPrecResVec();
317 Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
318 if (lhsMV == Teuchos::null || rhsMV == Teuchos::null) {
319 stateStorageInitialized_ =
false;
326 if (U_ == Teuchos::null) {
328 TEUCHOS_TEST_FOR_EXCEPTION(rhsMV == Teuchos::null, std::invalid_argument,
"LSQRIter::setStateSize(): linear problem does not specify right hand multivector to clone from.");
329 TEUCHOS_TEST_FOR_EXCEPTION(lhsMV == Teuchos::null, std::invalid_argument,
"LSQRIter::setStateSize(): linear problem does not specify left hand multivector to clone from.");
337 stateStorageInitialized_ =
true;
345 template <
class ScalarType,
class MV,
class OP>
351 if (!stateStorageInitialized_)
354 TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
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);
370 std::cout <<
"initializeLSQR lhsNorm " << lhsNorm[0] << std::endl;
374 RCP<const MV> rhsMV = lp_->getInitPrecResVec();
375 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
376 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
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);
387 std::cout <<
"initializeLSQR | U_ | : " << rhsNorm[0] << std::endl;
399 if ( M_left.is_null())
412 if (! M_right.is_null())
423 frob_mat_norm_ = zero;
434 template <
class ScalarType,
class MV,
class OP>
440 if (initialized_ ==
false) {
445 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
446 const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
447 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
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;
464 Teuchos::RCP<MV> AtU;
466 const bool debugSerialLSQR =
false;
469 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
474 "LSQRIter::iterate(): current linear system has more than one vector!" );
479 if( SCT::real(beta[0]) > zero )
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 )
502 if( beta[0] * alpha[0] > 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() )
545 if (! M_left.is_null())
552 if ( !( M_left.is_null() && M_right.is_null() )
553 && debugSerialLSQR && iter_ == 1)
559 if (! M_left.is_null())
569 if ( !( M_right.is_null() ) )
577 std::cout <<
"| V alpha - A' u |= " << xi[0] << std::endl;
579 Teuchos::SerialDenseMatrix<int,ScalarType> uotuo(1,1);
581 std::cout <<
"<U, U> = " << uotuo << std::endl;
583 std::cout <<
"alpha = " << alpha[0] << std::endl;
585 Teuchos::SerialDenseMatrix<int,ScalarType> utav(1,1);
587 std::cout <<
"<AV, U> = alpha = " << utav << std::endl;
593 anorm2 += alpha[0]*alpha[0] + beta[0]*beta[0] + lambda_*lambda_;
596 if ( SCT::real(beta[0]) > zero )
601 if (M_left.is_null())
611 if (! M_right.is_null())
625 if ( SCT::real(alpha[0]) > zero )
632 common = Teuchos::ScalarTraits< ScalarType >::squareroot(rhobar*rhobar + lambda_*lambda_);
633 cs1 = rhobar / common;
634 sn1 = lambda_ / common;
636 phibar = cs1 * phibar;
640 rho = Teuchos::ScalarTraits< ScalarType >::squareroot(rhobar*rhobar + lambda_*lambda_ + beta[0]*beta[0]);
643 theta = sn * alpha[0];
644 rhobar = -cs * alpha[0];
646 phibar = sn * phibar;
650 gammabar = -cs2 * rho;
651 zetabar = (phi - delta*zeta) / gammabar;
652 sol_norm_ = Teuchos::ScalarTraits< ScalarType >::squareroot(xxnorm + zetabar*zetabar);
653 gamma = Teuchos::ScalarTraits< ScalarType >::squareroot(gammabar*gammabar + theta*theta);
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);
668 MVT::MvAddMv( phi / rho, *tempInDomainOfA, one, *cur_soln_vec, *cur_soln_vec);
672 ddnorm += (one / rho)*(one / rho) * wnorm2[0]*wnorm2[0];
675 frob_mat_norm_ = Teuchos::ScalarTraits< ScalarType >::squareroot(anorm2);
676 mat_cond_num_ = frob_mat_norm_ * Teuchos::ScalarTraits< ScalarType >::squareroot(ddnorm);
678 resid_norm_ = Teuchos::ScalarTraits< ScalarType >::squareroot(phibar*phibar + res);
679 mat_resid_norm_ = alpha[0] * Teuchos::ScalarTraits< ScalarType >::magnitude(tau);
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...
Class which manages the output and verbosity of the Belos solvers.
ScalarType resid_norm
The current residual norm.
Teuchos::ScalarTraits< ScalarType >::magnitudeType lambda
The damping value.
void resetNumIters(int iter=0)
Reset the iteration count.
Pure virtual base class for defining the status testing capabilities of Belos.
Declaration of basic traits for the multivector type.
void initializeLSQR(LSQRIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, completing the initial state.
ScalarType sol_norm
An estimate of the norm of the solution.
ScalarType mat_cond_num
An approximation to the condition number of A.
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 .
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.
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 .
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.
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
ScalarType mat_resid_norm
An estimate of the norm of A^T*resid.
static void Apply(const OP &Op, const MV &x, MV &y, ETrans trans=NOTRANS)
Apply Op to x, putting the result into y.
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.
Class which describes the linear problem to be solved by the iterative solver.
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
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
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.
static void MvNorm(const MV &mv, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec, NormType type=TwoNorm)
Compute the 2-norm of each individual vector of mv. Upon return, normvec[i] holds the value of ...
Class which defines basic traits for the operator type.
Implementation of the LSQR iteration.
static void MvScale(MV &mv, const ScalarType alpha)
Scale each element of the vectors in mv with alpha.
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.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
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.