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];
221 Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms )
const;
229 Teuchos::RCP<MV> getCurrentUpdate()
const;
235 void updateLSQR(
int dim = -1 );
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.");
268 void setNumBlocks(
int numBlocks);
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>
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] );
Teuchos::ScalarTraits< ScalarType > SCT
std::vector< Teuchos::RCP< MV > > V_
bool isInitialized()
States whether the solver has been initialized or not.
void resetNumIters(int iter=0)
Reset the iteration count.
PseudoBlockGmresIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
const Teuchos::RCP< OrthoManager< ScalarType, MV > > ortho_
const Teuchos::RCP< OutputManager< ScalarType > > om_
Sacado::MP::Vector< Storage > ScalarType
KOKKOS_INLINE_FUNCTION MaskedAssign< scalar > mask_assign(bool b, scalar *s)
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Teuchos::RCP< MV > U_vec_
std::vector< Teuchos::RCP< Teuchos::SerialDenseVector< int, MagnitudeType > > > cs_
int getNumBlocks() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem...
SCT::magnitudeType MagnitudeType
std::vector< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > R_
std::vector< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > H_
OperatorTraits< ScalarType, MV, OP > OPT
void setBlockSize(int blockSize)
Set the blocksize.
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 initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
SVT::magnitudeType breakDownTol_
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
std::vector< Teuchos::RCP< Teuchos::SerialDenseVector< int, ScalarType > > > Z_
Teuchos::RCP< MV > cur_block_sol_
MultiVecTraits< ScalarType, MV > MVT
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
Mask< MagnitudeType > lucky_breakdown_
Teuchos::ScalarTraits< typename Storage::value_type > SVT
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
std::vector< Teuchos::RCP< Teuchos::SerialDenseVector< int, ScalarType > > > sn_
int getNumIters() const
Get the current iteration count.