Belos Version of the Day
Loading...
Searching...
No Matches
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
55#include "BelosStatusTest.hpp"
58#include "BelosCGIteration.hpp"
59
60#include "Teuchos_SerialDenseMatrix.hpp"
61#include "Teuchos_SerialDenseVector.hpp"
62#include "Teuchos_ScalarTraits.hpp"
63#include "Teuchos_ParameterList.hpp"
64#include "Teuchos_TimeMonitor.hpp"
65
77namespace Belos {
78
80
81
86 template <class ScalarType, class MV>
93 int curDim;
96
98 Teuchos::RCP<MV> R;
99
101 Teuchos::RCP<MV> Z;
102
104 Teuchos::RCP<MV> P;
105
107 Teuchos::RCP<MV> AP;
108
110 Teuchos::RCP<MV> U;
111
113 Teuchos::RCP<MV> C;
114
119 Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > D;
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 //
142 typedef Teuchos::ScalarTraits<ScalarType> SCT;
143 typedef typename SCT::magnitudeType MagnitudeType;
144
146
147
155 PCPGIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
156 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
157 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
158 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
159 Teuchos::ParameterList &params );
160
162 virtual ~PCPGIter() {};
164
165
167
168
187 void iterate();
188
208
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 //
318 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
319 const Teuchos::RCP<OutputManager<ScalarType> > om_;
320 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
321 const Teuchos::RCP<OrthoManager<ScalarType,MV> > ortho_;
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
379 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > D_;
380 };
381
383 // Constructor.
384 template<class ScalarType, class MV, class OP>
386 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
387 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
388 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
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>
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) {
518 D_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>() );
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.
637 Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
638 Teuchos::SerialDenseMatrix<int,ScalarType> pAp( 1, 1 );
639 Teuchos::SerialDenseMatrix<int,ScalarType> beta( 1, 1 );
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;
647 Teuchos::RCP<const MV> Uprev;
648 Teuchos::RCP<const MV> Cprev;
649 Teuchos::SerialDenseMatrix<int,ScalarType> CZ;
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 ) {
704 Teuchos::RCP<const MV> P;
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
733 TEUCHOS_TEST_FOR_EXCEPTION( pAp(0,0) <= zero, CGPositiveDefiniteFailure,
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
740 TEUCHOS_TEST_FOR_EXCEPTION( alpha(0,0) <= zero, CGPositiveDefiniteFailure,
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 */
Pure virtual base class which augments the basic interface for a conjugate gradient linear solver ite...
Belos header file which uses auto-configuration information to include necessary C++ headers.
Class which describes the linear problem to be solved by the iterative solver.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Declaration of basic traits for the multivector type.
Class which defines basic traits for the operator type.
Class which manages the output and verbosity of the Belos solvers.
Pure virtual base class for defining the status testing capabilities of Belos.
Collection of types and exceptions used within the Belos solvers.
CGIterationInitFailure is thrown when the CGIteration object is unable to generate an initial iterate...
CGPositiveDefiniteFailure is thrown when the the CG 'alpha = p^H*A*P' value is less than zero,...
A linear system to solve, and its associated information.
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
Traits class which defines basic operations on multivectors.
Class which defines basic traits for the operator type.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
This class implements the PCPG iteration, where a single-std::vector Krylov subspace is constructed....
PCPGIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
MultiVecTraits< ScalarType, MV > MVT
int getNumIters() const
Get the current iteration count.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system solution?.
void initialize()
Initialize the solver with the initial vectors from the linear problem. An exception is thrown if ini...
void setBlockSize(int blockSize)
Get the blocksize to be used by the iterative solver in solving this linear problem.
void resetNumIters(int iter=0)
Reset the iteration count.
virtual ~PCPGIter()
Destructor.
int getNumRecycledBlocks() const
Get the maximum number of recycled blocks used by the iterative solver in solving this linear problem...
void resetState()
tell the Iterator to "reset" itself; delete and rebuild the seed space.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
int getCurSubspaceDim() const
Get the current dimension of the whole seed subspace.
void setSize(int savedBlocks)
Set the maximum number of saved or recycled blocks used by the iterative solver.
SCT::magnitudeType MagnitudeType
OperatorTraits< ScalarType, MV, OP > OPT
int getBlockSize() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem.
int getPrevSubspaceDim() const
Get the dimension of the search subspace used to solve the current solution to the linear problem.
Teuchos::ScalarTraits< ScalarType > SCT
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.
void iterate()
PCPGIter iterates CG until the status test either requests a stop or detects an error....
bool isInitialized()
States whether the solver has been initialized or not.
A pure virtual class for defining the status tests for the Belos iterative solvers.
Structure to contain pointers to PCPGIter state variables.
Teuchos::RCP< MV > AP
The matrix A applied to current decent direction std::vector.
int prevUdim
Number of block columns in matrices C and U before current iteration.
Teuchos::RCP< MV > Z
The current preconditioned residual.
Teuchos::RCP< MV > U
The recycled subspace.
Teuchos::RCP< MV > C
C = AU, U spans recycled subspace.
int curDim
The current dimension of the reduction.
Teuchos::RCP< MV > P
The current decent direction std::vector.
Teuchos::RCP< MV > R
The current residual.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > D
The current Hessenberg matrix.

Generated for Belos by doxygen 1.10.0