42#ifndef BELOS_PSEUDO_BLOCK_GMRES_ITER_MP_VECTOR_HPP
43#define BELOS_PSEUDO_BLOCK_GMRES_ITER_MP_VECTOR_HPP
54#include "BelosPseudoBlockGmresIter.hpp"
71 template<
class Storage,
class MV,
class OP>
73 virtual public Iteration<Sacado::MP::Vector<Storage>, MV, OP> {
81 typedef MultiVecTraits<ScalarType,MV>
MVT;
82 typedef OperatorTraits<ScalarType,MV,OP>
OPT;
83 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
85 typedef Teuchos::ScalarTraits<typename Storage::value_type>
SVT;
99 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
100 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
101 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
102 Teuchos::ParameterList ¶ms );
156 void initialize(
const PseudoBlockGmresIterState<ScalarType,MV> & newstate);
163 PseudoBlockGmresIterState<ScalarType,MV> empty;
174 PseudoBlockGmresIterState<ScalarType,MV>
getState()
const {
175 PseudoBlockGmresIterState<ScalarType,MV> state;
176 state.curDim = curDim_;
177 state.V.resize(numRHS_);
178 state.H.resize(numRHS_);
179 state.Z.resize(numRHS_);
180 state.sn.resize(numRHS_);
181 state.cs.resize(numRHS_);
182 for (
int i=0; i<numRHS_; ++i) {
186 state.sn[i] = sn_[i];
187 state.cs[i] = cs_[i];
239 if (!initialized_)
return 0;
253 const LinearProblem<ScalarType,MV,OP>&
getProblem()
const {
return *lp_; }
260 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
261 "Belos::PseudoBlockGmresIter::setBlockSize(): Cannot use a block size that is not one.");
280 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >
lp_;
281 const Teuchos::RCP<OutputManager<ScalarType> >
om_;
282 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
stest_;
283 const Teuchos::RCP<OrthoManager<ScalarType,MV> >
ortho_;
297 std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > >
sn_;
298 std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,MagnitudeType> > >
cs_;
320 std::vector<Teuchos::RCP<MV> >
V_;
325 std::vector<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > >
H_;
330 std::vector<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > >
R_;
331 std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > >
Z_;
336#ifdef BELOS_TEUCHOS_TIME_MONITOR
337 Teuchos::RCP<Teuchos::Time> timerUpdateLSQR_, timerSolveLSQR_;
343 template<
class StorageType,
class MV,
class OP>
344 PseudoBlockGmresIter<Sacado::MP::Vector<StorageType>,MV,OP>::
345 PseudoBlockGmresIter(
const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
346 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
347 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
348 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
349 Teuchos::ParameterList ¶ms ):
359 breakDownTol_(params.get(
"Ensemble Breakdown Tolerance", 0.0))
362 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter(
"Num Blocks"), std::invalid_argument,
363 "Belos::PseudoBlockGmresIter::constructor: mandatory parameter 'Num Blocks' is not specified.");
364 int nb = Teuchos::getParameter<int>(params,
"Num Blocks");
371 template <
class StorageType,
class MV,
class OP>
372 void PseudoBlockGmresIter<Sacado::MP::Vector<StorageType>,MV,OP>::setNumBlocks (
int numBlocks)
377 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument,
"Belos::PseudoBlockGmresIter::setNumBlocks was passed a non-positive argument.");
379 numBlocks_ = numBlocks;
382 initialized_ =
false;
387 template <
class StorageType,
class MV,
class OP>
388 Teuchos::RCP<MV> PseudoBlockGmresIter<Sacado::MP::Vector<StorageType>,MV,OP>::getCurrentUpdate()
const
390#ifdef BELOS_TEUCHOS_TIME_MONITOR
391 Teuchos::TimeMonitor updateTimer( *timerSolveLSQR_ );
397 Teuchos::RCP<MV> currentUpdate = Teuchos::null;
399 return currentUpdate;
401 currentUpdate = MVT::Clone(*(V_[0]), numRHS_);
402 std::vector<int> index(1), index2(curDim_);
403 for (
int i=0; i<curDim_; ++i) {
406 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
407 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
408 Teuchos::BLAS<int,ScalarType> blas;
410 for (
int i=0; i<numRHS_; ++i) {
412 Teuchos::RCP<MV> cur_block_copy_vec = MVT::CloneViewNonConst( *currentUpdate, index );
416 Teuchos::SerialDenseVector<int,ScalarType>
y( Teuchos::Copy, Z_[i]->values(), curDim_ );
420 blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
421 Teuchos::NON_UNIT_DIAG, curDim_, 1, one,
422 H_[i]->values(), H_[i]->stride(),
y.values(),
y.stride() );
425 Teuchos::RCP<const MV> Vjp1 = MVT::CloneView( *V_[i], index2 );
426 MVT::MvTimesMatAddMv( one, *Vjp1, y, zero, *cur_block_copy_vec );
429 return currentUpdate;
436 template <
class StorageType,
class MV,
class OP>
437 Teuchos::RCP<const MV>
438 PseudoBlockGmresIter<Sacado::MP::Vector<StorageType>,MV,OP>::
439 getNativeResiduals (std::vector<MagnitudeType> *norms)
const
441 typedef typename Teuchos::ScalarTraits<ScalarType> STS;
447 if (
static_cast<int> (norms->size()) < numRHS_)
448 norms->resize (numRHS_);
450 Teuchos::BLAS<int, ScalarType> blas;
451 for (
int j = 0;
j < numRHS_; ++
j)
453 const ScalarType curNativeResid = (*Z_[
j])(curDim_);
454 (*norms)[
j] = STS::magnitude (curNativeResid);
457 return Teuchos::null;
461 template <
class StorageType,
class MV,
class OP>
463 PseudoBlockGmresIter<Sacado::MP::Vector<StorageType>,MV,OP>::
464 initialize (
const PseudoBlockGmresIterState<ScalarType,MV> & newstate)
470 this->numRHS_ = MVT::GetNumberVecs (*(lp_->getCurrLHSVec()));
472#ifdef BELOS_TEUCHOS_TIME_MONITOR
473 std::stringstream ss;
476 std::string updateLabel = ss.str() +
": Update LSQR";
477 timerUpdateLSQR_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
479 std::string solveLabel = ss.str() +
": Solve LSQR";
480 timerSolveLSQR_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
487 std::string errstr (
"Belos::PseudoBlockGmresIter::initialize(): "
488 "Specified multivectors must have a consistent "
489 "length and width.");
492 TEUCHOS_TEST_FOR_EXCEPTION((
int)newstate.V.size()==0 || (
int)newstate.Z.size()==0,
493 std::invalid_argument,
494 "Belos::PseudoBlockGmresIter::initialize(): "
495 "V and/or Z was not specified in the input state; "
496 "the V and/or Z arrays have length zero.");
507 RCP<const MV> lhsMV = lp_->getLHS();
508 RCP<const MV> rhsMV = lp_->getRHS();
512 RCP<const MV> vectorInBasisSpace = rhsMV.is_null() ? lhsMV : rhsMV;
515 TEUCHOS_TEST_FOR_EXCEPTION(vectorInBasisSpace.is_null(),
516 std::invalid_argument,
517 "Belos::PseudoBlockGmresIter::initialize(): "
518 "The linear problem to solve does not specify multi"
519 "vectors from which we can clone basis vectors. The "
520 "right-hand side(s), left-hand side(s), or both should "
525 TEUCHOS_TEST_FOR_EXCEPTION(newstate.curDim > numBlocks_+1,
526 std::invalid_argument,
528 curDim_ = newstate.curDim;
534 for (
int i=0; i<numRHS_; ++i) {
538 if (V_[i].is_null() || MVT::GetNumberVecs(*V_[i]) < numBlocks_ + 1) {
539 V_[i] = MVT::Clone (*vectorInBasisSpace, numBlocks_ + 1);
543 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.V[i]) != MVT::GetGlobalLength(*V_[i]),
544 std::invalid_argument, errstr );
545 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V[i]) < newstate.curDim,
546 std::invalid_argument, errstr );
551 int lclDim = MVT::GetNumberVecs(*newstate.V[i]);
552 if (newstate.V[i] != V_[i]) {
554 if (curDim_ == 0 && lclDim > 1) {
555 om_->stream(Warnings)
556 <<
"Belos::PseudoBlockGmresIter::initialize(): the solver was "
557 <<
"initialized with a kernel of " << lclDim
559 <<
"The block size however is only " << 1
561 <<
"The last " << lclDim - 1 <<
" vectors will be discarded."
564 std::vector<int> nevind (curDim_ + 1);
565 for (
int j = 0;
j < curDim_ + 1; ++
j)
568 RCP<const MV> newV = MVT::CloneView (*newstate.V[i], nevind);
569 RCP<MV> lclV = MVT::CloneViewNonConst( *V_[i], nevind );
570 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
571 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
572 MVT::MvAddMv (one, *newV, zero, *newV, *lclV);
575 lclV = Teuchos::null;
582 for (
int i=0; i<numRHS_; ++i) {
584 if (Z_[i] == Teuchos::null) {
585 Z_[i] = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,ScalarType>() );
587 if (Z_[i]->length() < numBlocks_+1) {
588 Z_[i]->shapeUninitialized(numBlocks_+1, 1);
592 TEUCHOS_TEST_FOR_EXCEPTION(newstate.Z[i]->numRows() < curDim_, std::invalid_argument, errstr);
595 if (newstate.Z[i] != Z_[i]) {
599 Teuchos::SerialDenseVector<int,ScalarType> newZ(Teuchos::View,newstate.Z[i]->values(),curDim_+1);
600 Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > lclZ;
601 lclZ = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,ScalarType>(Teuchos::View,Z_[i]->values(),curDim_+1) );
605 lclZ = Teuchos::null;
612 for (
int i=0; i<numRHS_; ++i) {
614 if (H_[i] == Teuchos::null) {
615 H_[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>() );
617 if (H_[i]->numRows() < numBlocks_+1 || H_[i]->numCols() < numBlocks_) {
618 H_[i]->shapeUninitialized(numBlocks_+1, numBlocks_);
622 if ((
int)newstate.H.size() == numRHS_) {
625 TEUCHOS_TEST_FOR_EXCEPTION((newstate.H[i]->numRows() < curDim_ || newstate.H[i]->numCols() < curDim_), std::invalid_argument,
626 "Belos::PseudoBlockGmresIter::initialize(): Specified Hessenberg matrices must have a consistent size to the current subspace dimension");
628 if (newstate.H[i] != H_[i]) {
631 Teuchos::SerialDenseMatrix<int,ScalarType> newH(Teuchos::View,*newstate.H[i],curDim_+1, curDim_);
632 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclH;
633 lclH = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*H_[i],curDim_+1, curDim_) );
637 lclH = Teuchos::null;
649 if ((
int)newstate.cs.size() == numRHS_ && (
int)newstate.sn.size() == numRHS_) {
650 for (
int i=0; i<numRHS_; ++i) {
651 if (cs_[i] != newstate.cs[i])
652 cs_[i] = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,MagnitudeType>(*newstate.cs[i]) );
653 if (sn_[i] != newstate.sn[i])
654 sn_[i] = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,ScalarType>(*newstate.sn[i]) );
659 for (
int i=0; i<numRHS_; ++i) {
660 if (cs_[i] == Teuchos::null)
661 cs_[i] = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,MagnitudeType>(numBlocks_+1) );
663 cs_[i]->resize(numBlocks_+1);
664 if (sn_[i] == Teuchos::null)
665 sn_[i] = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,ScalarType>(numBlocks_+1) );
667 sn_[i]->resize(numBlocks_+1);
689 template <
class StorageType,
class MV,
class OP>
690 void PseudoBlockGmresIter<Sacado::MP::Vector<StorageType>,MV,OP>::iterate()
695 if (initialized_ ==
false) {
699 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
700 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
703 int searchDim = numBlocks_;
708 std::vector<int> index(1);
709 std::vector<int> index2(1);
711 Teuchos::RCP<MV> U_vec = MVT::Clone( *V_[0], numRHS_ );
714 Teuchos::RCP<MV> AU_vec = MVT::Clone( *V_[0], numRHS_ );
716 for (
int i=0; i<numRHS_; ++i) {
718 Teuchos::RCP<const MV> tmp_vec = MVT::CloneView( *V_[i], index );
719 Teuchos::RCP<MV> U_vec_view = MVT::CloneViewNonConst( *U_vec, index2 );
720 MVT::MvAddMv( one, *tmp_vec, zero, *tmp_vec, *U_vec_view );
728 while (stest_->checkStatus(
this) != Passed && curDim_ < searchDim) {
734 lp_->apply( *U_vec, *AU_vec );
739 int num_prev = curDim_+1;
740 index.resize( num_prev );
741 for (
int i=0; i<num_prev; ++i) {
747 for (
int i=0; i<numRHS_; ++i) {
751 Teuchos::RCP<const MV> V_prev = MVT::CloneView( *V_[i], index );
752 Teuchos::Array< Teuchos::RCP<const MV> > V_array( 1, V_prev );
758 Teuchos::RCP<MV> V_new = MVT::CloneViewNonConst( *AU_vec, index2 );
763 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > h_new
764 = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>
765 ( Teuchos::View, *H_[i], num_prev, 1, 0, curDim_ ) );
766 Teuchos::Array< Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > h_array( 1, h_new );
768 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > r_new
769 = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>
770 ( Teuchos::View, *H_[i], 1, 1, num_prev, curDim_ ) );
776 ortho_->projectAndNormalize( *V_new, h_array, r_new, V_array );
782 index2[0] = curDim_+1;
783 Teuchos::RCP<MV> tmp_vec = MVT::CloneViewNonConst( *V_[i], index2 );
784 MVT::MvAddMv( one, *V_new, zero, *V_new, *tmp_vec );
790 Teuchos::RCP<MV> tmp_AU_vec = U_vec;
799#ifdef BELOS_TEUCHOS_TIME_MONITOR
800 Teuchos::TimeMonitor updateTimer( *timerUpdateLSQR_ );
831 template<
class StorageType,
class MV,
class OP>
832 void PseudoBlockGmresIter<Sacado::MP::Vector<StorageType>,MV,OP>::updateLSQR(
int dim )
837 int curDim = curDim_;
838 if (dim >= curDim_ && dim < getMaxSubspaceDim()) {
843 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
844 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
846 Teuchos::BLAS<int, ScalarType> blas;
848 for (i=0; i<numRHS_; ++i) {
854 for (
j=0;
j<curDim;
j++) {
858 blas.ROT( 1, &(*H_[i])(
j,curDim), 1, &(*H_[i])(
j+1, curDim), 1, &(*cs_[i])[
j], &(*sn_[i])[
j] );
866 lucky_breakdown_ = ((*Z_[i])(curDim) == zero);
868 auto lucky_breakdown_tmp = ((*H_[i])(curDim+1,curDim) == zero);
869 mask_assign(lucky_breakdown_,(*H_[i])(curDim,curDim)) = one;
870 lucky_breakdown_ = lucky_breakdown_ || lucky_breakdown_tmp;
872 blas.ROTG( &(*H_[i])(curDim,curDim), &(*H_[i])(curDim+1,curDim), &(*cs_[i])[curDim], &(*sn_[i])[curDim] );
873 (*H_[i])(curDim+1,curDim) = zero;
877 blas.ROT( 1, &(*Z_[i])(curDim), 1, &(*Z_[i])(curDim+1), 1, &(*cs_[i])[curDim], &(*sn_[i])[curDim] );
KOKKOS_INLINE_FUNCTION MaskedAssign< scalar > mask_assign(bool b, scalar *s)
std::vector< Teuchos::RCP< MV > > V_
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.
SVT::magnitudeType breakDownTol_
int getNumBlocks() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem.
PseudoBlockGmresIter(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)
PseudoBlockGmresIter constructor with linear problem, solver utilities, and parameter list of solver ...
void resetNumIters(int iter=0)
Reset the iteration count.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the "native" residual vectors.
std::vector< Teuchos::RCP< Teuchos::SerialDenseVector< int, ScalarType > > > Z_
Sacado::MP::Vector< Storage > ScalarType
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
std::vector< Teuchos::RCP< Teuchos::SerialDenseVector< int, ScalarType > > > sn_
void iterate()
This method performs block Gmres iterations until the status test indicates the need to stop or an er...
std::vector< Teuchos::RCP< Teuchos::SerialDenseVector< int, MagnitudeType > > > cs_
Teuchos::RCP< MV > cur_block_rhs_
SCT::magnitudeType MagnitudeType
std::vector< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > H_
std::vector< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > R_
PseudoBlockGmresIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
virtual ~PseudoBlockGmresIter()=default
Destructor.
const Teuchos::RCP< OrthoManager< ScalarType, MV > > ortho_
MultiVecTraits< ScalarType, MV > MVT
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
void setBlockSize(int blockSize)
Set the blocksize.
Mask< MagnitudeType > lucky_breakdown_
void setNumBlocks(int numBlocks)
Set the maximum number of blocks used by the iterative solver.
void initialize(const PseudoBlockGmresIterState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
Teuchos::ScalarTraits< ScalarType > SCT
void updateLSQR(int dim=-1)
Method for updating QR factorization of upper Hessenberg matrix.
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
int getNumIters() const
Get the current iteration count.
bool isInitialized()
States whether the solver has been initialized or not.
OperatorTraits< ScalarType, MV, OP > OPT
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
const Teuchos::RCP< OutputManager< ScalarType > > om_
Teuchos::RCP< MV > AU_vec_
Teuchos::ScalarTraits< typename Storage::value_type > SVT
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
This class implements the pseudo-block GMRES iteration, where a block Krylov subspace is constructed ...
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType ValueType * y