Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosPseudoBlockGmresIter.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Belos: Block Linear Solvers Package
4 //
5 // Copyright 2004-2016 NTESS and the Belos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef BELOS_PSEUDO_BLOCK_GMRES_ITER_HPP
11 #define BELOS_PSEUDO_BLOCK_GMRES_ITER_HPP
12 
17 #include "BelosConfigDefs.hpp"
18 #include "BelosTypes.hpp"
19 #include "BelosIteration.hpp"
20 
21 #include "BelosLinearProblem.hpp"
22 #include "BelosMatOrthoManager.hpp"
23 #include "BelosOutputManager.hpp"
24 #include "BelosStatusTest.hpp"
25 #include "BelosOperatorTraits.hpp"
26 #include "BelosMultiVecTraits.hpp"
27 
28 #include "Teuchos_BLAS.hpp"
31 #include "Teuchos_ScalarTraits.hpp"
33 #include "Teuchos_TimeMonitor.hpp"
34 
48 namespace Belos {
49 
51 
52 
57  template <class ScalarType, class MV>
59 
62 
67  int curDim;
69  std::vector<Teuchos::RCP<const MV> > V;
75  std::vector<Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > > H;
77  std::vector<Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > > R;
79  std::vector<Teuchos::RCP<const Teuchos::SerialDenseVector<int,ScalarType> > > Z;
81  std::vector<Teuchos::RCP<const Teuchos::SerialDenseVector<int,ScalarType> > > sn;
82  std::vector<Teuchos::RCP<const Teuchos::SerialDenseVector<int,MagnitudeType> > > cs;
83 
85  H(0), R(0), Z(0),
86  sn(0), cs(0)
87  {}
88  };
89 
91 
92 
100  PseudoBlockGmresIterOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
101  {}};
102 
104 
105 
106  template<class ScalarType, class MV, class OP>
107  class PseudoBlockGmresIter : virtual public Iteration<ScalarType,MV,OP> {
108 
109  public:
110 
111  //
112  // Convenience typedefs
113  //
118 
120 
121 
131  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
134  Teuchos::ParameterList &params );
135 
137  virtual ~PseudoBlockGmresIter() {};
139 
140 
142 
143 
165  void iterate();
166 
189 
193  void initialize()
194  {
196  initialize(empty);
197  }
198 
208  state.curDim = curDim_;
209  state.V.resize(numRHS_);
210  state.H.resize(numRHS_);
211  state.Z.resize(numRHS_);
212  state.sn.resize(numRHS_);
213  state.cs.resize(numRHS_);
214  for (int i=0; i<numRHS_; ++i) {
215  state.V[i] = V_[i];
216  state.H[i] = H_[i];
217  state.Z[i] = Z_[i];
218  state.sn[i] = sn_[i];
219  state.cs[i] = cs_[i];
220  }
221  return state;
222  }
223 
225 
226 
228 
229 
231  int getNumIters() const { return iter_; }
232 
234  void resetNumIters( int iter = 0 ) { iter_ = iter; }
235 
253  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
254 
256 
262 
264 
267  void updateLSQR( int dim = -1 );
268 
270  int getCurSubspaceDim() const {
271  if (!initialized_) return 0;
272  return curDim_;
273  };
274 
276  int getMaxSubspaceDim() const { return numBlocks_; }
277 
279 
280 
282 
283 
285  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
286 
288  int getBlockSize() const { return 1; }
289 
291  void setBlockSize(int blockSize) {
292  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
293  "Belos::PseudoBlockGmresIter::setBlockSize(): Cannot use a block size that is not one.");
294  }
295 
297  int getNumBlocks() const { return numBlocks_; }
298 
300  void setNumBlocks(int numBlocks);
301 
303  bool isInitialized() { return initialized_; }
304 
306 
307  private:
308 
309  //
310  // Classes inputed through constructor that define the linear problem to be solved.
311  //
316 
317  //
318  // Algorithmic parameters
319  //
320  // numRHS_ is the current number of linear systems being solved.
321  int numRHS_;
322  // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
323  int numBlocks_;
324 
325  // Storage for QR factorization of the least squares system.
326  std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > > sn_;
327  std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,MagnitudeType> > > cs_;
328 
329  // Pointers to a work vector used to improve aggregate performance.
330  Teuchos::RCP<MV> U_vec_, AU_vec_;
331 
332  // Pointers to the current right-hand side and solution multivecs being solved for.
333  Teuchos::RCP<MV> cur_block_rhs_, cur_block_sol_;
334 
335  //
336  // Current solver state
337  //
338  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
339  // is capable of running; _initialize is controlled by the initialize() member method
340  // For the implications of the state of initialized_, please see documentation for initialize()
341  bool initialized_;
342 
343  // Current subspace dimension, and number of iterations performed.
344  int curDim_, iter_;
345 
346  //
347  // State Storage
348  //
349  std::vector<Teuchos::RCP<MV> > V_;
350  //
351  // Projected matrices
352  // H_ : Projected matrix from the Krylov factorization AV = VH + FE^T
353  //
354  std::vector<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > H_;
355  //
356  // QR decomposition of Projected matrices for solving the least squares system HY = B.
357  // R_: Upper triangular reduction of H
358  // Z_: Q applied to right-hand side of the least squares system
359  std::vector<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > R_;
360  std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > > Z_;
361  };
362 
364  // Constructor.
365  template<class ScalarType, class MV, class OP>
367  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
370  Teuchos::ParameterList &params ):
371  lp_(problem),
372  om_(printer),
373  stest_(tester),
374  ortho_(ortho),
375  numRHS_(0),
376  numBlocks_(0),
377  initialized_(false),
378  curDim_(0),
379  iter_(0)
380  {
381  // Get the maximum number of blocks allowed for each Krylov subspace
382  TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Num Blocks"), std::invalid_argument,
383  "Belos::PseudoBlockGmresIter::constructor: mandatory parameter 'Num Blocks' is not specified.");
384  int nb = Teuchos::getParameter<int>(params, "Num Blocks");
385 
386  setNumBlocks( nb );
387  }
388 
390  // Set the block size and make necessary adjustments.
391  template <class ScalarType, class MV, class OP>
393  {
394  // This routine only allocates space; it doesn't not perform any computation
395  // any change in size will invalidate the state of the solver.
396 
397  TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument, "Belos::PseudoBlockGmresIter::setNumBlocks was passed a non-positive argument.");
398 
399  numBlocks_ = numBlocks;
400  curDim_ = 0;
401 
402  initialized_ = false;
403  }
404 
406  // Get the current update from this subspace.
407  template <class ScalarType, class MV, class OP>
409  {
410  //
411  // If this is the first iteration of the Arnoldi factorization,
412  // there is no update, so return Teuchos::null.
413  //
414  Teuchos::RCP<MV> currentUpdate = Teuchos::null;
415  if (curDim_==0) {
416  return currentUpdate;
417  } else {
418  currentUpdate = MVT::Clone(*(V_[0]), numRHS_);
419  std::vector<int> index(1), index2(curDim_);
420  for (int i=0; i<curDim_; ++i) {
421  index2[i] = i;
422  }
423  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
424  const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
426 
427  for (int i=0; i<numRHS_; ++i) {
428  index[0] = i;
429  Teuchos::RCP<MV> cur_block_copy_vec = MVT::CloneViewNonConst( *currentUpdate, index );
430  //
431  // Make a view and then copy the RHS of the least squares problem. DON'T OVERWRITE IT!
432  //
433  Teuchos::SerialDenseVector<int,ScalarType> y( Teuchos::Copy, Z_[i]->values(), curDim_ );
434  //
435  // Solve the least squares problem and compute current solutions.
436  //
438  Teuchos::NON_UNIT_DIAG, curDim_, 1, one,
439  H_[i]->values(), H_[i]->stride(), y.values(), y.stride() );
440 
441  Teuchos::RCP<const MV> Vjp1 = MVT::CloneView( *V_[i], index2 );
442  MVT::MvTimesMatAddMv( one, *Vjp1, y, zero, *cur_block_copy_vec );
443  }
444  }
445  return currentUpdate;
446  }
447 
448 
450  // Get the native residuals stored in this iteration.
451  // Note: No residual vector will be returned by Gmres.
452  template <class ScalarType, class MV, class OP>
455  getNativeResiduals (std::vector<MagnitudeType> *norms) const
456  {
457  typedef typename Teuchos::ScalarTraits<ScalarType> STS;
458 
459  if (norms)
460  { // Resize the incoming std::vector if necessary. The type
461  // cast avoids the compiler warning resulting from a signed /
462  // unsigned integer comparison.
463  if (static_cast<int> (norms->size()) < numRHS_)
464  norms->resize (numRHS_);
465 
467  for (int j = 0; j < numRHS_; ++j)
468  {
469  const ScalarType curNativeResid = (*Z_[j])(curDim_);
470  (*norms)[j] = STS::magnitude (curNativeResid);
471  }
472  }
473  return Teuchos::null;
474  }
475 
476 
477  template <class ScalarType, class MV, class OP>
478  void
481  {
482  using Teuchos::RCP;
483 
484  // (Re)set the number of right-hand sides, by interrogating the
485  // current LinearProblem to solve.
486  this->numRHS_ = MVT::GetNumberVecs (*(lp_->getCurrLHSVec()));
487 
488  // NOTE: In PseudoBlockGmresIter, V and Z are required!!!
489  // Inconsistent multivectors widths and lengths will not be tolerated, and
490  // will be treated with exceptions.
491  //
492  std::string errstr ("Belos::PseudoBlockGmresIter::initialize(): "
493  "Specified multivectors must have a consistent "
494  "length and width.");
495 
496  // Check that newstate has V and Z arrays with nonzero length.
497  TEUCHOS_TEST_FOR_EXCEPTION((int)newstate.V.size()==0 || (int)newstate.Z.size()==0,
498  std::invalid_argument,
499  "Belos::PseudoBlockGmresIter::initialize(): "
500  "V and/or Z was not specified in the input state; "
501  "the V and/or Z arrays have length zero.");
502 
503  // In order to create basis multivectors, we have to clone them
504  // from some already existing multivector. We require that at
505  // least one of the right-hand side B and left-hand side X in the
506  // LinearProblem be non-null. Thus, we can clone from either B or
507  // X. We prefer to close from B, since B is in the range of the
508  // operator A and the basis vectors should also be in the range of
509  // A (the first basis vector is a scaled residual vector).
510  // However, if B is not given, we will try our best by cloning
511  // from X.
512  RCP<const MV> lhsMV = lp_->getLHS();
513  RCP<const MV> rhsMV = lp_->getRHS();
514 
515  // If the right-hand side is null, we make do with the left-hand
516  // side, otherwise we use the right-hand side.
517  RCP<const MV> vectorInBasisSpace = rhsMV.is_null() ? lhsMV : rhsMV;
518  //RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
519 
520  TEUCHOS_TEST_FOR_EXCEPTION(vectorInBasisSpace.is_null(),
521  std::invalid_argument,
522  "Belos::PseudoBlockGmresIter::initialize(): "
523  "The linear problem to solve does not specify multi"
524  "vectors from which we can clone basis vectors. The "
525  "right-hand side(s), left-hand side(s), or both should "
526  "be nonnull.");
527 
528  // Check the new dimension is not more that the maximum number of
529  // allowable blocks.
530  TEUCHOS_TEST_FOR_EXCEPTION(newstate.curDim > numBlocks_+1,
531  std::invalid_argument,
532  errstr);
533  curDim_ = newstate.curDim;
534 
535  // Initialize the state storage. If the subspace has not be
536  // initialized before, generate it using the right-hand side or
537  // left-hand side from the LinearProblem lp_ to solve.
538  V_.resize(numRHS_);
539  for (int i=0; i<numRHS_; ++i) {
540  // Create a new vector if we need to. We "need to" if the
541  // current vector V_[i] is null, or if it doesn't have enough
542  // columns.
543  if (V_[i].is_null() || MVT::GetNumberVecs(*V_[i]) < numBlocks_ + 1) {
544  V_[i] = MVT::Clone (*vectorInBasisSpace, numBlocks_ + 1);
545  }
546  // Check that the newstate vector newstate.V[i] has dimensions
547  // consistent with those of V_[i].
548  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.V[i]) != MVT::GetGlobalLength(*V_[i]),
549  std::invalid_argument, errstr );
550  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V[i]) < newstate.curDim,
551  std::invalid_argument, errstr );
552  //
553  // If newstate.V[i] and V_[i] are not identically the same
554  // vector, then copy newstate.V[i] into V_[i].
555  //
556  int lclDim = MVT::GetNumberVecs(*newstate.V[i]);
557  if (newstate.V[i] != V_[i]) {
558  // Only copy over the first block and print a warning.
559  if (curDim_ == 0 && lclDim > 1) {
560  om_->stream(Warnings)
561  << "Belos::PseudoBlockGmresIter::initialize(): the solver was "
562  << "initialized with a kernel of " << lclDim
563  << std::endl
564  << "The block size however is only " << 1
565  << std::endl
566  << "The last " << lclDim - 1 << " vectors will be discarded."
567  << std::endl;
568  }
569  std::vector<int> nevind (curDim_ + 1);
570  for (int j = 0; j < curDim_ + 1; ++j)
571  nevind[j] = j;
572 
573  RCP<const MV> newV = MVT::CloneView (*newstate.V[i], nevind);
574  RCP<MV> lclV = MVT::CloneViewNonConst( *V_[i], nevind );
575  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
576  const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
577  MVT::MvAddMv (one, *newV, zero, *newV, *lclV);
578 
579  // Done with local pointers
580  lclV = Teuchos::null;
581  }
582  }
583 
584 
585  // Check size of Z
586  Z_.resize(numRHS_);
587  for (int i=0; i<numRHS_; ++i) {
588  // Create a vector if we need to.
589  if (Z_[i] == Teuchos::null) {
591  }
592  if (Z_[i]->length() < numBlocks_+1) {
593  Z_[i]->shapeUninitialized(numBlocks_+1, 1);
594  }
595 
596  // Check that the newstate vector is consistent.
597  TEUCHOS_TEST_FOR_EXCEPTION(newstate.Z[i]->numRows() < curDim_, std::invalid_argument, errstr);
598 
599  // Put data into Z_, make sure old information is not still hanging around.
600  if (newstate.Z[i] != Z_[i]) {
601  if (curDim_==0)
602  Z_[i]->putScalar();
603 
604  Teuchos::SerialDenseVector<int,ScalarType> newZ(Teuchos::View,newstate.Z[i]->values(),curDim_+1);
606  lclZ = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(Teuchos::View,Z_[i]->values(),curDim_+1) );
607  lclZ->assign(newZ);
608 
609  // Done with local pointers
610  lclZ = Teuchos::null;
611  }
612  }
613 
614 
615  // Check size of H
616  H_.resize(numRHS_);
617  for (int i=0; i<numRHS_; ++i) {
618  // Create a matrix if we need to.
619  if (H_[i] == Teuchos::null) {
621  }
622  if (H_[i]->numRows() < numBlocks_+1 || H_[i]->numCols() < numBlocks_) {
623  H_[i]->shapeUninitialized(numBlocks_+1, numBlocks_);
624  }
625 
626  // Put data into H_ if it exists, make sure old information is not still hanging around.
627  if ((int)newstate.H.size() == numRHS_) {
628 
629  // Check that the newstate matrix is consistent.
630  TEUCHOS_TEST_FOR_EXCEPTION((newstate.H[i]->numRows() < curDim_ || newstate.H[i]->numCols() < curDim_), std::invalid_argument,
631  "Belos::PseudoBlockGmresIter::initialize(): Specified Hessenberg matrices must have a consistent size to the current subspace dimension");
632 
633  if (newstate.H[i] != H_[i]) {
634  // H_[i]->putScalar();
635 
636  Teuchos::SerialDenseMatrix<int,ScalarType> newH(Teuchos::View,*newstate.H[i],curDim_+1, curDim_);
638  lclH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*H_[i],curDim_+1, curDim_) );
639  lclH->assign(newH);
640 
641  // Done with local pointers
642  lclH = Teuchos::null;
643  }
644  }
645  }
646 
648  // Reinitialize storage for least squares solve
649  //
650  cs_.resize(numRHS_);
651  sn_.resize(numRHS_);
652 
653  // Copy over rotation angles if they exist
654  if ((int)newstate.cs.size() == numRHS_ && (int)newstate.sn.size() == numRHS_) {
655  for (int i=0; i<numRHS_; ++i) {
656  if (cs_[i] != newstate.cs[i])
657  cs_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,MagnitudeType>(*newstate.cs[i]) );
658  if (sn_[i] != newstate.sn[i])
659  sn_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(*newstate.sn[i]) );
660  }
661  }
662 
663  // Resize or create the vectors as necessary
664  for (int i=0; i<numRHS_; ++i) {
665  if (cs_[i] == Teuchos::null)
666  cs_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,MagnitudeType>(numBlocks_+1) );
667  else
668  cs_[i]->resize(numBlocks_+1);
669  if (sn_[i] == Teuchos::null)
670  sn_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(numBlocks_+1) );
671  else
672  sn_[i]->resize(numBlocks_+1);
673  }
674 
675  // the solver is initialized
676  initialized_ = true;
677 
678  }
679 
680 
682  // Iterate until the status test informs us we should stop.
683  template <class ScalarType, class MV, class OP>
685  {
686  //
687  // Allocate/initialize data structures
688  //
689  if (initialized_ == false) {
690  initialize();
691  }
692 
693  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
694  const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
695 
696  // Compute the current search dimension.
697  int searchDim = numBlocks_;
698  //
699  // Associate each initial block of V_[i] with U_vec[i]
700  // Reset the index vector (this might have been changed if there was a restart)
701  //
702  std::vector<int> index(1);
703  std::vector<int> index2(1);
704  index[0] = curDim_;
705  Teuchos::RCP<MV> U_vec = MVT::Clone( *V_[0], numRHS_ );
706 
707  // Create AU_vec to hold A*U_vec.
708  Teuchos::RCP<MV> AU_vec = MVT::Clone( *V_[0], numRHS_ );
709 
710  for (int i=0; i<numRHS_; ++i) {
711  index2[0] = i;
712  Teuchos::RCP<const MV> tmp_vec = MVT::CloneView( *V_[i], index );
713  Teuchos::RCP<MV> U_vec_view = MVT::CloneViewNonConst( *U_vec, index2 );
714  MVT::MvAddMv( one, *tmp_vec, zero, *tmp_vec, *U_vec_view );
715  }
716 
718  // iterate until the status test tells us to stop.
719  //
720  // also break if our basis is full
721  //
722  while (stest_->checkStatus(this) != Passed && curDim_ < searchDim) {
723 
724  iter_++;
725  //
726  // Apply the operator to _work_vector
727  //
728  lp_->apply( *U_vec, *AU_vec );
729  //
730  //
731  // Resize index.
732  //
733  int num_prev = curDim_+1;
734  index.resize( num_prev );
735  for (int i=0; i<num_prev; ++i) {
736  index[i] = i;
737  }
738  //
739  // Orthogonalize next Krylov vector for each right-hand side.
740  //
741  for (int i=0; i<numRHS_; ++i) {
742  //
743  // Get previous Krylov vectors.
744  //
745  Teuchos::RCP<const MV> V_prev = MVT::CloneView( *V_[i], index );
746  Teuchos::Array< Teuchos::RCP<const MV> > V_array( 1, V_prev );
747  //
748  // Get a view of the new candidate std::vector.
749  //
750  index2[0] = i;
751  Teuchos::RCP<MV> V_new = MVT::CloneViewNonConst( *AU_vec, index2 );
752  //
753  // Get a view of the current part of the upper-hessenberg matrix.
754  //
757  ( Teuchos::View, *H_[i], num_prev, 1, 0, curDim_ ) );
759 
762  ( Teuchos::View, *H_[i], 1, 1, num_prev, curDim_ ) );
763  //
764  // Orthonormalize the new block of the Krylov expansion
765  // NOTE: Rank deficiencies are not checked because this is a single-std::vector Krylov method.
766  //
767  ortho_->projectAndNormalize( *V_new, h_array, r_new, V_array );
768  //
769  // NOTE: V_new is a copy of the iter+1 vector in V_[i], so the normalized vector has to be
770  // be copied back in when V_new is changed.
771  //
772  index2[0] = curDim_+1;
773  Teuchos::RCP<MV> tmp_vec = MVT::CloneViewNonConst( *V_[i], index2 );
774  MVT::MvAddMv( one, *V_new, zero, *V_new, *tmp_vec );
775  }
776  //
777  // Now _AU_vec is the new _U_vec, so swap these two vectors.
778  // NOTE: This alleviates the need for allocating a vector for AU_vec each iteration.
779  //
780  Teuchos::RCP<MV> tmp_AU_vec = U_vec;
781  U_vec = AU_vec;
782  AU_vec = tmp_AU_vec;
783  //
784  // V has been extended, and H has been extended.
785  //
786  // Update the QR factorization of the upper Hessenberg matrix
787  //
788  updateLSQR();
789  //
790  // Update basis dim and release all pointers.
791  //
792  curDim_ += 1;
793  //
794  } // end while (statusTest == false)
795 
796  }
797 
799  // Update the least squares solution for each right-hand side.
800  template<class ScalarType, class MV, class OP>
802  {
803  // Get correct dimension based on input "dim"
804  // Remember that ortho failures result in an exit before updateLSQR() is called.
805  // Therefore, it is possible that dim == curDim_.
806  int curDim = curDim_;
807  if (dim >= curDim_ && dim < getMaxSubspaceDim()) {
808  curDim = dim;
809  }
810 
811  int i, j;
812  const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
813 
815 
816  for (i=0; i<numRHS_; ++i) {
817  //
818  // Update the least-squares QR for each linear system.
819  //
820  // QR factorization of Least-Squares system with Givens rotations
821  //
822  for (j=0; j<curDim; j++) {
823  //
824  // Apply previous Givens rotations to new column of Hessenberg matrix
825  //
826  blas.ROT( 1, &(*H_[i])(j,curDim), 1, &(*H_[i])(j+1, curDim), 1, &(*cs_[i])[j], &(*sn_[i])[j] );
827  }
828  //
829  // Calculate new Givens rotation
830  //
831  blas.ROTG( &(*H_[i])(curDim,curDim), &(*H_[i])(curDim+1,curDim), &(*cs_[i])[curDim], &(*sn_[i])[curDim] );
832  (*H_[i])(curDim+1,curDim) = zero;
833  //
834  // Update RHS w/ new transformation
835  //
836  blas.ROT( 1, &(*Z_[i])(curDim), 1, &(*Z_[i])(curDim+1), 1, &(*cs_[i])[curDim], &(*sn_[i])[curDim] );
837  }
838 
839  } // end updateLSQR()
840 
841 } // end Belos namespace
842 
843 #endif /* BELOS_PSEUDO_BLOCK_GMRES_ITER_HPP */
ScalarType * values() const
std::vector< Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > > H
The current Hessenberg matrix.
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current solution to the linear problem...
Collection of types and exceptions used within the Belos solvers.
void TRSM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType &m, const OrdinalType &n, const alpha_type alpha, const A_type *A, const OrdinalType &lda, ScalarType *B, const OrdinalType &ldb) const
Belos&#39;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.
bool is_null(const boost::shared_ptr< T > &p)
PseudoBlockGmresIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
std::vector< Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > > R
The current upper-triangular matrix from the QR reduction of H.
Teuchos::ScalarTraits< ScalarType > SCT
void ROT(const OrdinalType &n, ScalarType *dx, const OrdinalType &incx, ScalarType *dy, const OrdinalType &incy, MagnitudeType *c, ScalarType *s) const
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.
int getNumBlocks() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem...
std::vector< Teuchos::RCP< const MV > > V
The current Krylov basis.
int resize(OrdinalType length_in)
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the &quot;native&quot; residual vectors.
Pure virtual base class which describes the basic interface to the linear solver iteration.
A pure virtual class for defining the status tests for the Belos iterative solvers.
Class which defines basic traits for the operator type.
int getNumIters() const
Get the current iteration count.
Structure to contain pointers to PseudoBlockGmresIter state variables.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
bool isInitialized()
States whether the solver has been initialized or not.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
Traits class which defines basic operations on multivectors.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
bool isParameter(const std::string &name) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void updateLSQR(int dim=-1)
Method for updating QR factorization of upper Hessenberg matrix.
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 &params)
PseudoBlockGmresIter constructor with linear problem, solver utilities, and parameter list of solver ...
PseudoBlockGmresIterOrthoFailure(const std::string &what_arg)
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, MagnitudeType > > > cs
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
void setNumBlocks(int numBlocks)
Set the maximum number of blocks used by the iterative solver.
Teuchos::ScalarTraits< ScalarType > SCT
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
void resetNumIters(int iter=0)
Reset the iteration count.
MultiVecTraits< ScalarType, MV > MVT
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, ScalarType > > > Z
The current right-hand side of the least squares system RY = Z.
OperatorTraits< ScalarType, MV, OP > OPT
void ROTG(ScalarType *da, ScalarType *db, rotg_c_type *c, ScalarType *s) const
void iterate()
This method performs block Gmres iterations until the status test indicates the need to stop or an er...
Class which defines basic traits for the operator type.
int curDim
The current dimension of the reduction.
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:28
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, ScalarType > > > sn
The current Given&#39;s rotation coefficients.
Belos&#39;s templated virtual class for providing routines for orthogonalization and orthonormzalition of...
Belos header file which uses auto-configuration information to include necessary C++ headers...
This class implements the pseudo-block GMRES iteration, where a block Krylov subspace is constructed ...
SerialDenseMatrix< OrdinalType, ScalarType > & assign(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
void setBlockSize(int blockSize)
Set the blocksize.
PseudoBlockGmresIterOrthoFailure is thrown when the orthogonalization manager is unable to generate o...
OrdinalType stride() const
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...

Generated on Thu Dec 5 2024 09:24:44 for Belos by doxygen 1.8.5