42 #ifndef BELOS_LSQR_ITER_HPP
43 #define BELOS_LSQR_ITER_HPP
74 template<
class ScalarType,
class MV,
class OP>
163 state.
bnorm = bnorm_;
204 "LSQRIter::setBlockSize(): Cannot use a block size that is not one.");
240 bool stateStorageInitialized_;
262 ScalarType resid_norm_;
265 ScalarType frob_mat_norm_;
268 ScalarType mat_cond_num_;
271 ScalarType mat_resid_norm_;
274 ScalarType sol_norm_;
286 template<
class ScalarType,
class MV,
class OP>
295 stateStorageInitialized_(false),
301 mat_resid_norm_(0.0),
309 template <
class ScalarType,
class MV,
class OP>
312 if (!stateStorageInitialized_) {
316 if (lhsMV == Teuchos::null || rhsMV == Teuchos::null) {
317 stateStorageInitialized_ =
false;
324 if (U_ == Teuchos::null) {
326 TEUCHOS_TEST_FOR_EXCEPTION(rhsMV == Teuchos::null, std::invalid_argument,
"LSQRIter::setStateSize(): linear problem does not specify right hand multivector to clone from.");
327 TEUCHOS_TEST_FOR_EXCEPTION(lhsMV == Teuchos::null, std::invalid_argument,
"LSQRIter::setStateSize(): linear problem does not specify left hand multivector to clone from.");
329 U_ = MVT::Clone( *rhsMV, 1 );
330 V_ = MVT::Clone( *lhsMV, 1 );
331 W_ = MVT::Clone( *lhsMV, 1 );
335 stateStorageInitialized_ =
true;
343 template <
class ScalarType,
class MV,
class OP>
349 if (!stateStorageInitialized_)
353 "LSQRIter::initialize(): Cannot initialize state storage!");
355 std::string errstr(
"LSQRIter::initialize(): Specified multivectors must have a consistent length and width.");
360 RCP<const MV> lhsMV = lp_->getLHS();
362 const bool debugSerialLSQR =
false;
364 if( debugSerialLSQR )
366 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> lhsNorm(1);
367 MVT::MvNorm( *lhsMV, lhsNorm );
368 std::cout <<
"initializeLSQR lhsNorm " << lhsNorm[0] << std::endl;
372 RCP<const MV> rhsMV = lp_->getInitPrecResVec();
375 MVT::MvAddMv( one, *rhsMV, zero, *U_, *U_);
377 RCP<const OP> M_left = lp_->getLeftPrec();
378 RCP<const OP> A = lp_->getOperator();
379 RCP<const OP> M_right = lp_->getRightPrec();
381 if( debugSerialLSQR )
383 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> rhsNorm(1);
384 MVT::MvNorm( *U_, rhsNorm );
385 std::cout <<
"initializeLSQR | U_ | : " << rhsNorm[0] << std::endl;
397 if ( M_left.is_null())
405 RCP<MV> tempInRangeOfA = MVT::CloneCopy (*U_);
406 OPT::Apply (*M_left, *U_, *tempInRangeOfA,
CONJTRANS);
407 OPT::Apply (*A, *tempInRangeOfA, *V_,
CONJTRANS);
410 if (! M_right.is_null())
412 RCP<MV> tempInDomainOfA = MVT::CloneCopy (*V_);
414 OPT::Apply (*M_right, *tempInDomainOfA, *V_,
CONJTRANS);
419 MVT::MvAddMv( one, *V_, zero, *V_, *W_);
421 frob_mat_norm_ = zero;
432 template <
class ScalarType,
class MV,
class OP>
438 if (initialized_ ==
false) {
448 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> alpha(1);
449 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> beta(1);
450 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> xi(1);
453 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> wnorm2(1);
454 ScalarType rhobar, phibar, cs1, phi, rho, cs, sn, theta, xxnorm = MTzero, common;
455 ScalarType zetabar, sn1, psi, res = zero, ddnorm = zero, gamma, tau;
456 ScalarType anorm2 = zero;
457 ScalarType cs2 = -one, sn2 = zero, gammabar, zeta = zero, delta;
461 AV = MVT::Clone( *U_, 1);
463 AtU = MVT::Clone( *V_, 1);
464 const bool debugSerialLSQR =
false;
472 "LSQRIter::iterate(): current linear system has more than one vector!" );
476 MVT::MvNorm( *U_, beta );
477 if( SCT::real(beta[0]) > zero )
479 MVT::MvScale( *U_, one / beta[0] );
484 MVT::MvScale( *V_, one / beta[0] );
489 MVT::MvNorm( *V_, alpha );
490 if( debugSerialLSQR )
494 std::cout << iter_ <<
" First alpha " << alpha[0] <<
" beta " << beta[0] <<
" lambda " << lambda_ << std::endl;
496 if( SCT::real(alpha[0]) > zero )
498 MVT::MvScale( *V_, one / alpha[0] );
500 if( beta[0] * alpha[0] > zero )
502 MVT::MvScale( *W_, one / (beta[0] * alpha[0]) );
506 MVT::MvScale( *W_, zero );
510 RCP<const OP> M_left = lp_->getLeftPrec();
511 RCP<const OP> A = lp_->getOperator();
512 RCP<const OP> M_right = lp_->getRightPrec();
517 resid_norm_ = beta[0];
518 mat_resid_norm_ = alpha[0] * beta[0];
532 if ( M_right.is_null() )
534 OPT::Apply(*A, *V_, *AV);
538 RCP<MV> tempInDomainOfA = MVT::CloneCopy (*V_);
539 OPT::Apply (*M_right, *V_, *tempInDomainOfA);
540 OPT::Apply(*A, *tempInDomainOfA, *AV);
543 if (! M_left.is_null())
545 RCP<MV> tempInRangeOfA = MVT::CloneCopy (*AV);
546 OPT::Apply (*M_left, *tempInRangeOfA, *AV);
550 if ( !( M_left.is_null() && M_right.is_null() )
551 && debugSerialLSQR && iter_ == 1)
557 if (! M_left.is_null())
559 RCP<MV> tempInRangeOfA = MVT::CloneCopy (*U_);
560 OPT::Apply (*M_left, *U_, *tempInRangeOfA,
CONJTRANS);
561 OPT::Apply (*A, *tempInRangeOfA, *AtU,
CONJTRANS);
567 if ( !( M_right.is_null() ) )
569 RCP<MV> tempInDomainOfA = MVT::CloneCopy (*AtU);
570 OPT::Apply (*M_right, *tempInDomainOfA, *AtU,
CONJTRANS);
573 MVT::MvAddMv( one, *AtU, -alpha[0], *V_, *AtU );
574 MVT::MvNorm( *AtU, xi );
575 std::cout <<
"| V alpha - A' u |= " << xi[0] << std::endl;
578 MVT::MvTransMv( one, *U_, *U_, uotuo );
579 std::cout <<
"<U, U> = " <<
printMat(uotuo) << std::endl;
581 std::cout <<
"alpha = " << alpha[0] << std::endl;
584 MVT::MvTransMv( one, *AV, *U_, utav );
585 std::cout <<
"<AV, U> = alpha = " <<
printMat(utav) << std::endl;
588 MVT::MvAddMv( one, *AV, -alpha[0], *U_, *U_ );
589 MVT::MvNorm( *U_, beta);
591 anorm2 += alpha[0]*alpha[0] + beta[0]*beta[0] + lambda_*lambda_;
594 if ( SCT::real(beta[0]) > zero )
597 MVT::MvScale( *U_, one / beta[0] );
599 if (M_left.is_null())
605 RCP<MV> tempInRangeOfA = MVT::CloneCopy (*U_);
606 OPT::Apply (*M_left, *U_, *tempInRangeOfA,
CONJTRANS);
607 OPT::Apply(*A, *tempInRangeOfA, *AtU,
CONJTRANS);
609 if (! M_right.is_null())
611 RCP<MV> tempInDomainOfA = MVT::CloneCopy (*AtU);
612 OPT::Apply (*M_right, *tempInDomainOfA, *AtU,
CONJTRANS);
615 MVT::MvAddMv( one, *AtU, -beta[0], *V_, *V_ );
616 MVT::MvNorm( *V_, alpha );
623 if ( SCT::real(alpha[0]) > zero )
625 MVT::MvScale( *V_, one / alpha[0] );
631 cs1 = rhobar / common;
632 sn1 = lambda_ / common;
634 phibar = cs1 * phibar;
641 theta = sn * alpha[0];
642 rhobar = -cs * alpha[0];
644 phibar = sn * phibar;
648 gammabar = -cs2 * rho;
649 zetabar = (phi - delta*zeta) / gammabar;
652 cs2 = gammabar / gamma;
654 zeta = (phi - delta*zeta) / gamma;
658 if ( M_right.is_null())
660 MVT::MvAddMv( phi / rho, *W_, one, *cur_soln_vec, *cur_soln_vec);
664 RCP<MV> tempInDomainOfA = MVT::CloneCopy (*W_);
665 OPT::Apply (*M_right, *W_, *tempInDomainOfA);
666 MVT::MvAddMv( phi / rho, *tempInDomainOfA, one, *cur_soln_vec, *cur_soln_vec);
669 MVT::MvNorm( *W_, wnorm2 );
670 ddnorm += (one / rho)*(one / rho) * wnorm2[0]*wnorm2[0];
671 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...
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.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
SerialBandDenseMatrixPrinter< OrdinalType, ScalarType > printMat(const SerialBandDenseMatrix< OrdinalType, ScalarType > &obj)
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.
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 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.
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
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.