541 TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
542 "Belos::PCPGIter::initialize(): Cannot initialize state storage!");
546 std::string errstr(
"Belos::PCPGIter::initialize(): Specified multivectors must have a consistent length and width.");
548 if (newstate.
R != Teuchos::null){
551 if (newstate.
U == Teuchos::null){
557 prevUdim_ = newstate.
curDim;
558 if (newstate.
C == Teuchos::null){
559 std::vector<int> index(prevUdim_);
560 for (
int i=0; i< prevUdim_; ++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 );
567 curDim_ = prevUdim_ ;
571 if (!stateStorageInitialized_)
578 newstate.
curDim = curDim_;
582 std::vector<int> zero_index(1);
584 if ( lp_->getLeftPrec() != Teuchos::null ) {
585 lp_->applyLeftPrec( *R_, *Z_ );
586 MVT::SetBlock( *Z_, zero_index , *P_ );
589 MVT::SetBlock( *R_, zero_index, *P_ );
592 std::vector<int> nextind(1);
593 nextind[0] = curDim_;
595 MVT::SetBlock( *P_, nextind, *newstate.
U );
598 newstate.
curDim = curDim_;
600 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
U) != savedBlocks_ ,
601 std::invalid_argument, errstr );
602 if (newstate.
U != U_) {
606 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.
C) != savedBlocks_ ,
607 std::invalid_argument, errstr );
608 if (newstate.
C != C_) {
614 TEUCHOS_TEST_FOR_EXCEPTION(newstate.
R == Teuchos::null,std::invalid_argument,
615 "Belos::PCPGIter::initialize(): PCPGStateIterState does not have initial kernel R_0.");
631 if (initialized_ ==
false) {
634 const bool debug =
false;
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 );
643 std::cout <<
" Iterate Warning: begin from nonzero iter_ ?" << std::endl;
646 std::vector<int> prevInd;
647 Teuchos::RCP<const MV> Uprev;
648 Teuchos::RCP<const MV> Cprev;
649 Teuchos::SerialDenseMatrix<int,ScalarType> CZ;
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);
660 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
664 "Belos::PCPGIter::iterate(): current linear system has more than one std::vector!" );
668 "Belos::PCPGIter::iterate(): mistake in initialization !" );
671 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
672 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
675 std::vector<int> curind(1);
676 std::vector<ScalarType> rnorm(MVT::GetNumberVecs(*cur_soln_vec));
679 curind[0] = curDim_ - 1;
680 P = MVT::CloneViewNonConst(*U_,curind);
681 MVT::MvTransMv( one, *Cprev, *P, CZ );
682 MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *P );
685 MVT::MvTransMv( one, *Cprev, *P, CZ );
686 std::cout <<
" Input CZ post ortho " << std::endl;
687 CZ.print( std::cout );
689 if( curDim_ == savedBlocks_ ){
690 std::vector<int> zero_index(1);
692 MVT::SetBlock( *P, zero_index, *P_ );
698 MVT::MvTransMv( one, *R_, *Z_, rHz );
703 while (stest_->checkStatus(
this) !=
Passed ) {
704 Teuchos::RCP<const MV> P;
708 curind[0] = curDim_ - 1;
710 MVT::MvNorm(*R_, rnorm);
711 std::cout << iter_ <<
" " << curDim_ <<
" " << rnorm[0] << std::endl;
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 );
719 if( prevUdim_ + iter_ == savedBlocks_ ){
720 AP = MVT::CloneViewNonConst(*C_,curind);
721 lp_->applyOp( *P_, *AP );
722 MVT::MvTransMv( one, *P_, *AP, pAp );
724 lp_->applyOp( *P_, *AP_ );
725 MVT::MvTransMv( one, *P_, *AP_, pAp );
729 if( keepDiagonal_ && prevUdim_ + iter_ <= savedBlocks_ )
730 (*D_)(iter_ -1 ,iter_ -1 ) = pAp(0,0);
734 "Belos::PCPGIter::iterate(): non-positive value for p^H*A*p encountered!" );
737 alpha(0,0) = rHz(0,0) / pAp(0,0);
741 "Belos::PCPGIter::iterate(): non-positive value for alpha encountered!" );
744 if( curDim_ < savedBlocks_ ){
745 MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P, *cur_soln_vec );
747 MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P_, *cur_soln_vec );
753 rHz_old(0,0) = rHz(0,0);
757 if( prevUdim_ + iter_ <= savedBlocks_ ){
758 MVT::MvAddMv( one, *R_, -alpha(0,0), *AP, *R_ );
761 MVT::MvAddMv( one, *R_, -alpha(0,0), *AP_, *R_ );
766 if ( lp_->getLeftPrec() != Teuchos::null ) {
767 lp_->applyLeftPrec( *R_, *Z_ );
772 MVT::MvTransMv( one, *R_, *Z_, rHz );
774 beta(0,0) = rHz(0,0) / rHz_old(0,0);
776 if( curDim_ < savedBlocks_ ){
778 curind[0] = curDim_ - 1;
779 Teuchos::RCP<MV> Pnext = MVT::CloneViewNonConst(*U_,curind);
780 MVT::MvAddMv( one, *Z_, beta(0,0), *P, *Pnext );
782 MVT::MvTransMv( one, *Cprev, *Z_, CZ );
783 MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *Pnext );
785 std::cout <<
" Check CZ " << std::endl;
786 MVT::MvTransMv( one, *Cprev, *Pnext, CZ );
787 CZ.print( std::cout );
791 if( curDim_ == savedBlocks_ ){
792 std::vector<int> zero_index(1);
794 MVT::SetBlock( *Pnext, zero_index, *P_ );
796 Pnext = Teuchos::null;
798 MVT::MvAddMv( one, *Z_, beta(0,0), *P_, *P_ );
800 MVT::MvTransMv( one, *Cprev, *Z_, CZ );
801 MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *P_ );
804 std::cout <<
" Check CZ " << std::endl;
805 MVT::MvTransMv( one, *Cprev, *P_, CZ );
806 CZ.print( std::cout );
815 TEUCHOS_TEST_FOR_EXCEPTION( AP != Teuchos::null || P != Teuchos::null, std::logic_error,
"Loop recurrence violated. Please contact Belos team.");
817 if( prevUdim_ + iter_ < savedBlocks_ ) --curDim_;