Belos Package Browser (Single Doxygen Collection)  Development
BelosGmresPolyOp.hpp
Go to the documentation of this file.
1 //@HEADER
2 // ************************************************************************
3 //
4 // Belos: Block Linear Solvers Package
5 // Copyright 2004 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 //@HEADER
41 
42 #ifndef BELOS_GMRESPOLYOP_HPP
43 #define BELOS_GMRESPOLYOP_HPP
44 
49 #include "BelosConfigDefs.hpp"
50 #include "BelosTypes.hpp"
51 
52 #include "BelosOperator.hpp"
53 #include "BelosMultiVec.hpp"
54 #include "BelosOperatorTraits.hpp"
55 #include "BelosMultiVecTraits.hpp"
56 #include "BelosLinearProblem.hpp"
57 
58 #include "BelosGmresIteration.hpp"
59 #include "BelosBlockGmresIter.hpp"
61 
65 #include "BelosStatusTestCombo.hpp"
67 
68 #include "BelosOutputManager.hpp"
69 
70 #include "Teuchos_BLAS.hpp"
71 #include "Teuchos_LAPACK.hpp"
72 #include "Teuchos_as.hpp"
73 #include "Teuchos_RCP.hpp"
74 #include "Teuchos_SerialDenseMatrix.hpp"
75 #include "Teuchos_SerialDenseVector.hpp"
76 #include "Teuchos_SerialDenseSolver.hpp"
77 #include "Teuchos_ParameterList.hpp"
78 
79 #ifdef BELOS_TEUCHOS_TIME_MONITOR
80  #include "Teuchos_TimeMonitor.hpp"
81 #endif // BELOS_TEUCHOS_TIME_MONITOR
82 
83 namespace Belos {
84 
91  class GmresPolyOpOrthoFailure : public BelosError {public:
92  GmresPolyOpOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
93  {}};
94 
95  // Create a shell class for the MV, inherited off MultiVec<> that will operate with the GmresPolyOp.
96  template <class ScalarType, class MV>
97  class GmresPolyMv : public MultiVec< ScalarType >
98  {
99  public:
100 
101  GmresPolyMv ( const Teuchos::RCP<MV>& mv_in )
102  : mv_(mv_in)
103  {}
104  GmresPolyMv ( const Teuchos::RCP<const MV>& mv_in )
105  {
106  mv_ = Teuchos::rcp_const_cast<MV>( mv_in );
107  }
108  Teuchos::RCP<MV> getMV() { return mv_; }
109  Teuchos::RCP<const MV> getConstMV() const { return mv_; }
110 
111  GmresPolyMv * Clone ( const int numvecs ) const
112  {
113  GmresPolyMv * newMV = new GmresPolyMv( MVT::Clone( *mv_, numvecs ) );
114  return newMV;
115  }
116  GmresPolyMv * CloneCopy () const
117  {
118  GmresPolyMv * newMV = new GmresPolyMv( MVT::CloneCopy( *mv_ ) );
119  return newMV;
120  }
121  GmresPolyMv * CloneCopy ( const std::vector<int>& index ) const
122  {
123  GmresPolyMv * newMV = new GmresPolyMv( MVT::CloneCopy( *mv_, index ) );
124  return newMV;
125  }
126  GmresPolyMv * CloneViewNonConst ( const std::vector<int>& index )
127  {
128  GmresPolyMv * newMV = new GmresPolyMv( MVT::CloneViewNonConst( *mv_, index ) );
129  return newMV;
130  }
131  const GmresPolyMv * CloneView ( const std::vector<int>& index ) const
132  {
133  const GmresPolyMv * newMV = new GmresPolyMv( MVT::CloneView( *mv_, index ) );
134  return newMV;
135  }
136  ptrdiff_t GetGlobalLength () const { return MVT::GetGlobalLength( *mv_ ); }
137  int GetNumberVecs () const { return MVT::GetNumberVecs( *mv_ ); }
138  void MvTimesMatAddMv (const ScalarType alpha,
139  const MultiVec<ScalarType>& A,
140  const Teuchos::SerialDenseMatrix<int,ScalarType>& B, const ScalarType beta)
141  {
142  const GmresPolyMv<ScalarType,MV>& A_in = dynamic_cast<const GmresPolyMv<ScalarType,MV>&>(A);
143  MVT::MvTimesMatAddMv( alpha, *(A_in.getConstMV()), B, beta, *mv_ );
144  }
145  void MvAddMv ( const ScalarType alpha, const MultiVec<ScalarType>& A, const ScalarType beta, const MultiVec<ScalarType>& B )
146  {
147  const GmresPolyMv<ScalarType,MV>& A_in = dynamic_cast<const GmresPolyMv<ScalarType,MV>&>(A);
148  const GmresPolyMv<ScalarType,MV>& B_in = dynamic_cast<const GmresPolyMv<ScalarType,MV>&>(B);
149  MVT::MvAddMv( alpha, *(A_in.getConstMV()), beta, *(B_in.getConstMV()), *mv_ );
150  }
151  void MvScale ( const ScalarType alpha ) { MVT::MvScale( *mv_, alpha ); }
152  void MvScale ( const std::vector<ScalarType>& alpha ) { MVT::MvScale( *mv_, alpha ); }
153  void MvTransMv ( const ScalarType alpha, const MultiVec<ScalarType>& A, Teuchos::SerialDenseMatrix<int,ScalarType>& B) const
154  {
155  const GmresPolyMv<ScalarType,MV>& A_in = dynamic_cast<const GmresPolyMv<ScalarType,MV>&>(A);
156  MVT::MvTransMv( alpha, *(A_in.getConstMV()), *mv_, B );
157  }
158  void MvDot ( const MultiVec<ScalarType>& A, std::vector<ScalarType>& b ) const
159  {
160  const GmresPolyMv<ScalarType,MV>& A_in = dynamic_cast<const GmresPolyMv<ScalarType,MV>&>(A);
161  MVT::MvDot( *(A_in.getConstMV()), *mv_, b );
162  }
163  void MvNorm ( std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& normvec, NormType type = TwoNorm ) const
164  {
165  MVT::MvNorm( *mv_, normvec, type );
166  }
167  void SetBlock ( const MultiVec<ScalarType>& A, const std::vector<int>& index )
168  {
169  const GmresPolyMv<ScalarType,MV>& A_in = dynamic_cast<const GmresPolyMv<ScalarType,MV>&>(A);
170  MVT::SetBlock( *(A_in.getConstMV()), index, *mv_ );
171  }
172  void MvRandom () { MVT::MvRandom( *mv_ ); }
173  void MvInit ( const ScalarType alpha ) { MVT::MvInit( *mv_, alpha ); }
174  void MvPrint ( std::ostream& os ) const { MVT::MvPrint( *mv_, os ); }
175 
176  private:
177 
179 
180  Teuchos::RCP<MV> mv_;
181 
182  };
183 
194  template <class ScalarType, class MV, class OP>
195  class GmresPolyOp : public Operator<ScalarType> {
196  public:
197 
199 
200 
202  GmresPolyOp( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >& problem_in,
203  const Teuchos::RCP<Teuchos::ParameterList>& params_in
204  )
205  : problem_(problem_in),
206  params_(params_in),
207  LP_(problem_in->getLeftPrec()),
208  RP_(problem_in->getRightPrec())
209  {
211 
212  polyUpdateLabel_ = label_ + ": Hybrid Gmres: Vector Update";
213 #ifdef BELOS_TEUCHOS_TIME_MONITOR
214  timerPolyUpdate_ = Teuchos::TimeMonitor::getNewCounter(polyUpdateLabel_);
215 #endif // BELOS_TEUCHOS_TIME_MONITOR
216 
217  if (polyType_ == "Arnoldi" || polyType_=="Roots")
219  else if (polyType_ == "Gmres")
221  else
222  TEUCHOS_TEST_FOR_EXCEPTION(polyType_!="Arnoldi"&&polyType_!="Gmres"&&polyType_!="Roots",std::invalid_argument,
223  "Belos::GmresPolyOp: \"Polynomial Type\" must be either \"Arnoldi\", \"Gmres\", or \"Roots\".");
224  }
225 
227  GmresPolyOp( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >& problem_in )
228  : problem_(problem_in)
229  {
230  // If dimension is zero, it will just apply the operator from problem_in in the Apply method.
231  dim_ = 0;
232  }
233 
235  virtual ~GmresPolyOp() {};
237 
239 
240 
242  void setParameters( const Teuchos::RCP<Teuchos::ParameterList>& params_in );
244 
246 
247 
251  void generateArnoldiPoly();
252 
256  void generateGmresPoly();
257 
259 
261 
262 
268  void ApplyPoly ( const MV& x, MV& y ) const;
269  void ApplyArnoldiPoly ( const MV& x, MV& y ) const;
270  void ApplyGmresPoly ( const MV& x, MV& y ) const;
271  void ApplyRootsPoly ( const MV& x, MV& y ) const;
272 
276  void Apply ( const MultiVec<ScalarType>& x, MultiVec<ScalarType>& y, ETrans /* trans */=NOTRANS ) const
277  {
278  const GmresPolyMv<ScalarType,MV>& x_in = dynamic_cast<const GmresPolyMv<ScalarType,MV>&>(x);
279  GmresPolyMv<ScalarType,MV>& y_in = dynamic_cast<GmresPolyMv<ScalarType,MV>&>(y);
280  ApplyPoly( *(x_in.getConstMV()), *(y_in.getMV()) );
281  }
282 
283  int polyDegree() const { return dim_; }
284 
285  private:
286 
287 #ifdef BELOS_TEUCHOS_TIME_MONITOR
288  Teuchos::RCP<Teuchos::Time> timerPolyUpdate_;
289 #endif // BELOS_TEUCHOS_TIME_MONITOR
290  std::string polyUpdateLabel_;
291 
292  typedef int OT; //Ordinal type
294  typedef Teuchos::ScalarTraits<ScalarType> SCT ;
295  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
296  typedef Teuchos::ScalarTraits<MagnitudeType> MCT ;
297 
298  // Default polynomial parameters
299  static constexpr int maxDegree_default_ = 25;
300  static constexpr int verbosity_default_ = Belos::Errors;
301  static constexpr bool randomRHS_default_ = true;
302  static constexpr const char * label_default_ = "Belos";
303  static constexpr const char * polyType_default_ = "Roots";
304  static constexpr const char * orthoType_default_ = "DGKS";
305 // https://stackoverflow.com/questions/24398102/constexpr-and-initialization-of-a-static-const-void-pointer-with-reinterpret-cas
306 #if defined(_WIN32) && defined(__clang__)
307  static constexpr std::ostream * outputStream_default_ =
308  __builtin_constant_p(reinterpret_cast<const std::ostream*>(&std::cout));
309 #else
310  static constexpr std::ostream * outputStream_default_ = &std::cout;
311 #endif
312  static constexpr bool damp_default_ = false;
313  static constexpr bool addRoots_default_ = true;
314 
315  // Variables for generating the polynomial
316  Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
317  Teuchos::RCP<Teuchos::ParameterList> params_;
318  Teuchos::RCP<const OP> LP_, RP_;
319 
320  // Output manager.
321  Teuchos::RCP<OutputManager<ScalarType> > printer_;
322  Teuchos::RCP<std::ostream> outputStream_ = Teuchos::rcp(outputStream_default_,false);
323 
324  // Orthogonalization manager.
325  Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
326 
327  // Current polynomial parameters
332  std::string label_ = label_default_;
335  int dim_ = 0;
338 
339  // Variables for Arnoldi polynomial
340  mutable Teuchos::RCP<MV> V_, wL_, wR_;
341  Teuchos::SerialDenseMatrix<OT,ScalarType> H_, y_;
342  Teuchos::SerialDenseVector<OT,ScalarType> r0_;
343 
344  // Variables for Gmres polynomial;
345  bool autoDeg = false;
346  Teuchos::SerialDenseMatrix< OT, ScalarType > pCoeff_;
347 
348  // Variables for Roots polynomial:
349  Teuchos::SerialDenseMatrix< OT, MagnitudeType > theta_;
350 
351  // Modified Leja sorting function. Takes a serial dense matrix of M harmonic Ritz values and an index
352  // of values from 0 to M. Returns the sorted values and sorted index, similar to Matlab.
353  void SortModLeja(Teuchos::SerialDenseMatrix< OT, MagnitudeType > &thetaN, std::vector<int> &index) const ;
354 
355  //Function determines whether added roots are needed and adds them if option is turned on.
356  void ComputeAddedRoots();
357  };
358 
359  template <class ScalarType, class MV, class OP>
360  void GmresPolyOp<ScalarType, MV, OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList>& params_in )
361  {
362  // Check which Gmres polynomial to use
363  if (params_in->isParameter("Polynomial Type")) {
364  polyType_ = params_in->get("Polynomial Type", polyType_default_);
365  }
366 
367  // Check for polynomial convergence tolerance
368  if (params_in->isParameter("Polynomial Tolerance")) {
369  if (params_in->isType<MagnitudeType> ("Polynomial Tolerance")) {
370  polyTol_ = params_in->get ("Polynomial Tolerance",
371  static_cast<MagnitudeType> (DefaultSolverParameters::polyTol));
372  }
373  else {
374  polyTol_ = params_in->get ("Polynomial Tolerance", DefaultSolverParameters::polyTol);
375  }
376  }
377 
378  // Check for maximum polynomial degree
379  if (params_in->isParameter("Maximum Degree")) {
380  maxDegree_ = params_in->get("Maximum Degree", maxDegree_default_);
381  }
382 
383  // Check for maximum polynomial degree
384  if (params_in->isParameter("Random RHS")) {
385  randomRHS_ = params_in->get("Random RHS", randomRHS_default_);
386  }
387 
388  // Check for a change in verbosity level
389  if (params_in->isParameter("Verbosity")) {
390  if (Teuchos::isParameterType<int>(*params_in,"Verbosity")) {
391  verbosity_ = params_in->get("Verbosity", verbosity_default_);
392  }
393  else {
394  verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params_in,"Verbosity");
395  }
396  }
397 
398  if (params_in->isParameter("Orthogonalization")) {
399  orthoType_ = params_in->get("Orthogonalization",orthoType_default_);
400  }
401 
402  // Check for timer label
403  if (params_in->isParameter("Timer Label")) {
404  label_ = params_in->get("Timer Label", label_default_);
405  }
406 
407  // Output stream
408  if (params_in->isParameter("Output Stream")) {
409  outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params_in,"Output Stream");
410  }
411 
412  // Check for damped polynomial
413  if (params_in->isParameter("Damped Poly")) {
414  damp_ = params_in->get("Damped Poly", damp_default_);
415  }
416 
417  // Check for root-adding
418  if (params_in->isParameter("Add Roots")) {
419  addRoots_ = params_in->get("Add Roots", addRoots_default_);
420  }
421  }
422 
423  template <class ScalarType, class MV, class OP>
425  {
426  Teuchos::RCP< MV > V = MVT::Clone( *problem_->getRHS(), maxDegree_+1 );
427 
428  //Make power basis:
429  std::vector<int> index(1,0);
430  Teuchos::RCP< MV > V0 = MVT::CloneViewNonConst(*V, index);
431  if (randomRHS_)
432  MVT::MvRandom( *V0 );
433  else
434  MVT::Assign( *problem_->getRHS(), *V0 );
435 
436  if ( !LP_.is_null() ) {
437  Teuchos::RCP< MV > Vtemp = MVT::CloneCopy(*V0);
438  problem_->applyLeftPrec( *Vtemp, *V0);
439  }
440  if ( damp_ ) {
441  Teuchos::RCP< MV > Vtemp = MVT::CloneCopy(*V0);
442  problem_->apply( *Vtemp, *V0);
443  }
444 
445  for(int i=0; i< maxDegree_; i++)
446  {
447  index[0] = i;
448  Teuchos::RCP< const MV > Vi = MVT::CloneView(*V, index);
449  index[0] = i+1;
450  Teuchos::RCP< MV > Vip1 = MVT::CloneViewNonConst(*V, index);
451  problem_->apply( *Vi, *Vip1);
452  }
453 
454  //Consider AV:
455  Teuchos::Range1D range( 1, maxDegree_);
456  Teuchos::RCP< const MV > AV = MVT::CloneView( *V, range);
457 
458  //Make lhs (AV)^T(AV)
459  Teuchos::SerialDenseMatrix< OT, ScalarType > AVtransAV( maxDegree_, maxDegree_);
460  MVT::MvTransMv( SCT::one(), *AV, *AV, AVtransAV);
461  //This process adds pDeg*pDeg + pDeg inner products that aren't in the final count.
462 
463  Teuchos::LAPACK< OT, ScalarType > lapack;
464  int infoInt;
465  bool status = true; //Keep adjusting poly deg when true.
466 
467  dim_ = maxDegree_;
468  Teuchos::SerialDenseMatrix< OT, ScalarType > lhs;
469  while( status && dim_ >= 1)
470  {
471  Teuchos::SerialDenseMatrix< OT, ScalarType > lhstemp(Teuchos::Copy, AVtransAV, dim_, dim_);
472  lapack.POTRF( 'U', dim_, lhstemp.values(), lhstemp.stride(), &infoInt);
473 
474  if(autoDeg == false)
475  {
476  status = false;
477  if(infoInt != 0)
478  {
479  std::cout << "BelosGmresPolyOp.hpp: LAPACK POTRF was not successful!!" << std::endl;
480  std::cout << "Error code: " << infoInt << std::endl;
481  }
482  }
483  else
484  {
485  if(infoInt != 0)
486  {//Had bad factor. Reduce poly degree.
487  dim_--;
488  }
489  else
490  {
491  status = false;
492  }
493  }
494  if(status == false)
495  {
496  lhs = lhstemp;
497  }
498  }
499  if(dim_ == 0)
500  {
501  pCoeff_.shape( 1, 1);
502  pCoeff_(0,0) = SCT::one();
503  std::cout << "Poly Degree is zero. No preconditioner created." << std::endl;
504  }
505  else
506  {
507  pCoeff_.shape( dim_, 1);
508  //Get correct submatrix of AV:
509  Teuchos::Range1D rangeSub( 1, dim_);
510  Teuchos::RCP< const MV > AVsub = MVT::CloneView( *V, rangeSub);
511 
512  //Compute rhs (AV)^T V0
513  MVT::MvTransMv( SCT::one(), *AVsub, *V0, pCoeff_);
514  lapack.POTRS( 'U', dim_, 1, lhs.values(), lhs.stride(), pCoeff_.values(), pCoeff_.stride(), &infoInt);
515  if(infoInt != 0)
516  {
517  std::cout << "BelosGmresPolyOp.hpp: LAPACK POTRS was not successful!!" << std::endl;
518  std::cout << "Error code: " << infoInt << std::endl;
519  }
520  }
521  }
522 
523  template <class ScalarType, class MV, class OP>
525  {
526  std::string polyLabel = label_ + ": GmresPolyOp creation";
527 
528  // Create a copy of the linear problem that has a zero initial guess and random RHS.
529  std::vector<int> idx(1,0);
530  Teuchos::RCP<MV> newX = MVT::Clone( *(problem_->getLHS()), 1 );
531  Teuchos::RCP<MV> newB = MVT::Clone( *(problem_->getRHS()), 1 );
532  MVT::MvInit( *newX, SCT::zero() );
533  if (randomRHS_) {
534  MVT::MvRandom( *newB );
535  }
536  else {
537  MVT::Assign( *(MVT::CloneView(*(problem_->getRHS()), idx)), *newB );
538  }
539  Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > newProblem =
540  Teuchos::rcp( new LinearProblem<ScalarType,MV,OP>( problem_->getOperator(), newX, newB ) );
541  newProblem->setInitResVec( newB );
542  newProblem->setLeftPrec( problem_->getLeftPrec() );
543  newProblem->setRightPrec( problem_->getRightPrec() );
544  newProblem->setLabel(polyLabel);
545  newProblem->setProblem();
546  newProblem->setLSIndex( idx );
547 
548  // Create a parameter list for the GMRES iteration.
549  Teuchos::ParameterList polyList;
550 
551  // Tell the block solver that the block size is one.
552  polyList.set("Num Blocks",maxDegree_);
553  polyList.set("Block Size",1);
554  polyList.set("Keep Hessenberg", true);
555 
556  // Create output manager.
557  printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
558 
559  // Create orthogonalization manager if we need to.
560  if (ortho_.is_null()) {
561  params_->set("Orthogonalization", orthoType_);
563  Teuchos::RCP<Teuchos::ParameterList> paramsOrtho; // can be null
564 
565  ortho_ = factory.makeMatOrthoManager (orthoType_, Teuchos::null, printer_, polyLabel, paramsOrtho);
566  }
567 
568  // Create a simple status test that either reaches the relative residual tolerance or maximum polynomial size.
569  Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxItrTst =
570  Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxDegree_ ) );
571 
572  // Implicit residual test, using the native residual to determine if convergence was achieved.
573  Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTst =
574  Teuchos::rcp( new StatusTestGenResNorm<ScalarType,MV,OP>( polyTol_ ) );
575  convTst->defineScaleForm( convertStringToScaleType("Norm of RHS"), Belos::TwoNorm );
576 
577  // Convergence test that stops the iteration when either are satisfied.
578  Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > polyTest =
579  Teuchos::rcp( new StatusTestCombo<ScalarType,MV,OP>( StatusTestCombo<ScalarType,MV,OP>::OR, maxItrTst, convTst ) );
580 
581  // Create Gmres iteration object to perform one cycle of Gmres.
582  Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > gmres_iter;
583  gmres_iter = Teuchos::rcp( new BlockGmresIter<ScalarType,MV,OP>(newProblem,printer_,polyTest,ortho_,polyList) );
584 
585  // Create the first block in the current Krylov basis (residual).
586  Teuchos::RCP<MV> V_0 = MVT::CloneCopy( *newB );
587  if ( !LP_.is_null() )
588  newProblem->applyLeftPrec( *newB, *V_0 );
589  if ( damp_ )
590  {
591  Teuchos::RCP< MV > Vtemp = MVT::CloneCopy(*V_0);
592  newProblem->apply( *Vtemp, *V_0 );
593  }
594 
595  // Get a matrix to hold the orthonormalization coefficients.
596  r0_.resize(1);
597 
598  // Orthonormalize the new V_0
599  int rank = ortho_->normalize( *V_0, Teuchos::rcpFromRef(r0_) );
600  TEUCHOS_TEST_FOR_EXCEPTION(rank != 1,GmresPolyOpOrthoFailure,
601  "Belos::GmresPolyOp::generateArnoldiPoly(): Failed to compute initial block of orthonormal vectors for polynomial generation.");
602 
603  // Set the new state and initialize the solver.
605  newstate.V = V_0;
606  newstate.z = Teuchos::rcpFromRef( r0_);
607  newstate.curDim = 0;
608  gmres_iter->initializeGmres(newstate);
609 
610  // Perform Gmres iteration
611  try {
612  gmres_iter->iterate();
613  }
614  catch (GmresIterationOrthoFailure& e) {
615  // Try to recover the most recent least-squares solution
616  gmres_iter->updateLSQR( gmres_iter->getCurSubspaceDim() );
617  }
618  catch (std::exception& e) {
619  using std::endl;
620  printer_->stream(Errors) << "Error! Caught exception in BlockGmresIter::iterate() at iteration "
621  << gmres_iter->getNumIters() << endl << e.what () << endl;
622  throw;
623  }
624 
625  // Get the solution for this polynomial, use in comparison below
626  Teuchos::RCP<MV> currX = gmres_iter->getCurrentUpdate();
627 
628  // Record polynomial info, get current GMRES state
629  GmresIterationState<ScalarType,MV> gmresState = gmres_iter->getState();
630 
631  // If the polynomial has no dimension, the tolerance is too low, return false
632  dim_ = gmresState.curDim;
633  if (dim_ == 0) {
634  return;
635  }
636  if(polyType_ == "Arnoldi"){
637  // Make a view and then copy the RHS of the least squares problem.
638  //
639  y_ = Teuchos::SerialDenseMatrix<OT,ScalarType>( Teuchos::Copy, *gmresState.z, dim_, 1 );
640  H_ = *gmresState.H;
641 
642  //
643  // Solve the least squares problem.
644  //
645  Teuchos::BLAS<OT,ScalarType> blas;
646  blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
647  Teuchos::NON_UNIT_DIAG, dim_, 1, SCT::one(),
648  gmresState.R->values(), gmresState.R->stride(),
649  y_.values(), y_.stride() );
650  }
651  else{ //Generate Roots Poly
652  //Find Harmonic Ritz Values to use as polynomial roots:
653 
654  //Copy of square H used to find poly roots:
655  H_ = Teuchos::SerialDenseMatrix<OT,ScalarType>(Teuchos::Copy, *gmresState.H, dim_, dim_);
656  //Zero out below subdiagonal of H:
657  for(int i=0; i <= dim_-3; i++) {
658  for(int k=i+2; k <= dim_-1; k++) {
659  H_(k,i) = SCT::zero();
660  }
661  }
662  //Extra copy of H because equilibrate changes the matrix:
663  Teuchos::SerialDenseMatrix<OT,ScalarType> Htemp (Teuchos::Copy, H_, dim_, dim_);
664 
665  //View the m+1,m element and last col of H:
666  ScalarType Hlast = (*gmresState.H)(dim_,dim_-1);
667  Teuchos::SerialDenseMatrix<OT,ScalarType> HlastCol (Teuchos::View, H_, dim_, 1, 0, dim_-1);
668 
669  //Set up linear system for H^{-*}e_m:
670  Teuchos::SerialDenseMatrix< OT, ScalarType > F(dim_,1), E(dim_,1);
671  E.putScalar(SCT::zero());
672  E(dim_-1,0) = SCT::one();
673 
674  Teuchos::SerialDenseSolver< OT, ScalarType > HSolver;
675  HSolver.setMatrix( Teuchos::rcpFromRef(Htemp));
676  HSolver.solveWithTransposeFlag( Teuchos::CONJ_TRANS );
677  HSolver.setVectors( Teuchos::rcpFromRef(F), Teuchos::rcpFromRef(E));
678  HSolver.factorWithEquilibration( true );
679 
680  //Factor matrix and solve for F = H^{-*}e_m:
681  int info = 0;
682  info = HSolver.factor();
683  if(info != 0){
684  std::cout << "Hsolver factor: info = " << info << std::endl;
685  }
686  info = HSolver.solve();
687  if(info != 0){
688  std::cout << "Hsolver solve : info = " << info << std::endl;
689  }
690 
691  //Scale F and adjust H for Harmonic Ritz value eigenproblem:
692  F.scale(Hlast*Hlast);
693  HlastCol += F;
694 
695  //Set up for eigenvalue problem to get Harmonic Ritz Values:
696  Teuchos::LAPACK< OT, ScalarType > lapack;
697  theta_.shape(dim_,2);//1st col for real part, 2nd col for imaginary
698 
699  const int ldv = 1;
700  ScalarType* vlr = 0;
701 
702  // Size of workspace and workspace for DGEEV
703  int lwork = -1;
704  std::vector<ScalarType> work(1);
705  std::vector<MagnitudeType> rwork(2*dim_);
706 
707  //Find workspace size for DGEEV:
708  lapack.GEEV('N','N',dim_,H_.values(),H_.stride(),theta_[0],theta_[1],vlr, ldv, vlr, ldv, &work[0], lwork, &rwork[0], &info);
709  lwork = std::abs (static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work[0])));
710  work.resize( lwork );
711  // Solve for Harmonic Ritz Values:
712  lapack.GEEV('N','N',dim_,H_.values(),H_.stride(),theta_[0],theta_[1],vlr, ldv, vlr, ldv, &work[0], lwork, &rwork[0], &info);
713 
714  if(info != 0){
715  std::cout << "GEEV solve : info = " << info << std::endl;
716  }
717 
718  // Set index for sort function, verify roots are non-zero,
719  // and sort Harmonic Ritz Values:
720  const MagnitudeType tol = 10.0 * Teuchos::ScalarTraits<MagnitudeType>::eps();
721  std::vector<int> index(dim_);
722  for(int i=0; i<dim_; ++i){
723  index[i] = i;
724  // Check if real + imag parts of roots < tol.
725  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.");
726  }
727  SortModLeja(theta_,index);
728 
729  //Add roots if neded.
730  ComputeAddedRoots();
731 
732  }
733  }
734 
735  //Function determines whether added roots are needed and adds them if option is turned on.
736  template <class ScalarType, class MV, class OP>
738  {
739  // Store theta (with cols for real and imag parts of Harmonic Ritz Vals)
740  // as one vector of complex numbers to perform arithmetic:
741  std::vector<std::complex<MagnitudeType>> cmplxHRitz (dim_);
742  for(unsigned int i=0; i<cmplxHRitz.size(); ++i){
743  cmplxHRitz[i] = std::complex<MagnitudeType>( theta_(i,0), theta_(i,1) );
744  }
745 
746  // Compute product of factors (pof) to determine added roots:
747  const MagnitudeType one(1.0);
748  std::vector<MagnitudeType> pof (dim_,one);
749  for(int j=0; j<dim_; ++j) {
750  for(int i=0; i<dim_; ++i) {
751  if(i!=j) {
752  pof[j] = std::abs(pof[j]*(one-(cmplxHRitz[j]/cmplxHRitz[i])));
753  }
754  }
755  }
756 
757  // Compute number of extra roots needed:
758  std::vector<int> extra (dim_);
759  int totalExtra = 0;
760  for(int i=0; i<dim_; ++i){
761  if (pof[i] > MCT::zero())
762  extra[i] = ceil((log10(pof[i])-MagnitudeType(4.0))/MagnitudeType(14.0));
763  else
764  extra[i] = 0;
765  if(extra[i] > 0){
766  totalExtra += extra[i];
767  }
768  }
769  if (totalExtra){
770  printer_->stream(Warnings) << "Warning: Need to add " << totalExtra << " extra roots." << std::endl;}
771 
772  // If requested to add roots, append them to the theta matrix:
773  if(addRoots_ && totalExtra>0)
774  {
775  theta_.reshape(dim_+totalExtra,2);
776  // Make a matrix copy for perturbed roots:
777  Teuchos::SerialDenseMatrix<OT,MagnitudeType> thetaPert (Teuchos::Copy, theta_, dim_+totalExtra, 2);
778 
779  //Add extra eigenvalues to matrix and perturb for sort:
780  int count = dim_;
781  for(int i=0; i<dim_; ++i){
782  for(int j=0; j< extra[i]; ++j){
783  theta_(count,0) = theta_(i,0);
784  theta_(count,1) = theta_(i,1);
785  thetaPert(count,0) = theta_(i,0)+(j+MCT::one())*MagnitudeType(5e-8);
786  thetaPert(count,1) = theta_(i,1);
787  ++count;
788  }
789  }
790 
791  // Update polynomial degree:
792  dim_ += totalExtra;
793  if (totalExtra){
794  printer_->stream(Warnings) << "New poly degree is: " << dim_ << std::endl;}
795 
796  // Create a new index and sort perturbed roots:
797  std::vector<int> index2(dim_);
798  for(int i=0; i<dim_; ++i){
799  index2[i] = i;
800  }
801  SortModLeja(thetaPert,index2);
802  //Apply sorting to non-perturbed roots:
803  for(int i=0; i<dim_; ++i)
804  {
805  thetaPert(i,0) = theta_(index2[i],0);
806  thetaPert(i,1) = theta_(index2[i],1);
807  }
808  theta_ = thetaPert;
809 
810  }
811  }
812 
813  // Modified Leja sorting function. Takes a serial dense matrix of M harmonic Ritz values and an index
814  // of values from 0 to M. Returns the sorted values and sorted index, similar to Matlab.
815  template <class ScalarType, class MV, class OP>
816  void GmresPolyOp<ScalarType, MV, OP>::SortModLeja(Teuchos::SerialDenseMatrix< OT, MagnitudeType > &thetaN, std::vector<int> &index) const
817  {
818  //Sort theta values via Modified Leja Ordering:
819 
820  // Set up blank matrices to track sorting:
821  int dimN = index.size();
822  std::vector<int> newIndex(dimN);
823  Teuchos::SerialDenseMatrix< OT, MagnitudeType > sorted (thetaN.numRows(), thetaN.numCols());
824  Teuchos::SerialDenseVector< OT, MagnitudeType > absVal (thetaN.numRows());
825  Teuchos::SerialDenseVector< OT, MagnitudeType > prod (thetaN.numRows());
826 
827  //Compute all absolute values and find maximum:
828  for(int i = 0; i < dimN; i++){
829  absVal(i) = hypot(thetaN(i,0), thetaN(i,1));
830  }
831  MagnitudeType * maxPointer = std::max_element(absVal.values(), (absVal.values()+dimN));
832  int maxIndex = int (maxPointer- absVal.values());
833 
834  //Put largest abs value first in the list:
835  sorted(0,0) = thetaN(maxIndex,0);
836  sorted(0,1) = thetaN(maxIndex,1);
837  newIndex[0] = index[maxIndex];
838 
839  int j;
840  // If largest value was complex (for real scalar type) put its conjugate in the next slot.
841  if(sorted(0,1)!= SCT::zero() && !SCT::isComplex)
842  {
843  sorted(1,0) = thetaN(maxIndex,0);
844  sorted(1,1) = -thetaN(maxIndex,1);
845  newIndex[1] = index[maxIndex+1];
846  j = 2;
847  }
848  else
849  {
850  j = 1;
851  }
852 
853  //Sort remaining values:
854  MagnitudeType a, b;
855  while( j < dimN )
856  {
857  //For each value, compute (a log of) a product of differences:
858  for(int i = 0; i < dimN; i++)
859  {
860  prod(i) = MCT::one();
861  for(int k = 0; k < j; k++)
862  {
863  a = thetaN(i,0) - sorted(k,0);
864  b = thetaN(i,1) - sorted(k,1);
865  if (a*a + b*b > MCT::zero())
866  prod(i) = prod(i) + log10(hypot(a,b));
867  else {
868  prod(i) = -std::numeric_limits<MagnitudeType>::infinity();
869  break;
870  }
871  }
872  }
873 
874  //Value with largest product goes in the next slot:
875  MagnitudeType * maxPointer = std::max_element(prod.values(), (prod.values()+dimN));
876  int maxIndex = int (maxPointer- prod.values());
877  sorted(j,0) = thetaN(maxIndex,0);
878  sorted(j,1) = thetaN(maxIndex,1);
879  newIndex[j] = index[maxIndex];
880 
881  //If it was complex (and scalar type real) put its conjugate in next slot:
882  if(sorted(j,1)!= SCT::zero() && !SCT::isComplex)
883  {
884  j++;
885  sorted(j,0) = thetaN(maxIndex,0);
886  sorted(j,1) = -thetaN(maxIndex,1);
887  newIndex[j] = index[maxIndex+1];
888  }
889  j++;
890  }
891 
892  //Return sorted values and sorted indices:
893  thetaN = sorted;
894  index = newIndex;
895  } //End Modified Leja ordering
896 
897  template <class ScalarType, class MV, class OP>
898  void GmresPolyOp<ScalarType, MV, OP>::ApplyPoly( const MV& x, MV& y ) const
899  {
900  if (dim_) {
901  if (polyType_ == "Arnoldi")
902  ApplyArnoldiPoly(x, y);
903  else if (polyType_ == "Gmres")
904  ApplyGmresPoly(x, y);
905  else if (polyType_ == "Roots")
906  ApplyRootsPoly(x, y);
907  }
908  else {
909  // Just apply the operator in problem_ to x and return y.
910  problem_->applyOp( x, y );
911  }
912  }
913 
914  template <class ScalarType, class MV, class OP>
915  void GmresPolyOp<ScalarType, MV, OP>::ApplyGmresPoly( const MV& x, MV& y ) const
916  {
917  Teuchos::RCP<MV> AX = MVT::CloneCopy(x);
918  Teuchos::RCP<MV> AX2 = MVT::Clone( x, MVT::GetNumberVecs(x) );
919 
920  // Apply left preconditioner.
921  if (!LP_.is_null()) {
922  Teuchos::RCP<MV> Xtmp = MVT::Clone( x, MVT::GetNumberVecs(x) );
923  problem_->applyLeftPrec( *AX, *Xtmp ); // Left precondition x into the first vector
924  AX = Xtmp;
925  }
926 
927  {
928 #ifdef BELOS_TEUCHOS_TIME_MONITOR
929  Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
930 #endif
931  MVT::MvAddMv(pCoeff_(0,0), *AX, SCT::zero(), y, y); //y= coeff_i(A^ix)
932  }
933  for( int i=1; i < dim_; i++)
934  {
935  Teuchos::RCP<MV> X, Y;
936  if ( i%2 )
937  {
938  X = AX;
939  Y = AX2;
940  }
941  else
942  {
943  X = AX2;
944  Y = AX;
945  }
946  problem_->apply(*X, *Y);
947  {
948 #ifdef BELOS_TEUCHOS_TIME_MONITOR
949  Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
950 #endif
951  MVT::MvAddMv(pCoeff_(i,0), *Y, SCT::one(), y, y); //y= coeff_i(A^ix) +y
952  }
953  }
954 
955  // Apply right preconditioner.
956  if (!RP_.is_null()) {
957  Teuchos::RCP<MV> Ytmp = MVT::CloneCopy(y);
958  problem_->applyRightPrec( *Ytmp, y );
959  }
960  }
961 
962  template <class ScalarType, class MV, class OP>
963  void GmresPolyOp<ScalarType, MV, OP>::ApplyRootsPoly( const MV& x, MV& y ) const
964  {
965  MVT::MvInit( y, SCT::zero() ); //Zero out y to take the vector with poly applied.
966  Teuchos::RCP<MV> prod = MVT::CloneCopy(x);
967  Teuchos::RCP<MV> Xtmp = MVT::Clone( x, MVT::GetNumberVecs(x) );
968  Teuchos::RCP<MV> Xtmp2 = MVT::Clone( x, MVT::GetNumberVecs(x) );
969 
970  // Apply left preconditioner.
971  if (!LP_.is_null()) {
972  problem_->applyLeftPrec( *prod, *Xtmp ); // Left precondition x into the first vector
973  prod = Xtmp;
974  }
975 
976  int i=0;
977  while(i < dim_-1)
978  {
979  if(theta_(i,1)== SCT::zero() || SCT::isComplex) //Real Harmonic Ritz value or complex scalars
980  {
981  {
982 #ifdef BELOS_TEUCHOS_TIME_MONITOR
983  Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
984 #endif
985  MVT::MvAddMv(SCT::one(), y, SCT::one()/theta_(i,0), *prod, y); //poly = poly + 1/theta_i * prod
986  }
987  problem_->apply(*prod, *Xtmp); // temp = A*prod
988  {
989 #ifdef BELOS_TEUCHOS_TIME_MONITOR
990  Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
991 #endif
992  MVT::MvAddMv(SCT::one(), *prod, -SCT::one()/theta_(i,0), *Xtmp, *prod); //prod = prod - 1/theta_i * temp
993  }
994  i++;
995  }
996  else //Current theta is complex and has a conjugate; combine to preserve real arithmetic
997  {
998  MagnitudeType mod = theta_(i,0)*theta_(i,0) + theta_(i,1)*theta_(i,1); //mod = a^2 + b^2
999  problem_->apply(*prod, *Xtmp); // temp = A*prod
1000  {
1001 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1002  Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
1003 #endif
1004  MVT::MvAddMv(2*theta_(i,0), *prod, -SCT::one(), *Xtmp, *Xtmp); //temp = 2a*prod-temp
1005  MVT::MvAddMv(SCT::one(), y, SCT::one()/mod, *Xtmp, y); //poly = poly + 1/mod*temp
1006  }
1007  if( i < dim_-2 )
1008  {
1009  problem_->apply(*Xtmp, *Xtmp2); // temp2 = A*temp
1010  {
1011 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1012  Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
1013 #endif
1014  MVT::MvAddMv(SCT::one(), *prod, -SCT::one()/mod, *Xtmp2, *prod); //prod = prod - 1/mod * temp2
1015  }
1016  }
1017  i = i + 2;
1018  }
1019  }
1020  if(theta_(dim_-1,1)== SCT::zero() || SCT::isComplex)
1021  {
1022 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1023  Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
1024 #endif
1025  MVT::MvAddMv(SCT::one(), y, SCT::one()/theta_(dim_-1,0), *prod, y); //poly = poly + 1/theta_i * prod
1026  }
1027 
1028  // Apply right preconditioner.
1029  if (!RP_.is_null()) {
1030  Teuchos::RCP<MV> Ytmp = MVT::CloneCopy(y);
1031  problem_->applyRightPrec( *Ytmp, y );
1032  }
1033  }
1034 
1035  template <class ScalarType, class MV, class OP>
1036  void GmresPolyOp<ScalarType, MV, OP>::ApplyArnoldiPoly( const MV& x, MV& y ) const
1037  {
1038  // Initialize vector storage.
1039  if (V_.is_null()) {
1040  V_ = MVT::Clone( x, dim_ );
1041  if (!LP_.is_null()) {
1042  wL_ = MVT::Clone( y, 1 );
1043  }
1044  if (!RP_.is_null()) {
1045  wR_ = MVT::Clone( y, 1 );
1046  }
1047  }
1048  //
1049  // Apply polynomial to x.
1050  //
1051  int n = MVT::GetNumberVecs( x );
1052  std::vector<int> idxi(1), idxi2, idxj(1);
1053 
1054  // Select vector x[j].
1055  for (int j=0; j<n; ++j) {
1056 
1057  idxi[0] = 0;
1058  idxj[0] = j;
1059  Teuchos::RCP<const MV> x_view = MVT::CloneView( x, idxj );
1060  Teuchos::RCP<MV> y_view = MVT::CloneViewNonConst( y, idxj );
1061  if (!LP_.is_null()) {
1062  Teuchos::RCP<MV> v_curr = MVT::CloneViewNonConst( *V_, idxi );
1063  problem_->applyLeftPrec( *x_view, *v_curr ); // Left precondition x into the first vector of V
1064  } else {
1065  MVT::SetBlock( *x_view, idxi, *V_ ); // Set x as the first vector of V
1066  }
1067 
1068  for (int i=0; i<dim_-1; ++i) {
1069 
1070  // Get views into the current and next vectors
1071  idxi2.resize(i+1);
1072  for (int ii=0; ii<i+1; ++ii) { idxi2[ii] = ii; }
1073  Teuchos::RCP<const MV> v_prev = MVT::CloneView( *V_, idxi2 );
1074  // the tricks below with wR_ and wL_ (potentially set to v_curr and v_next) unfortunately imply that
1075  // v_curr and v_next must be non-const views.
1076  Teuchos::RCP<MV> v_curr = MVT::CloneViewNonConst( *V_, idxi );
1077  idxi[0] = i+1;
1078  Teuchos::RCP<MV> v_next = MVT::CloneViewNonConst( *V_, idxi );
1079 
1080  //---------------------------------------------
1081  // Apply operator to next vector
1082  //---------------------------------------------
1083  // 1) Apply right preconditioner, if we have one.
1084  if (!RP_.is_null()) {
1085  problem_->applyRightPrec( *v_curr, *wR_ );
1086  } else {
1087  wR_ = v_curr;
1088  }
1089  // 2) Check for left preconditioner, if none exists, point at the next vector.
1090  if (LP_.is_null()) {
1091  wL_ = v_next;
1092  }
1093  // 3) Apply operator A.
1094  problem_->applyOp( *wR_, *wL_ );
1095  // 4) Apply left preconditioner, if we have one.
1096  if (!LP_.is_null()) {
1097  problem_->applyLeftPrec( *wL_, *v_next );
1098  }
1099 
1100  // Compute A*v_curr - v_prev*H(1:i,i)
1101  Teuchos::SerialDenseMatrix<OT,ScalarType> h(Teuchos::View,H_,i+1,1,0,i);
1102  {
1103 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1104  Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
1105 #endif
1106  MVT::MvTimesMatAddMv( -SCT::one(), *v_prev, h, SCT::one(), *v_next );
1107  }
1108 
1109  // Scale by H(i+1,i)
1110  MVT::MvScale( *v_next, SCT::one()/H_(i+1,i) );
1111  }
1112 
1113  // Compute output y = V*y_./r0_
1114  if (!RP_.is_null()) {
1115  {
1116 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1117  Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
1118 #endif
1119  MVT::MvTimesMatAddMv( SCT::one()/r0_(0), *V_, y_, SCT::zero(), *wR_ );
1120  }
1121  problem_->applyRightPrec( *wR_, *y_view );
1122  }
1123  else {
1124 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1125  Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
1126 #endif
1127  MVT::MvTimesMatAddMv( SCT::one()/r0_(0), *V_, y_, SCT::zero(), *y_view );
1128  }
1129  } // (int j=0; j<n; ++j)
1130  } // end Apply()
1131 } // end Belos namespace
1132 
1133 #endif
1134 
1135 // end of file BelosGmresPolyOp.hpp
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
Definition: BelosTypes.cpp:106
GmresPolyMv(const Teuchos::RCP< const MV > &mv_in)
Collection of types and exceptions used within the Belos solvers.
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
static constexpr bool randomRHS_default_
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.
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...
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 ...
Teuchos::RCP< MV > wR_
void MvAddMv(const ScalarType alpha, const MultiVec< ScalarType > &A, const ScalarType beta, const MultiVec< ScalarType > &B)
Replace *this with alpha * A + beta * B.
void ApplyGmresPoly(const MV &x, MV &y) const
void MvScale(const ScalarType alpha)
Scale each element of the vectors in *this with alpha.
Teuchos::RCP< OutputManager< ScalarType > > printer_
static constexpr const char * label_default_
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()
Teuchos::RCP< MV > wL_
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...
MultiVecTraits< ScalarType, MV > MVT
Declaration of basic traits for the multivector type.
GmresPolyMv * CloneCopy() const
Create a new MultiVec and copy contents of *this into it (deep copy).
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...
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 * Clone(const int numvecs) const
Create a new MultiVec with numvecs columns.
Teuchos::RCP< const OP > RP_
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...
MultiVecTraits< ScalarType, MV > MVT
void MvRandom()
Fill all the vectors in *this with random numbers.
Structure to contain pointers to GmresIteration state variables.
static constexpr int verbosity_default_
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
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:296
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.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > ortho_
A factory class for generating StatusTestOutput objects.
Teuchos::RCP< const MV > getConstMV() const
static constexpr const char * polyType_default_
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
void MvNorm(std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec, NormType type=TwoNorm) const
Compute the norm of each vector in *this.
Teuchos::ScalarTraits< MagnitudeType > MCT
ETrans
Whether to apply the (conjugate) transpose of an operator.
Definition: BelosTypes.hpp:81
void ApplyRootsPoly(const MV &x, MV &y) const
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
int GetNumberVecs() const
The number of vectors (i.e., columns) in the multivector.
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...
void SortModLeja(Teuchos::SerialDenseMatrix< OT, MagnitudeType > &thetaN, std::vector< int > &index) const
static constexpr std::ostream * outputStream_default_
A Belos::StatusTest class for specifying a maximum number of iterations.
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.
static constexpr int maxDegree_default_
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::SerialDenseMatrix< OT, ScalarType > y_
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 curDim
The current dimension of the reduction.
void MvPrint(std::ostream &os) const
Print *this multivector to the os output stream.
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...
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > R
The current upper-triangular matrix from the QR reduction of H.
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...
A linear system to solve, and its associated information.
Teuchos::RCP< const OP > LP_
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.
Teuchos::RCP< MV > V_
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
Teuchos::SerialDenseMatrix< OT, MagnitudeType > theta_
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
void ApplyArnoldiPoly(const MV &x, MV &y) const
Belos&#39;s class for applying the GMRES polynomial operator that is used by the hybrid-GMRES linear solv...
Teuchos::SerialDenseMatrix< OT, ScalarType > pCoeff_
Belos concrete class for performing the block GMRES iteration.
Teuchos::RCP< MV > mv_
std::string polyUpdateLabel_
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
Teuchos::ScalarTraits< ScalarType > SCT
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:97
void MvInit(const ScalarType alpha)
Replace each element of the vectors in *this with alpha.
Alternative run-time polymorphic interface for operators.
virtual ~GmresPolyOp()
Destructor.
Teuchos::RCP< Teuchos::ParameterList > params_
static constexpr bool damp_default_
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Belos::StatusTest for specifying an implicit residual norm stopping criteria that checks for loss of ...
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:60
static constexpr bool addRoots_default_
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.
static constexpr const char * orthoType_default_
GmresPolyOpOrthoFailure is thrown when the orthogonalization manager is unable to generate orthonorma...
Teuchos::SerialDenseMatrix< OT, ScalarType > H_
Teuchos::RCP< std::ostream > outputStream_
Teuchos::SerialDenseVector< OT, ScalarType > r0_
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > &params_in)
Process the passed in parameters.
ptrdiff_t GetGlobalLength() const
The number of rows in the multivector.
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.