Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosGmresPolyOp.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_GMRESPOLYOP_HPP
11 #define BELOS_GMRESPOLYOP_HPP
12 
17 #include "BelosConfigDefs.hpp"
18 #include "BelosTypes.hpp"
19 
20 #include "BelosOperator.hpp"
21 #include "BelosMultiVec.hpp"
22 #include "BelosOperatorTraits.hpp"
23 #include "BelosMultiVecTraits.hpp"
24 #include "BelosLinearProblem.hpp"
25 
26 #include "BelosGmresIteration.hpp"
27 #include "BelosBlockGmresIter.hpp"
29 
33 #include "BelosStatusTestCombo.hpp"
35 
36 #include "BelosOutputManager.hpp"
37 
38 #include "Teuchos_BLAS.hpp"
39 #include "Teuchos_LAPACK.hpp"
40 #include "Teuchos_as.hpp"
41 #include "Teuchos_RCP.hpp"
46 
47 #ifdef BELOS_TEUCHOS_TIME_MONITOR
48  #include "Teuchos_TimeMonitor.hpp"
49 #endif // BELOS_TEUCHOS_TIME_MONITOR
50 
51 namespace Belos {
52 
59  class GmresPolyOpOrthoFailure : public BelosError {public:
60  GmresPolyOpOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
61  {}};
62 
63  // Create a shell class for the MV, inherited off MultiVec<> that will operate with the GmresPolyOp.
64  template <class ScalarType, class MV>
65  class GmresPolyMv : public MultiVec< ScalarType >
66  {
67  public:
68 
69  GmresPolyMv ( const Teuchos::RCP<MV>& mv_in )
70  : mv_(mv_in)
71  {}
73  {
74  mv_ = Teuchos::rcp_const_cast<MV>( mv_in );
75  }
76  Teuchos::RCP<MV> getMV() { return mv_; }
77  Teuchos::RCP<const MV> getConstMV() const { return mv_; }
78 
79  GmresPolyMv * Clone ( const int numvecs ) const
80  {
81  GmresPolyMv * newMV = new GmresPolyMv( MVT::Clone( *mv_, numvecs ) );
82  return newMV;
83  }
84  GmresPolyMv * CloneCopy () const
85  {
86  GmresPolyMv * newMV = new GmresPolyMv( MVT::CloneCopy( *mv_ ) );
87  return newMV;
88  }
89  GmresPolyMv * CloneCopy ( const std::vector<int>& index ) const
90  {
91  GmresPolyMv * newMV = new GmresPolyMv( MVT::CloneCopy( *mv_, index ) );
92  return newMV;
93  }
94  GmresPolyMv * CloneViewNonConst ( const std::vector<int>& index )
95  {
96  GmresPolyMv * newMV = new GmresPolyMv( MVT::CloneViewNonConst( *mv_, index ) );
97  return newMV;
98  }
99  const GmresPolyMv * CloneView ( const std::vector<int>& index ) const
100  {
101  const GmresPolyMv * newMV = new GmresPolyMv( MVT::CloneView( *mv_, index ) );
102  return newMV;
103  }
104  ptrdiff_t GetGlobalLength () const { return MVT::GetGlobalLength( *mv_ ); }
105  int GetNumberVecs () const { return MVT::GetNumberVecs( *mv_ ); }
106  void MvTimesMatAddMv (const ScalarType alpha,
107  const MultiVec<ScalarType>& A,
108  const Teuchos::SerialDenseMatrix<int,ScalarType>& B, const ScalarType beta)
109  {
110  const GmresPolyMv<ScalarType,MV>& A_in = dynamic_cast<const GmresPolyMv<ScalarType,MV>&>(A);
111  MVT::MvTimesMatAddMv( alpha, *(A_in.getConstMV()), B, beta, *mv_ );
112  }
113  void MvAddMv ( const ScalarType alpha, const MultiVec<ScalarType>& A, const ScalarType beta, const MultiVec<ScalarType>& B )
114  {
115  const GmresPolyMv<ScalarType,MV>& A_in = dynamic_cast<const GmresPolyMv<ScalarType,MV>&>(A);
116  const GmresPolyMv<ScalarType,MV>& B_in = dynamic_cast<const GmresPolyMv<ScalarType,MV>&>(B);
117  MVT::MvAddMv( alpha, *(A_in.getConstMV()), beta, *(B_in.getConstMV()), *mv_ );
118  }
119  void MvScale ( const ScalarType alpha ) { MVT::MvScale( *mv_, alpha ); }
120  void MvScale ( const std::vector<ScalarType>& alpha ) { MVT::MvScale( *mv_, alpha ); }
121  void MvTransMv ( const ScalarType alpha, const MultiVec<ScalarType>& A, Teuchos::SerialDenseMatrix<int,ScalarType>& B) const
122  {
123  const GmresPolyMv<ScalarType,MV>& A_in = dynamic_cast<const GmresPolyMv<ScalarType,MV>&>(A);
124  MVT::MvTransMv( alpha, *(A_in.getConstMV()), *mv_, B );
125  }
126  void MvDot ( const MultiVec<ScalarType>& A, std::vector<ScalarType>& b ) const
127  {
128  const GmresPolyMv<ScalarType,MV>& A_in = dynamic_cast<const GmresPolyMv<ScalarType,MV>&>(A);
129  MVT::MvDot( *(A_in.getConstMV()), *mv_, b );
130  }
131  void MvNorm ( std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& normvec, NormType type = TwoNorm ) const
132  {
133  MVT::MvNorm( *mv_, normvec, type );
134  }
135  void SetBlock ( const MultiVec<ScalarType>& A, const std::vector<int>& index )
136  {
137  const GmresPolyMv<ScalarType,MV>& A_in = dynamic_cast<const GmresPolyMv<ScalarType,MV>&>(A);
138  MVT::SetBlock( *(A_in.getConstMV()), index, *mv_ );
139  }
140  void MvRandom () { MVT::MvRandom( *mv_ ); }
141  void MvInit ( const ScalarType alpha ) { MVT::MvInit( *mv_, alpha ); }
142  void MvPrint ( std::ostream& os ) const { MVT::MvPrint( *mv_, os ); }
143 
144  private:
145 
146  typedef MultiVecTraits<ScalarType,MV> MVT;
147 
148  Teuchos::RCP<MV> mv_;
149 
150  };
151 
162  template <class ScalarType, class MV, class OP>
163  class GmresPolyOp : public Operator<ScalarType> {
164  public:
165 
167 
168 
171  const Teuchos::RCP<Teuchos::ParameterList>& params_in
172  )
173  : problem_(problem_in),
174  params_(params_in),
175  LP_(problem_in->getLeftPrec()),
176  RP_(problem_in->getRightPrec())
177  {
178  setParameters( params_ );
179 
180  polyUpdateLabel_ = label_ + ": Hybrid Gmres: Vector Update";
181 #ifdef BELOS_TEUCHOS_TIME_MONITOR
182  timerPolyUpdate_ = Teuchos::TimeMonitor::getNewCounter(polyUpdateLabel_);
183 #endif // BELOS_TEUCHOS_TIME_MONITOR
184 
185  if (polyType_ == "Arnoldi" || polyType_=="Roots")
187  else if (polyType_ == "Gmres")
189  else
190  TEUCHOS_TEST_FOR_EXCEPTION(polyType_!="Arnoldi"&&polyType_!="Gmres"&&polyType_!="Roots",std::invalid_argument,
191  "Belos::GmresPolyOp: \"Polynomial Type\" must be either \"Arnoldi\", \"Gmres\", or \"Roots\".");
192  }
193 
196  : problem_(problem_in)
197  {
198  // If dimension is zero, it will just apply the operator from problem_in in the Apply method.
199  dim_ = 0;
200  }
201 
203  virtual ~GmresPolyOp() {};
205 
207 
208 
210  void setParameters( const Teuchos::RCP<Teuchos::ParameterList>& params_in );
212 
214 
215 
219  void generateArnoldiPoly();
220 
224  void generateGmresPoly();
225 
227 
229 
230 
236  void ApplyPoly ( const MV& x, MV& y ) const;
237  void ApplyArnoldiPoly ( const MV& x, MV& y ) const;
238  void ApplyGmresPoly ( const MV& x, MV& y ) const;
239  void ApplyRootsPoly ( const MV& x, MV& y ) const;
240 
244  void Apply ( const MultiVec<ScalarType>& x, MultiVec<ScalarType>& y, ETrans /* trans */=NOTRANS ) const
245  {
246  const GmresPolyMv<ScalarType,MV>& x_in = dynamic_cast<const GmresPolyMv<ScalarType,MV>&>(x);
247  GmresPolyMv<ScalarType,MV>& y_in = dynamic_cast<GmresPolyMv<ScalarType,MV>&>(y);
248  ApplyPoly( *(x_in.getConstMV()), *(y_in.getMV()) );
249  }
250 
251  int polyDegree() const { return dim_; }
252 
253  private:
254 
255 #ifdef BELOS_TEUCHOS_TIME_MONITOR
256  Teuchos::RCP<Teuchos::Time> timerPolyUpdate_;
257 #endif // BELOS_TEUCHOS_TIME_MONITOR
258  std::string polyUpdateLabel_;
259 
260  typedef int OT; //Ordinal type
261  typedef MultiVecTraits<ScalarType,MV> MVT;
263  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
265 
266  // Default polynomial parameters
267  static constexpr int maxDegree_default_ = 25;
268  static constexpr int verbosity_default_ = Belos::Errors;
269  static constexpr bool randomRHS_default_ = true;
270  static constexpr const char * label_default_ = "Belos";
271  static constexpr const char * polyType_default_ = "Roots";
272  static constexpr const char * orthoType_default_ = "DGKS";
273  static constexpr bool damp_default_ = false;
274  static constexpr bool addRoots_default_ = true;
275 
276  // Variables for generating the polynomial
279  Teuchos::RCP<const OP> LP_, RP_;
280 
281  // Output manager.
283  Teuchos::RCP<std::ostream> outputStream_ = Teuchos::rcpFromRef(std::cout);
284 
285  // Orthogonalization manager.
287 
288  // Current polynomial parameters
289  MagnitudeType polyTol_ = DefaultSolverParameters::polyTol;
290  int maxDegree_ = maxDegree_default_;
291  int verbosity_ = verbosity_default_;
292  bool randomRHS_ = randomRHS_default_;
293  std::string label_ = label_default_;
294  std::string polyType_ = polyType_default_;
295  std::string orthoType_ = orthoType_default_;
296  int dim_ = 0;
297  bool damp_ = damp_default_;
298  bool addRoots_ = addRoots_default_;
299 
300  // Variables for Arnoldi polynomial
301  mutable Teuchos::RCP<MV> V_, wL_, wR_;
304 
305  // Variables for Gmres polynomial;
306  bool autoDeg = false;
308 
309  // Variables for Roots polynomial:
311 
312  // Modified Leja sorting function. Takes a serial dense matrix of M harmonic Ritz values and an index
313  // of values from 0 to M. Returns the sorted values and sorted index, similar to Matlab.
314  void SortModLeja(Teuchos::SerialDenseMatrix< OT, MagnitudeType > &thetaN, std::vector<int> &index) const ;
315 
316  //Function determines whether added roots are needed and adds them if option is turned on.
317  void ComputeAddedRoots();
318  };
319 
320  template <class ScalarType, class MV, class OP>
322  {
323  // Check which Gmres polynomial to use
324  if (params_in->isParameter("Polynomial Type")) {
325  polyType_ = params_in->get("Polynomial Type", polyType_default_);
326  }
327 
328  // Check for polynomial convergence tolerance
329  if (params_in->isParameter("Polynomial Tolerance")) {
330  if (params_in->isType<MagnitudeType> ("Polynomial Tolerance")) {
331  polyTol_ = params_in->get ("Polynomial Tolerance",
332  static_cast<MagnitudeType> (DefaultSolverParameters::polyTol));
333  }
334  else {
335  polyTol_ = params_in->get ("Polynomial Tolerance", DefaultSolverParameters::polyTol);
336  }
337  }
338 
339  // Check for maximum polynomial degree
340  if (params_in->isParameter("Maximum Degree")) {
341  maxDegree_ = params_in->get("Maximum Degree", maxDegree_default_);
342  }
343 
344  // Check for maximum polynomial degree
345  if (params_in->isParameter("Random RHS")) {
346  randomRHS_ = params_in->get("Random RHS", randomRHS_default_);
347  }
348 
349  // Check for a change in verbosity level
350  if (params_in->isParameter("Verbosity")) {
351  if (Teuchos::isParameterType<int>(*params_in,"Verbosity")) {
352  verbosity_ = params_in->get("Verbosity", verbosity_default_);
353  }
354  else {
355  verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params_in,"Verbosity");
356  }
357  }
358 
359  if (params_in->isParameter("Orthogonalization")) {
360  orthoType_ = params_in->get("Orthogonalization",orthoType_default_);
361  }
362 
363  // Check for timer label
364  if (params_in->isParameter("Timer Label")) {
365  label_ = params_in->get("Timer Label", label_default_);
366  }
367 
368  // Output stream
369  if (params_in->isParameter("Output Stream")) {
370  outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params_in,"Output Stream");
371  }
372 
373  // Check for damped polynomial
374  if (params_in->isParameter("Damped Poly")) {
375  damp_ = params_in->get("Damped Poly", damp_default_);
376  }
377 
378  // Check for root-adding
379  if (params_in->isParameter("Add Roots")) {
380  addRoots_ = params_in->get("Add Roots", addRoots_default_);
381  }
382  }
383 
384  template <class ScalarType, class MV, class OP>
386  {
387  Teuchos::RCP< MV > V = MVT::Clone( *problem_->getRHS(), maxDegree_+1 );
388 
389  //Make power basis:
390  std::vector<int> index(1,0);
391  Teuchos::RCP< MV > V0 = MVT::CloneViewNonConst(*V, index);
392  if (randomRHS_)
393  MVT::MvRandom( *V0 );
394  else
395  MVT::Assign( *problem_->getRHS(), *V0 );
396 
397  if ( !LP_.is_null() ) {
398  Teuchos::RCP< MV > Vtemp = MVT::CloneCopy(*V0);
399  problem_->applyLeftPrec( *Vtemp, *V0);
400  }
401  if ( damp_ ) {
402  Teuchos::RCP< MV > Vtemp = MVT::CloneCopy(*V0);
403  problem_->apply( *Vtemp, *V0);
404  }
405 
406  for(int i=0; i< maxDegree_; i++)
407  {
408  index[0] = i;
409  Teuchos::RCP< const MV > Vi = MVT::CloneView(*V, index);
410  index[0] = i+1;
411  Teuchos::RCP< MV > Vip1 = MVT::CloneViewNonConst(*V, index);
412  problem_->apply( *Vi, *Vip1);
413  }
414 
415  //Consider AV:
416  Teuchos::Range1D range( 1, maxDegree_);
417  Teuchos::RCP< const MV > AV = MVT::CloneView( *V, range);
418 
419  //Make lhs (AV)^T(AV)
420  Teuchos::SerialDenseMatrix< OT, ScalarType > AVtransAV( maxDegree_, maxDegree_);
421  MVT::MvTransMv( SCT::one(), *AV, *AV, AVtransAV);
422  //This process adds pDeg*pDeg + pDeg inner products that aren't in the final count.
423 
425  int infoInt;
426  bool status = true; //Keep adjusting poly deg when true.
427 
428  dim_ = maxDegree_;
430  while( status && dim_ >= 1)
431  {
432  Teuchos::SerialDenseMatrix< OT, ScalarType > lhstemp(Teuchos::Copy, AVtransAV, dim_, dim_);
433  lapack.POTRF( 'U', dim_, lhstemp.values(), lhstemp.stride(), &infoInt);
434 
435  if(autoDeg == false)
436  {
437  status = false;
438  if(infoInt != 0)
439  {
440  std::cout << "BelosGmresPolyOp.hpp: LAPACK POTRF was not successful!!" << std::endl;
441  std::cout << "Error code: " << infoInt << std::endl;
442  }
443  }
444  else
445  {
446  if(infoInt != 0)
447  {//Had bad factor. Reduce poly degree.
448  dim_--;
449  }
450  else
451  {
452  status = false;
453  }
454  }
455  if(status == false)
456  {
457  lhs = lhstemp;
458  }
459  }
460  if(dim_ == 0)
461  {
462  pCoeff_.shape( 1, 1);
463  pCoeff_(0,0) = SCT::one();
464  std::cout << "Poly Degree is zero. No preconditioner created." << std::endl;
465  }
466  else
467  {
468  pCoeff_.shape( dim_, 1);
469  //Get correct submatrix of AV:
470  Teuchos::Range1D rangeSub( 1, dim_);
471  Teuchos::RCP< const MV > AVsub = MVT::CloneView( *V, rangeSub);
472 
473  //Compute rhs (AV)^T V0
474  MVT::MvTransMv( SCT::one(), *AVsub, *V0, pCoeff_);
475  lapack.POTRS( 'U', dim_, 1, lhs.values(), lhs.stride(), pCoeff_.values(), pCoeff_.stride(), &infoInt);
476  if(infoInt != 0)
477  {
478  std::cout << "BelosGmresPolyOp.hpp: LAPACK POTRS was not successful!!" << std::endl;
479  std::cout << "Error code: " << infoInt << std::endl;
480  }
481  }
482  }
483 
484  template <class ScalarType, class MV, class OP>
486  {
487  std::string polyLabel = label_ + ": GmresPolyOp creation";
488 
489  // Create a copy of the linear problem that has a zero initial guess and random RHS.
490  std::vector<int> idx(1,0);
491  Teuchos::RCP<MV> newX = MVT::Clone( *(problem_->getLHS()), 1 );
492  Teuchos::RCP<MV> newB = MVT::Clone( *(problem_->getRHS()), 1 );
493  MVT::MvInit( *newX, SCT::zero() );
494  if (randomRHS_) {
495  MVT::MvRandom( *newB );
496  }
497  else {
498  MVT::Assign( *(MVT::CloneView(*(problem_->getRHS()), idx)), *newB );
499  }
501  Teuchos::rcp( new LinearProblem<ScalarType,MV,OP>( problem_->getOperator(), newX, newB ) );
502  newProblem->setInitResVec( newB );
503  newProblem->setLeftPrec( problem_->getLeftPrec() );
504  newProblem->setRightPrec( problem_->getRightPrec() );
505  newProblem->setLabel(polyLabel);
506  newProblem->setProblem();
507  newProblem->setLSIndex( idx );
508 
509  // Create a parameter list for the GMRES iteration.
510  Teuchos::ParameterList polyList;
511 
512  // Tell the block solver that the block size is one.
513  polyList.set("Num Blocks",maxDegree_);
514  polyList.set("Block Size",1);
515  polyList.set("Keep Hessenberg", true);
516 
517  // Create output manager.
518  printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
519 
520  // Create orthogonalization manager if we need to.
521  if (ortho_.is_null()) {
522  params_->set("Orthogonalization", orthoType_);
524  Teuchos::RCP<Teuchos::ParameterList> paramsOrtho; // can be null
525 
526  ortho_ = factory.makeMatOrthoManager (orthoType_, Teuchos::null, printer_, polyLabel, paramsOrtho);
527  }
528 
529  // Create a simple status test that either reaches the relative residual tolerance or maximum polynomial size.
532 
533  // Implicit residual test, using the native residual to determine if convergence was achieved.
536  convTst->defineScaleForm( convertStringToScaleType("Norm of RHS"), Belos::TwoNorm );
537 
538  // Convergence test that stops the iteration when either are satisfied.
541 
542  // Create Gmres iteration object to perform one cycle of Gmres.
544  gmres_iter = Teuchos::rcp( new BlockGmresIter<ScalarType,MV,OP>(newProblem,printer_,polyTest,ortho_,polyList) );
545 
546  // Create the first block in the current Krylov basis (residual).
547  Teuchos::RCP<MV> V_0 = MVT::CloneCopy( *newB );
548  if ( !LP_.is_null() )
549  newProblem->applyLeftPrec( *newB, *V_0 );
550  if ( damp_ )
551  {
552  Teuchos::RCP< MV > Vtemp = MVT::CloneCopy(*V_0);
553  newProblem->apply( *Vtemp, *V_0 );
554  }
555 
556  // Get a matrix to hold the orthonormalization coefficients.
557  r0_.resize(1);
558 
559  // Orthonormalize the new V_0
560  int rank = ortho_->normalize( *V_0, Teuchos::rcpFromRef(r0_) );
562  "Belos::GmresPolyOp::generateArnoldiPoly(): Failed to compute initial block of orthonormal vectors for polynomial generation.");
563 
564  // Set the new state and initialize the solver.
566  newstate.V = V_0;
567  newstate.z = Teuchos::rcpFromRef( r0_);
568  newstate.curDim = 0;
569  gmres_iter->initializeGmres(newstate);
570 
571  // Perform Gmres iteration
572  try {
573  gmres_iter->iterate();
574  }
575  catch (GmresIterationOrthoFailure& e) {
576  // Try to recover the most recent least-squares solution
577  gmres_iter->updateLSQR( gmres_iter->getCurSubspaceDim() );
578  }
579  catch (std::exception& e) {
580  using std::endl;
581  printer_->stream(Errors) << "Error! Caught exception in BlockGmresIter::iterate() at iteration "
582  << gmres_iter->getNumIters() << endl << e.what () << endl;
583  throw;
584  }
585 
586  // Get the solution for this polynomial, use in comparison below
587  Teuchos::RCP<MV> currX = gmres_iter->getCurrentUpdate();
588 
589  // Record polynomial info, get current GMRES state
590  GmresIterationState<ScalarType,MV> gmresState = gmres_iter->getState();
591 
592  // If the polynomial has no dimension, the tolerance is too low, return false
593  dim_ = gmresState.curDim;
594  if (dim_ == 0) {
595  return;
596  }
597  if(polyType_ == "Arnoldi"){
598  // Make a view and then copy the RHS of the least squares problem.
599  //
600  y_ = Teuchos::SerialDenseMatrix<OT,ScalarType>( Teuchos::Copy, *gmresState.z, dim_, 1 );
601  H_ = *gmresState.H;
602 
603  //
604  // Solve the least squares problem.
605  //
608  Teuchos::NON_UNIT_DIAG, dim_, 1, SCT::one(),
609  gmresState.R->values(), gmresState.R->stride(),
610  y_.values(), y_.stride() );
611  }
612  else{ //Generate Roots Poly
613  //Find Harmonic Ritz Values to use as polynomial roots:
614 
615  //Copy of square H used to find poly roots:
616  H_ = Teuchos::SerialDenseMatrix<OT,ScalarType>(Teuchos::Copy, *gmresState.H, dim_, dim_);
617  //Zero out below subdiagonal of H:
618  for(int i=0; i <= dim_-3; i++) {
619  for(int k=i+2; k <= dim_-1; k++) {
620  H_(k,i) = SCT::zero();
621  }
622  }
623  //Extra copy of H because equilibrate changes the matrix:
624  Teuchos::SerialDenseMatrix<OT,ScalarType> Htemp (Teuchos::Copy, H_, dim_, dim_);
625 
626  //View the m+1,m element and last col of H:
627  ScalarType Hlast = (*gmresState.H)(dim_,dim_-1);
628  Teuchos::SerialDenseMatrix<OT,ScalarType> HlastCol (Teuchos::View, H_, dim_, 1, 0, dim_-1);
629 
630  //Set up linear system for H^{-*}e_m:
632  E.putScalar(SCT::zero());
633  E(dim_-1,0) = SCT::one();
634 
636  HSolver.setMatrix( Teuchos::rcpFromRef(Htemp));
638  HSolver.setVectors( Teuchos::rcpFromRef(F), Teuchos::rcpFromRef(E));
639  HSolver.factorWithEquilibration( true );
640 
641  //Factor matrix and solve for F = H^{-*}e_m:
642  int info = 0;
643  info = HSolver.factor();
644  if(info != 0){
645  std::cout << "Hsolver factor: info = " << info << std::endl;
646  }
647  info = HSolver.solve();
648  if(info != 0){
649  std::cout << "Hsolver solve : info = " << info << std::endl;
650  }
651 
652  //Scale F and adjust H for Harmonic Ritz value eigenproblem:
653  F.scale(Hlast*Hlast);
654  HlastCol += F;
655 
656  //Set up for eigenvalue problem to get Harmonic Ritz Values:
658  theta_.shape(dim_,2);//1st col for real part, 2nd col for imaginary
659 
660  const int ldv = 1;
661  ScalarType* vlr = 0;
662 
663  // Size of workspace and workspace for DGEEV
664  int lwork = -1;
665  std::vector<ScalarType> work(1);
666  std::vector<MagnitudeType> rwork(2*dim_);
667 
668  //Find workspace size for DGEEV:
669  lapack.GEEV('N','N',dim_,H_.values(),H_.stride(),theta_[0],theta_[1],vlr, ldv, vlr, ldv, &work[0], lwork, &rwork[0], &info);
670  lwork = std::abs (static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work[0])));
671  work.resize( lwork );
672  // Solve for Harmonic Ritz Values:
673  lapack.GEEV('N','N',dim_,H_.values(),H_.stride(),theta_[0],theta_[1],vlr, ldv, vlr, ldv, &work[0], lwork, &rwork[0], &info);
674 
675  if(info != 0){
676  std::cout << "GEEV solve : info = " << info << std::endl;
677  }
678 
679  // Set index for sort function, verify roots are non-zero,
680  // and sort Harmonic Ritz Values:
681  const MagnitudeType tol = 10.0 * Teuchos::ScalarTraits<MagnitudeType>::eps();
682  std::vector<int> index(dim_);
683  for(int i=0; i<dim_; ++i){
684  index[i] = i;
685  // Check if real + imag parts of roots < tol.
686  TEUCHOS_TEST_FOR_EXCEPTION(hypot(theta_(i,0),theta_(i,1)) < tol, std::runtime_error, "BelosGmresPolyOp Error: One of the computed polynomial roots is approximately zero. This will cause a divide by zero error! Your matrix may be close to singular. Please select a lower polynomial degree or give a shifted matrix.");
687  }
688  SortModLeja(theta_,index);
689 
690  //Add roots if neded.
691  ComputeAddedRoots();
692 
693  }
694  }
695 
696  //Function determines whether added roots are needed and adds them if option is turned on.
697  template <class ScalarType, class MV, class OP>
699  {
700  // Store theta (with cols for real and imag parts of Harmonic Ritz Vals)
701  // as one vector of complex numbers to perform arithmetic:
702  std::vector<std::complex<MagnitudeType>> cmplxHRitz (dim_);
703  for(unsigned int i=0; i<cmplxHRitz.size(); ++i){
704  cmplxHRitz[i] = std::complex<MagnitudeType>( theta_(i,0), theta_(i,1) );
705  }
706 
707  // Compute product of factors (pof) to determine added roots:
708  const MagnitudeType one(1.0);
709  std::vector<MagnitudeType> pof (dim_,one);
710  for(int j=0; j<dim_; ++j) {
711  for(int i=0; i<dim_; ++i) {
712  if(i!=j) {
713  pof[j] = std::abs(pof[j]*(one-(cmplxHRitz[j]/cmplxHRitz[i])));
714  }
715  }
716  }
717 
718  // Compute number of extra roots needed:
719  std::vector<int> extra (dim_);
720  int totalExtra = 0;
721  for(int i=0; i<dim_; ++i){
722  if (pof[i] > MCT::zero())
723  extra[i] = ceil((log10(pof[i])-MagnitudeType(4.0))/MagnitudeType(14.0));
724  else
725  extra[i] = 0;
726  if(extra[i] > 0){
727  totalExtra += extra[i];
728  }
729  }
730  if (totalExtra){
731  printer_->stream(Warnings) << "Warning: Need to add " << totalExtra << " extra roots." << std::endl;}
732 
733  // If requested to add roots, append them to the theta matrix:
734  if(addRoots_ && totalExtra>0)
735  {
736  theta_.reshape(dim_+totalExtra,2);
737  // Make a matrix copy for perturbed roots:
738  Teuchos::SerialDenseMatrix<OT,MagnitudeType> thetaPert (Teuchos::Copy, theta_, dim_+totalExtra, 2);
739 
740  //Add extra eigenvalues to matrix and perturb for sort:
741  int count = dim_;
742  for(int i=0; i<dim_; ++i){
743  for(int j=0; j< extra[i]; ++j){
744  theta_(count,0) = theta_(i,0);
745  theta_(count,1) = theta_(i,1);
746  thetaPert(count,0) = theta_(i,0)+(j+MCT::one())*MagnitudeType(5e-8);
747  thetaPert(count,1) = theta_(i,1);
748  ++count;
749  }
750  }
751 
752  // Update polynomial degree:
753  dim_ += totalExtra;
754  if (totalExtra){
755  printer_->stream(Warnings) << "New poly degree is: " << dim_ << std::endl;}
756 
757  // Create a new index and sort perturbed roots:
758  std::vector<int> index2(dim_);
759  for(int i=0; i<dim_; ++i){
760  index2[i] = i;
761  }
762  SortModLeja(thetaPert,index2);
763  //Apply sorting to non-perturbed roots:
764  for(int i=0; i<dim_; ++i)
765  {
766  thetaPert(i,0) = theta_(index2[i],0);
767  thetaPert(i,1) = theta_(index2[i],1);
768  }
769  theta_ = thetaPert;
770 
771  }
772  }
773 
774  // Modified Leja sorting function. Takes a serial dense matrix of M harmonic Ritz values and an index
775  // of values from 0 to M. Returns the sorted values and sorted index, similar to Matlab.
776  template <class ScalarType, class MV, class OP>
777  void GmresPolyOp<ScalarType, MV, OP>::SortModLeja(Teuchos::SerialDenseMatrix< OT, MagnitudeType > &thetaN, std::vector<int> &index) const
778  {
779  //Sort theta values via Modified Leja Ordering:
780 
781  // Set up blank matrices to track sorting:
782  int dimN = index.size();
783  std::vector<int> newIndex(dimN);
787 
788  //Compute all absolute values and find maximum:
789  for(int i = 0; i < dimN; i++){
790  absVal(i) = hypot(thetaN(i,0), thetaN(i,1));
791  }
792  MagnitudeType * maxPointer = std::max_element(absVal.values(), (absVal.values()+dimN));
793  int maxIndex = int (maxPointer- absVal.values());
794 
795  //Put largest abs value first in the list:
796  sorted(0,0) = thetaN(maxIndex,0);
797  sorted(0,1) = thetaN(maxIndex,1);
798  newIndex[0] = index[maxIndex];
799 
800  int j;
801  // If largest value was complex (for real scalar type) put its conjugate in the next slot.
802  if(sorted(0,1)!= SCT::zero() && !SCT::isComplex)
803  {
804  sorted(1,0) = thetaN(maxIndex,0);
805  sorted(1,1) = -thetaN(maxIndex,1);
806  newIndex[1] = index[maxIndex+1];
807  j = 2;
808  }
809  else
810  {
811  j = 1;
812  }
813 
814  //Sort remaining values:
815  MagnitudeType a, b;
816  while( j < dimN )
817  {
818  //For each value, compute (a log of) a product of differences:
819  for(int i = 0; i < dimN; i++)
820  {
821  prod(i) = MCT::one();
822  for(int k = 0; k < j; k++)
823  {
824  a = thetaN(i,0) - sorted(k,0);
825  b = thetaN(i,1) - sorted(k,1);
826  if (a*a + b*b > MCT::zero())
827  prod(i) = prod(i) + log10(hypot(a,b));
828  else {
829  prod(i) = -std::numeric_limits<MagnitudeType>::infinity();
830  break;
831  }
832  }
833  }
834 
835  //Value with largest product goes in the next slot:
836  maxPointer = std::max_element(prod.values(), (prod.values()+dimN));
837  maxIndex = int (maxPointer- prod.values());
838  sorted(j,0) = thetaN(maxIndex,0);
839  sorted(j,1) = thetaN(maxIndex,1);
840  newIndex[j] = index[maxIndex];
841 
842  //If it was complex (and scalar type real) put its conjugate in next slot:
843  if(sorted(j,1)!= SCT::zero() && !SCT::isComplex)
844  {
845  j++;
846  sorted(j,0) = thetaN(maxIndex,0);
847  sorted(j,1) = -thetaN(maxIndex,1);
848  newIndex[j] = index[maxIndex+1];
849  }
850  j++;
851  }
852 
853  //Return sorted values and sorted indices:
854  thetaN = sorted;
855  index = newIndex;
856  } //End Modified Leja ordering
857 
858  template <class ScalarType, class MV, class OP>
859  void GmresPolyOp<ScalarType, MV, OP>::ApplyPoly( const MV& x, MV& y ) const
860  {
861  if (dim_) {
862  if (polyType_ == "Arnoldi")
863  ApplyArnoldiPoly(x, y);
864  else if (polyType_ == "Gmres")
865  ApplyGmresPoly(x, y);
866  else if (polyType_ == "Roots")
867  ApplyRootsPoly(x, y);
868  }
869  else {
870  // Just apply the operator in problem_ to x and return y.
871  problem_->applyOp( x, y );
872  }
873  }
874 
875  template <class ScalarType, class MV, class OP>
876  void GmresPolyOp<ScalarType, MV, OP>::ApplyGmresPoly( const MV& x, MV& y ) const
877  {
878  Teuchos::RCP<MV> AX = MVT::CloneCopy(x);
879  Teuchos::RCP<MV> AX2 = MVT::Clone( x, MVT::GetNumberVecs(x) );
880 
881  // Apply left preconditioner.
882  if (!LP_.is_null()) {
883  Teuchos::RCP<MV> Xtmp = MVT::Clone( x, MVT::GetNumberVecs(x) );
884  problem_->applyLeftPrec( *AX, *Xtmp ); // Left precondition x into the first vector
885  AX = Xtmp;
886  }
887 
888  {
889 #ifdef BELOS_TEUCHOS_TIME_MONITOR
890  Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
891 #endif
892  MVT::MvAddMv(pCoeff_(0,0), *AX, SCT::zero(), y, y); //y= coeff_i(A^ix)
893  }
894  for( int i=1; i < dim_; i++)
895  {
896  Teuchos::RCP<MV> X, Y;
897  if ( i%2 )
898  {
899  X = AX;
900  Y = AX2;
901  }
902  else
903  {
904  X = AX2;
905  Y = AX;
906  }
907  problem_->apply(*X, *Y);
908  {
909 #ifdef BELOS_TEUCHOS_TIME_MONITOR
910  Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
911 #endif
912  MVT::MvAddMv(pCoeff_(i,0), *Y, SCT::one(), y, y); //y= coeff_i(A^ix) +y
913  }
914  }
915 
916  // Apply right preconditioner.
917  if (!RP_.is_null()) {
918  Teuchos::RCP<MV> Ytmp = MVT::CloneCopy(y);
919  problem_->applyRightPrec( *Ytmp, y );
920  }
921  }
922 
923  template <class ScalarType, class MV, class OP>
924  void GmresPolyOp<ScalarType, MV, OP>::ApplyRootsPoly( const MV& x, MV& y ) const
925  {
926  MVT::MvInit( y, SCT::zero() ); //Zero out y to take the vector with poly applied.
927  Teuchos::RCP<MV> prod = MVT::CloneCopy(x);
928  Teuchos::RCP<MV> Xtmp = MVT::Clone( x, MVT::GetNumberVecs(x) );
929  Teuchos::RCP<MV> Xtmp2 = MVT::Clone( x, MVT::GetNumberVecs(x) );
930 
931  // Apply left preconditioner.
932  if (!LP_.is_null()) {
933  problem_->applyLeftPrec( *prod, *Xtmp ); // Left precondition x into the first vector
934  prod = Xtmp;
935  }
936 
937  int i=0;
938  while(i < dim_-1)
939  {
940  if(theta_(i,1)== SCT::zero() || SCT::isComplex) //Real Harmonic Ritz value or complex scalars
941  {
942  {
943 #ifdef BELOS_TEUCHOS_TIME_MONITOR
944  Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
945 #endif
946  MVT::MvAddMv(SCT::one(), y, SCT::one()/theta_(i,0), *prod, y); //poly = poly + 1/theta_i * prod
947  }
948  problem_->apply(*prod, *Xtmp); // temp = A*prod
949  {
950 #ifdef BELOS_TEUCHOS_TIME_MONITOR
951  Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
952 #endif
953  MVT::MvAddMv(SCT::one(), *prod, -SCT::one()/theta_(i,0), *Xtmp, *prod); //prod = prod - 1/theta_i * temp
954  }
955  i++;
956  }
957  else //Current theta is complex and has a conjugate; combine to preserve real arithmetic
958  {
959  MagnitudeType mod = theta_(i,0)*theta_(i,0) + theta_(i,1)*theta_(i,1); //mod = a^2 + b^2
960  problem_->apply(*prod, *Xtmp); // temp = A*prod
961  {
962 #ifdef BELOS_TEUCHOS_TIME_MONITOR
963  Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
964 #endif
965  MVT::MvAddMv(2*theta_(i,0), *prod, -SCT::one(), *Xtmp, *Xtmp); //temp = 2a*prod-temp
966  MVT::MvAddMv(SCT::one(), y, SCT::one()/mod, *Xtmp, y); //poly = poly + 1/mod*temp
967  }
968  if( i < dim_-2 )
969  {
970  problem_->apply(*Xtmp, *Xtmp2); // temp2 = A*temp
971  {
972 #ifdef BELOS_TEUCHOS_TIME_MONITOR
973  Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
974 #endif
975  MVT::MvAddMv(SCT::one(), *prod, -SCT::one()/mod, *Xtmp2, *prod); //prod = prod - 1/mod * temp2
976  }
977  }
978  i = i + 2;
979  }
980  }
981  if(theta_(dim_-1,1)== SCT::zero() || SCT::isComplex)
982  {
983 #ifdef BELOS_TEUCHOS_TIME_MONITOR
984  Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
985 #endif
986  MVT::MvAddMv(SCT::one(), y, SCT::one()/theta_(dim_-1,0), *prod, y); //poly = poly + 1/theta_i * prod
987  }
988 
989  // Apply right preconditioner.
990  if (!RP_.is_null()) {
991  Teuchos::RCP<MV> Ytmp = MVT::CloneCopy(y);
992  problem_->applyRightPrec( *Ytmp, y );
993  }
994  }
995 
996  template <class ScalarType, class MV, class OP>
997  void GmresPolyOp<ScalarType, MV, OP>::ApplyArnoldiPoly( const MV& x, MV& y ) const
998  {
999  // Initialize vector storage.
1000  if (V_.is_null()) {
1001  V_ = MVT::Clone( x, dim_ );
1002  if (!LP_.is_null()) {
1003  wL_ = MVT::Clone( y, 1 );
1004  }
1005  if (!RP_.is_null()) {
1006  wR_ = MVT::Clone( y, 1 );
1007  }
1008  }
1009  //
1010  // Apply polynomial to x.
1011  //
1012  int n = MVT::GetNumberVecs( x );
1013  std::vector<int> idxi(1), idxi2, idxj(1);
1014 
1015  // Select vector x[j].
1016  for (int j=0; j<n; ++j) {
1017 
1018  idxi[0] = 0;
1019  idxj[0] = j;
1020  Teuchos::RCP<const MV> x_view = MVT::CloneView( x, idxj );
1021  Teuchos::RCP<MV> y_view = MVT::CloneViewNonConst( y, idxj );
1022  if (!LP_.is_null()) {
1023  Teuchos::RCP<MV> v_curr = MVT::CloneViewNonConst( *V_, idxi );
1024  problem_->applyLeftPrec( *x_view, *v_curr ); // Left precondition x into the first vector of V
1025  } else {
1026  MVT::SetBlock( *x_view, idxi, *V_ ); // Set x as the first vector of V
1027  }
1028 
1029  for (int i=0; i<dim_-1; ++i) {
1030 
1031  // Get views into the current and next vectors
1032  idxi2.resize(i+1);
1033  for (int ii=0; ii<i+1; ++ii) { idxi2[ii] = ii; }
1034  Teuchos::RCP<const MV> v_prev = MVT::CloneView( *V_, idxi2 );
1035  // the tricks below with wR_ and wL_ (potentially set to v_curr and v_next) unfortunately imply that
1036  // v_curr and v_next must be non-const views.
1037  Teuchos::RCP<MV> v_curr = MVT::CloneViewNonConst( *V_, idxi );
1038  idxi[0] = i+1;
1039  Teuchos::RCP<MV> v_next = MVT::CloneViewNonConst( *V_, idxi );
1040 
1041  //---------------------------------------------
1042  // Apply operator to next vector
1043  //---------------------------------------------
1044  // 1) Apply right preconditioner, if we have one.
1045  if (!RP_.is_null()) {
1046  problem_->applyRightPrec( *v_curr, *wR_ );
1047  } else {
1048  wR_ = v_curr;
1049  }
1050  // 2) Check for left preconditioner, if none exists, point at the next vector.
1051  if (LP_.is_null()) {
1052  wL_ = v_next;
1053  }
1054  // 3) Apply operator A.
1055  problem_->applyOp( *wR_, *wL_ );
1056  // 4) Apply left preconditioner, if we have one.
1057  if (!LP_.is_null()) {
1058  problem_->applyLeftPrec( *wL_, *v_next );
1059  }
1060 
1061  // Compute A*v_curr - v_prev*H(1:i,i)
1063  {
1064 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1065  Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
1066 #endif
1067  MVT::MvTimesMatAddMv( -SCT::one(), *v_prev, h, SCT::one(), *v_next );
1068  }
1069 
1070  // Scale by H(i+1,i)
1071  MVT::MvScale( *v_next, SCT::one()/H_(i+1,i) );
1072  }
1073 
1074  // Compute output y = V*y_./r0_
1075  if (!RP_.is_null()) {
1076  {
1077 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1078  Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
1079 #endif
1080  MVT::MvTimesMatAddMv( SCT::one()/r0_(0), *V_, y_, SCT::zero(), *wR_ );
1081  }
1082  problem_->applyRightPrec( *wR_, *y_view );
1083  }
1084  else {
1085 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1086  Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
1087 #endif
1088  MVT::MvTimesMatAddMv( SCT::one()/r0_(0), *V_, y_, SCT::zero(), *y_view );
1089  }
1090  } // (int j=0; j<n; ++j)
1091  } // end Apply()
1092 } // end Belos namespace
1093 
1094 #endif
1095 
1096 // end of file BelosGmresPolyOp.hpp
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
Definition: BelosTypes.cpp:74
ScalarType * values() const
GmresPolyMv(const Teuchos::RCP< const MV > &mv_in)
Collection of types and exceptions used within the Belos solvers.
const GmresPolyMv * CloneView(const std::vector< int > &index) const
Creates a new Belos::MultiVec that shares the selected contents of *this. The index of the numvecs ve...
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...
void MvScale(const std::vector< ScalarType > &alpha)
Scale each element of the i-th vector in *this with alpha[i].
Class which manages the output and verbosity of the Belos solvers.
static void MvDot(const MV &mv, const MV &A, std::vector< ScalarType > &b)
Compute a vector b where the components are the individual dot-products of the i-th columns of A and ...
void MvAddMv(const ScalarType alpha, const MultiVec< ScalarType > &A, const ScalarType beta, const MultiVec< ScalarType > &B)
Replace *this with alpha * A + beta * B.
GmresPolyMv * CloneCopy() const
Create a new MultiVec and copy contents of *this into it (deep copy).
void MvScale(const ScalarType alpha)
Scale each element of the vectors in *this with alpha.
static magnitudeType eps()
ptrdiff_t GetGlobalLength() const
The number of rows in the multivector.
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
This class implements the block GMRES iteration, where a block Krylov subspace is constructed...
Teuchos::RCP< const MV > V
The current Krylov basis.
Teuchos::RCP< MV > getMV()
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
static RCP< Time > getNewCounter(const std::string &name)
void ApplyGmresPoly(const MV &x, MV &y) const
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
void generateGmresPoly()
This routine takes the matrix, preconditioner, and vectors from the linear problem as well as the par...
Declaration of basic traits for the multivector type.
void POTRS(const char &UPLO, const OrdinalType &n, const OrdinalType &nrhs, const ScalarType *A, const OrdinalType &lda, ScalarType *B, const OrdinalType &ldb, OrdinalType *info) const
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > z
The current right-hand side of the least squares system RY = Z.
An implementation of StatusTestResNorm using a family of residual norms.
GmresPolyMv * CloneViewNonConst(const std::vector< int > &index)
Creates a new Belos::MultiVec that shares the selected contents of *this. The index of the numvecs ve...
GmresPolyMv * CloneCopy(const std::vector< int > &index) const
Creates a new Belos::MultiVec and copies the selected contents of *this into the new multivector (dee...
void MvRandom()
Fill all the vectors in *this with random numbers.
Structure to contain pointers to GmresIteration state variables.
GmresPolyOp(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem_in, const Teuchos::RCP< Teuchos::ParameterList > &params_in)
Basic contstructor.
Belos::StatusTest class for specifying a maximum number of iterations.
static void MvInit(MV &mv, const ScalarType alpha=Teuchos::ScalarTraits< ScalarType >::zero())
Replace each element of the vectors in mv with alpha.
static const double polyTol
Relative residual tolerance for matrix polynomial construction.
Definition: BelosTypes.hpp:264
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &mv, Teuchos::SerialDenseMatrix< int, ScalarType > &B)
Compute a dense matrix B through the matrix-matrix multiply .
Class which defines basic traits for the operator type.
int GetNumberVecs() const
The number of vectors (i.e., columns) in the multivector.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
A factory class for generating StatusTestOutput objects.
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
ETrans
Whether to apply the (conjugate) transpose of an operator.
Definition: BelosTypes.hpp:49
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
bool isParameter(const std::string &name) const
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index...
A Belos::StatusTest class for specifying a maximum number of iterations.
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
void MvTimesMatAddMv(const ScalarType alpha, const MultiVec< ScalarType > &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta)
Update *this with alpha * A * B + beta * (*this).
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
Alternative run-time polymorphic interface for operators.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
int curDim
The current dimension of the reduction.
void ApplyPoly(const MV &x, MV &y) const
This routine takes the MV x and applies the polynomial operator phi(OP) to it resulting in the MV y...
void POTRF(const char &UPLO, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, OrdinalType *info) const
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > R
The current upper-triangular matrix from the QR reduction of H.
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
void factorWithEquilibration(bool flag)
void ApplyArnoldiPoly(const MV &x, MV &y) const
void MvNorm(std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec, NormType type=TwoNorm) const
Compute the norm of each vector in *this.
void solveWithTransposeFlag(Teuchos::ETransp trans)
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
GmresPolyOpOrthoFailure(const std::string &what_arg)
void SetBlock(const MultiVec< ScalarType > &A, const std::vector< int > &index)
Copy the vectors in A to a set of vectors in *this.
void ApplyRootsPoly(const MV &x, MV &y) const
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
Teuchos::RCP< const MV > getConstMV() const
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
Belos&#39;s class for applying the GMRES polynomial operator that is used by the hybrid-GMRES linear solv...
Belos concrete class for performing the block GMRES iteration.
OrdinalType numCols() const
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
NormType
The type of vector norm to compute.
Definition: BelosTypes.hpp:65
void MvInit(const ScalarType alpha)
Replace each element of the vectors in *this with alpha.
void MvDot(const MultiVec< ScalarType > &A, std::vector< ScalarType > &b) const
Compute the dot product of each column of *this with the corresponding column of A.
void MvTransMv(const ScalarType alpha, const MultiVec< ScalarType > &A, Teuchos::SerialDenseMatrix< int, ScalarType > &B) const
Compute a dense matrix B through the matrix-matrix multiply alpha * A^T * (*this).
int setVectors(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &X, const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &B)
Alternative run-time polymorphic interface for operators.
virtual ~GmresPolyOp()
Destructor.
void GEEV(const char &JOBVL, const char &JOBVR, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, MagnitudeType *WR, MagnitudeType *WI, ScalarType *VL, const OrdinalType &ldvl, ScalarType *VR, const OrdinalType &ldvr, ScalarType *WORK, const OrdinalType &lwork, MagnitudeType *RWORK, OrdinalType *info) const
GmresPolyMv * Clone(const int numvecs) const
Create a new MultiVec with numvecs columns.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Belos::StatusTest for specifying an implicit residual norm stopping criteria that checks for loss of ...
bool isType(const std::string &name) const
A class for extending the status testing capabilities of Belos via logical combinations.
Interface for multivectors used by Belos&#39; linear solvers.
GmresPolyOp(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem_in)
Given no ParameterList, constructor creates no polynomial and only applies the given operator...
static void MvNorm(const MV &mv, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec, NormType type=TwoNorm)
Compute the norm of each individual vector of mv. Upon return, normvec[i] holds the value of ...
void generateArnoldiPoly()
This routine takes the matrix, preconditioner, and vectors from the linear problem as well as the par...
static void MvScale(MV &mv, const ScalarType alpha)
Scale each element of the vectors in mv with alpha.
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:28
int shape(OrdinalType numRows, OrdinalType numCols)
Pure virtual base class which augments the basic interface for a Gmres linear solver iteration...
GmresIterationOrthoFailure is thrown when the GmresIteration object is unable to compute independent ...
Belos header file which uses auto-configuration information to include necessary C++ headers...
static void MvPrint(const MV &mv, std::ostream &os)
Print the mv multi-vector to the os output stream.
GmresPolyOpOrthoFailure is thrown when the orthogonalization manager is unable to generate orthonorma...
void MvPrint(std::ostream &os) const
Print *this multivector to the os output stream.
OrdinalType stride() const
void Apply(const MultiVec< ScalarType > &x, MultiVec< ScalarType > &y, ETrans=NOTRANS) const
This routine casts the MultiVec to GmresPolyMv to retrieve the MV. Then the above apply method is cal...
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > &params_in)
Process the passed in parameters.
OrdinalType numRows() const
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)
Interface for multivectors used by Belos&#39; linear solvers.
GmresPolyMv(const Teuchos::RCP< MV > &mv_in)
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
Teuchos::RCP< Belos::MatOrthoManager< Scalar, MV, OP > > makeMatOrthoManager(const std::string &ortho, const Teuchos::RCP< const OP > &M, const Teuchos::RCP< OutputManager< Scalar > > &, const std::string &label, const Teuchos::RCP< Teuchos::ParameterList > &params)
Return an instance of the specified MatOrthoManager subclass.

Generated on Mon Jul 15 2024 09:24:23 for Belos by doxygen 1.8.5