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 //
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"
60 
64 
68 #include "BelosStatusTestCombo.hpp"
70 
71 #include "BelosOutputManager.hpp"
72 
73 #include "Teuchos_BLAS.hpp"
74 #include "Teuchos_LAPACK.hpp"
75 #include "Teuchos_as.hpp"
76 #include "Teuchos_RCP.hpp"
80 
81 namespace Belos {
82 
89  class GmresPolyOpOrthoFailure : public BelosError {public:
90  GmresPolyOpOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
91  {}};
92 
93  // Create a shell class for the MV, inherited off MultiVec<> that will operate with the GmresPolyOp.
94  template <class ScalarType, class MV>
95  class GmresPolyMv : public MultiVec< ScalarType >
96  {
97  public:
98 
99  GmresPolyMv ( const Teuchos::RCP<MV>& mv_in )
100  : mv_(mv_in)
101  {}
103  {
104  mv_ = Teuchos::rcp_const_cast<MV>( mv_in );
105  }
106  Teuchos::RCP<MV> getMV() { return mv_; }
107  Teuchos::RCP<const MV> getConstMV() const { return mv_; }
108 
109  GmresPolyMv * Clone ( const int numvecs ) const
110  {
111  GmresPolyMv * newMV = new GmresPolyMv( MVT::Clone( *mv_, numvecs ) );
112  return newMV;
113  }
114  GmresPolyMv * CloneCopy () const
115  {
116  GmresPolyMv * newMV = new GmresPolyMv( MVT::CloneCopy( *mv_ ) );
117  return newMV;
118  }
119  GmresPolyMv * CloneCopy ( const std::vector<int>& index ) const
120  {
121  GmresPolyMv * newMV = new GmresPolyMv( MVT::CloneCopy( *mv_, index ) );
122  return newMV;
123  }
124  GmresPolyMv * CloneViewNonConst ( const std::vector<int>& index )
125  {
126  GmresPolyMv * newMV = new GmresPolyMv( MVT::CloneViewNonConst( *mv_, index ) );
127  return newMV;
128  }
129  const GmresPolyMv * CloneView ( const std::vector<int>& index ) const
130  {
131  const GmresPolyMv * newMV = new GmresPolyMv( MVT::CloneView( *mv_, index ) );
132  return newMV;
133  }
134  ptrdiff_t GetGlobalLength () const { return MVT::GetGlobalLength( *mv_ ); }
135  int GetNumberVecs () const { return MVT::GetNumberVecs( *mv_ ); }
136  void MvTimesMatAddMv (const ScalarType alpha,
137  const MultiVec<ScalarType>& A,
138  const Teuchos::SerialDenseMatrix<int,ScalarType>& B, const ScalarType beta)
139  {
140  const GmresPolyMv<ScalarType,MV>& A_in = dynamic_cast<const GmresPolyMv<ScalarType,MV>&>(A);
141  MVT::MvTimesMatAddMv( alpha, *(A_in.getConstMV()), B, beta, *mv_ );
142  }
143  void MvAddMv ( const ScalarType alpha, const MultiVec<ScalarType>& A, const ScalarType beta, const MultiVec<ScalarType>& B )
144  {
145  const GmresPolyMv<ScalarType,MV>& A_in = dynamic_cast<const GmresPolyMv<ScalarType,MV>&>(A);
146  const GmresPolyMv<ScalarType,MV>& B_in = dynamic_cast<const GmresPolyMv<ScalarType,MV>&>(B);
147  MVT::MvAddMv( alpha, *(A_in.getConstMV()), beta, *(B_in.getConstMV()), *mv_ );
148  }
149  void MvScale ( const ScalarType alpha ) { MVT::MvScale( *mv_, alpha ); }
150  void MvScale ( const std::vector<ScalarType>& alpha ) { MVT::MvScale( *mv_, alpha ); }
151  void MvTransMv ( const ScalarType alpha, const MultiVec<ScalarType>& A, Teuchos::SerialDenseMatrix<int,ScalarType>& B) const
152  {
153  const GmresPolyMv<ScalarType,MV>& A_in = dynamic_cast<const GmresPolyMv<ScalarType,MV>&>(A);
154  MVT::MvTransMv( alpha, *(A_in.getConstMV()), *mv_, B );
155  }
156  void MvDot ( const MultiVec<ScalarType>& A, std::vector<ScalarType>& b ) const
157  {
158  const GmresPolyMv<ScalarType,MV>& A_in = dynamic_cast<const GmresPolyMv<ScalarType,MV>&>(A);
159  MVT::MvDot( *(A_in.getConstMV()), *mv_, b );
160  }
161  void MvNorm ( std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& normvec, NormType type = TwoNorm ) const
162  {
163  MVT::MvNorm( *mv_, normvec, type );
164  }
165  void SetBlock ( const MultiVec<ScalarType>& A, const std::vector<int>& index )
166  {
167  const GmresPolyMv<ScalarType,MV>& A_in = dynamic_cast<const GmresPolyMv<ScalarType,MV>&>(A);
168  MVT::SetBlock( *(A_in.getConstMV()), index, *mv_ );
169  }
170  void MvRandom () { MVT::MvRandom( *mv_ ); }
171  void MvInit ( const ScalarType alpha ) { MVT::MvInit( *mv_, alpha ); }
172  void MvPrint ( std::ostream& os ) const { MVT::MvPrint( *mv_, os ); }
173 
174  private:
175 
176  typedef MultiVecTraits<ScalarType,MV> MVT;
177 
178  Teuchos::RCP<MV> mv_;
179 
180  };
181 
192  template <class ScalarType, class MV, class OP>
193  class GmresPolyOp : public Operator<ScalarType> {
194  public:
195 
197 
198 
201  const Teuchos::RCP<Teuchos::ParameterList>& params_in
202  )
203  : problem_(problem_in),
204  params_(params_in),
205  LP_(problem_in->getLeftPrec()),
206  RP_(problem_in->getRightPrec()),
207  outputStream_ (Teuchos::rcp(outputStream_default_,false)),
208  polyTol_ (DefaultSolverParameters::polyTol),
209  maxDegree_ (maxDegree_default_),
210  verbosity_ (verbosity_default_),
211  randomRHS_ (randomRHS_default_),
212  label_ (label_default_),
213  polyType_ (polyType_default_),
214  orthoType_ (orthoType_default_),
215  dim_(0)
216  {
217  setParameters( params_ );
218 
219  if (polyType_ == "Arnoldi")
221  else if (polyType_ == "Gmres")
223  else
224  TEUCHOS_TEST_FOR_EXCEPTION(polyType_!="Arnoldi"&&polyType_!="Gmres",std::logic_error,
225  "Belos::GmresPolyOp(): Invalid polynomial type.");
226  }
227 
230  : problem_(problem_in)
231  {
232  // If dimension is zero, it will just apply the operator from problem_in in the Apply method.
233  dim_ = 0;
234  }
235 
237  virtual ~GmresPolyOp() {};
239 
241 
242 
244  void setParameters( const Teuchos::RCP<Teuchos::ParameterList>& params_in );
246 
248 
249 
253  void generateArnoldiPoly();
254 
258  void generateGmresPoly();
259 
261 
263 
264 
270  void ApplyPoly ( const MV& x, MV& y ) const;
271  void ApplyArnoldiPoly ( const MV& x, MV& y ) const;
272  void ApplyGmresPoly ( const MV& x, MV& y ) const;
273 
277  void Apply ( const MultiVec<ScalarType>& x, MultiVec<ScalarType>& y, ETrans /* trans */=NOTRANS ) const
278  {
279  const GmresPolyMv<ScalarType,MV>& x_in = dynamic_cast<const GmresPolyMv<ScalarType,MV>&>(x);
280  GmresPolyMv<ScalarType,MV>& y_in = dynamic_cast<GmresPolyMv<ScalarType,MV>&>(y);
281  ApplyPoly( *(x_in.getConstMV()), *(y_in.getMV()) );
282  }
283 
284  int polyDegree() const { return dim_; }
285 
286  private:
287 
289  GmresPolyOp() {}
290 
291  typedef int OT; //Ordinal type
292  typedef MultiVecTraits<ScalarType,MV> MVT;
294  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
295 
296  // Default polynomial parameters
297  static constexpr int maxDegree_default_ = 25;
298  static constexpr int verbosity_default_ = Belos::Errors;
299  static constexpr bool randomRHS_default_ = true;
300  static constexpr const char * label_default_ = "Belos";
301  static constexpr const char * polyType_default_ = "Arnoldi";
302  static constexpr const char * orthoType_default_ = "DGKS";
303  static constexpr std::ostream * outputStream_default_ = &std::cout;
304 
305  // Variables for generating the polynomial
308  Teuchos::RCP<const OP> LP_, RP_;
309 
310  // Output manager.
312  Teuchos::RCP<std::ostream> outputStream_;
313 
314  // Orthogonalization manager.
316 
317  // Current polynomial parameters
318  MagnitudeType polyTol_;
319  int maxDegree_;
320  int verbosity_;
321  bool randomRHS_;
322  std::string label_;
323  std::string polyType_;
324  std::string orthoType_;
325  int dim_;
326 
327  // Variables for Arnoldi polynomial
328  mutable Teuchos::RCP<MV> V_, wL_, wR_;
331 
332  // Variables for Gmres polynomial;
333  bool autoDeg = false;
335 
336  };
337 
338  template <class ScalarType, class MV, class OP>
340  {
341  // Check which Gmres polynomial to use
342  if (params_in->isParameter("Polynomial Type")) {
343  polyType_ = params_in->get("Polynomial Type", polyType_default_);
344  }
345 
346  // Check for polynomial convergence tolerance
347  if (params_in->isParameter("Polynomial Tolerance")) {
348  if (params_in->isType<MagnitudeType> ("Polynomial Tolerance")) {
349  polyTol_ = params_in->get ("Polynomial Tolerance",
350  static_cast<MagnitudeType> (DefaultSolverParameters::polyTol));
351  }
352  else {
353  polyTol_ = params_in->get ("Polynomial Tolerance", DefaultSolverParameters::polyTol);
354  }
355  }
356 
357  // Check for maximum polynomial degree
358  if (params_in->isParameter("Maximum Degree")) {
359  maxDegree_ = params_in->get("Maximum Degree", maxDegree_default_);
360  }
361 
362  // Check for maximum polynomial degree
363  if (params_in->isParameter("Random RHS")) {
364  randomRHS_ = params_in->get("Random RHS", randomRHS_default_);
365  }
366 
367  // Check for a change in verbosity level
368  if (params_in->isParameter("Verbosity")) {
369  if (Teuchos::isParameterType<int>(*params_in,"Verbosity")) {
370  verbosity_ = params_in->get("Verbosity", verbosity_default_);
371  }
372  else {
373  verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params_in,"Verbosity");
374  }
375  }
376 
377  if (params_in->isParameter("Orthogonalization")) {
378  orthoType_ = params_in->get("Orthogonalization",orthoType_default_);
379  TEUCHOS_TEST_FOR_EXCEPTION( orthoType_ != "DGKS" && orthoType_ != "ICGS" && orthoType_ != "IMGS",
380  std::invalid_argument,
381  "Belos::GmresPolyOp: \"Orthogonalization\" must be either \"DGKS\", \"ICGS\", or \"IMGS\".");
382  }
383 
384  // Check for timer label
385  if (params_in->isParameter("Timer Label")) {
386  label_ = params_in->get("Timer Label", label_default_);
387  }
388 
389  // Output stream
390  if (params_in->isParameter("Output Stream")) {
391  outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params_in,"Output Stream");
392  }
393  }
394 
395  template <class ScalarType, class MV, class OP>
397  {
398  Teuchos::RCP< MV > V = MVT::Clone( *problem_->getRHS(), maxDegree_+2 );
399 
400  //Make power basis:
401  std::vector<int> index(1,0);
402  Teuchos::RCP< MV > V0 = MVT::CloneViewNonConst(*V, index);
403  if (randomRHS_)
404  MVT::MvRandom( *V0 );
405  else
406  MVT::Assign( *problem_->getRHS(), *V0 );
407 
408  if ( !LP_.is_null() )
409  {
410  Teuchos::RCP< MV > Vtemp = MVT::CloneCopy(*V0);
411  problem_->applyLeftPrec( *Vtemp, *V0);
412  }
413 
414  for(int i=0; i< maxDegree_+1; i++)
415  {
416  index[0] = i;
417  Teuchos::RCP< const MV > Vi = MVT::CloneView(*V, index);
418  index[0] = i+1;
419  Teuchos::RCP< MV > Vip1 = MVT::CloneViewNonConst(*V, index);
420  problem_->apply( *Vi, *Vip1);
421  }
422 
423  //Consider AV:
424  Teuchos::Range1D range( 1, maxDegree_+1);
425  Teuchos::RCP< const MV > AV = MVT::CloneView( *V, range);
426 
427  //Make lhs (AV)^T(AV)
428  Teuchos::SerialDenseMatrix< OT, ScalarType > AVtransAV( maxDegree_+1, maxDegree_+1);
429  MVT::MvTransMv( SCT::one(), *AV, *AV, AVtransAV);
430  //This process adds pDeg*pDeg + pDeg inner products that aren't in the final count.
431 
433  int infoInt;
434  bool status = true; //Keep adjusting poly deg when true.
435 
436  dim_ = maxDegree_;
438  while( status && dim_ >= 1)
439  {
440  Teuchos::SerialDenseMatrix< OT, ScalarType > lhstemp(Teuchos::Copy, AVtransAV, dim_+1, dim_+1);
441  lapack.POTRF( 'U', dim_+1, lhstemp.values(), lhstemp.stride(), &infoInt);
442 
443  if(autoDeg == false)
444  {
445  status = false;
446  if(infoInt != 0)
447  {
448  std::cout << "BelosGmresPolyOp.hpp: LAPACK POTRF was not successful!!" << std::endl;
449  std::cout << "Error code: " << infoInt << std::endl;
450  }
451  }
452  else
453  {
454  if(infoInt != 0)
455  {//Had bad factor. Reduce poly degree.
456  dim_--;
457  }
458  else
459  {
460  status = false;
461  }
462  }
463  if(status == false)
464  {
465  lhs = lhstemp;
466  }
467  }
468  if(dim_ == 0)
469  {
470  pCoeff_.shape( 1, 1);
471  pCoeff_(0,0) = SCT::one();
472  std::cout << "Poly Degree is zero. No preconditioner created." << std::endl;
473  }
474  else
475  {
476  pCoeff_.shape( dim_+1, 1);
477  //Get correct submatrix of AV:
478  Teuchos::Range1D rangeSub( 1, dim_+1);
479  Teuchos::RCP< const MV > AVsub = MVT::CloneView( *V, rangeSub);
480 
481  //Compute rhs (AV)^T V0
482  MVT::MvTransMv( SCT::one(), *AVsub, *V0, pCoeff_);
483  lapack.POTRS( 'U', dim_+1, 1, lhs.values(), lhs.stride(), pCoeff_.values(), pCoeff_.stride(), &infoInt);
484  if(infoInt != 0)
485  {
486  std::cout << "BelosGmresPolyOp.hpp: LAPACK POTRS was not successful!!" << std::endl;
487  std::cout << "Error code: " << infoInt << std::endl;
488  }
489  }
490  }
491 
492  template <class ScalarType, class MV, class OP>
494  {
495  std::string polyLabel = label_ + ": GmresPolyOp creation";
496 
497  // Create a copy of the linear problem that has a zero initial guess and random RHS.
498  std::vector<int> idx(1,0);
499  Teuchos::RCP<MV> newX = MVT::Clone( *(problem_->getLHS()), 1 );
500  Teuchos::RCP<MV> newB = MVT::Clone( *(problem_->getRHS()), 1 );
501  MVT::MvInit( *newX, SCT::zero() );
502  if (randomRHS_) {
503  MVT::MvRandom( *newB );
504  }
505  else {
506  MVT::Assign( *(MVT::CloneView(*(problem_->getRHS()), idx)), *newB );
507  }
509  Teuchos::rcp( new LinearProblem<ScalarType,MV,OP>( problem_->getOperator(), newX, newB ) );
510  newProblem->setInitResVec( newB );
511  newProblem->setLeftPrec( problem_->getLeftPrec() );
512  newProblem->setRightPrec( problem_->getRightPrec() );
513  newProblem->setLabel(polyLabel);
514  newProblem->setProblem();
515  newProblem->setLSIndex( idx );
516 
517  // Create a parameter list for the GMRES iteration.
518  Teuchos::ParameterList polyList;
519 
520  // Tell the block solver that the block size is one.
521  polyList.set("Num Blocks",maxDegree_);
522  polyList.set("Block Size",1);
523  polyList.set("Keep Hessenberg", true);
524 
525  // Create orthogonalization manager if we need to.
526  if (ortho_.is_null()) {
527  params_->set("Orthogonalization", orthoType_);
528  if (orthoType_=="DGKS") {
529  ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( polyLabel ) );
530  }
531  else if (orthoType_=="ICGS") {
532  ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( polyLabel ) );
533  }
534  else if (orthoType_=="IMGS") {
535  ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( polyLabel ) );
536  }
537  else {
538  TEUCHOS_TEST_FOR_EXCEPTION(orthoType_!="ICGS"&&orthoType_!="DGKS"&&orthoType_!="IMGS",std::logic_error,
539  "Belos::GmresPolyOp(): Invalid orthogonalization type.");
540  }
541  }
542 
543  // Create output manager.
544  printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
545 
546  // Create a simple status test that either reaches the relative residual tolerance or maximum polynomial size.
549 
550  // Implicit residual test, using the native residual to determine if convergence was achieved.
553  convTst->defineScaleForm( convertStringToScaleType("Norm of RHS"), Belos::TwoNorm );
554 
555  // Convergence test that stops the iteration when either are satisfied.
558 
559  // Create Gmres iteration object to perform one cycle of Gmres.
561  gmres_iter = Teuchos::rcp( new BlockGmresIter<ScalarType,MV,OP>(newProblem,printer_,polyTest,ortho_,polyList) );
562 
563  // Create the first block in the current Krylov basis (residual).
564  Teuchos::RCP<MV> V_0 = MVT::CloneCopy( *newB );
565  if ( !LP_.is_null() )
566  problem_->applyLeftPrec( *newB, *V_0 );
567 
568  // Get a matrix to hold the orthonormalization coefficients.
569  r0_.resize(1);
570 
571  // Orthonormalize the new V_0
572  int rank = ortho_->normalize( *V_0, Teuchos::rcp( &r0_, false ) );
574  "Belos::GmresPolyOp::generateArnoldiPoly(): Failed to compute initial block of orthonormal vectors for polynomial generation.");
575 
576  // Set the new state and initialize the solver.
578  newstate.V = V_0;
579  newstate.z = Teuchos::rcp( &r0_,false );
580  newstate.curDim = 0;
581  gmres_iter->initializeGmres(newstate);
582 
583  // Perform Gmres iteration
584  try {
585  gmres_iter->iterate();
586  }
587  catch (GmresIterationOrthoFailure e) {
588  // Try to recover the most recent least-squares solution
589  gmres_iter->updateLSQR( gmres_iter->getCurSubspaceDim() );
590  }
591  catch (std::exception e) {
592  using std::endl;
593  printer_->stream(Errors) << "Error! Caught exception in BlockGmresIter::iterate() at iteration "
594  << gmres_iter->getNumIters() << endl << e.what () << endl;
595  throw;
596  }
597 
598  // Get the solution for this polynomial, use in comparison below
599  Teuchos::RCP<MV> currX = gmres_iter->getCurrentUpdate();
600 
601  // Record polynomial info, get current GMRES state
602  GmresIterationState<ScalarType,MV> gmresState = gmres_iter->getState();
603 
604  // If the polynomial has no dimension, the tolerance is too low, return false
605  dim_ = gmresState.curDim;
606  if (dim_ == 0) {
607  return;
608  }
609 
610  // Make a view and then copy the RHS of the least squares problem.
611  //
612  y_ = Teuchos::SerialDenseMatrix<OT,ScalarType>( Teuchos::Copy, *gmresState.z, dim_, 1 );
613  H_ = *gmresState.H;
614 
615  //
616  // Solve the least squares problem.
617  //
620  Teuchos::NON_UNIT_DIAG, dim_, 1, SCT::one(),
621  gmresState.R->values(), gmresState.R->stride(),
622  y_.values(), y_.stride() );
623 
624  }
625 
626  template <class ScalarType, class MV, class OP>
627  void GmresPolyOp<ScalarType, MV, OP>::ApplyPoly( const MV& x, MV& y ) const
628  {
629  if (dim_) {
630  if (polyType_ == "Arnoldi")
631  ApplyArnoldiPoly(x, y);
632  else if (polyType_ == "Gmres")
633  ApplyGmresPoly(x, y);
634  }
635  else {
636  // Just apply the operator in problem_ to x and return y.
637  problem_->applyOp( x, y );
638  }
639  }
640 
641  template <class ScalarType, class MV, class OP>
642  void GmresPolyOp<ScalarType, MV, OP>::ApplyGmresPoly( const MV& x, MV& y ) const
643  {
644  Teuchos::RCP<MV> AX = MVT::CloneCopy(x);
645  Teuchos::RCP<MV> AX2 = MVT::Clone( x, MVT::GetNumberVecs(x) );
646 
647  // Apply left preconditioner.
648  if (!LP_.is_null()) {
649  Teuchos::RCP<MV> Xtmp = MVT::Clone( x, MVT::GetNumberVecs(x) );
650  problem_->applyLeftPrec( *AX, *Xtmp ); // Left precondition x into the first vector
651  AX = Xtmp;
652  }
653 
654  MVT::MvAddMv(pCoeff_(0,0), *AX, SCT::zero(), y, y); //y= coeff_i(A^ix)
655  for( int i=1; i < dim_+1; i++)
656  {
657  Teuchos::RCP<MV> X, Y;
658  if ( i%2 )
659  {
660  X = AX;
661  Y = AX2;
662  }
663  else
664  {
665  X = AX2;
666  Y = AX;
667  }
668  problem_->apply(*X, *Y);
669  MVT::MvAddMv(pCoeff_(i,0), *Y, SCT::one(), y, y); //y= coeff_i(A^ix) +y
670  }
671 
672  // Apply right preconditioner.
673  if (!RP_.is_null()) {
674  Teuchos::RCP<MV> Ytmp = MVT::CloneCopy(y);
675  problem_->applyRightPrec( *Ytmp, y );
676  }
677  }
678 
679  template <class ScalarType, class MV, class OP>
680  void GmresPolyOp<ScalarType, MV, OP>::ApplyArnoldiPoly( const MV& x, MV& y ) const
681  {
682  // Initialize vector storage.
683  if (V_.is_null()) {
684  V_ = MVT::Clone( x, dim_ );
685  if (!LP_.is_null()) {
686  wL_ = MVT::Clone( y, 1 );
687  }
688  if (!RP_.is_null()) {
689  wR_ = MVT::Clone( y, 1 );
690  }
691  }
692  //
693  // Apply polynomial to x.
694  //
695  int n = MVT::GetNumberVecs( x );
696  std::vector<int> idxi(1), idxi2, idxj(1);
697 
698  // Select vector x[j].
699  for (int j=0; j<n; ++j) {
700 
701  idxi[0] = 0;
702  idxj[0] = j;
703  Teuchos::RCP<const MV> x_view = MVT::CloneView( x, idxj );
704  Teuchos::RCP<MV> y_view = MVT::CloneViewNonConst( y, idxj );
705  if (!LP_.is_null()) {
706  Teuchos::RCP<MV> v_curr = MVT::CloneViewNonConst( *V_, idxi );
707  problem_->applyLeftPrec( *x_view, *v_curr ); // Left precondition x into the first vector of V
708  } else {
709  MVT::SetBlock( *x_view, idxi, *V_ ); // Set x as the first vector of V
710  }
711 
712  for (int i=0; i<dim_-1; ++i) {
713 
714  // Get views into the current and next vectors
715  idxi2.resize(i+1);
716  for (int ii=0; ii<i+1; ++ii) { idxi2[ii] = ii; }
717  Teuchos::RCP<const MV> v_prev = MVT::CloneView( *V_, idxi2 );
718  // the tricks below with wR_ and wL_ (potentially set to v_curr and v_next) unfortunately imply that
719  // v_curr and v_next must be non-const views.
720  Teuchos::RCP<MV> v_curr = MVT::CloneViewNonConst( *V_, idxi );
721  idxi[0] = i+1;
722  Teuchos::RCP<MV> v_next = MVT::CloneViewNonConst( *V_, idxi );
723 
724  //---------------------------------------------
725  // Apply operator to next vector
726  //---------------------------------------------
727  // 1) Apply right preconditioner, if we have one.
728  if (!RP_.is_null()) {
729  problem_->applyRightPrec( *v_curr, *wR_ );
730  } else {
731  wR_ = v_curr;
732  }
733  // 2) Check for left preconditioner, if none exists, point at the next vector.
734  if (LP_.is_null()) {
735  wL_ = v_next;
736  }
737  // 3) Apply operator A.
738  problem_->applyOp( *wR_, *wL_ );
739  // 4) Apply left preconditioner, if we have one.
740  if (!LP_.is_null()) {
741  problem_->applyLeftPrec( *wL_, *v_next );
742  }
743 
744  // Compute A*v_curr - v_prev*H(1:i,i)
746  MVT::MvTimesMatAddMv( -SCT::one(), *v_prev, h, SCT::one(), *v_next );
747 
748  // Scale by H(i+1,i)
749  MVT::MvScale( *v_next, SCT::one()/H_(i+1,i) );
750  }
751 
752  // Compute output y = V*y_./r0_
753  if (!RP_.is_null()) {
754  MVT::MvTimesMatAddMv( SCT::one()/r0_(0), *V_, y_, SCT::zero(), *wR_ );
755  problem_->applyRightPrec( *wR_, *y_view );
756  }
757  else {
758  MVT::MvTimesMatAddMv( SCT::one()/r0_(0), *V_, y_, SCT::zero(), *y_view );
759  }
760  } // (int j=0; j<n; ++j)
761  } // end Apply()
762 } // end Belos namespace
763 
764 #endif
765 
766 // end of file BelosGmresPolyOp.hpp
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
Definition: BelosTypes.cpp:106
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.
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)
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:304
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.
Iterated Modified Gram-Schmidt (IMGS) implementation of the Belos::OrthoManager class.
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:81
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...
Classical Gram-Schmidt (with DGKS correction) implementation of the Belos::OrthoManager class...
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.
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.
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.
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.
Iterated Classical Gram-Schmidt (ICGS) implementation of the Belos::OrthoManager class.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
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.
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).
Alternative run-time polymorphic interface for operators.
virtual ~GmresPolyOp()
Destructor.
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
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
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)
Create a simple polynomial of dimension 1.
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
int shape(OrdinalType numRows, OrdinalType numCols)
Default parameters common to most Belos solvers.
Definition: BelosTypes.hpp:291
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.
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.

Generated on Fri Jun 5 2020 10:20:40 for Belos by doxygen 1.8.5