Belos  Version of the Day
BelosPCPGIter.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_PCPG_ITER_HPP
43 #define BELOS_PCPG_ITER_HPP
44 
49 #include "BelosConfigDefs.hpp"
50 #include "BelosTypes.hpp"
51 
52 #include "BelosLinearProblem.hpp"
53 #include "BelosMatOrthoManager.hpp"
54 #include "BelosOutputManager.hpp"
55 #include "BelosStatusTest.hpp"
56 #include "BelosOperatorTraits.hpp"
57 #include "BelosMultiVecTraits.hpp"
58 
59 #include "Teuchos_SerialDenseMatrix.hpp"
60 #include "Teuchos_SerialDenseVector.hpp"
61 #include "Teuchos_ScalarTraits.hpp"
62 #include "Teuchos_ParameterList.hpp"
63 #include "Teuchos_TimeMonitor.hpp"
64 
76 namespace Belos {
77 
79 
80 
85  template <class ScalarType, class MV>
86  struct PCPGIterState {
92  int curDim;
94  int prevUdim;
95 
97  Teuchos::RCP<MV> R;
98 
100  Teuchos::RCP<MV> Z;
101 
103  Teuchos::RCP<MV> P;
104 
106  Teuchos::RCP<MV> AP;
107 
109  Teuchos::RCP<MV> U;
110 
112  Teuchos::RCP<MV> C;
113 
118  Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > D;
119 
121  prevUdim(0),
122  R(Teuchos::null), Z(Teuchos::null),
123  P(Teuchos::null), AP(Teuchos::null),
124  U(Teuchos::null), C(Teuchos::null),
125  D(Teuchos::null)
126  {}
127  };
128 
130 
132 
133 
145  class PCPGIterInitFailure : public BelosError {public:
146  PCPGIterInitFailure(const std::string& what_arg) : BelosError(what_arg)
147  {}};
148 
153  class PCPGIterateFailure : public BelosError {public:
154  PCPGIterateFailure(const std::string& what_arg) : BelosError(what_arg)
155  {}};
156 
163  class PCPGIterOrthoFailure : public BelosError {public:
164  PCPGIterOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
165  {}};
166 
167  template<class ScalarType, class MV, class OP>
168  class PCPGIter : virtual public Iteration<ScalarType,MV,OP> {
169 
170  public:
171 
172  //
173  // Convenience typedefs
174  //
177  typedef Teuchos::ScalarTraits<ScalarType> SCT;
178  typedef typename SCT::magnitudeType MagnitudeType;
179 
181 
182 
190  PCPGIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
191  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
192  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
193  const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
194  Teuchos::ParameterList &params );
195 
197  virtual ~PCPGIter() {};
199 
200 
202 
203 
222  void iterate();
223 
243 
247  void initialize()
248  {
250  initialize(empty);
251  }
252 
262  state.Z = Z_; // CG state
263  state.P = P_;
264  state.AP = AP_;
265  state.R = R_;
266  state.U = U_; // seed state
267  state.C = C_;
268  state.D = D_;
269  state.curDim = curDim_;
270  state.prevUdim = prevUdim_;
271  return state;
272  }
273 
275 
276 
278 
279 
281  int getNumIters() const { return iter_; }
282 
284  void resetNumIters( int iter = 0 ) { iter_ = iter; }
285 
288  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> * /* norms */ ) const { return R_; }
289 
291 
294  Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
295 
297  int getCurSubspaceDim() const {
298  if (!initialized_) return 0;
299  return curDim_;
300  };
301 
303  int getPrevSubspaceDim() const {
304  if (!initialized_) return 0;
305  return prevUdim_;
306  };
307 
309 
310 
312 
313 
315  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
316 
318  int getBlockSize() const { return 1; }
319 
321  int getNumRecycledBlocks() const { return savedBlocks_; }
322 
324 
326  void setBlockSize(int blockSize) {
327  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
328  "Belos::PCPGIter::setBlockSize(): Cannot use a block size that is not one.");
329  }
330 
332  void setSize( int savedBlocks );
333 
335  bool isInitialized() { return initialized_; }
336 
338  void resetState();
339 
341 
342  private:
343 
344  //
345  // Internal methods
346  //
348  void setStateSize();
349 
350  //
351  // Classes inputed through constructor that define the linear problem to be solved.
352  //
353  const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
354  const Teuchos::RCP<OutputManager<ScalarType> > om_;
355  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
356  const Teuchos::RCP<OrthoManager<ScalarType,MV> > ortho_;
357 
358  //
359  // Algorithmic parameters
360  // savedBlocks_ is the number of blocks allocated for the reused subspace
361  int savedBlocks_;
362  //
363  //
364  // Current solver state
365  //
366  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
367  // is capable of running; _initialize is controlled by the initialize() member method
368  // For the implications of the state of initialized_, please see documentation for initialize()
369  bool initialized_;
370 
371  // stateStorageInitialized_ indicates that the state storage has be initialized to the current
372  // savedBlocks_. State storage initialization may be postponed if the linear problem was
373  // generated without either the right-hand side or solution vectors.
374  bool stateStorageInitialized_;
375 
376  // keepDiagonal_ specifies that the iteration must keep the diagonal matrix of pivots
377  bool keepDiagonal_;
378 
379  // initDiagonal_ specifies that the iteration will reinitialize the diagonal matrix by zeroing
380  // out all entries before an iteration is started.
381  bool initDiagonal_;
382 
383  // Current subspace dimension
384  int curDim_;
385 
386  // Dimension of seed space used to solve current linear system
387  int prevUdim_;
388 
389  // Number of iterations performed
390  int iter_;
391  //
392  // State Storage ... of course this part is different for CG
393  //
394  // Residual
395  Teuchos::RCP<MV> R_;
396  //
397  // Preconditioned residual
398  Teuchos::RCP<MV> Z_;
399  //
400  // Direction std::vector
401  Teuchos::RCP<MV> P_;
402  //
403  // Operator applied to direction std::vector
404  Teuchos::RCP<MV> AP_;
405  //
406  // Recycled subspace vectors.
407  Teuchos::RCP<MV> U_;
408  //
409  // C = A * U, linear system is Ax=b
410  Teuchos::RCP<MV> C_;
411  //
412  // Projected matrices
413  // D_ : Diagonal matrix of pivots D = P'AP
414  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > D_;
415  };
416 
418  // Constructor.
419  template<class ScalarType, class MV, class OP>
421  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
422  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
423  const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
424  Teuchos::ParameterList &params ):
425  lp_(problem),
426  om_(printer),
427  stest_(tester),
428  ortho_(ortho),
429  savedBlocks_(0),
430  initialized_(false),
431  stateStorageInitialized_(false),
432  keepDiagonal_(false),
433  initDiagonal_(false),
434  curDim_(0),
435  prevUdim_(0),
436  iter_(0)
437  {
438  // Get the maximum number of blocks allowed for this Krylov subspace
439 
440  TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Saved Blocks"), std::invalid_argument,
441  "Belos::PCPGIter::constructor: mandatory parameter \"Saved Blocks\" is not specified.");
442  int rb = Teuchos::getParameter<int>(params, "Saved Blocks");
443 
444  // Find out whether we are saving the Diagonal matrix.
445  keepDiagonal_ = params.get("Keep Diagonal", false);
446 
447  // Find out whether we are initializing the Diagonal matrix.
448  initDiagonal_ = params.get("Initialize Diagonal", false);
449 
450  // Set the number of blocks and allocate data
451  setSize( rb );
452  }
453 
455  // Set the block size and adjust as necessary
456  template <class ScalarType, class MV, class OP>
457  void PCPGIter<ScalarType,MV,OP>::setSize( int savedBlocks )
458  {
459  // allocate space only; perform no computation
460  // Any change in size invalidates the state of the solver as implemented here.
461 
462  TEUCHOS_TEST_FOR_EXCEPTION(savedBlocks <= 0, std::invalid_argument, "Belos::PCPGIter::setSize() was passed a non-positive argument for \"Num Saved Blocks\".");
463 
464  if ( savedBlocks_ != savedBlocks) {
465  stateStorageInitialized_ = false;
466  savedBlocks_ = savedBlocks;
467  initialized_ = false;
468  curDim_ = 0;
469  prevUdim_ = 0;
470  setStateSize(); // Use the current savedBlocks_ to initialize the state storage.
471  }
472  }
473 
475  // Enable the reuse of a single solver object for completely different linear systems
476  template <class ScalarType, class MV, class OP>
478  {
479  stateStorageInitialized_ = false;
480  initialized_ = false;
481  curDim_ = 0;
482  prevUdim_ = 0;
483  setStateSize();
484  }
485 
487  // Setup the state storage. Called by either initialize or, if savedBlocks_ changes, setSize.
488  template <class ScalarType, class MV, class OP>
490  {
491  if (!stateStorageInitialized_) {
492 
493  // Check if there is any multivector to clone from.
494  Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
495  Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
496  if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
497  return; // postpone exception
498  }
499  else {
500 
502  // blockSize*recycledBlocks dependent
503  int newsd = savedBlocks_ ; //int newsd = blockSize_* savedBlocks_ ;
504  //
505  // Initialize the CG state storage
506  // If the subspace is not initialized, generate it using the LHS or RHS from lp_.
507  // Generate CG state only if it does not exist, otherwise resize it.
508  if (Z_ == Teuchos::null) {
509  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
510  Z_ = MVT::Clone( *tmp, 1 );
511  }
512  if (P_ == Teuchos::null) {
513  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
514  P_ = MVT::Clone( *tmp, 1 );
515  }
516  if (AP_ == Teuchos::null) {
517  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
518  AP_ = MVT::Clone( *tmp, 1 );
519  }
520 
521  if (C_ == Teuchos::null) {
522 
523  // Get the multivector that is not null.
524  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
525  TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
526  "Belos::PCPGIter::setStateSize(): linear problem does not specify multivectors to clone from.");
527  TEUCHOS_TEST_FOR_EXCEPTION( 0 != prevUdim_,std::invalid_argument,
528  "Belos::PCPGIter::setStateSize(): prevUdim not zero and C is null.");
529  C_ = MVT::Clone( *tmp, savedBlocks_ );
530  }
531  else {
532  // Generate C_ by cloning itself ONLY if more space is needed.
533  if (MVT::GetNumberVecs(*C_) < savedBlocks_ ) {
534  Teuchos::RCP<const MV> tmp = C_;
535  C_ = MVT::Clone( *tmp, savedBlocks_ );
536  }
537  }
538  if (U_ == Teuchos::null) {
539  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
540  TEUCHOS_TEST_FOR_EXCEPTION( 0 != prevUdim_,std::invalid_argument,
541  "Belos::PCPGIter::setStateSize(): prevUdim not zero and U is null.");
542  U_ = MVT::Clone( *tmp, savedBlocks_ );
543  }
544  else {
545  // Generate U_ by cloning itself ONLY if more space is needed.
546  if (MVT::GetNumberVecs(*U_) < savedBlocks_ ) {
547  Teuchos::RCP<const MV> tmp = U_;
548  U_ = MVT::Clone( *tmp, savedBlocks_ );
549  }
550  }
551  if (keepDiagonal_) {
552  if (D_ == Teuchos::null) {
553  D_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>() );
554  }
555  if (initDiagonal_) {
556  D_->shape( newsd, newsd );
557  }
558  else {
559  if (D_->numRows() < newsd || D_->numCols() < newsd) {
560  D_->shapeUninitialized( newsd, newsd );
561  }
562  }
563  }
564  // State storage has now been initialized.
565  stateStorageInitialized_ = true;
566  } // if there is a vector to clone from
567  } // if !stateStorageInitialized_
568  } // end of setStateSize
569 
571  // Initialize the iteration object
572  template <class ScalarType, class MV, class OP>
574  {
575 
576  TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
577  "Belos::PCPGIter::initialize(): Cannot initialize state storage!");
578 
579  // Requirements: R_ and consistent multivectors widths and lengths
580  //
581  std::string errstr("Belos::PCPGIter::initialize(): Specified multivectors must have a consistent length and width.");
582 
583  if (newstate.R != Teuchos::null){
584 
585  R_ = newstate.R; // SolverManager::R_ == newstate.R == Iterator::R_
586  if (newstate.U == Teuchos::null){
587  prevUdim_ = 0;
588  newstate.U = U_;
589  newstate.C = C_;
590  }
591  else {
592  prevUdim_ = newstate.curDim;
593  if (newstate.C == Teuchos::null){ // Stub for new feature
594  std::vector<int> index(prevUdim_);
595  for (int i=0; i< prevUdim_; ++i)
596  index[i] = i;
597  Teuchos::RCP<const MV> Ukeff = MVT::CloneView( *newstate.U, index );
598  newstate.C = MVT::Clone( *newstate.U, prevUdim_ );
599  Teuchos::RCP<MV> Ckeff = MVT::CloneViewNonConst( *newstate.C, index );
600  lp_->apply( *Ukeff, *Ckeff );
601  }
602  curDim_ = prevUdim_ ;
603  }
604 
605  // Initialize the state storage if not already allocated in the constructor
606  if (!stateStorageInitialized_)
607  setStateSize();
608 
609  //TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.V) != MVT::GetGlobalLength(*V_), std::invalid_argument, errstr );
610  //TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < 1, std::invalid_argument, errstr );
611 
612  newstate.prevUdim = prevUdim_; // big change in functionality from GCRODR
613  newstate.curDim = curDim_;
614 
615  //TEUCHOS_TEST_FOR_EXCEPTION(newstate.z->numRows() < curDim_ || newstate.z->numCols() < 1, std::invalid_argument, errstr);
616 
617  std::vector<int> zero_index(1);
618  zero_index[0] = 0;
619  if ( lp_->getLeftPrec() != Teuchos::null ) { // Compute the initial search direction
620  lp_->applyLeftPrec( *R_, *Z_ );
621  MVT::SetBlock( *Z_, zero_index , *P_ ); // P(:,zero_index) := Z
622  } else {
623  Z_ = R_;
624  MVT::SetBlock( *R_, zero_index, *P_ );
625  }
626 
627  std::vector<int> nextind(1);
628  nextind[0] = curDim_;
629 
630  MVT::SetBlock( *P_, nextind, *newstate.U ); // U(:,curDim_ ) := P_
631 
632  ++curDim_;
633  newstate.curDim = curDim_;
634 
635  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.U) != savedBlocks_ ,
636  std::invalid_argument, errstr );
637  if (newstate.U != U_) { // Why this is needed?
638  U_ = newstate.U;
639  }
640 
641  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.C) != savedBlocks_ ,
642  std::invalid_argument, errstr );
643  if (newstate.C != C_) {
644  C_ = newstate.C;
645  }
646  }
647  else {
648 
649  TEUCHOS_TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
650  "Belos::PCPGIter::initialize(): PCPGStateIterState does not have initial kernel R_0.");
651  }
652 
653  // the solver is initialized
654  initialized_ = true;
655  }
656 
657 
659  // Iterate until the status test informs us we should stop.
660  template <class ScalarType, class MV, class OP>
662  {
663  //
664  // Allocate/initialize data structures
665  //
666  if (initialized_ == false) {
667  initialize();
668  }
669  const bool debug = false;
670 
671  // Allocate memory for scalars.
672  Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
673  Teuchos::SerialDenseMatrix<int,ScalarType> pAp( 1, 1 );
674  Teuchos::SerialDenseMatrix<int,ScalarType> beta( 1, 1 );
675  Teuchos::SerialDenseMatrix<int,ScalarType> rHz( 1, 1 ), rHz_old( 1, 1 );
676 
677  if( iter_ != 0 )
678  std::cout << " Iterate Warning: begin from nonzero iter_ ?" << std::endl; //DMD
679 
680  // GenOrtho Project Stubs
681  std::vector<int> prevInd;
682  Teuchos::RCP<const MV> Uprev;
683  Teuchos::RCP<const MV> Cprev;
684  Teuchos::SerialDenseMatrix<int,ScalarType> CZ;
685 
686  if( prevUdim_ ){
687  prevInd.resize( prevUdim_ );
688  for( int i=0; i<prevUdim_ ; i++) prevInd[i] = i;
689  CZ.reshape( prevUdim_ , 1 );
690  Uprev = MVT::CloneView(*U_, prevInd);
691  Cprev = MVT::CloneView(*C_, prevInd);
692  }
693 
694  // Get the current solution std::vector.
695  Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
696 
697  // Check that the current solution std::vector only has one column.
698  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, PCPGIterInitFailure,
699  "Belos::CGIter::iterate(): current linear system has more than one std::vector!" );
700 
701  //Check that the input is correctly set up
702  TEUCHOS_TEST_FOR_EXCEPTION( curDim_ != prevUdim_ + 1, PCPGIterInitFailure,
703  "Belos::CGIter::iterate(): mistake in initialization !" );
704 
705 
706  const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
707  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
708 
709 
710  std::vector<int> curind(1);
711  std::vector<ScalarType> rnorm(MVT::GetNumberVecs(*cur_soln_vec));
712  if (prevUdim_ > 0){ // A-orthonalize P=Z to Uprev
713  Teuchos::RCP<MV> P;
714  curind[0] = curDim_ - 1; // column = dimension - 1
715  P = MVT::CloneViewNonConst(*U_,curind);
716  MVT::MvTransMv( one, *Cprev, *P, CZ );
717  MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *P ); // P -= U*(C'Z)
718 
719  if( debug ){
720  MVT::MvTransMv( one, *Cprev, *P, CZ );
721  std::cout << " Input CZ post ortho " << std::endl;
722  CZ.print( std::cout );
723  }
724  if( curDim_ == savedBlocks_ ){
725  std::vector<int> zero_index(1);
726  zero_index[0] = 0;
727  MVT::SetBlock( *P, zero_index, *P_ );
728  }
729  P = Teuchos::null;
730  }
731 
732  // Compute first <r,z> a.k.a. rHz
733  MVT::MvTransMv( one, *R_, *Z_, rHz );
734 
736  // iterate until the status test is satisfied
737  //
738  while (stest_->checkStatus(this) != Passed ) {
739  Teuchos::RCP<const MV> P;
740  Teuchos::RCP<MV> AP;
741  iter_++; // The next iteration begins.
742  //std::vector<int> curind(1);
743  curind[0] = curDim_ - 1; // column = dimension - 1
744  if( debug ){
745  MVT::MvNorm(*R_, rnorm);
746  std::cout << iter_ << " " << curDim_ << " " << rnorm[0] << std::endl;
747  }
748  if( prevUdim_ + iter_ < savedBlocks_ ){
749  P = MVT::CloneView(*U_,curind);
750  AP = MVT::CloneViewNonConst(*C_,curind);
751  lp_->applyOp( *P, *AP );
752  MVT::MvTransMv( one, *P, *AP, pAp );
753  }else{
754  if( prevUdim_ + iter_ == savedBlocks_ ){
755  AP = MVT::CloneViewNonConst(*C_,curind);
756  lp_->applyOp( *P_, *AP );
757  MVT::MvTransMv( one, *P_, *AP, pAp );
758  }else{
759  lp_->applyOp( *P_, *AP_ );
760  MVT::MvTransMv( one, *P_, *AP_, pAp );
761  }
762  }
763 
764  if( keepDiagonal_ && prevUdim_ + iter_ <= savedBlocks_ )
765  (*D_)(iter_ -1 ,iter_ -1 ) = pAp(0,0);
766 
767  // positive pAp required
768  TEUCHOS_TEST_FOR_EXCEPTION( pAp(0,0) <= zero, PCPGIterateFailure,
769  "Belos::CGIter::iterate(): non-positive value for p^H*A*p encountered!" );
770 
771  // alpha := <R_,Z_> / <P,AP>
772  alpha(0,0) = rHz(0,0) / pAp(0,0);
773 
774  // positive alpha required
775  TEUCHOS_TEST_FOR_EXCEPTION( alpha(0,0) <= zero, PCPGIterateFailure,
776  "Belos::CGIter::iterate(): non-positive value for alpha encountered!" );
777 
778  // solution update x += alpha * P
779  if( curDim_ < savedBlocks_ ){
780  MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P, *cur_soln_vec );
781  }else{
782  MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P_, *cur_soln_vec );
783  }
784  //lp_->updateSolution(); ... does nothing.
785  //
786  // The denominator of beta is saved before residual is updated [ old <R_, Z_> ].
787  //
788  rHz_old(0,0) = rHz(0,0);
789  //
790  // residual update R_ := R_ - alpha * AP
791  //
792  if( prevUdim_ + iter_ <= savedBlocks_ ){
793  MVT::MvAddMv( one, *R_, -alpha(0,0), *AP, *R_ );
794  AP = Teuchos::null;
795  }else{
796  MVT::MvAddMv( one, *R_, -alpha(0,0), *AP_, *R_ );
797  }
798  //
799  // update beta := [ new <R_, Z_> ] / [ old <R_, Z_> ] and the search direction p.
800  //
801  if ( lp_->getLeftPrec() != Teuchos::null ) {
802  lp_->applyLeftPrec( *R_, *Z_ );
803  } else {
804  Z_ = R_;
805  }
806  //
807  MVT::MvTransMv( one, *R_, *Z_, rHz );
808  //
809  beta(0,0) = rHz(0,0) / rHz_old(0,0);
810  //
811  if( curDim_ < savedBlocks_ ){
812  curDim_++; // update basis dim
813  curind[0] = curDim_ - 1;
814  Teuchos::RCP<MV> Pnext = MVT::CloneViewNonConst(*U_,curind);
815  MVT::MvAddMv( one, *Z_, beta(0,0), *P, *Pnext );
816  if( prevUdim_ ){ // Deflate seed space
817  MVT::MvTransMv( one, *Cprev, *Z_, CZ );
818  MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *Pnext ); // Pnext -= U*(C'Z)
819  if( debug ){
820  std::cout << " Check CZ " << std::endl;
821  MVT::MvTransMv( one, *Cprev, *Pnext, CZ );
822  CZ.print( std::cout );
823  }
824  }
825  P = Teuchos::null;
826  if( curDim_ == savedBlocks_ ){
827  std::vector<int> zero_index(1);
828  zero_index[0] = 0;
829  MVT::SetBlock( *Pnext, zero_index, *P_ );
830  }
831  Pnext = Teuchos::null;
832  }else{
833  MVT::MvAddMv( one, *Z_, beta(0,0), *P_, *P_ );
834  if( prevUdim_ ){ // Deflate seed space
835  MVT::MvTransMv( one, *Cprev, *Z_, CZ );
836  MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *P_ ); // P_ -= U*(C'Z)
837 
838  if( debug ){
839  std::cout << " Check CZ " << std::endl;
840  MVT::MvTransMv( one, *Cprev, *P_, CZ );
841  CZ.print( std::cout );
842  }
843  }
844  }
845  // CGB: 5/26/2010
846  // this RCP<const MV> P was previously a variable outside the loop. however, it didn't appear to be see any use between
847  // loop iterations. therefore, I moved it inside to avoid scoping errors with previously used variables named P.
848  // to ensure that this wasn't a bug, I verify below that we have set P == null, i.e., that we are not going to use it again
849  // same for AP
850  TEUCHOS_TEST_FOR_EXCEPTION( AP != Teuchos::null || P != Teuchos::null, std::logic_error, "Loop recurrence violated. Please contact Belos team.");
851  } // end coupled two-term recursion
852  if( prevUdim_ + iter_ < savedBlocks_ ) --curDim_; // discard negligible search direction
853  }
854 
855 } // end Belos namespace
856 
857 #endif /* BELOS_PCPG_ITER_HPP */
PCPGIterateFailure is thrown when the PCPGIter object breaks down.
Collection of types and exceptions used within the Belos solvers.
int prevUdim
Number of block columns in matrices C and U before current iteration.
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
bool isInitialized()
States whether the solver has been initialized or not.
Class which manages the output and verbosity of the Belos solvers.
Teuchos::RCP< MV > R
The current residual.
Pure virtual base class for defining the status testing capabilities of Belos.
Teuchos::RCP< MV > C
C = AU, U spans recycled subspace.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system solution?.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
Declaration of basic traits for the multivector type.
void setSize(int savedBlocks)
Set the maximum number of saved or recycled blocks used by the iterative solver.
int curDim
The current dimension of the reduction.
void resetNumIters(int iter=0)
Reset the iteration count.
Teuchos::RCP< MV > P
The current decent direction std::vector.
Teuchos::RCP< MV > AP
The matrix A applied to current decent direction std::vector.
Structure to contain pointers to PCPGIter state variables.
OperatorTraits< ScalarType, MV, OP > OPT
A pure virtual class for defining the status tests for the Belos iterative solvers.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Class which defines basic traits for the operator type.
Traits class which defines basic operations on multivectors.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > D
The current Hessenberg matrix.
virtual ~PCPGIter()
Destructor.
void setBlockSize(int blockSize)
Get the blocksize to be used by the iterative solver in solving this linear problem.
int getNumRecycledBlocks() const
Get the maximum number of recycled blocks used by the iterative solver in solving this linear problem...
void initialize()
Initialize the solver with the initial vectors from the linear problem. An exception is thrown if ini...
int getNumIters() const
Get the current iteration count.
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
PCPGIterOrthoFailure(const std::string &what_arg)
SCT::magnitudeType MagnitudeType
PCPGIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Teuchos::RCP< MV > U
The recycled subspace.
int getPrevSubspaceDim() const
Get the dimension of the search subspace used to solve the current solution to the linear problem...
void resetState()
tell the Iterator to "reset" itself; delete and rebuild the seed space.
PCPGIterOrthoFailure is thrown when the PCPGIter object is unable to compute independent direction ve...
void iterate()
PCPGIter iterates CG until the status test either requests a stop or detects an error. In the latter case, std::exception is thrown.
PCPGIterateFailure(const std::string &what_arg)
Teuchos::RCP< MV > Z
The current preconditioned residual.
Class which defines basic traits for the operator type.
PCPGIter(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< OutputManager< ScalarType > > &printer, const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &tester, const Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > &ortho, Teuchos::ParameterList &params)
PCPGIter constructor with linear problem, solver utilities, and parameter list of solver options...
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:60
Belos&#39;s templated virtual class for providing routines for orthogonalization and orthonormzalition of...
This class implements the PCPG iteration, where a single-std::vector Krylov subspace is constructed...
Belos header file which uses auto-configuration information to include necessary C++ headers...
MultiVecTraits< ScalarType, MV > MVT
int getBlockSize() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem...
int getCurSubspaceDim() const
Get the current dimension of the whole seed subspace.
PCPGIterInitFailure is thrown when the PCPGIter object is unable to generate an initial iterate in th...
PCPGIterInitFailure(const std::string &what_arg)
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Teuchos::ScalarTraits< ScalarType > SCT

Generated for Belos by doxygen 1.8.14