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 #include "BelosCGIteration.hpp"
59 
62 #include "Teuchos_ScalarTraits.hpp"
64 #include "Teuchos_TimeMonitor.hpp"
65 
77 namespace Belos {
78 
80 
81 
86  template <class ScalarType, class MV>
87  struct PCPGIterState {
93  int curDim;
95  int prevUdim;
96 
99 
102 
105 
108 
111 
114 
120 
122  prevUdim(0),
123  R(Teuchos::null), Z(Teuchos::null),
124  P(Teuchos::null), AP(Teuchos::null),
125  U(Teuchos::null), C(Teuchos::null),
126  D(Teuchos::null)
127  {}
128  };
129 
131 
132  template<class ScalarType, class MV, class OP>
133  class PCPGIter : virtual public Iteration<ScalarType,MV,OP> {
134 
135  public:
136 
137  //
138  // Convenience typedefs
139  //
144 
146 
147 
156  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
159  Teuchos::ParameterList &params );
160 
162  virtual ~PCPGIter() {};
164 
165 
167 
168 
187  void iterate();
188 
208 
212  void initialize()
213  {
215  initialize(empty);
216  }
217 
227  state.Z = Z_; // CG state
228  state.P = P_;
229  state.AP = AP_;
230  state.R = R_;
231  state.U = U_; // seed state
232  state.C = C_;
233  state.D = D_;
234  state.curDim = curDim_;
235  state.prevUdim = prevUdim_;
236  return state;
237  }
238 
240 
241 
243 
244 
246  int getNumIters() const { return iter_; }
247 
249  void resetNumIters( int iter = 0 ) { iter_ = iter; }
250 
253  Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> * /* norms */ ) const { return R_; }
254 
256 
259  Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
260 
262  int getCurSubspaceDim() const {
263  if (!initialized_) return 0;
264  return curDim_;
265  };
266 
268  int getPrevSubspaceDim() const {
269  if (!initialized_) return 0;
270  return prevUdim_;
271  };
272 
274 
275 
277 
278 
280  const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
281 
283  int getBlockSize() const { return 1; }
284 
286  int getNumRecycledBlocks() const { return savedBlocks_; }
287 
289 
291  void setBlockSize(int blockSize) {
292  TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
293  "Belos::PCPGIter::setBlockSize(): Cannot use a block size that is not one.");
294  }
295 
297  void setSize( int savedBlocks );
298 
300  bool isInitialized() { return initialized_; }
301 
303  void resetState();
304 
306 
307  private:
308 
309  //
310  // Internal methods
311  //
313  void setStateSize();
314 
315  //
316  // Classes inputed through constructor that define the linear problem to be solved.
317  //
322 
323  //
324  // Algorithmic parameters
325  // savedBlocks_ is the number of blocks allocated for the reused subspace
326  int savedBlocks_;
327  //
328  //
329  // Current solver state
330  //
331  // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
332  // is capable of running; _initialize is controlled by the initialize() member method
333  // For the implications of the state of initialized_, please see documentation for initialize()
334  bool initialized_;
335 
336  // stateStorageInitialized_ indicates that the state storage has be initialized to the current
337  // savedBlocks_. State storage initialization may be postponed if the linear problem was
338  // generated without either the right-hand side or solution vectors.
339  bool stateStorageInitialized_;
340 
341  // keepDiagonal_ specifies that the iteration must keep the diagonal matrix of pivots
342  bool keepDiagonal_;
343 
344  // initDiagonal_ specifies that the iteration will reinitialize the diagonal matrix by zeroing
345  // out all entries before an iteration is started.
346  bool initDiagonal_;
347 
348  // Current subspace dimension
349  int curDim_;
350 
351  // Dimension of seed space used to solve current linear system
352  int prevUdim_;
353 
354  // Number of iterations performed
355  int iter_;
356  //
357  // State Storage ... of course this part is different for CG
358  //
359  // Residual
360  Teuchos::RCP<MV> R_;
361  //
362  // Preconditioned residual
363  Teuchos::RCP<MV> Z_;
364  //
365  // Direction std::vector
366  Teuchos::RCP<MV> P_;
367  //
368  // Operator applied to direction std::vector
369  Teuchos::RCP<MV> AP_;
370  //
371  // Recycled subspace vectors.
372  Teuchos::RCP<MV> U_;
373  //
374  // C = A * U, linear system is Ax=b
375  Teuchos::RCP<MV> C_;
376  //
377  // Projected matrices
378  // D_ : Diagonal matrix of pivots D = P'AP
380  };
381 
383  // Constructor.
384  template<class ScalarType, class MV, class OP>
386  const Teuchos::RCP<OutputManager<ScalarType> > &printer,
389  Teuchos::ParameterList &params ):
390  lp_(problem),
391  om_(printer),
392  stest_(tester),
393  ortho_(ortho),
394  savedBlocks_(0),
395  initialized_(false),
396  stateStorageInitialized_(false),
397  keepDiagonal_(false),
398  initDiagonal_(false),
399  curDim_(0),
400  prevUdim_(0),
401  iter_(0)
402  {
403  // Get the maximum number of blocks allowed for this Krylov subspace
404 
405  TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Saved Blocks"), std::invalid_argument,
406  "Belos::PCPGIter::constructor: mandatory parameter \"Saved Blocks\" is not specified.");
407  int rb = Teuchos::getParameter<int>(params, "Saved Blocks");
408 
409  // Find out whether we are saving the Diagonal matrix.
410  keepDiagonal_ = params.get("Keep Diagonal", false);
411 
412  // Find out whether we are initializing the Diagonal matrix.
413  initDiagonal_ = params.get("Initialize Diagonal", false);
414 
415  // Set the number of blocks and allocate data
416  setSize( rb );
417  }
418 
420  // Set the block size and adjust as necessary
421  template <class ScalarType, class MV, class OP>
422  void PCPGIter<ScalarType,MV,OP>::setSize( int savedBlocks )
423  {
424  // allocate space only; perform no computation
425  // Any change in size invalidates the state of the solver as implemented here.
426 
427  TEUCHOS_TEST_FOR_EXCEPTION(savedBlocks <= 0, std::invalid_argument, "Belos::PCPGIter::setSize() was passed a non-positive argument for \"Num Saved Blocks\".");
428 
429  if ( savedBlocks_ != savedBlocks) {
430  stateStorageInitialized_ = false;
431  savedBlocks_ = savedBlocks;
432  initialized_ = false;
433  curDim_ = 0;
434  prevUdim_ = 0;
435  setStateSize(); // Use the current savedBlocks_ to initialize the state storage.
436  }
437  }
438 
440  // Enable the reuse of a single solver object for completely different linear systems
441  template <class ScalarType, class MV, class OP>
443  {
444  stateStorageInitialized_ = false;
445  initialized_ = false;
446  curDim_ = 0;
447  prevUdim_ = 0;
448  setStateSize();
449  }
450 
452  // Setup the state storage. Called by either initialize or, if savedBlocks_ changes, setSize.
453  template <class ScalarType, class MV, class OP>
455  {
456  if (!stateStorageInitialized_) {
457 
458  // Check if there is any multivector to clone from.
459  Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
460  Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
461  if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
462  return; // postpone exception
463  }
464  else {
465 
467  // blockSize*recycledBlocks dependent
468  int newsd = savedBlocks_ ; //int newsd = blockSize_* savedBlocks_ ;
469  //
470  // Initialize the CG state storage
471  // If the subspace is not initialized, generate it using the LHS or RHS from lp_.
472  // Generate CG state only if it does not exist, otherwise resize it.
473  if (Z_ == Teuchos::null) {
474  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
475  Z_ = MVT::Clone( *tmp, 1 );
476  }
477  if (P_ == Teuchos::null) {
478  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
479  P_ = MVT::Clone( *tmp, 1 );
480  }
481  if (AP_ == Teuchos::null) {
482  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
483  AP_ = MVT::Clone( *tmp, 1 );
484  }
485 
486  if (C_ == Teuchos::null) {
487 
488  // Get the multivector that is not null.
489  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
490  TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
491  "Belos::PCPGIter::setStateSize(): linear problem does not specify multivectors to clone from.");
492  TEUCHOS_TEST_FOR_EXCEPTION( 0 != prevUdim_,std::invalid_argument,
493  "Belos::PCPGIter::setStateSize(): prevUdim not zero and C is null.");
494  C_ = MVT::Clone( *tmp, savedBlocks_ );
495  }
496  else {
497  // Generate C_ by cloning itself ONLY if more space is needed.
498  if (MVT::GetNumberVecs(*C_) < savedBlocks_ ) {
499  Teuchos::RCP<const MV> tmp = C_;
500  C_ = MVT::Clone( *tmp, savedBlocks_ );
501  }
502  }
503  if (U_ == Teuchos::null) {
504  Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
505  TEUCHOS_TEST_FOR_EXCEPTION( 0 != prevUdim_,std::invalid_argument,
506  "Belos::PCPGIter::setStateSize(): prevUdim not zero and U is null.");
507  U_ = MVT::Clone( *tmp, savedBlocks_ );
508  }
509  else {
510  // Generate U_ by cloning itself ONLY if more space is needed.
511  if (MVT::GetNumberVecs(*U_) < savedBlocks_ ) {
512  Teuchos::RCP<const MV> tmp = U_;
513  U_ = MVT::Clone( *tmp, savedBlocks_ );
514  }
515  }
516  if (keepDiagonal_) {
517  if (D_ == Teuchos::null) {
519  }
520  if (initDiagonal_) {
521  D_->shape( newsd, newsd );
522  }
523  else {
524  if (D_->numRows() < newsd || D_->numCols() < newsd) {
525  D_->shapeUninitialized( newsd, newsd );
526  }
527  }
528  }
529  // State storage has now been initialized.
530  stateStorageInitialized_ = true;
531  } // if there is a vector to clone from
532  } // if !stateStorageInitialized_
533  } // end of setStateSize
534 
536  // Initialize the iteration object
537  template <class ScalarType, class MV, class OP>
539  {
540 
541  TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
542  "Belos::PCPGIter::initialize(): Cannot initialize state storage!");
543 
544  // Requirements: R_ and consistent multivectors widths and lengths
545  //
546  std::string errstr("Belos::PCPGIter::initialize(): Specified multivectors must have a consistent length and width.");
547 
548  if (newstate.R != Teuchos::null){
549 
550  R_ = newstate.R; // SolverManager::R_ == newstate.R == Iterator::R_
551  if (newstate.U == Teuchos::null){
552  prevUdim_ = 0;
553  newstate.U = U_;
554  newstate.C = C_;
555  }
556  else {
557  prevUdim_ = newstate.curDim;
558  if (newstate.C == Teuchos::null){ // Stub for new feature
559  std::vector<int> index(prevUdim_);
560  for (int i=0; i< prevUdim_; ++i)
561  index[i] = i;
562  Teuchos::RCP<const MV> Ukeff = MVT::CloneView( *newstate.U, index );
563  newstate.C = MVT::Clone( *newstate.U, prevUdim_ );
564  Teuchos::RCP<MV> Ckeff = MVT::CloneViewNonConst( *newstate.C, index );
565  lp_->apply( *Ukeff, *Ckeff );
566  }
567  curDim_ = prevUdim_ ;
568  }
569 
570  // Initialize the state storage if not already allocated in the constructor
571  if (!stateStorageInitialized_)
572  setStateSize();
573 
574  //TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.V) != MVT::GetGlobalLength(*V_), std::invalid_argument, errstr );
575  //TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < 1, std::invalid_argument, errstr );
576 
577  newstate.prevUdim = prevUdim_; // big change in functionality from GCRODR
578  newstate.curDim = curDim_;
579 
580  //TEUCHOS_TEST_FOR_EXCEPTION(newstate.z->numRows() < curDim_ || newstate.z->numCols() < 1, std::invalid_argument, errstr);
581 
582  std::vector<int> zero_index(1);
583  zero_index[0] = 0;
584  if ( lp_->getLeftPrec() != Teuchos::null ) { // Compute the initial search direction
585  lp_->applyLeftPrec( *R_, *Z_ );
586  MVT::SetBlock( *Z_, zero_index , *P_ ); // P(:,zero_index) := Z
587  } else {
588  Z_ = R_;
589  MVT::SetBlock( *R_, zero_index, *P_ );
590  }
591 
592  std::vector<int> nextind(1);
593  nextind[0] = curDim_;
594 
595  MVT::SetBlock( *P_, nextind, *newstate.U ); // U(:,curDim_ ) := P_
596 
597  ++curDim_;
598  newstate.curDim = curDim_;
599 
600  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.U) != savedBlocks_ ,
601  std::invalid_argument, errstr );
602  if (newstate.U != U_) { // Why this is needed?
603  U_ = newstate.U;
604  }
605 
606  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.C) != savedBlocks_ ,
607  std::invalid_argument, errstr );
608  if (newstate.C != C_) {
609  C_ = newstate.C;
610  }
611  }
612  else {
613 
614  TEUCHOS_TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
615  "Belos::PCPGIter::initialize(): PCPGStateIterState does not have initial kernel R_0.");
616  }
617 
618  // the solver is initialized
619  initialized_ = true;
620  }
621 
622 
624  // Iterate until the status test informs us we should stop.
625  template <class ScalarType, class MV, class OP>
627  {
628  //
629  // Allocate/initialize data structures
630  //
631  if (initialized_ == false) {
632  initialize();
633  }
634  const bool debug = false;
635 
636  // Allocate memory for scalars.
640  Teuchos::SerialDenseMatrix<int,ScalarType> rHz( 1, 1 ), rHz_old( 1, 1 );
641 
642  if( iter_ != 0 )
643  std::cout << " Iterate Warning: begin from nonzero iter_ ?" << std::endl; //DMD
644 
645  // GenOrtho Project Stubs
646  std::vector<int> prevInd;
650 
651  if( prevUdim_ ){
652  prevInd.resize( prevUdim_ );
653  for( int i=0; i<prevUdim_ ; i++) prevInd[i] = i;
654  CZ.reshape( prevUdim_ , 1 );
655  Uprev = MVT::CloneView(*U_, prevInd);
656  Cprev = MVT::CloneView(*C_, prevInd);
657  }
658 
659  // Get the current solution std::vector.
660  Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
661 
662  // Check that the current solution std::vector only has one column.
663  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, CGIterationInitFailure,
664  "Belos::PCPGIter::iterate(): current linear system has more than one std::vector!" );
665 
666  //Check that the input is correctly set up
667  TEUCHOS_TEST_FOR_EXCEPTION( curDim_ != prevUdim_ + 1, CGIterationInitFailure,
668  "Belos::PCPGIter::iterate(): mistake in initialization !" );
669 
670 
671  const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
672  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
673 
674 
675  std::vector<int> curind(1);
676  std::vector<ScalarType> rnorm(MVT::GetNumberVecs(*cur_soln_vec));
677  if (prevUdim_ > 0){ // A-orthonalize P=Z to Uprev
678  Teuchos::RCP<MV> P;
679  curind[0] = curDim_ - 1; // column = dimension - 1
680  P = MVT::CloneViewNonConst(*U_,curind);
681  MVT::MvTransMv( one, *Cprev, *P, CZ );
682  MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *P ); // P -= U*(C'Z)
683 
684  if( debug ){
685  MVT::MvTransMv( one, *Cprev, *P, CZ );
686  std::cout << " Input CZ post ortho " << std::endl;
687  CZ.print( std::cout );
688  }
689  if( curDim_ == savedBlocks_ ){
690  std::vector<int> zero_index(1);
691  zero_index[0] = 0;
692  MVT::SetBlock( *P, zero_index, *P_ );
693  }
694  P = Teuchos::null;
695  }
696 
697  // Compute first <r,z> a.k.a. rHz
698  MVT::MvTransMv( one, *R_, *Z_, rHz );
699 
701  // iterate until the status test is satisfied
702  //
703  while (stest_->checkStatus(this) != Passed ) {
705  Teuchos::RCP<MV> AP;
706  iter_++; // The next iteration begins.
707  //std::vector<int> curind(1);
708  curind[0] = curDim_ - 1; // column = dimension - 1
709  if( debug ){
710  MVT::MvNorm(*R_, rnorm);
711  std::cout << iter_ << " " << curDim_ << " " << rnorm[0] << std::endl;
712  }
713  if( prevUdim_ + iter_ < savedBlocks_ ){
714  P = MVT::CloneView(*U_,curind);
715  AP = MVT::CloneViewNonConst(*C_,curind);
716  lp_->applyOp( *P, *AP );
717  MVT::MvTransMv( one, *P, *AP, pAp );
718  }else{
719  if( prevUdim_ + iter_ == savedBlocks_ ){
720  AP = MVT::CloneViewNonConst(*C_,curind);
721  lp_->applyOp( *P_, *AP );
722  MVT::MvTransMv( one, *P_, *AP, pAp );
723  }else{
724  lp_->applyOp( *P_, *AP_ );
725  MVT::MvTransMv( one, *P_, *AP_, pAp );
726  }
727  }
728 
729  if( keepDiagonal_ && prevUdim_ + iter_ <= savedBlocks_ )
730  (*D_)(iter_ -1 ,iter_ -1 ) = pAp(0,0);
731 
732  // positive pAp required
734  "Belos::PCPGIter::iterate(): non-positive value for p^H*A*p encountered!" );
735 
736  // alpha := <R_,Z_> / <P,AP>
737  alpha(0,0) = rHz(0,0) / pAp(0,0);
738 
739  // positive alpha required
741  "Belos::PCPGIter::iterate(): non-positive value for alpha encountered!" );
742 
743  // solution update x += alpha * P
744  if( curDim_ < savedBlocks_ ){
745  MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P, *cur_soln_vec );
746  }else{
747  MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P_, *cur_soln_vec );
748  }
749  //lp_->updateSolution(); ... does nothing.
750  //
751  // The denominator of beta is saved before residual is updated [ old <R_, Z_> ].
752  //
753  rHz_old(0,0) = rHz(0,0);
754  //
755  // residual update R_ := R_ - alpha * AP
756  //
757  if( prevUdim_ + iter_ <= savedBlocks_ ){
758  MVT::MvAddMv( one, *R_, -alpha(0,0), *AP, *R_ );
759  AP = Teuchos::null;
760  }else{
761  MVT::MvAddMv( one, *R_, -alpha(0,0), *AP_, *R_ );
762  }
763  //
764  // update beta := [ new <R_, Z_> ] / [ old <R_, Z_> ] and the search direction p.
765  //
766  if ( lp_->getLeftPrec() != Teuchos::null ) {
767  lp_->applyLeftPrec( *R_, *Z_ );
768  } else {
769  Z_ = R_;
770  }
771  //
772  MVT::MvTransMv( one, *R_, *Z_, rHz );
773  //
774  beta(0,0) = rHz(0,0) / rHz_old(0,0);
775  //
776  if( curDim_ < savedBlocks_ ){
777  curDim_++; // update basis dim
778  curind[0] = curDim_ - 1;
779  Teuchos::RCP<MV> Pnext = MVT::CloneViewNonConst(*U_,curind);
780  MVT::MvAddMv( one, *Z_, beta(0,0), *P, *Pnext );
781  if( prevUdim_ ){ // Deflate seed space
782  MVT::MvTransMv( one, *Cprev, *Z_, CZ );
783  MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *Pnext ); // Pnext -= U*(C'Z)
784  if( debug ){
785  std::cout << " Check CZ " << std::endl;
786  MVT::MvTransMv( one, *Cprev, *Pnext, CZ );
787  CZ.print( std::cout );
788  }
789  }
790  P = Teuchos::null;
791  if( curDim_ == savedBlocks_ ){
792  std::vector<int> zero_index(1);
793  zero_index[0] = 0;
794  MVT::SetBlock( *Pnext, zero_index, *P_ );
795  }
796  Pnext = Teuchos::null;
797  }else{
798  MVT::MvAddMv( one, *Z_, beta(0,0), *P_, *P_ );
799  if( prevUdim_ ){ // Deflate seed space
800  MVT::MvTransMv( one, *Cprev, *Z_, CZ );
801  MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *P_ ); // P_ -= U*(C'Z)
802 
803  if( debug ){
804  std::cout << " Check CZ " << std::endl;
805  MVT::MvTransMv( one, *Cprev, *P_, CZ );
806  CZ.print( std::cout );
807  }
808  }
809  }
810  // CGB: 5/26/2010
811  // this RCP<const MV> P was previously a variable outside the loop. however, it didn't appear to be see any use between
812  // loop iterations. therefore, I moved it inside to avoid scoping errors with previously used variables named P.
813  // 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
814  // same for AP
815  TEUCHOS_TEST_FOR_EXCEPTION( AP != Teuchos::null || P != Teuchos::null, std::logic_error, "Loop recurrence violated. Please contact Belos team.");
816  } // end coupled two-term recursion
817  if( prevUdim_ + iter_ < savedBlocks_ ) --curDim_; // discard negligible search direction
818  }
819 
820 } // end Belos namespace
821 
822 #endif /* BELOS_PCPG_ITER_HPP */
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 which augments the basic interface for a conjugate gradient linear solver ite...
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.
virtual std::ostream & print(std::ostream &os) const
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.
Traits class which defines basic operations on multivectors.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > D
The current Hessenberg matrix.
CGPositiveDefiniteFailure is thrown when the the CG &#39;alpha = p^H*A*P&#39; value is less than zero...
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.
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.
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.
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.
CGIterationInitFailure is thrown when the CGIteration object is unable to generate an initial iterate...
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...
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?.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Teuchos::ScalarTraits< ScalarType > SCT

Generated on Thu Mar 28 2024 09:24:29 for Belos by doxygen 1.8.5