42#ifndef BELOS_GCRODR_ITER_HPP
43#define BELOS_GCRODR_ITER_HPP
59#include "Teuchos_BLAS.hpp"
60#include "Teuchos_SerialDenseMatrix.hpp"
61#include "Teuchos_SerialDenseVector.hpp"
62#include "Teuchos_ScalarTraits.hpp"
63#include "Teuchos_ParameterList.hpp"
64#include "Teuchos_TimeMonitor.hpp"
87 template <
class ScalarType,
class MV>
99 Teuchos::RCP<MV>
U,
C;
106 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
H;
110 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
B;
113 U(Teuchos::null),
C(Teuchos::null),
114 H(Teuchos::null),
B(Teuchos::null)
153 template<
class ScalarType,
class MV,
class OP>
163 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
181 Teuchos::ParameterList ¶ms );
328 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
"Belos::GCRODRIter::setBlockSize(): Cannot use a block size that is not one.");
332 void setSize(
int recycledBlocks,
int numBlocks ) {
358 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >
lp_;
359 const Teuchos::RCP<OutputManager<ScalarType> >
om_;
360 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
stest_;
361 const Teuchos::RCP<OrthoManager<ScalarType,MV> >
ortho_;
373 Teuchos::SerialDenseVector<int,ScalarType>
sn_;
374 Teuchos::SerialDenseVector<int,MagnitudeType>
cs_;
398 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
H_;
401 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
B_;
406 Teuchos::SerialDenseMatrix<int,ScalarType>
R_;
407 Teuchos::SerialDenseVector<int,ScalarType>
z_;
412 template<
class ScalarType,
class MV,
class OP>
417 Teuchos::ParameterList ¶ms ):
418 lp_(problem), om_(printer), stest_(tester), ortho_(ortho) {
431 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter(
"Num Blocks"), std::invalid_argument,
"Belos::GCRODRIter::constructor: mandatory parameter \"Num Blocks\" is not specified.");
432 int nb = Teuchos::getParameter<int>(params,
"Num Blocks");
434 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter(
"Recycled Blocks"), std::invalid_argument,
"Belos::GCRODRIter::constructor: mandatory parameter \"Recycled Blocks\" is not specified.");
435 int rb = Teuchos::getParameter<int>(params,
"Recycled Blocks");
437 TEUCHOS_TEST_FOR_EXCEPTION(nb <= 0, std::invalid_argument,
"Belos::GCRODRIter() was passed a non-positive argument for \"Num Blocks\".");
438 TEUCHOS_TEST_FOR_EXCEPTION(rb >= nb, std::invalid_argument,
"Belos::GCRODRIter() the number of recycled blocks is larger than the allowable subspace.");
452 template <
class ScalarType,
class MV,
class OP>
458 Teuchos::RCP<MV> currentUpdate = Teuchos::null;
460 return currentUpdate;
462 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
463 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
464 Teuchos::BLAS<int,ScalarType> blas;
465 currentUpdate = MVT::Clone( *V_, 1 );
469 Teuchos::SerialDenseMatrix<int,ScalarType> y( Teuchos::Copy, z_, curDim_, 1 );
473 blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
474 Teuchos::NON_UNIT_DIAG, curDim_, 1, one,
475 R_.values(), R_.stride(), y.values(), y.stride() );
479 std::vector<int> index(curDim_);
480 for (
int i=0; i<curDim_; i++ ) index[i] = i;
481 Teuchos::RCP<const MV> Vjp1 = MVT::CloneView( *V_, index );
482 MVT::MvTimesMatAddMv( one, *Vjp1, y, zero, *currentUpdate );
486 if (U_ != Teuchos::null) {
487 Teuchos::SerialDenseMatrix<int,ScalarType> z(recycledBlocks_,1);
488 Teuchos::SerialDenseMatrix<int,ScalarType> subB( Teuchos::View, *B_, recycledBlocks_, curDim_ );
489 z.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, subB, y, zero );
490 MVT::MvTimesMatAddMv( -one, *U_, z, one, *currentUpdate );
493 return currentUpdate;
500 template <
class ScalarType,
class MV,
class OP>
505 if ( norms && (
int)norms->size()==0 )
509 Teuchos::BLAS<int,ScalarType> blas;
510 (*norms)[0] = blas.NRM2( 1, &z_(curDim_), 1);
512 return Teuchos::null;
519 template <
class ScalarType,
class MV,
class OP>
522 if (newstate.
V != Teuchos::null && newstate.
H != Teuchos::null) {
523 curDim_ = newstate.
curDim;
531 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
V == Teuchos::null,std::invalid_argument,
"Belos::GCRODRIter::initialize(): GCRODRIterState does not have V initialized.");
532 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
H == Teuchos::null,std::invalid_argument,
"Belos::GCRODRIter::initialize(): GCRODRIterState does not have H initialized.");
543 template <
class ScalarType,
class MV,
class OP>
546 TEUCHOS_TEST_FOR_EXCEPTION( initialized_ ==
false,
GCRODRIterInitFailure,
"Belos::GCRODRIter::iterate(): GCRODRIter class not initialized." );
549 setSize( recycledBlocks_, numBlocks_ );
551 Teuchos::RCP<MV> Vnext;
552 Teuchos::RCP<const MV> Vprev;
553 std::vector<int> curind(1);
560 Vnext = MVT::CloneViewNonConst(*V_,curind);
563 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > z0 =
564 Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(1,1) );
565 int rank = ortho_->normalize( *Vnext, z0 );
566 TEUCHOS_TEST_FOR_EXCEPTION(rank != 1,
GCRODRIterOrthoFailure,
"Belos::GCRODRIter::iterate(): couldn't generate basis of full rank.");
570 std::vector<int> prevind(numBlocks_+1);
577 if (U_ == Teuchos::null) {
578 while (stest_->checkStatus(
this) !=
Passed && curDim_+1 <= numBlocks_) {
580 int lclDim = curDim_ + 1;
584 Vnext = MVT::CloneViewNonConst(*V_,curind);
588 Vprev = MVT::CloneView(*V_,curind);
591 lp_->apply(*Vprev,*Vnext);
596 prevind.resize(lclDim);
597 for (
int i=0; i<lclDim; i++) { prevind[i] = i; }
598 Vprev = MVT::CloneView(*V_,prevind);
599 Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
602 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
603 subH = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View,*H_,lclDim,1,0,curDim_ ) );
604 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubH;
605 AsubH.append( subH );
608 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
609 subR = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View,*H_,1,1,lclDim,curDim_ ) );
612 rank = ortho_->projectAndNormalize(*Vnext,AsubH,subR,AVprev);
615 Teuchos::SerialDenseMatrix<int,ScalarType> subR2( Teuchos::View,R_,lclDim+1,1,0,curDim_ );
616 Teuchos::SerialDenseMatrix<int,ScalarType> subH2( Teuchos::View,*H_,lclDim+1,1,0,curDim_ );
619 TEUCHOS_TEST_FOR_EXCEPTION(rank != 1,
GCRODRIterOrthoFailure,
"Belos::GCRODRIter::iterate(): couldn't generate basis of full rank.");
629 while (stest_->checkStatus(
this) !=
Passed && curDim_+1 <= numBlocks_) {
631 int lclDim = curDim_ + 1;
635 Vnext = MVT::CloneViewNonConst(*V_,curind);
640 Vprev = MVT::CloneView(*V_,curind);
643 lp_->apply(*Vprev,*Vnext);
644 Vprev = Teuchos::null;
647 Teuchos::Array<Teuchos::RCP<const MV> > C(1, C_);
648 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
649 subB = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View,*B_,recycledBlocks_,1,0,curDim_ ) );
650 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubB;
651 AsubB.append( subB );
654 ortho_->project( *Vnext, AsubB, C );
658 prevind.resize(lclDim);
659 for (
int i=0; i<lclDim; i++) { prevind[i] = i; }
660 Vprev = MVT::CloneView(*V_,prevind);
661 Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
664 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > subH = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View,*H_,lclDim,1,0,curDim_ ) );
665 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubH;
666 AsubH.append( subH );
669 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > subR = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View,*H_,1,1,lclDim,curDim_ ) );
672 rank = ortho_->projectAndNormalize(*Vnext,AsubH,subR,AVprev);
675 Teuchos::SerialDenseMatrix<int,ScalarType> subR2( Teuchos::View,R_,lclDim+1,1,0,curDim_ );
676 Teuchos::SerialDenseMatrix<int,ScalarType> subH2( Teuchos::View,*H_,lclDim+1,1,0,curDim_ );
679 TEUCHOS_TEST_FOR_EXCEPTION(rank != 1,
GCRODRIterOrthoFailure,
"Belos::GCRODRIter::iterate(): couldn't generate basis of full rank.");
693 template<
class ScalarType,
class MV,
class OP>
697 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
702 int curDim = curDim_;
703 if ( (dim >= curDim_) && (dim < getMaxSubspaceDim()) )
706 Teuchos::BLAS<int, ScalarType> blas;
713 for (i=0; i<curDim; i++) {
717 blas.ROT( 1, &R_(i,curDim), 1, &R_(i+1, curDim), 1, &cs_[i], &sn_[i] );
722 blas.ROTG( &R_(curDim,curDim), &R_(curDim+1,curDim), &cs_[curDim], &sn_[curDim] );
723 R_(curDim+1,curDim) = zero;
727 blas.ROT( 1, &z_(curDim), 1, &z_(curDim+1), 1, &cs_[curDim], &sn_[curDim] );
Belos header file which uses auto-configuration information to include necessary C++ headers.
Class which describes the linear problem to be solved by the iterative solver.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Declaration of basic traits for the multivector type.
Class which defines basic traits for the operator type.
Class which manages the output and verbosity of the Belos solvers.
Pure virtual base class for defining the status testing capabilities of Belos.
Collection of types and exceptions used within the Belos solvers.
Parent class to all Belos exceptions.
GCRODRIterInitFailure is thrown when the GCRODRIter object is unable to generate an initial iterate i...
GCRODRIterInitFailure(const std::string &what_arg)
GCRODRIterOrthoFailure is thrown when the GCRODRIter object is unable to compute independent directio...
GCRODRIterOrthoFailure(const std::string &what_arg)
This class implements the GCRODR iteration, where a single-std::vector Krylov subspace is constructed...
Teuchos::SerialDenseMatrix< int, ScalarType > R_
Teuchos::SerialDenseVector< int, ScalarType > sn_
int getNumBlocks() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem.
virtual ~GCRODRIter()
Destructor.
void setRecycledBlocks(int recycledBlocks)
Set the maximum number of recycled blocks used by the iterative solver.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
void setBlockSize(int blockSize)
Set the blocksize.
const Teuchos::RCP< OrthoManager< ScalarType, MV > > ortho_
Teuchos::SerialDenseVector< int, ScalarType > z_
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current solution to the linear problem.
void setNumBlocks(int numBlocks)
Set the maximum number of blocks used by the iterative solver.
OperatorTraits< ScalarType, MV, OP > OPT
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
void resetNumIters(int iter=0)
Reset the iteration count.
GCRODRIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
int getRecycledBlocks() const
Get the maximum number of recycled blocks used by the iterative solver in solving this linear problem...
Teuchos::ScalarTraits< ScalarType > SCT
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
bool isInitialized()
States whether the solver has been initialized or not.
MultiVecTraits< ScalarType, MV > MVT
int getNumIters() const
Get the current iteration count.
GCRODRIter(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< OutputManager< ScalarType > > &printer, const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &tester, const Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > &ortho, Teuchos::ParameterList ¶ms)
GCRODRIter constructor with linear problem, solver utilities, and parameter list of solver options.
void initialize()
Initialize the solver with empty data. Calling this method will result in error, as GCRODRIter must b...
void iterate()
This method performs block Gmres iterations until the status test indicates the need to stop or an er...
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
void setSize(int recycledBlocks, int numBlocks)
Set the maximum number of blocks used by the iterative solver and the number of recycled vectors.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
void updateLSQR(int dim=-1)
Method for updating QR factorization of upper Hessenberg matrix.
SCT::magnitudeType MagnitudeType
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > H_
Teuchos::SerialDenseVector< int, MagnitudeType > cs_
const Teuchos::RCP< OutputManager< ScalarType > > om_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B_
A linear system to solve, and its associated information.
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
Traits class which defines basic operations on multivectors.
Class which defines basic traits for the operator type.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
A pure virtual class for defining the status tests for the Belos iterative solvers.
Structure to contain pointers to GCRODRIter state variables.
Teuchos::RCP< MV > U
The recycled subspace and its projection.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B
The projection of the Krylov subspace against the recycled subspace.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
Teuchos::RCP< MV > V
The current Krylov basis.
int curDim
The current dimension of the reduction.