Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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_BLAS.hpp"
60 #include "Teuchos_LAPACK.hpp"
63 #include "Teuchos_ScalarTraits.hpp"
65 #include "Teuchos_TimeMonitor.hpp"
66 
78 namespace Belos {
79 
81 
82 
87  template <class ScalarType, class MV>
88  struct PCPGIterState {
94  int curDim;
96  int prevUdim;
97 
100 
103 
106 
109 
112 
115 
121 
123  prevUdim(0),
124  R(Teuchos::null), Z(Teuchos::null),
125  P(Teuchos::null), AP(Teuchos::null),
126  U(Teuchos::null), C(Teuchos::null),
127  D(Teuchos::null)
128  {}
129  };
130 
132 
134 
135 
147  class PCPGIterInitFailure : public BelosError {public:
148  PCPGIterInitFailure(const std::string& what_arg) : BelosError(what_arg)
149  {}};
150 
155  class PCPGIterateFailure : public BelosError {public:
156  PCPGIterateFailure(const std::string& what_arg) : BelosError(what_arg)
157  {}};
158 
159 
160 
167  class PCPGIterOrthoFailure : public BelosError {public:
168  PCPGIterOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
169  {}};
170 
177  class PCPGIterLAPACKFailure : public BelosError {public:
178  PCPGIterLAPACKFailure(const std::string& what_arg) : BelosError(what_arg)
179  {}};
180 
182 
183 
184  template<class ScalarType, class MV, class OP>
185  class PCPGIter : virtual public Iteration<ScalarType,MV,OP> {
186 
187  public:
188 
189  //
190  // Convenience typedefs
191  //
196 
198 
199 
208  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
211  Teuchos::ParameterList &params );
212 
214  virtual ~PCPGIter() {};
216 
217 
219 
220 
239  void iterate();
240 
260 
264  void initialize()
265  {
267  initialize(empty);
268  }
269 
279  state.Z = Z_; // CG state
280  state.P = P_;
281  state.AP = AP_;
282  state.R = R_;
283  state.U = U_; // seed state
284  state.C = C_;
285  state.D = D_;
286  state.curDim = curDim_;
287  state.prevUdim = prevUdim_;
288  return state;
289  }
290 
292 
293 
295 
296 
298  int getNumIters() const { return iter_; }
299 
301  void resetNumIters( int iter = 0 ) { iter_ = iter; }
302 
305  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> * /* norms */ ) const { return R_; }
306 
308 
311  Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
312 
314  int getCurSubspaceDim() const {
315  if (!initialized_) return 0;
316  return curDim_;
317  };
318 
320  int getPrevSubspaceDim() const {
321  if (!initialized_) return 0;
322  return prevUdim_;
323  };
324 
326 
327 
329 
330 
332  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
333 
335  int getBlockSize() const { return 1; }
336 
338  int getNumRecycledBlocks() const { return savedBlocks_; }
339 
341 
343  void setBlockSize(int blockSize) {
344  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
345  "Belos::PCPGIter::setBlockSize(): Cannot use a block size that is not one.");
346  }
347 
349  void setSize( int savedBlocks );
350 
352  bool isInitialized() { return initialized_; }
353 
355  void resetState();
356 
358 
359  private:
360 
361  //
362  // Internal methods
363  //
365  void setStateSize();
366 
367  //
368  // Classes inputed through constructor that define the linear problem to be solved.
369  //
374 
375  //
376  // Algorithmic parameters
377  // savedBlocks_ is the number of blocks allocated for the reused subspace
378  int savedBlocks_;
379  //
380  //
381  // Current solver state
382  //
383  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
384  // is capable of running; _initialize is controlled by the initialize() member method
385  // For the implications of the state of initialized_, please see documentation for initialize()
386  bool initialized_;
387 
388  // stateStorageInitialized_ indicates that the state storage has be initialized to the current
389  // savedBlocks_. State storage initialization may be postponed if the linear problem was
390  // generated without either the right-hand side or solution vectors.
391  bool stateStorageInitialized_;
392 
393  // keepDiagonal_ specifies that the iteration must keep the diagonal matrix of pivots
394  bool keepDiagonal_;
395 
396  // initDiagonal_ specifies that the iteration will reinitialize the diagonal matrix by zeroing
397  // out all entries before an iteration is started.
398  bool initDiagonal_;
399 
400  // Current subspace dimension
401  int curDim_;
402 
403  // Dimension of seed space used to solve current linear system
404  int prevUdim_;
405 
406  // Number of iterations performed
407  int iter_;
408  //
409  // State Storage ... of course this part is different for CG
410  //
411  // Residual
412  Teuchos::RCP<MV> R_;
413  //
414  // Preconditioned residual
415  Teuchos::RCP<MV> Z_;
416  //
417  // Direction std::vector
418  Teuchos::RCP<MV> P_;
419  //
420  // Operator applied to direction std::vector
421  Teuchos::RCP<MV> AP_;
422  //
423  // Recycled subspace vectors.
424  Teuchos::RCP<MV> U_;
425  //
426  // C = A * U, linear system is Ax=b
427  Teuchos::RCP<MV> C_;
428  //
429  // Projected matrices
430  // D_ : Diagonal matrix of pivots D = P'AP
432  };
433 
435  // Constructor.
436  template<class ScalarType, class MV, class OP>
438  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
441  Teuchos::ParameterList &params ):
442  lp_(problem),
443  om_(printer),
444  stest_(tester),
445  ortho_(ortho),
446  savedBlocks_(0),
447  initialized_(false),
448  stateStorageInitialized_(false),
449  keepDiagonal_(false),
450  initDiagonal_(false),
451  curDim_(0),
452  prevUdim_(0),
453  iter_(0)
454  {
455  // Get the maximum number of blocks allowed for this Krylov subspace
456 
457  TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Saved Blocks"), std::invalid_argument,
458  "Belos::PCPGIter::constructor: mandatory parameter \"Saved Blocks\" is not specified.");
459  int rb = Teuchos::getParameter<int>(params, "Saved Blocks");
460 
461  // Find out whether we are saving the Diagonal matrix.
462  keepDiagonal_ = params.get("Keep Diagonal", false);
463 
464  // Find out whether we are initializing the Diagonal matrix.
465  initDiagonal_ = params.get("Initialize Diagonal", false);
466 
467  // Set the number of blocks and allocate data
468  setSize( rb );
469  }
470 
472  // Set the block size and adjust as necessary
473  template <class ScalarType, class MV, class OP>
474  void PCPGIter<ScalarType,MV,OP>::setSize( int savedBlocks )
475  {
476  // allocate space only; perform no computation
477  // Any change in size invalidates the state of the solver as implemented here.
478 
479  TEUCHOS_TEST_FOR_EXCEPTION(savedBlocks <= 0, std::invalid_argument, "Belos::PCPGIter::setSize() was passed a non-positive argument for \"Num Saved Blocks\".");
480 
481  if ( savedBlocks_ != savedBlocks) {
482  stateStorageInitialized_ = false;
483  savedBlocks_ = savedBlocks;
484  initialized_ = false;
485  curDim_ = 0;
486  prevUdim_ = 0;
487  setStateSize(); // Use the current savedBlocks_ to initialize the state storage.
488  }
489  }
490 
492  // Enable the reuse of a single solver object for completely different linear systems
493  template <class ScalarType, class MV, class OP>
495  {
496  stateStorageInitialized_ = false;
497  initialized_ = false;
498  curDim_ = 0;
499  prevUdim_ = 0;
500  setStateSize();
501  }
502 
504  // Setup the state storage. Called by either initialize or, if savedBlocks_ changes, setSize.
505  template <class ScalarType, class MV, class OP>
507  {
508  if (!stateStorageInitialized_) {
509 
510  // Check if there is any multivector to clone from.
511  Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
512  Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
513  if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
514  return; // postpone exception
515  }
516  else {
517 
519  // blockSize*recycledBlocks dependent
520  int newsd = savedBlocks_ ; //int newsd = blockSize_* savedBlocks_ ;
521  //
522  // Initialize the CG state storage
523  // If the subspace is not initialized, generate it using the LHS or RHS from lp_.
524  // Generate CG state only if it does not exist, otherwise resize it.
525  if (Z_ == Teuchos::null) {
526  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
527  Z_ = MVT::Clone( *tmp, 1 );
528  }
529  if (P_ == Teuchos::null) {
530  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
531  P_ = MVT::Clone( *tmp, 1 );
532  }
533  if (AP_ == Teuchos::null) {
534  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
535  AP_ = MVT::Clone( *tmp, 1 );
536  }
537 
538  if (C_ == Teuchos::null) {
539 
540  // Get the multivector that is not null.
541  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
542  TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
543  "Belos::PCPGIter::setStateSize(): linear problem does not specify multivectors to clone from.");
544  TEUCHOS_TEST_FOR_EXCEPTION( 0 != prevUdim_,std::invalid_argument,
545  "Belos::PCPGIter::setStateSize(): prevUdim not zero and C is null.");
546  C_ = MVT::Clone( *tmp, savedBlocks_ );
547  }
548  else {
549  // Generate C_ by cloning itself ONLY if more space is needed.
550  if (MVT::GetNumberVecs(*C_) < savedBlocks_ ) {
551  Teuchos::RCP<const MV> tmp = C_;
552  C_ = MVT::Clone( *tmp, savedBlocks_ );
553  }
554  }
555  if (U_ == Teuchos::null) {
556  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
557  TEUCHOS_TEST_FOR_EXCEPTION( 0 != prevUdim_,std::invalid_argument,
558  "Belos::PCPGIter::setStateSize(): prevUdim not zero and U is null.");
559  U_ = MVT::Clone( *tmp, savedBlocks_ );
560  }
561  else {
562  // Generate U_ by cloning itself ONLY if more space is needed.
563  if (MVT::GetNumberVecs(*U_) < savedBlocks_ ) {
564  Teuchos::RCP<const MV> tmp = U_;
565  U_ = MVT::Clone( *tmp, savedBlocks_ );
566  }
567  }
568  if (keepDiagonal_) {
569  if (D_ == Teuchos::null) {
571  }
572  if (initDiagonal_) {
573  D_->shape( newsd, newsd );
574  }
575  else {
576  if (D_->numRows() < newsd || D_->numCols() < newsd) {
577  D_->shapeUninitialized( newsd, newsd );
578  }
579  }
580  }
581  // State storage has now been initialized.
582  stateStorageInitialized_ = true;
583  } // if there is a vector to clone from
584  } // if !stateStorageInitialized_
585  } // end of setStateSize
586 
588  // Initialize the iteration object
589  template <class ScalarType, class MV, class OP>
591  {
592 
593  TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
594  "Belos::PCPGIter::initialize(): Cannot initialize state storage!");
595 
596  // Requirements: R_ and consistent multivectors widths and lengths
597  //
598  std::string errstr("Belos::PCPGIter::initialize(): Specified multivectors must have a consistent length and width.");
599 
600  if (newstate.R != Teuchos::null){
601 
602  R_ = newstate.R; // SolverManager::R_ == newstate.R == Iterator::R_
603  if (newstate.U == Teuchos::null){
604  prevUdim_ = 0;
605  newstate.U = U_;
606  newstate.C = C_;
607  }
608  else {
609  prevUdim_ = newstate.curDim;
610  if (newstate.C == Teuchos::null){ // Stub for new feature
611  std::vector<int> index(prevUdim_);
612  for (int i=0; i< prevUdim_; ++i)
613  index[i] = i;
614  Teuchos::RCP<const MV> Ukeff = MVT::CloneView( *newstate.U, index );
615  newstate.C = MVT::Clone( *newstate.U, prevUdim_ );
616  Teuchos::RCP<MV> Ckeff = MVT::CloneViewNonConst( *newstate.C, index );
617  lp_->apply( *Ukeff, *Ckeff );
618  }
619  curDim_ = prevUdim_ ;
620  }
621 
622  // Initialize the state storage if not already allocated in the constructor
623  if (!stateStorageInitialized_)
624  setStateSize();
625 
626  //TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.V) != MVT::GetGlobalLength(*V_), std::invalid_argument, errstr );
627  //TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < 1, std::invalid_argument, errstr );
628 
629  newstate.prevUdim = prevUdim_; // big change in functionality from GCRODR
630  newstate.curDim = curDim_;
631 
632  //TEUCHOS_TEST_FOR_EXCEPTION(newstate.z->numRows() < curDim_ || newstate.z->numCols() < 1, std::invalid_argument, errstr);
633 
634  std::vector<int> zero_index(1);
635  zero_index[0] = 0;
636  if ( lp_->getLeftPrec() != Teuchos::null ) { // Compute the initial search direction
637  lp_->applyLeftPrec( *R_, *Z_ );
638  MVT::SetBlock( *Z_, zero_index , *P_ ); // P(:,zero_index) := Z
639  } else {
640  Z_ = R_;
641  MVT::SetBlock( *R_, zero_index, *P_ );
642  }
643 
644  std::vector<int> nextind(1);
645  nextind[0] = curDim_;
646 
647  MVT::SetBlock( *P_, nextind, *newstate.U ); // U(:,curDim_ ) := P_
648 
649  ++curDim_;
650  newstate.curDim = curDim_;
651 
652  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.U) != savedBlocks_ ,
653  std::invalid_argument, errstr );
654  if (newstate.U != U_) { // Why this is needed?
655  U_ = newstate.U;
656  }
657 
658  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.C) != savedBlocks_ ,
659  std::invalid_argument, errstr );
660  if (newstate.C != C_) {
661  C_ = newstate.C;
662  }
663  }
664  else {
665 
666  TEUCHOS_TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
667  "Belos::PCPGIter::initialize(): PCPGStateIterState does not have initial kernel R_0.");
668  }
669 
670  // the solver is initialized
671  initialized_ = true;
672  }
673 
674 
676  // Iterate until the status test informs us we should stop.
677  template <class ScalarType, class MV, class OP>
679  {
680  //
681  // Allocate/initialize data structures
682  //
683  if (initialized_ == false) {
684  initialize();
685  }
686  const bool debug = false;
687 
688  // Allocate memory for scalars.
692  Teuchos::SerialDenseMatrix<int,ScalarType> rHz( 1, 1 ), rHz_old( 1, 1 );
693 
694  if( iter_ != 0 )
695  std::cout << " Iterate Warning: begin from nonzero iter_ ?" << std::endl; //DMD
696 
697  // GenOrtho Project Stubs
698  std::vector<int> prevInd;
702 
703  if( prevUdim_ ){
704  prevInd.resize( prevUdim_ );
705  for( int i=0; i<prevUdim_ ; i++) prevInd[i] = i;
706  CZ.reshape( prevUdim_ , 1 );
707  Uprev = MVT::CloneView(*U_, prevInd);
708  Cprev = MVT::CloneView(*C_, prevInd);
709  }
710 
711  // Get the current solution std::vector.
712  Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
713 
714  // Check that the current solution std::vector only has one column.
715  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, PCPGIterInitFailure,
716  "Belos::CGIter::iterate(): current linear system has more than one std::vector!" );
717 
718  //Check that the input is correctly set up
719  TEUCHOS_TEST_FOR_EXCEPTION( curDim_ != prevUdim_ + 1, PCPGIterInitFailure,
720  "Belos::CGIter::iterate(): mistake in initialization !" );
721 
722 
723  const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
724  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
725 
726 
727  std::vector<int> curind(1);
728  std::vector<ScalarType> rnorm(MVT::GetNumberVecs(*cur_soln_vec));
729  if (prevUdim_ > 0){ // A-orthonalize P=Z to Uprev
730  Teuchos::RCP<MV> P;
731  curind[0] = curDim_ - 1; // column = dimension - 1
732  P = MVT::CloneViewNonConst(*U_,curind);
733  MVT::MvTransMv( one, *Cprev, *P, CZ );
734  MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *P ); // P -= U*(C'Z)
735 
736  if( debug ){
737  MVT::MvTransMv( one, *Cprev, *P, CZ );
738  std::cout << " Input CZ post ortho " << std::endl;
739  CZ.print( std::cout );
740  }
741  if( curDim_ == savedBlocks_ ){
742  std::vector<int> zero_index(1);
743  zero_index[0] = 0;
744  MVT::SetBlock( *P, zero_index, *P_ );
745  }
746  P = Teuchos::null;
747  }
748 
749  // Compute first <r,z> a.k.a. rHz
750  MVT::MvTransMv( one, *R_, *Z_, rHz );
751 
753  // iterate until the status test is satisfied
754  //
755  while (stest_->checkStatus(this) != Passed ) {
757  Teuchos::RCP<MV> AP;
758  iter_++; // The next iteration begins.
759  //std::vector<int> curind(1);
760  curind[0] = curDim_ - 1; // column = dimension - 1
761  if( debug ){
762  MVT::MvNorm(*R_, rnorm);
763  std::cout << iter_ << " " << curDim_ << " " << rnorm[0] << std::endl;
764  }
765  if( prevUdim_ + iter_ < savedBlocks_ ){
766  P = MVT::CloneView(*U_,curind);
767  AP = MVT::CloneViewNonConst(*C_,curind);
768  lp_->applyOp( *P, *AP );
769  MVT::MvTransMv( one, *P, *AP, pAp );
770  }else{
771  if( prevUdim_ + iter_ == savedBlocks_ ){
772  AP = MVT::CloneViewNonConst(*C_,curind);
773  lp_->applyOp( *P_, *AP );
774  MVT::MvTransMv( one, *P_, *AP, pAp );
775  }else{
776  lp_->applyOp( *P_, *AP_ );
777  MVT::MvTransMv( one, *P_, *AP_, pAp );
778  }
779  }
780 
781  if( keepDiagonal_ && prevUdim_ + iter_ <= savedBlocks_ )
782  (*D_)(iter_ -1 ,iter_ -1 ) = pAp(0,0);
783 
784  // positive pAp required
786  "Belos::CGIter::iterate(): non-positive value for p^H*A*p encountered!" );
787 
788  // alpha := <R_,Z_> / <P,AP>
789  alpha(0,0) = rHz(0,0) / pAp(0,0);
790 
791  // positive alpha required
792  TEUCHOS_TEST_FOR_EXCEPTION( alpha(0,0) <= zero, PCPGIterateFailure,
793  "Belos::CGIter::iterate(): non-positive value for alpha encountered!" );
794 
795  // solution update x += alpha * P
796  if( curDim_ < savedBlocks_ ){
797  MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P, *cur_soln_vec );
798  }else{
799  MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P_, *cur_soln_vec );
800  }
801  //lp_->updateSolution(); ... does nothing.
802  //
803  // The denominator of beta is saved before residual is updated [ old <R_, Z_> ].
804  //
805  rHz_old(0,0) = rHz(0,0);
806  //
807  // residual update R_ := R_ - alpha * AP
808  //
809  if( prevUdim_ + iter_ <= savedBlocks_ ){
810  MVT::MvAddMv( one, *R_, -alpha(0,0), *AP, *R_ );
811  AP = Teuchos::null;
812  }else{
813  MVT::MvAddMv( one, *R_, -alpha(0,0), *AP_, *R_ );
814  }
815  //
816  // update beta := [ new <R_, Z_> ] / [ old <R_, Z_> ] and the search direction p.
817  //
818  if ( lp_->getLeftPrec() != Teuchos::null ) {
819  lp_->applyLeftPrec( *R_, *Z_ );
820  } else {
821  Z_ = R_;
822  }
823  //
824  MVT::MvTransMv( one, *R_, *Z_, rHz );
825  //
826  beta(0,0) = rHz(0,0) / rHz_old(0,0);
827  //
828  if( curDim_ < savedBlocks_ ){
829  curDim_++; // update basis dim
830  curind[0] = curDim_ - 1;
831  Teuchos::RCP<MV> Pnext = MVT::CloneViewNonConst(*U_,curind);
832  MVT::MvAddMv( one, *Z_, beta(0,0), *P, *Pnext );
833  if( prevUdim_ ){ // Deflate seed space
834  MVT::MvTransMv( one, *Cprev, *Z_, CZ );
835  MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *Pnext ); // Pnext -= U*(C'Z)
836  if( debug ){
837  std::cout << " Check CZ " << std::endl;
838  MVT::MvTransMv( one, *Cprev, *Pnext, CZ );
839  CZ.print( std::cout );
840  }
841  }
842  P = Teuchos::null;
843  if( curDim_ == savedBlocks_ ){
844  std::vector<int> zero_index(1);
845  zero_index[0] = 0;
846  MVT::SetBlock( *Pnext, zero_index, *P_ );
847  }
848  Pnext = Teuchos::null;
849  }else{
850  MVT::MvAddMv( one, *Z_, beta(0,0), *P_, *P_ );
851  if( prevUdim_ ){ // Deflate seed space
852  MVT::MvTransMv( one, *Cprev, *Z_, CZ );
853  MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *P_ ); // P_ -= U*(C'Z)
854 
855  if( debug ){
856  std::cout << " Check CZ " << std::endl;
857  MVT::MvTransMv( one, *Cprev, *P_, CZ );
858  CZ.print( std::cout );
859  }
860  }
861  }
862  // CGB: 5/26/2010
863  // this RCP<const MV> P was previously a variable outside the loop. however, it didn't appear to be see any use between
864  // loop iterations. therefore, I moved it inside to avoid scoping errors with previously used variables named P.
865  // 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
866  // same for AP
867  TEUCHOS_TEST_FOR_EXCEPTION( AP != Teuchos::null || P != Teuchos::null, std::logic_error, "Loop recurrence violated. Please contact Belos team.");
868  } // end coupled two-term recursion
869  if( prevUdim_ + iter_ < savedBlocks_ ) --curDim_; // discard negligible search direction
870  }
871 
872 } // end Belos namespace
873 
874 #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.
virtual void print(std::ostream &os) const
Teuchos::RCP< MV > R
The current residual.
PCPGIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Pure virtual base class for defining the status testing capabilities of Belos.
Teuchos::RCP< MV > C
C = AU, U spans recycled subspace.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
int getCurSubspaceDim() const
Get the current dimension of the whole seed subspace.
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
int getPrevSubspaceDim() const
Get the dimension of the search subspace used to solve the current solution to the linear problem...
A pure virtual class for defining the status tests for the Belos iterative solvers.
int getBlockSize() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem...
Class which defines basic traits for the operator type.
PCPGIterLAPACKFailure is thrown when a nonzero return value is passed back from an LAPACK routine...
PCPGIterLAPACKFailure(const std::string &what_arg)
Traits class which defines basic operations on multivectors.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > D
The current Hessenberg matrix.
bool isParameter(const std::string &name) const
virtual ~PCPGIter()
Destructor.
void setBlockSize(int blockSize)
Get the blocksize to be used by the iterative solver in solving this linear problem.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
int getNumIters() const
Get the current iteration count.
void initialize()
Initialize the solver with the initial vectors from the linear problem. An exception is thrown if ini...
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)
int getNumRecycledBlocks() const
Get the maximum number of recycled blocks used by the iterative solver in solving this linear problem...
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
SCT::magnitudeType MagnitudeType
Teuchos::RCP< MV > U
The recycled subspace.
void resetState()
tell the Iterator to &quot;reset&quot; itself; delete and rebuild the seed space.
PCPGIterOrthoFailure is thrown when the PCPGIter object is unable to compute independent direction ve...
int reshape(OrdinalType numRows, OrdinalType numCols)
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.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
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
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system solution?.
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 on Fri Aug 14 2020 10:48:34 for Belos by doxygen 1.8.5