Belos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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"
57
61
67
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
83namespace 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 }
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);
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 static constexpr bool damp_default_ = false;
306 static constexpr bool addRoots_default_ = true;
307
308 // Variables for generating the polynomial
309 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
310 Teuchos::RCP<Teuchos::ParameterList> params_;
311 Teuchos::RCP<const OP> LP_, RP_;
312
313 // Output manager.
314 Teuchos::RCP<OutputManager<ScalarType> > printer_;
315 Teuchos::RCP<std::ostream> outputStream_ = Teuchos::rcpFromRef(std::cout);
316
317 // Orthogonalization manager.
318 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
319
320 // Current polynomial parameters
325 std::string label_ = label_default_;
328 int dim_ = 0;
331
332 // Variables for Arnoldi polynomial
333 mutable Teuchos::RCP<MV> V_, wL_, wR_;
334 Teuchos::SerialDenseMatrix<OT,ScalarType> H_, y_;
335 Teuchos::SerialDenseVector<OT,ScalarType> r0_;
336
337 // Variables for Gmres polynomial;
338 bool autoDeg = false;
339 Teuchos::SerialDenseMatrix< OT, ScalarType > pCoeff_;
340
341 // Variables for Roots polynomial:
342 Teuchos::SerialDenseMatrix< OT, MagnitudeType > theta_;
343
344 // Modified Leja sorting function. Takes a serial dense matrix of M harmonic Ritz values and an index
345 // of values from 0 to M. Returns the sorted values and sorted index, similar to Matlab.
346 void SortModLeja(Teuchos::SerialDenseMatrix< OT, MagnitudeType > &thetaN, std::vector<int> &index) const ;
347
348 //Function determines whether added roots are needed and adds them if option is turned on.
349 void ComputeAddedRoots();
350 };
351
352 template <class ScalarType, class MV, class OP>
353 void GmresPolyOp<ScalarType, MV, OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList>& params_in )
354 {
355 // Check which Gmres polynomial to use
356 if (params_in->isParameter("Polynomial Type")) {
357 polyType_ = params_in->get("Polynomial Type", polyType_default_);
358 }
359
360 // Check for polynomial convergence tolerance
361 if (params_in->isParameter("Polynomial Tolerance")) {
362 if (params_in->isType<MagnitudeType> ("Polynomial Tolerance")) {
363 polyTol_ = params_in->get ("Polynomial Tolerance",
365 }
366 else {
367 polyTol_ = params_in->get ("Polynomial Tolerance", DefaultSolverParameters::polyTol);
368 }
369 }
370
371 // Check for maximum polynomial degree
372 if (params_in->isParameter("Maximum Degree")) {
373 maxDegree_ = params_in->get("Maximum Degree", maxDegree_default_);
374 }
375
376 // Check for maximum polynomial degree
377 if (params_in->isParameter("Random RHS")) {
378 randomRHS_ = params_in->get("Random RHS", randomRHS_default_);
379 }
380
381 // Check for a change in verbosity level
382 if (params_in->isParameter("Verbosity")) {
383 if (Teuchos::isParameterType<int>(*params_in,"Verbosity")) {
384 verbosity_ = params_in->get("Verbosity", verbosity_default_);
385 }
386 else {
387 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params_in,"Verbosity");
388 }
389 }
390
391 if (params_in->isParameter("Orthogonalization")) {
392 orthoType_ = params_in->get("Orthogonalization",orthoType_default_);
393 }
394
395 // Check for timer label
396 if (params_in->isParameter("Timer Label")) {
397 label_ = params_in->get("Timer Label", label_default_);
398 }
399
400 // Output stream
401 if (params_in->isParameter("Output Stream")) {
402 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params_in,"Output Stream");
403 }
404
405 // Check for damped polynomial
406 if (params_in->isParameter("Damped Poly")) {
407 damp_ = params_in->get("Damped Poly", damp_default_);
408 }
409
410 // Check for root-adding
411 if (params_in->isParameter("Add Roots")) {
412 addRoots_ = params_in->get("Add Roots", addRoots_default_);
413 }
414 }
415
416 template <class ScalarType, class MV, class OP>
418 {
419 Teuchos::RCP< MV > V = MVT::Clone( *problem_->getRHS(), maxDegree_+1 );
420
421 //Make power basis:
422 std::vector<int> index(1,0);
423 Teuchos::RCP< MV > V0 = MVT::CloneViewNonConst(*V, index);
424 if (randomRHS_)
425 MVT::MvRandom( *V0 );
426 else
427 MVT::Assign( *problem_->getRHS(), *V0 );
428
429 if ( !LP_.is_null() ) {
430 Teuchos::RCP< MV > Vtemp = MVT::CloneCopy(*V0);
431 problem_->applyLeftPrec( *Vtemp, *V0);
432 }
433 if ( damp_ ) {
434 Teuchos::RCP< MV > Vtemp = MVT::CloneCopy(*V0);
435 problem_->apply( *Vtemp, *V0);
436 }
437
438 for(int i=0; i< maxDegree_; i++)
439 {
440 index[0] = i;
441 Teuchos::RCP< const MV > Vi = MVT::CloneView(*V, index);
442 index[0] = i+1;
443 Teuchos::RCP< MV > Vip1 = MVT::CloneViewNonConst(*V, index);
444 problem_->apply( *Vi, *Vip1);
445 }
446
447 //Consider AV:
448 Teuchos::Range1D range( 1, maxDegree_);
449 Teuchos::RCP< const MV > AV = MVT::CloneView( *V, range);
450
451 //Make lhs (AV)^T(AV)
452 Teuchos::SerialDenseMatrix< OT, ScalarType > AVtransAV( maxDegree_, maxDegree_);
453 MVT::MvTransMv( SCT::one(), *AV, *AV, AVtransAV);
454 //This process adds pDeg*pDeg + pDeg inner products that aren't in the final count.
455
456 Teuchos::LAPACK< OT, ScalarType > lapack;
457 int infoInt;
458 bool status = true; //Keep adjusting poly deg when true.
459
460 dim_ = maxDegree_;
461 Teuchos::SerialDenseMatrix< OT, ScalarType > lhs;
462 while( status && dim_ >= 1)
463 {
464 Teuchos::SerialDenseMatrix< OT, ScalarType > lhstemp(Teuchos::Copy, AVtransAV, dim_, dim_);
465 lapack.POTRF( 'U', dim_, lhstemp.values(), lhstemp.stride(), &infoInt);
466
467 if(autoDeg == false)
468 {
469 status = false;
470 if(infoInt != 0)
471 {
472 std::cout << "BelosGmresPolyOp.hpp: LAPACK POTRF was not successful!!" << std::endl;
473 std::cout << "Error code: " << infoInt << std::endl;
474 }
475 }
476 else
477 {
478 if(infoInt != 0)
479 {//Had bad factor. Reduce poly degree.
480 dim_--;
481 }
482 else
483 {
484 status = false;
485 }
486 }
487 if(status == false)
488 {
489 lhs = lhstemp;
490 }
491 }
492 if(dim_ == 0)
493 {
494 pCoeff_.shape( 1, 1);
495 pCoeff_(0,0) = SCT::one();
496 std::cout << "Poly Degree is zero. No preconditioner created." << std::endl;
497 }
498 else
499 {
500 pCoeff_.shape( dim_, 1);
501 //Get correct submatrix of AV:
502 Teuchos::Range1D rangeSub( 1, dim_);
503 Teuchos::RCP< const MV > AVsub = MVT::CloneView( *V, rangeSub);
504
505 //Compute rhs (AV)^T V0
506 MVT::MvTransMv( SCT::one(), *AVsub, *V0, pCoeff_);
507 lapack.POTRS( 'U', dim_, 1, lhs.values(), lhs.stride(), pCoeff_.values(), pCoeff_.stride(), &infoInt);
508 if(infoInt != 0)
509 {
510 std::cout << "BelosGmresPolyOp.hpp: LAPACK POTRS was not successful!!" << std::endl;
511 std::cout << "Error code: " << infoInt << std::endl;
512 }
513 }
514 }
515
516 template <class ScalarType, class MV, class OP>
518 {
519 std::string polyLabel = label_ + ": GmresPolyOp creation";
520
521 // Create a copy of the linear problem that has a zero initial guess and random RHS.
522 std::vector<int> idx(1,0);
523 Teuchos::RCP<MV> newX = MVT::Clone( *(problem_->getLHS()), 1 );
524 Teuchos::RCP<MV> newB = MVT::Clone( *(problem_->getRHS()), 1 );
525 MVT::MvInit( *newX, SCT::zero() );
526 if (randomRHS_) {
527 MVT::MvRandom( *newB );
528 }
529 else {
530 MVT::Assign( *(MVT::CloneView(*(problem_->getRHS()), idx)), *newB );
531 }
532 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > newProblem =
533 Teuchos::rcp( new LinearProblem<ScalarType,MV,OP>( problem_->getOperator(), newX, newB ) );
534 newProblem->setInitResVec( newB );
535 newProblem->setLeftPrec( problem_->getLeftPrec() );
536 newProblem->setRightPrec( problem_->getRightPrec() );
537 newProblem->setLabel(polyLabel);
538 newProblem->setProblem();
539 newProblem->setLSIndex( idx );
540
541 // Create a parameter list for the GMRES iteration.
542 Teuchos::ParameterList polyList;
543
544 // Tell the block solver that the block size is one.
545 polyList.set("Num Blocks",maxDegree_);
546 polyList.set("Block Size",1);
547 polyList.set("Keep Hessenberg", true);
548
549 // Create output manager.
550 printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
551
552 // Create orthogonalization manager if we need to.
553 if (ortho_.is_null()) {
554 params_->set("Orthogonalization", orthoType_);
556 Teuchos::RCP<Teuchos::ParameterList> paramsOrtho; // can be null
557
558 ortho_ = factory.makeMatOrthoManager (orthoType_, Teuchos::null, printer_, polyLabel, paramsOrtho);
559 }
560
561 // Create a simple status test that either reaches the relative residual tolerance or maximum polynomial size.
562 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxItrTst =
563 Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxDegree_ ) );
564
565 // Implicit residual test, using the native residual to determine if convergence was achieved.
566 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTst =
567 Teuchos::rcp( new StatusTestGenResNorm<ScalarType,MV,OP>( polyTol_ ) );
568 convTst->defineScaleForm( convertStringToScaleType("Norm of RHS"), Belos::TwoNorm );
569
570 // Convergence test that stops the iteration when either are satisfied.
571 Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > polyTest =
572 Teuchos::rcp( new StatusTestCombo<ScalarType,MV,OP>( StatusTestCombo<ScalarType,MV,OP>::OR, maxItrTst, convTst ) );
573
574 // Create Gmres iteration object to perform one cycle of Gmres.
575 Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > gmres_iter;
576 gmres_iter = Teuchos::rcp( new BlockGmresIter<ScalarType,MV,OP>(newProblem,printer_,polyTest,ortho_,polyList) );
577
578 // Create the first block in the current Krylov basis (residual).
579 Teuchos::RCP<MV> V_0 = MVT::CloneCopy( *newB );
580 if ( !LP_.is_null() )
581 newProblem->applyLeftPrec( *newB, *V_0 );
582 if ( damp_ )
583 {
584 Teuchos::RCP< MV > Vtemp = MVT::CloneCopy(*V_0);
585 newProblem->apply( *Vtemp, *V_0 );
586 }
587
588 // Get a matrix to hold the orthonormalization coefficients.
589 r0_.resize(1);
590
591 // Orthonormalize the new V_0
592 int rank = ortho_->normalize( *V_0, Teuchos::rcpFromRef(r0_) );
593 TEUCHOS_TEST_FOR_EXCEPTION(rank != 1,GmresPolyOpOrthoFailure,
594 "Belos::GmresPolyOp::generateArnoldiPoly(): Failed to compute initial block of orthonormal vectors for polynomial generation.");
595
596 // Set the new state and initialize the solver.
598 newstate.V = V_0;
599 newstate.z = Teuchos::rcpFromRef( r0_);
600 newstate.curDim = 0;
601 gmres_iter->initializeGmres(newstate);
602
603 // Perform Gmres iteration
604 try {
605 gmres_iter->iterate();
606 }
607 catch (GmresIterationOrthoFailure& e) {
608 // Try to recover the most recent least-squares solution
609 gmres_iter->updateLSQR( gmres_iter->getCurSubspaceDim() );
610 }
611 catch (std::exception& e) {
612 using std::endl;
613 printer_->stream(Errors) << "Error! Caught exception in BlockGmresIter::iterate() at iteration "
614 << gmres_iter->getNumIters() << endl << e.what () << endl;
615 throw;
616 }
617
618 // Get the solution for this polynomial, use in comparison below
619 Teuchos::RCP<MV> currX = gmres_iter->getCurrentUpdate();
620
621 // Record polynomial info, get current GMRES state
622 GmresIterationState<ScalarType,MV> gmresState = gmres_iter->getState();
623
624 // If the polynomial has no dimension, the tolerance is too low, return false
625 dim_ = gmresState.curDim;
626 if (dim_ == 0) {
627 return;
628 }
629 if(polyType_ == "Arnoldi"){
630 // Make a view and then copy the RHS of the least squares problem.
631 //
632 y_ = Teuchos::SerialDenseMatrix<OT,ScalarType>( Teuchos::Copy, *gmresState.z, dim_, 1 );
633 H_ = *gmresState.H;
634
635 //
636 // Solve the least squares problem.
637 //
638 Teuchos::BLAS<OT,ScalarType> blas;
639 blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
640 Teuchos::NON_UNIT_DIAG, dim_, 1, SCT::one(),
641 gmresState.R->values(), gmresState.R->stride(),
642 y_.values(), y_.stride() );
643 }
644 else{ //Generate Roots Poly
645 //Find Harmonic Ritz Values to use as polynomial roots:
646
647 //Copy of square H used to find poly roots:
648 H_ = Teuchos::SerialDenseMatrix<OT,ScalarType>(Teuchos::Copy, *gmresState.H, dim_, dim_);
649 //Zero out below subdiagonal of H:
650 for(int i=0; i <= dim_-3; i++) {
651 for(int k=i+2; k <= dim_-1; k++) {
652 H_(k,i) = SCT::zero();
653 }
654 }
655 //Extra copy of H because equilibrate changes the matrix:
656 Teuchos::SerialDenseMatrix<OT,ScalarType> Htemp (Teuchos::Copy, H_, dim_, dim_);
657
658 //View the m+1,m element and last col of H:
659 ScalarType Hlast = (*gmresState.H)(dim_,dim_-1);
660 Teuchos::SerialDenseMatrix<OT,ScalarType> HlastCol (Teuchos::View, H_, dim_, 1, 0, dim_-1);
661
662 //Set up linear system for H^{-*}e_m:
663 Teuchos::SerialDenseMatrix< OT, ScalarType > F(dim_,1), E(dim_,1);
664 E.putScalar(SCT::zero());
665 E(dim_-1,0) = SCT::one();
666
667 Teuchos::SerialDenseSolver< OT, ScalarType > HSolver;
668 HSolver.setMatrix( Teuchos::rcpFromRef(Htemp));
669 HSolver.solveWithTransposeFlag( Teuchos::CONJ_TRANS );
670 HSolver.setVectors( Teuchos::rcpFromRef(F), Teuchos::rcpFromRef(E));
671 HSolver.factorWithEquilibration( true );
672
673 //Factor matrix and solve for F = H^{-*}e_m:
674 int info = 0;
675 info = HSolver.factor();
676 if(info != 0){
677 std::cout << "Hsolver factor: info = " << info << std::endl;
678 }
679 info = HSolver.solve();
680 if(info != 0){
681 std::cout << "Hsolver solve : info = " << info << std::endl;
682 }
683
684 //Scale F and adjust H for Harmonic Ritz value eigenproblem:
685 F.scale(Hlast*Hlast);
686 HlastCol += F;
687
688 //Set up for eigenvalue problem to get Harmonic Ritz Values:
689 Teuchos::LAPACK< OT, ScalarType > lapack;
690 theta_.shape(dim_,2);//1st col for real part, 2nd col for imaginary
691
692 const int ldv = 1;
693 ScalarType* vlr = 0;
694
695 // Size of workspace and workspace for DGEEV
696 int lwork = -1;
697 std::vector<ScalarType> work(1);
698 std::vector<MagnitudeType> rwork(2*dim_);
699
700 //Find workspace size for DGEEV:
701 lapack.GEEV('N','N',dim_,H_.values(),H_.stride(),theta_[0],theta_[1],vlr, ldv, vlr, ldv, &work[0], lwork, &rwork[0], &info);
702 lwork = std::abs (static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work[0])));
703 work.resize( lwork );
704 // Solve for Harmonic Ritz Values:
705 lapack.GEEV('N','N',dim_,H_.values(),H_.stride(),theta_[0],theta_[1],vlr, ldv, vlr, ldv, &work[0], lwork, &rwork[0], &info);
706
707 if(info != 0){
708 std::cout << "GEEV solve : info = " << info << std::endl;
709 }
710
711 // Set index for sort function, verify roots are non-zero,
712 // and sort Harmonic Ritz Values:
713 const MagnitudeType tol = 10.0 * Teuchos::ScalarTraits<MagnitudeType>::eps();
714 std::vector<int> index(dim_);
715 for(int i=0; i<dim_; ++i){
716 index[i] = i;
717 // Check if real + imag parts of roots < tol.
718 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.");
719 }
720 SortModLeja(theta_,index);
721
722 //Add roots if neded.
723 ComputeAddedRoots();
724
725 }
726 }
727
728 //Function determines whether added roots are needed and adds them if option is turned on.
729 template <class ScalarType, class MV, class OP>
731 {
732 // Store theta (with cols for real and imag parts of Harmonic Ritz Vals)
733 // as one vector of complex numbers to perform arithmetic:
734 std::vector<std::complex<MagnitudeType>> cmplxHRitz (dim_);
735 for(unsigned int i=0; i<cmplxHRitz.size(); ++i){
736 cmplxHRitz[i] = std::complex<MagnitudeType>( theta_(i,0), theta_(i,1) );
737 }
738
739 // Compute product of factors (pof) to determine added roots:
740 const MagnitudeType one(1.0);
741 std::vector<MagnitudeType> pof (dim_,one);
742 for(int j=0; j<dim_; ++j) {
743 for(int i=0; i<dim_; ++i) {
744 if(i!=j) {
745 pof[j] = std::abs(pof[j]*(one-(cmplxHRitz[j]/cmplxHRitz[i])));
746 }
747 }
748 }
749
750 // Compute number of extra roots needed:
751 std::vector<int> extra (dim_);
752 int totalExtra = 0;
753 for(int i=0; i<dim_; ++i){
754 if (pof[i] > MCT::zero())
755 extra[i] = ceil((log10(pof[i])-MagnitudeType(4.0))/MagnitudeType(14.0));
756 else
757 extra[i] = 0;
758 if(extra[i] > 0){
759 totalExtra += extra[i];
760 }
761 }
762 if (totalExtra){
763 printer_->stream(Warnings) << "Warning: Need to add " << totalExtra << " extra roots." << std::endl;}
764
765 // If requested to add roots, append them to the theta matrix:
766 if(addRoots_ && totalExtra>0)
767 {
768 theta_.reshape(dim_+totalExtra,2);
769 // Make a matrix copy for perturbed roots:
770 Teuchos::SerialDenseMatrix<OT,MagnitudeType> thetaPert (Teuchos::Copy, theta_, dim_+totalExtra, 2);
771
772 //Add extra eigenvalues to matrix and perturb for sort:
773 int count = dim_;
774 for(int i=0; i<dim_; ++i){
775 for(int j=0; j< extra[i]; ++j){
776 theta_(count,0) = theta_(i,0);
777 theta_(count,1) = theta_(i,1);
778 thetaPert(count,0) = theta_(i,0)+(j+MCT::one())*MagnitudeType(5e-8);
779 thetaPert(count,1) = theta_(i,1);
780 ++count;
781 }
782 }
783
784 // Update polynomial degree:
785 dim_ += totalExtra;
786 if (totalExtra){
787 printer_->stream(Warnings) << "New poly degree is: " << dim_ << std::endl;}
788
789 // Create a new index and sort perturbed roots:
790 std::vector<int> index2(dim_);
791 for(int i=0; i<dim_; ++i){
792 index2[i] = i;
793 }
794 SortModLeja(thetaPert,index2);
795 //Apply sorting to non-perturbed roots:
796 for(int i=0; i<dim_; ++i)
797 {
798 thetaPert(i,0) = theta_(index2[i],0);
799 thetaPert(i,1) = theta_(index2[i],1);
800 }
801 theta_ = thetaPert;
802
803 }
804 }
805
806 // Modified Leja sorting function. Takes a serial dense matrix of M harmonic Ritz values and an index
807 // of values from 0 to M. Returns the sorted values and sorted index, similar to Matlab.
808 template <class ScalarType, class MV, class OP>
809 void GmresPolyOp<ScalarType, MV, OP>::SortModLeja(Teuchos::SerialDenseMatrix< OT, MagnitudeType > &thetaN, std::vector<int> &index) const
810 {
811 //Sort theta values via Modified Leja Ordering:
812
813 // Set up blank matrices to track sorting:
814 int dimN = index.size();
815 std::vector<int> newIndex(dimN);
816 Teuchos::SerialDenseMatrix< OT, MagnitudeType > sorted (thetaN.numRows(), thetaN.numCols());
817 Teuchos::SerialDenseVector< OT, MagnitudeType > absVal (thetaN.numRows());
818 Teuchos::SerialDenseVector< OT, MagnitudeType > prod (thetaN.numRows());
819
820 //Compute all absolute values and find maximum:
821 for(int i = 0; i < dimN; i++){
822 absVal(i) = hypot(thetaN(i,0), thetaN(i,1));
823 }
824 MagnitudeType * maxPointer = std::max_element(absVal.values(), (absVal.values()+dimN));
825 int maxIndex = int (maxPointer- absVal.values());
826
827 //Put largest abs value first in the list:
828 sorted(0,0) = thetaN(maxIndex,0);
829 sorted(0,1) = thetaN(maxIndex,1);
830 newIndex[0] = index[maxIndex];
831
832 int j;
833 // If largest value was complex (for real scalar type) put its conjugate in the next slot.
834 if(sorted(0,1)!= SCT::zero() && !SCT::isComplex)
835 {
836 sorted(1,0) = thetaN(maxIndex,0);
837 sorted(1,1) = -thetaN(maxIndex,1);
838 newIndex[1] = index[maxIndex+1];
839 j = 2;
840 }
841 else
842 {
843 j = 1;
844 }
845
846 //Sort remaining values:
847 MagnitudeType a, b;
848 while( j < dimN )
849 {
850 //For each value, compute (a log of) a product of differences:
851 for(int i = 0; i < dimN; i++)
852 {
853 prod(i) = MCT::one();
854 for(int k = 0; k < j; k++)
855 {
856 a = thetaN(i,0) - sorted(k,0);
857 b = thetaN(i,1) - sorted(k,1);
858 if (a*a + b*b > MCT::zero())
859 prod(i) = prod(i) + log10(hypot(a,b));
860 else {
861 prod(i) = -std::numeric_limits<MagnitudeType>::infinity();
862 break;
863 }
864 }
865 }
866
867 //Value with largest product goes in the next slot:
868 maxPointer = std::max_element(prod.values(), (prod.values()+dimN));
869 maxIndex = int (maxPointer- prod.values());
870 sorted(j,0) = thetaN(maxIndex,0);
871 sorted(j,1) = thetaN(maxIndex,1);
872 newIndex[j] = index[maxIndex];
873
874 //If it was complex (and scalar type real) put its conjugate in next slot:
875 if(sorted(j,1)!= SCT::zero() && !SCT::isComplex)
876 {
877 j++;
878 sorted(j,0) = thetaN(maxIndex,0);
879 sorted(j,1) = -thetaN(maxIndex,1);
880 newIndex[j] = index[maxIndex+1];
881 }
882 j++;
883 }
884
885 //Return sorted values and sorted indices:
886 thetaN = sorted;
887 index = newIndex;
888 } //End Modified Leja ordering
889
890 template <class ScalarType, class MV, class OP>
891 void GmresPolyOp<ScalarType, MV, OP>::ApplyPoly( const MV& x, MV& y ) const
892 {
893 if (dim_) {
894 if (polyType_ == "Arnoldi")
895 ApplyArnoldiPoly(x, y);
896 else if (polyType_ == "Gmres")
897 ApplyGmresPoly(x, y);
898 else if (polyType_ == "Roots")
899 ApplyRootsPoly(x, y);
900 }
901 else {
902 // Just apply the operator in problem_ to x and return y.
903 problem_->applyOp( x, y );
904 }
905 }
906
907 template <class ScalarType, class MV, class OP>
908 void GmresPolyOp<ScalarType, MV, OP>::ApplyGmresPoly( const MV& x, MV& y ) const
909 {
910 Teuchos::RCP<MV> AX = MVT::CloneCopy(x);
911 Teuchos::RCP<MV> AX2 = MVT::Clone( x, MVT::GetNumberVecs(x) );
912
913 // Apply left preconditioner.
914 if (!LP_.is_null()) {
915 Teuchos::RCP<MV> Xtmp = MVT::Clone( x, MVT::GetNumberVecs(x) );
916 problem_->applyLeftPrec( *AX, *Xtmp ); // Left precondition x into the first vector
917 AX = Xtmp;
918 }
919
920 {
921#ifdef BELOS_TEUCHOS_TIME_MONITOR
922 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
923#endif
924 MVT::MvAddMv(pCoeff_(0,0), *AX, SCT::zero(), y, y); //y= coeff_i(A^ix)
925 }
926 for( int i=1; i < dim_; i++)
927 {
928 Teuchos::RCP<MV> X, Y;
929 if ( i%2 )
930 {
931 X = AX;
932 Y = AX2;
933 }
934 else
935 {
936 X = AX2;
937 Y = AX;
938 }
939 problem_->apply(*X, *Y);
940 {
941#ifdef BELOS_TEUCHOS_TIME_MONITOR
942 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
943#endif
944 MVT::MvAddMv(pCoeff_(i,0), *Y, SCT::one(), y, y); //y= coeff_i(A^ix) +y
945 }
946 }
947
948 // Apply right preconditioner.
949 if (!RP_.is_null()) {
950 Teuchos::RCP<MV> Ytmp = MVT::CloneCopy(y);
951 problem_->applyRightPrec( *Ytmp, y );
952 }
953 }
954
955 template <class ScalarType, class MV, class OP>
956 void GmresPolyOp<ScalarType, MV, OP>::ApplyRootsPoly( const MV& x, MV& y ) const
957 {
958 MVT::MvInit( y, SCT::zero() ); //Zero out y to take the vector with poly applied.
959 Teuchos::RCP<MV> prod = MVT::CloneCopy(x);
960 Teuchos::RCP<MV> Xtmp = MVT::Clone( x, MVT::GetNumberVecs(x) );
961 Teuchos::RCP<MV> Xtmp2 = MVT::Clone( x, MVT::GetNumberVecs(x) );
962
963 // Apply left preconditioner.
964 if (!LP_.is_null()) {
965 problem_->applyLeftPrec( *prod, *Xtmp ); // Left precondition x into the first vector
966 prod = Xtmp;
967 }
968
969 int i=0;
970 while(i < dim_-1)
971 {
972 if(theta_(i,1)== SCT::zero() || SCT::isComplex) //Real Harmonic Ritz value or complex scalars
973 {
974 {
975#ifdef BELOS_TEUCHOS_TIME_MONITOR
976 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
977#endif
978 MVT::MvAddMv(SCT::one(), y, SCT::one()/theta_(i,0), *prod, y); //poly = poly + 1/theta_i * prod
979 }
980 problem_->apply(*prod, *Xtmp); // temp = A*prod
981 {
982#ifdef BELOS_TEUCHOS_TIME_MONITOR
983 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
984#endif
985 MVT::MvAddMv(SCT::one(), *prod, -SCT::one()/theta_(i,0), *Xtmp, *prod); //prod = prod - 1/theta_i * temp
986 }
987 i++;
988 }
989 else //Current theta is complex and has a conjugate; combine to preserve real arithmetic
990 {
991 MagnitudeType mod = theta_(i,0)*theta_(i,0) + theta_(i,1)*theta_(i,1); //mod = a^2 + b^2
992 problem_->apply(*prod, *Xtmp); // temp = A*prod
993 {
994#ifdef BELOS_TEUCHOS_TIME_MONITOR
995 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
996#endif
997 MVT::MvAddMv(2*theta_(i,0), *prod, -SCT::one(), *Xtmp, *Xtmp); //temp = 2a*prod-temp
998 MVT::MvAddMv(SCT::one(), y, SCT::one()/mod, *Xtmp, y); //poly = poly + 1/mod*temp
999 }
1000 if( i < dim_-2 )
1001 {
1002 problem_->apply(*Xtmp, *Xtmp2); // temp2 = A*temp
1003 {
1004#ifdef BELOS_TEUCHOS_TIME_MONITOR
1005 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
1006#endif
1007 MVT::MvAddMv(SCT::one(), *prod, -SCT::one()/mod, *Xtmp2, *prod); //prod = prod - 1/mod * temp2
1008 }
1009 }
1010 i = i + 2;
1011 }
1012 }
1013 if(theta_(dim_-1,1)== SCT::zero() || SCT::isComplex)
1014 {
1015#ifdef BELOS_TEUCHOS_TIME_MONITOR
1016 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
1017#endif
1018 MVT::MvAddMv(SCT::one(), y, SCT::one()/theta_(dim_-1,0), *prod, y); //poly = poly + 1/theta_i * prod
1019 }
1020
1021 // Apply right preconditioner.
1022 if (!RP_.is_null()) {
1023 Teuchos::RCP<MV> Ytmp = MVT::CloneCopy(y);
1024 problem_->applyRightPrec( *Ytmp, y );
1025 }
1026 }
1027
1028 template <class ScalarType, class MV, class OP>
1030 {
1031 // Initialize vector storage.
1032 if (V_.is_null()) {
1033 V_ = MVT::Clone( x, dim_ );
1034 if (!LP_.is_null()) {
1035 wL_ = MVT::Clone( y, 1 );
1036 }
1037 if (!RP_.is_null()) {
1038 wR_ = MVT::Clone( y, 1 );
1039 }
1040 }
1041 //
1042 // Apply polynomial to x.
1043 //
1044 int n = MVT::GetNumberVecs( x );
1045 std::vector<int> idxi(1), idxi2, idxj(1);
1046
1047 // Select vector x[j].
1048 for (int j=0; j<n; ++j) {
1049
1050 idxi[0] = 0;
1051 idxj[0] = j;
1052 Teuchos::RCP<const MV> x_view = MVT::CloneView( x, idxj );
1053 Teuchos::RCP<MV> y_view = MVT::CloneViewNonConst( y, idxj );
1054 if (!LP_.is_null()) {
1055 Teuchos::RCP<MV> v_curr = MVT::CloneViewNonConst( *V_, idxi );
1056 problem_->applyLeftPrec( *x_view, *v_curr ); // Left precondition x into the first vector of V
1057 } else {
1058 MVT::SetBlock( *x_view, idxi, *V_ ); // Set x as the first vector of V
1059 }
1060
1061 for (int i=0; i<dim_-1; ++i) {
1062
1063 // Get views into the current and next vectors
1064 idxi2.resize(i+1);
1065 for (int ii=0; ii<i+1; ++ii) { idxi2[ii] = ii; }
1066 Teuchos::RCP<const MV> v_prev = MVT::CloneView( *V_, idxi2 );
1067 // the tricks below with wR_ and wL_ (potentially set to v_curr and v_next) unfortunately imply that
1068 // v_curr and v_next must be non-const views.
1069 Teuchos::RCP<MV> v_curr = MVT::CloneViewNonConst( *V_, idxi );
1070 idxi[0] = i+1;
1071 Teuchos::RCP<MV> v_next = MVT::CloneViewNonConst( *V_, idxi );
1072
1073 //---------------------------------------------
1074 // Apply operator to next vector
1075 //---------------------------------------------
1076 // 1) Apply right preconditioner, if we have one.
1077 if (!RP_.is_null()) {
1078 problem_->applyRightPrec( *v_curr, *wR_ );
1079 } else {
1080 wR_ = v_curr;
1081 }
1082 // 2) Check for left preconditioner, if none exists, point at the next vector.
1083 if (LP_.is_null()) {
1084 wL_ = v_next;
1085 }
1086 // 3) Apply operator A.
1087 problem_->applyOp( *wR_, *wL_ );
1088 // 4) Apply left preconditioner, if we have one.
1089 if (!LP_.is_null()) {
1090 problem_->applyLeftPrec( *wL_, *v_next );
1091 }
1092
1093 // Compute A*v_curr - v_prev*H(1:i,i)
1094 Teuchos::SerialDenseMatrix<OT,ScalarType> h(Teuchos::View,H_,i+1,1,0,i);
1095 {
1096#ifdef BELOS_TEUCHOS_TIME_MONITOR
1097 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
1098#endif
1099 MVT::MvTimesMatAddMv( -SCT::one(), *v_prev, h, SCT::one(), *v_next );
1100 }
1101
1102 // Scale by H(i+1,i)
1103 MVT::MvScale( *v_next, SCT::one()/H_(i+1,i) );
1104 }
1105
1106 // Compute output y = V*y_./r0_
1107 if (!RP_.is_null()) {
1108 {
1109#ifdef BELOS_TEUCHOS_TIME_MONITOR
1110 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
1111#endif
1112 MVT::MvTimesMatAddMv( SCT::one()/r0_(0), *V_, y_, SCT::zero(), *wR_ );
1113 }
1114 problem_->applyRightPrec( *wR_, *y_view );
1115 }
1116 else {
1117#ifdef BELOS_TEUCHOS_TIME_MONITOR
1118 Teuchos::TimeMonitor updateTimer( *timerPolyUpdate_ );
1119#endif
1120 MVT::MvTimesMatAddMv( SCT::one()/r0_(0), *V_, y_, SCT::zero(), *y_view );
1121 }
1122 } // (int j=0; j<n; ++j)
1123 } // end Apply()
1124} // end Belos namespace
1125
1126#endif
1127
1128// end of file BelosGmresPolyOp.hpp
Belos concrete class for performing the block GMRES iteration.
Belos header file which uses auto-configuration information to include necessary C++ headers.
Pure virtual base class which augments the basic interface for a Gmres linear solver iteration.
Class which describes the linear problem to be solved by the iterative solver.
Declaration of basic traits for the multivector type.
Interface for multivectors used by Belos' linear solvers.
Class which defines basic traits for the operator type.
Alternative run-time polymorphic interface for operators.
Class which manages the output and verbosity of the Belos solvers.
Belos::StatusTest for logically combining several status tests.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Belos::StatusTest for specifying an implicit residual norm stopping criteria that checks for loss of ...
Belos::StatusTest class for specifying a maximum number of iterations.
A factory class for generating StatusTestOutput objects.
Collection of types and exceptions used within the Belos solvers.
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:60
This class implements the block GMRES iteration, where a block Krylov subspace is constructed....
GmresIterationOrthoFailure is thrown when the GmresIteration object is unable to compute independent ...
int GetNumberVecs() const
The number of vectors (i.e., columns) in the multivector.
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 MvScale(const ScalarType alpha)
Scale each element of the vectors in *this with alpha.
void MvAddMv(const ScalarType alpha, const MultiVec< ScalarType > &A, const ScalarType beta, const MultiVec< ScalarType > &B)
Replace *this with alpha * A + beta * B.
Teuchos::RCP< const MV > getConstMV() const
void MvPrint(std::ostream &os) const
Print *this multivector to the os output stream.
void MvScale(const std::vector< ScalarType > &alpha)
Scale each element of the i-th vector in *this with alpha[i].
void MvRandom()
Fill all the vectors in *this with random numbers.
void SetBlock(const MultiVec< ScalarType > &A, const std::vector< int > &index)
Copy the vectors in A to a set of vectors in *this.
GmresPolyMv(const Teuchos::RCP< const MV > &mv_in)
MultiVecTraits< ScalarType, MV > MVT
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...
ptrdiff_t GetGlobalLength() const
The number of rows in the multivector.
GmresPolyMv(const Teuchos::RCP< MV > &mv_in)
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 MvNorm(std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec, NormType type=TwoNorm) const
Compute the norm of each vector in *this.
Teuchos::RCP< MV > getMV()
Teuchos::RCP< MV > mv_
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).
GmresPolyMv * Clone(const int numvecs) const
Create a new MultiVec with numvecs columns.
GmresPolyMv * CloneCopy() const
Create a new MultiVec and copy contents of *this into it (deep 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).
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 MvInit(const ScalarType alpha)
Replace each element of the vectors in *this with alpha.
GmresPolyOpOrthoFailure is thrown when the orthogonalization manager is unable to generate orthonorma...
GmresPolyOpOrthoFailure(const std::string &what_arg)
Belos's class for applying the GMRES polynomial operator that is used by the hybrid-GMRES linear solv...
static constexpr bool addRoots_default_
Teuchos::SerialDenseMatrix< OT, MagnitudeType > theta_
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > &params_in)
Process the passed in parameters.
Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > ortho_
Teuchos::ScalarTraits< MagnitudeType > MCT
void ApplyRootsPoly(const MV &x, MV &y) const
static constexpr const char * polyType_default_
void generateGmresPoly()
This routine takes the matrix, preconditioner, and vectors from the linear problem as well as the par...
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
void SortModLeja(Teuchos::SerialDenseMatrix< OT, MagnitudeType > &thetaN, std::vector< int > &index) const
Teuchos::RCP< OutputManager< ScalarType > > printer_
GmresPolyOp(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem_in)
Given no ParameterList, constructor creates no polynomial and only applies the given operator.
MultiVecTraits< ScalarType, MV > MVT
static constexpr const char * label_default_
Teuchos::RCP< MV > V_
static constexpr bool damp_default_
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 ApplyGmresPoly(const MV &x, MV &y) const
Teuchos::RCP< Teuchos::ParameterList > params_
void generateArnoldiPoly()
This routine takes the matrix, preconditioner, and vectors from the linear problem as well as the par...
static constexpr int maxDegree_default_
Teuchos::SerialDenseMatrix< OT, ScalarType > y_
Teuchos::ScalarTraits< ScalarType > SCT
void ApplyArnoldiPoly(const MV &x, MV &y) const
Teuchos::RCP< const OP > LP_
Teuchos::RCP< MV > wR_
static constexpr const char * orthoType_default_
std::string polyUpdateLabel_
static constexpr int verbosity_default_
Teuchos::RCP< MV > wL_
Teuchos::RCP< const OP > RP_
Teuchos::SerialDenseMatrix< OT, ScalarType > pCoeff_
static constexpr bool randomRHS_default_
virtual ~GmresPolyOp()
Destructor.
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...
GmresPolyOp(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem_in, const Teuchos::RCP< Teuchos::ParameterList > &params_in)
Basic contstructor.
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
Teuchos::SerialDenseVector< OT, ScalarType > r0_
Teuchos::RCP< std::ostream > outputStream_
Teuchos::SerialDenseMatrix< OT, ScalarType > H_
A linear system to solve, and its associated information.
Traits class which defines basic operations on multivectors.
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 .
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
static void MvScale(MV &mv, const ScalarType alpha)
Scale each element of the vectors in mv with alpha.
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
static void MvPrint(const MV &mv, std::ostream &os)
Print the mv multi-vector to the os output stream.
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 ,...
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 ...
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
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).
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).
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
static void MvInit(MV &mv, const ScalarType alpha=Teuchos::ScalarTraits< ScalarType >::zero())
Replace each element of the vectors in mv with alpha.
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.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
Interface for multivectors used by Belos' linear solvers.
Alternative run-time polymorphic interface for operators.
Enumeration of all valid Belos (Mat)OrthoManager classes.
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.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
A class for extending the status testing capabilities of Belos via logical combinations.
An implementation of StatusTestResNorm using a family of residual norms.
A Belos::StatusTest class for specifying a maximum number of iterations.
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
Definition: BelosTypes.cpp:106
NormType
The type of vector norm to compute.
Definition: BelosTypes.hpp:97
@ TwoNorm
Definition: BelosTypes.hpp:98
@ Warnings
Definition: BelosTypes.hpp:256
ETrans
Whether to apply the (conjugate) transpose of an operator.
Definition: BelosTypes.hpp:81
@ NOTRANS
Definition: BelosTypes.hpp:81
static const double polyTol
Relative residual tolerance for matrix polynomial construction.
Definition: BelosTypes.hpp:296
Structure to contain pointers to GmresIteration state variables.
Teuchos::RCP< const MV > V
The current Krylov basis.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > z
The current right-hand side of the least squares system RY = Z.
int curDim
The current dimension of the reduction.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > R
The current upper-triangular matrix from the QR reduction of H.