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 );
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_;
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.");
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_)
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);
369 MVT::MvNorm( *lhsMV, lhsNorm );
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();
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) {
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;
463 AV = MVT::Clone( *U_, 1);
464 Teuchos::RCP<MV> AtU;
465 AtU = MVT::Clone( *V_, 1);
466 const bool debugSerialLSQR =
false;
469 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
473 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1,
LSQRIterateFailure,
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;
579 Teuchos::SerialDenseMatrix<int,ScalarType> uotuo(1,1);
580 MVT::MvTransMv( one, *U_, *U_, uotuo );
581 std::cout <<
"<U, U> = " << uotuo << std::endl;
583 std::cout <<
"alpha = " << alpha[0] << std::endl;
585 Teuchos::SerialDenseMatrix<int,ScalarType> utav(1,1);
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] );
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);
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_ );
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...
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.
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.
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.
const Belos::LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
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...
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Teuchos::ScalarTraits< ScalarType > SCT
const Teuchos::RCP< Belos::LinearProblem< ScalarType, MV, OP > > lp_
int getNumIters() const
Get the current iteration count.
Belos::MultiVecTraits< ScalarType, MV > MVT
SCT::magnitudeType MagnitudeType
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.
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...
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.
LSQRIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.