42#ifndef BELOS_BLOCK_GCRODR_ITER_HPP
43#define BELOS_BLOCK_GCRODR_ITER_HPP
60#include "Teuchos_BLAS.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"
94 template <
class ScalarType,
class MV>
107 Teuchos::RCP<MV>
U,
C;
114 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
H;
118 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
B;
121 U(Teuchos::null),
C(Teuchos::null),
122 H(Teuchos::null),
B(Teuchos::null)
166 template<
class ScalarType,
class MV,
class OP>
175 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
177 typedef Teuchos::SerialDenseMatrix<int,ScalarType>
SDM;
178 typedef Teuchos::SerialDenseVector<int,ScalarType>
SDV;
195 Teuchos::ParameterList ¶ms );
321 if (!initialized_)
return 0;
349 void setSize(
int recycledBlocks,
int numBlocks ) {
351 if ( (recycledBlocks_ != recycledBlocks) || (numBlocks_ != numBlocks) ) {
352 recycledBlocks_ = recycledBlocks;
353 numBlocks_ = numBlocks;
354 cs_.sizeUninitialized( numBlocks_ );
355 sn_.sizeUninitialized( numBlocks_ );
356 Z_.shapeUninitialized( numBlocks_*blockSize_, blockSize_ );
372 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
373 const Teuchos::RCP<OutputManager<ScalarType> > om_;
374 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
375 const Teuchos::RCP<OrthoManager<ScalarType,MV> > ortho_;
383 int numBlocks_, blockSize_;
388 std::vector<bool> trueRHSIndices_;
395 Teuchos::SerialDenseVector<int,MagnitudeType> cs_;
402 std::vector< SDM >House_;
414 int curDim_, iter_, lclIter_;
425 Teuchos::RCP<MV> U_, C_;
434 Teuchos::RCP<SDM > H_;
439 Teuchos::RCP<SDM > B_;
447 Teuchos::RCP<SDM> R_;
462 template<
class ScalarType,
class MV,
class OP>
467 Teuchos::ParameterList ¶ms ):lp_(problem),
468 om_(printer), stest_(tester), ortho_(ortho) {
472 initialized_ =
false;
483 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter(
"Num Blocks"), std::invalid_argument,
"Belos::BlockGCRODRIter::constructor: mandatory parameter \"Num Blocks\" is not specified.");
484 int nb = Teuchos::getParameter<int>(params,
"Num Blocks");
486 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter(
"Recycled Blocks"), std::invalid_argument,
"Belos::BlockGCRODRIter::constructor: mandatory parameter \"Recycled Blocks\" is not specified.");
487 int rb = Teuchos::getParameter<int>(params,
"Recycled Blocks");
489 TEUCHOS_TEST_FOR_EXCEPTION(nb <= 0, std::invalid_argument,
"Belos::BlockGCRODRIter() was passed a non-positive argument for \"Num Blocks\".");
490 TEUCHOS_TEST_FOR_EXCEPTION(rb >= nb, std::invalid_argument,
"Belos::BlockGCRODRIter() the number of recycled blocks is larger than the allowable subspace.");
493 int bs = Teuchos::getParameter<int>(params,
"Block Size");
495 TEUCHOS_TEST_FOR_EXCEPTION(bs <= 0, std::invalid_argument,
"Belos::BlockGCRODRIter() the block size was passed a non-postitive argument.");
499 recycledBlocks_ = rb;
504 trueRHSIndices_.resize(blockSize_);
506 for(i=0; i<blockSize_; i++){
507 trueRHSIndices_[i] =
true;
512 cs_.sizeUninitialized( numBlocks_+1 );
513 sn_.sizeUninitialized( numBlocks_+1 );
514 Z_.shapeUninitialized( (numBlocks_+1)*blockSize_,blockSize_ );
516 House_.resize(numBlocks_);
518 for(i=0; i<numBlocks_;i++){
519 House_[i].shapeUninitialized(2*blockSize_, 2*blockSize_);
525 template <
class ScalarType,
class MV,
class OP>
527 TEUCHOS_TEST_FOR_EXCEPTION( initialized_ ==
false,
BlockGCRODRIterInitFailure,
"Belos::BlockGCRODRIter::iterate(): GCRODRIter class not initialized." );
533 setSize( recycledBlocks_, numBlocks_ );
535 Teuchos::RCP<MV> Vnext;
536 Teuchos::RCP<const MV> Vprev;
537 std::vector<int> curind(blockSize_);
543 for(
int i = 0; i<blockSize_; i++){curind[i] = i;};
548 Vnext = MVT::CloneViewNonConst(*V_,curind);
551 Teuchos::RCP<SDM > Z0 =
552 Teuchos::rcp(
new SDM(blockSize_,blockSize_) );
553 int rank = ortho_->normalize(*Vnext,Z0);
558 TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_,
BlockGCRODRIterOrthoFailure,
"Belos::BlockGCRODRIter::iterate(): couldn't generate basis of full rank at the initial step.");
560 Teuchos::RCP<SDM > Z_block = Teuchos::rcp(
new SDM(Teuchos::View, Z_, blockSize_,blockSize_) );
561 Z_block->assign(*Z0);
563 std::vector<int> prevind(blockSize_*(numBlocks_ + 1));
570 while( (stest_->checkStatus(
this) !=
Passed) && (curDim_+blockSize_-1) < (numBlocks_*blockSize_)) {
576 int HFirstCol = curDim_-blockSize_;
577 int HLastCol = HFirstCol + blockSize_-1 ;
578 int HLastOrthRow = HLastCol;
579 int HFirstNormRow = HLastOrthRow + 1;
583 for(
int i = 0; i< blockSize_; i++){
584 curind[i] = curDim_ + i;
586 Vnext = MVT::CloneViewNonConst(*V_,curind);
591 for(
int i = 0; i< blockSize_; i++){
592 curind[blockSize_ - 1 - i] = curDim_ - i - 1;
594 Vprev = MVT::CloneView(*V_,curind);
596 lp_->apply(*Vprev,*Vnext);
597 Vprev = Teuchos::null;
603 Teuchos::Array<Teuchos::RCP<const MV> > C(1, C_);
606 subB = Teuchos::rcp(
new SDM ( Teuchos::View,*B_,recycledBlocks_,blockSize_,0, HFirstCol ) );
608 Teuchos::Array<Teuchos::RCP<SDM > > AsubB;
609 AsubB.append( subB );
611 ortho_->project( *Vnext, AsubB, C );
615 prevind.resize(curDim_);
616 for (
int i=0; i<curDim_; i++) { prevind[i] = i; }
617 Vprev = MVT::CloneView(*V_,prevind);
618 Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
621 Teuchos::RCP<SDM> subH = Teuchos::rcp(
new SDM ( Teuchos::View,*H_,curDim_,blockSize_,0,HFirstCol ) );
622 Teuchos::Array<Teuchos::RCP<SDM > > AsubH;
623 AsubH.append( subH );
625 Teuchos::RCP<SDM > subR = Teuchos::rcp(
new SDM ( Teuchos::View,*H_,blockSize_,blockSize_,HFirstNormRow,HFirstCol ) );
627 int rank = ortho_->projectAndNormalize(*Vnext,AsubH,subR,AVprev);
630 SDM subR2( Teuchos::View,*R_,(lclIter_+1)*blockSize_,blockSize_,0,HFirstCol);
631 SDM subH2( Teuchos::View,*H_,(lclIter_+1)*blockSize_,blockSize_,0,HFirstCol);
634 TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_,
BlockGCRODRIterOrthoFailure,
"Belos::BlockGCRODRIter::iterate(): couldn't generate basis of full rank.");
640 curDim_ = curDim_ + blockSize_;
650 template <
class ScalarType,
class MV,
class OP>
652 if (newstate.
V != Teuchos::null && newstate.
H != Teuchos::null) {
653 curDim_ = newstate.
curDim;
660 R_ = Teuchos::rcp(
new SDM(H_->numRows(), H_->numCols() ));
664 SDM Identity(2*blockSize_, 2*blockSize_);
665 for(
int i=0;i<2*blockSize_; i++){
668 for(
int i=0; i<numBlocks_;i++){
669 House_[i].assign(Identity);
673 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
V == Teuchos::null,std::invalid_argument,
"Belos::GCRODRIter::initialize(): BlockGCRODRIterState does not have V initialized.");
674 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
H == Teuchos::null,std::invalid_argument,
"Belos::GCRODRIter::initialize(): BlockGCRODRIterState does not have H initialized.");
686 template <
class ScalarType,
class MV,
class OP>
687 Teuchos::RCP<const MV>
694 if (
static_cast<int> (norms->size()) < blockSize_) {
695 norms->resize( blockSize_ );
697 Teuchos::BLAS<int,ScalarType> blas;
698 for (
int j=0; j<blockSize_; j++) {
699 if(trueRHSIndices_[j]){
700 (*norms)[j] = blas.NRM2( blockSize_, &Z_(curDim_-blockSize_+j, j), 1);
706 return Teuchos::null;
709 return Teuchos::null;
715 template <
class ScalarType,
class MV,
class OP>
721 Teuchos::RCP<MV> currentUpdate = Teuchos::null;
723 if(curDim_<=blockSize_) {
724 return currentUpdate;
727 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
728 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
729 Teuchos::BLAS<int,ScalarType> blas;
730 currentUpdate = MVT::Clone( *V_, blockSize_ );
734 SDM Y( Teuchos::Copy, Z_, curDim_-blockSize_, blockSize_ );
735 Teuchos::RCP<SDM> Rtmp = Teuchos::rcp(
new SDM(Teuchos::View, *R_, curDim_, curDim_-blockSize_));
752 blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
753 Teuchos::NON_UNIT_DIAG, curDim_-blockSize_, blockSize_, one,
754 Rtmp->values(), Rtmp->stride(), Y.values(), Y.stride() );
761 std::vector<int> index(curDim_-blockSize_);
762 for (
int i=0; i<curDim_-blockSize_; i++ ) index[i] = i;
763 Teuchos::RCP<const MV> Vjp1 = MVT::CloneView( *V_, index );
764 MVT::MvTimesMatAddMv( one, *Vjp1, Y, zero, *currentUpdate );
771 if (U_ != Teuchos::null) {
772 SDM z(recycledBlocks_,blockSize_);
773 SDM subB( Teuchos::View, *B_, recycledBlocks_, curDim_-blockSize_ );
774 z.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, subB, Y, zero );
777 MVT::MvTimesMatAddMv( -one, *U_, z, one, *currentUpdate );
783 return currentUpdate;
786 template<
class ScalarType,
class MV,
class OP>
790 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
791 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
798 int curDim = curDim_;
799 if ( (dim >= curDim_) && (dim < getMaxSubspaceDim()) ){
803 Teuchos::BLAS<int, ScalarType> blas;
812 for (i=0; i<curDim-1; i++) {
816 blas.ROT( 1, &(*R_)(i,curDim-1), 1, &(*R_)(i+1, curDim-1), 1, &cs_[i], &sn_[i] );
821 blas.ROTG( &(*R_)(curDim-1,curDim-1), &(*R_)(curDim,curDim-1), &cs_[curDim-1], &sn_[curDim-1] );
822 (*R_)(curDim,curDim-1) = zero;
826 blas.ROT( 1, &Z_(curDim-1,0), 1, &Z_(curDim,0), 1, &cs_[curDim-1], &sn_[curDim-1] );
841 Teuchos::RCP< SDM > workmatrix = Teuchos::null;
842 Teuchos::RCP< SDV > workvec = Teuchos::null;
843 Teuchos::RCP<SDV> v_refl = Teuchos::null;
844 int R_colStart = curDim_-blockSize_;
845 Teuchos::RCP< SDM >Rblock = Teuchos::null;
850 for(i=0; i<lclIter_-1; i++){
851 int R_rowStart = i*blockSize_;
853 Teuchos::RCP< SDM > RblockCopy = rcp(
new SDM (Teuchos::Copy, *R_, 2*blockSize_,blockSize_, R_rowStart, R_colStart));
854 Teuchos::RCP< SDM > RblockView = rcp(
new SDM (Teuchos::View, *R_, 2*blockSize_,blockSize_, R_rowStart, R_colStart));
855 blas.GEMM(Teuchos::NO_TRANS,Teuchos::NO_TRANS, 2*blockSize_,blockSize_,2*blockSize_,one,House_[i].values(),House_[i].stride(), RblockCopy->values(),RblockCopy -> stride(), zero, RblockView->values(),RblockView -> stride());
862 Rblock = rcp(
new SDM (Teuchos::View, *R_, 2*blockSize_,blockSize_, curDim_-blockSize_, curDim_-blockSize_));
865 for(i=0; i<blockSize_; i++){
869 int curcol = (lclIter_ - 1)*blockSize_ + i;
871 ScalarType signDiag = (*R_)(curcol,curcol) / Teuchos::ScalarTraits<ScalarType>::magnitude((*R_)(curcol,curcol));
875 ScalarType nvs = blas.NRM2(blockSize_+1,&((*R_)[curcol][curcol]),1);
876 ScalarType alpha = -signDiag*nvs;
894 v_refl = rcp(
new SDV(Teuchos::Copy, &((*R_)(curcol,curcol)), blockSize_ + 1 ));
895 (*v_refl)[0] -= alpha;
896 nvs = blas.NRM2(blockSize_+1,v_refl -> values(),1);
912 if(i < blockSize_ - 1){
913 workvec = Teuchos::rcp(
new SDV(blockSize_ - i -1));
915 workmatrix = Teuchos::rcp(
new SDM (Teuchos::View, *Rblock, blockSize_+1, blockSize_ - i -1, lclCurcol, lclCurcol +1 ) );
916 blas.GEMV(Teuchos::TRANS, workmatrix->numRows(), workmatrix->numCols(), one, workmatrix->values(), workmatrix->stride(), v_refl->values(), 1, zero, workvec->values(), 1);
917 blas.GER(workmatrix->numRows(),workmatrix->numCols(), -2.0*one/nvs, v_refl->values(),1,workvec->values(),1,workmatrix->values(),workmatrix->stride());
924 workvec = Teuchos::rcp(
new SDV(2*blockSize_));
925 workmatrix = Teuchos::rcp(
new SDM (Teuchos::View, House_[lclIter_ -1], blockSize_+1, 2*blockSize_, i, 0 ) );
926 blas.GEMV(Teuchos::TRANS,workmatrix->numRows(),workmatrix->numCols(),one,workmatrix->values(),workmatrix->stride(), v_refl->values(), 1,zero,workvec->values(),1);
927 blas.GER(workmatrix->numRows(),workmatrix->numCols(), -2.0*one/nvs, v_refl -> values(),1,workvec->values(),1,workmatrix->values(),(*workmatrix).stride());
932 workvec = Teuchos::rcp(
new SDV(blockSize_));
933 workmatrix = Teuchos::rcp(
new SDM (Teuchos::View, Z_, blockSize_+1, blockSize_, curcol, 0 ) );
934 blas.GEMV(Teuchos::TRANS, workmatrix->numRows(), workmatrix->numCols(), one, workmatrix-> values(), workmatrix->stride(), v_refl -> values(), 1, zero, workvec->values(), 1);
935 blas.GER((*workmatrix).numRows(),(*workmatrix).numCols(), -2.0*one/nvs,v_refl -> values(), 1,&((*workvec)[0]),1,(*workmatrix)[0],(*workmatrix).stride());
940 (*R_)[curcol][curcol] = alpha;
941 for(
int ii=1; ii<= blockSize_; ii++){
942 (*R_)[curcol][curcol+ii] = 0;
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.
BlockGCRODRIterInitFailure is thrown when the BlockGCRODRIter object is unable to generate an initial...
BlockGCRODRIterInitFailure(const std::string &what_arg)
BlockGCRODRIterOrthoFailure is thrown when the BlockGCRODRIter object is unable to compute independen...
BlockGCRODRIterOrthoFailure(const std::string &what_arg)
Implementation of the Block GCRO-DR (Block Recycling GMRES) iteration.
BlockGCRODRIter(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)
BlockGCRODRIter constructor with linear problem, solver utilities, and parameter list of solver optio...
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
void setSize(int recycledBlocks, int numBlocks)
Set the maximum number of blocks used by the iterative solver and the number of recycled vectors.
bool isInitialized()
States whether the solver has been initialized or not.
int getNumBlocks() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem.
void initialize()
Initialize the solver to an iterate, providing a complete state.
Teuchos::SerialDenseVector< int, ScalarType > SDV
void setRecycledBlocks(int recycledBlocks)
Set the maximum number of recycled blocks used by the iterative solver.
int getRecycledBlocks() const
Set the maximum number of recycled blocks used by the iterative solver.
void setBlockSize(int blockSize)
Set the blocksize.
Teuchos::SerialDenseMatrix< int, ScalarType > SDM
MultiVecTraits< ScalarType, MV > MVT
BlockGCRODRIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
void setNumBlocks(int numBlocks)
Set the maximum number of blocks used by the iterative solver.
SCT::magnitudeType MagnitudeType
Teuchos::ScalarTraits< ScalarType > SCT
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
virtual ~BlockGCRODRIter()
Destructor.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current solution to the linear problem.
void iterate()
This method performs block GCRODR iterations until the status test indicates the need to stop or an e...
void resetNumIters(int iter=0)
Reset the iteration count.
void updateLSQR(int dim=-1)
OperatorTraits< ScalarType, MV, OP > OPT
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
int getNumIters() const
Get the current iteration count.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
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 BlockGCRODRIter state variables.
Teuchos::RCP< MV > V
The current Krylov basis.
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.
int curDim
The current dimension of the reduction.