42 #ifndef BELOS_RCG_SOLMGR_HPP 43 #define BELOS_RCG_SOLMGR_HPP 64 #include "Teuchos_BLAS.hpp" 65 #include "Teuchos_LAPACK.hpp" 66 #include "Teuchos_as.hpp" 67 #ifdef BELOS_TEUCHOS_TIME_MONITOR 68 #include "Teuchos_TimeMonitor.hpp" 150 template<
class ScalarType,
class MV,
class OP,
151 const bool supportsScalarType =
153 ! Teuchos::ScalarTraits<ScalarType>::isComplex>
156 Belos::Details::LapackSupportsScalar<ScalarType>::value &&
157 ! Teuchos::ScalarTraits<ScalarType>::isComplex>
159 static const bool scalarTypeIsSupported =
161 ! Teuchos::ScalarTraits<ScalarType>::isComplex;
170 const Teuchos::RCP<Teuchos::ParameterList> &pl) :
176 Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const override {
184 template<
class ScalarType,
class MV,
class OP>
190 typedef Teuchos::ScalarTraits<ScalarType> SCT;
191 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
192 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
228 const Teuchos::RCP<Teuchos::ParameterList> &pl );
234 Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const override {
247 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters()
const override;
257 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
258 return Teuchos::tuple(timerSolve_);
286 void setParameters(
const Teuchos::RCP<Teuchos::ParameterList> ¶ms )
override;
299 if ((type &
Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem();
332 std::string description()
const override;
343 void getHarmonicVecs(
const Teuchos::SerialDenseMatrix<int,ScalarType> &F,
344 const Teuchos::SerialDenseMatrix<int,ScalarType> &G,
345 Teuchos::SerialDenseMatrix<int,ScalarType>& Y);
348 void sort(std::vector<ScalarType>& dlist,
int n, std::vector<int>& iperm);
351 void initializeStateStorage();
354 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
357 Teuchos::RCP<OutputManager<ScalarType> > printer_;
358 Teuchos::RCP<std::ostream> outputStream_;
361 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
362 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
363 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
364 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
367 Teuchos::RCP<Teuchos::ParameterList> params_;
370 static constexpr
int maxIters_default_ = 1000;
371 static constexpr
int blockSize_default_ = 1;
372 static constexpr
int numBlocks_default_ = 25;
373 static constexpr
int recycleBlocks_default_ = 3;
374 static constexpr
bool showMaxResNormOnly_default_ =
false;
377 static constexpr
int outputFreq_default_ = -1;
378 static constexpr
const char * label_default_ =
"Belos";
379 static constexpr std::ostream * outputStream_default_ = &std::cout;
386 MagnitudeType convtol_;
392 MagnitudeType achievedTol_;
400 int numBlocks_, recycleBlocks_;
401 bool showMaxResNormOnly_;
402 int verbosity_, outputStyle_, outputFreq_;
411 Teuchos::RCP<MV> Ap_;
426 Teuchos::RCP<MV> U_, AU_;
429 Teuchos::RCP<MV> U1_;
432 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Alpha_;
433 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Beta_;
434 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > D_;
437 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Delta_;
440 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > UTAU_;
443 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > LUUTAU_;
446 Teuchos::RCP<std::vector<int> > ipiv_;
449 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AUTAU_;
452 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > rTz_old_;
455 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > F_,G_,Y_;
458 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > L2_,DeltaL2_,AU1TUDeltaL2_;
459 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AU1TAU1_, AU1TU1_, AU1TAP_;
460 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > FY_,GY_;
461 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > APTAP_;
462 Teuchos::RCP<MV> U1Y1_, PY2_;
463 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AUTAP_, AU1TU_;
469 Teuchos::RCP<Teuchos::Time> timerSolve_;
477 template<
class ScalarType,
class MV,
class OP>
486 template<
class ScalarType,
class MV,
class OP>
489 const Teuchos::RCP<Teuchos::ParameterList> &pl ) :
495 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument,
"Problem not given to solver manager.");
498 if ( !is_null(pl) ) {
504 template<
class ScalarType,
class MV,
class OP>
507 outputStream_ = Teuchos::rcp(outputStream_default_,
false);
509 maxIters_ = maxIters_default_;
510 numBlocks_ = numBlocks_default_;
511 recycleBlocks_ = recycleBlocks_default_;
512 verbosity_ = verbosity_default_;
513 outputStyle_= outputStyle_default_;
514 outputFreq_= outputFreq_default_;
515 showMaxResNormOnly_ = showMaxResNormOnly_default_;
516 label_ = label_default_;
527 Alpha_ = Teuchos::null;
528 Beta_ = Teuchos::null;
530 Delta_ = Teuchos::null;
531 UTAU_ = Teuchos::null;
532 LUUTAU_ = Teuchos::null;
533 ipiv_ = Teuchos::null;
534 AUTAU_ = Teuchos::null;
535 rTz_old_ = Teuchos::null;
540 DeltaL2_ = Teuchos::null;
541 AU1TUDeltaL2_ = Teuchos::null;
542 AU1TAU1_ = Teuchos::null;
543 AU1TU1_ = Teuchos::null;
544 AU1TAP_ = Teuchos::null;
547 APTAP_ = Teuchos::null;
548 U1Y1_ = Teuchos::null;
549 PY2_ = Teuchos::null;
550 AUTAP_ = Teuchos::null;
551 AU1TU_ = Teuchos::null;
555 template<
class ScalarType,
class MV,
class OP>
559 if (params_ == Teuchos::null) {
560 params_ = Teuchos::rcp(
new Teuchos::ParameterList(*getValidParameters()) );
563 params->validateParameters(*getValidParameters());
567 if (params->isParameter(
"Maximum Iterations")) {
568 maxIters_ = params->get(
"Maximum Iterations",maxIters_default_);
571 params_->set(
"Maximum Iterations", maxIters_);
572 if (maxIterTest_!=Teuchos::null)
573 maxIterTest_->setMaxIters( maxIters_ );
577 if (params->isParameter(
"Num Blocks")) {
578 numBlocks_ = params->get(
"Num Blocks",numBlocks_default_);
579 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
580 "Belos::RCGSolMgr: \"Num Blocks\" must be strictly positive.");
583 params_->set(
"Num Blocks", numBlocks_);
587 if (params->isParameter(
"Num Recycled Blocks")) {
588 recycleBlocks_ = params->get(
"Num Recycled Blocks",recycleBlocks_default_);
589 TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks_ <= 0, std::invalid_argument,
590 "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be strictly positive.");
592 TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks_ >= numBlocks_, std::invalid_argument,
593 "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be less than \"Num Blocks\".");
596 params_->set(
"Num Recycled Blocks", recycleBlocks_);
600 if (params->isParameter(
"Timer Label")) {
601 std::string tempLabel = params->get(
"Timer Label", label_default_);
604 if (tempLabel != label_) {
606 params_->set(
"Timer Label", label_);
607 std::string solveLabel = label_ +
": RCGSolMgr total solve time";
608 #ifdef BELOS_TEUCHOS_TIME_MONITOR 609 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
615 if (params->isParameter(
"Verbosity")) {
616 if (Teuchos::isParameterType<int>(*params,
"Verbosity")) {
617 verbosity_ = params->get(
"Verbosity", verbosity_default_);
619 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,
"Verbosity");
623 params_->set(
"Verbosity", verbosity_);
624 if (printer_ != Teuchos::null)
625 printer_->setVerbosity(verbosity_);
629 if (params->isParameter(
"Output Style")) {
630 if (Teuchos::isParameterType<int>(*params,
"Output Style")) {
631 outputStyle_ = params->get(
"Output Style", outputStyle_default_);
633 outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,
"Output Style");
637 params_->set(
"Output Style", outputStyle_);
638 outputTest_ = Teuchos::null;
642 if (params->isParameter(
"Output Stream")) {
643 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,
"Output Stream");
646 params_->set(
"Output Stream", outputStream_);
647 if (printer_ != Teuchos::null)
648 printer_->setOStream( outputStream_ );
653 if (params->isParameter(
"Output Frequency")) {
654 outputFreq_ = params->get(
"Output Frequency", outputFreq_default_);
658 params_->set(
"Output Frequency", outputFreq_);
659 if (outputTest_ != Teuchos::null)
660 outputTest_->setOutputFrequency( outputFreq_ );
664 if (printer_ == Teuchos::null) {
673 if (params->isParameter(
"Convergence Tolerance")) {
674 if (params->isType<MagnitudeType> (
"Convergence Tolerance")) {
675 convtol_ = params->get (
"Convergence Tolerance",
683 params_->set(
"Convergence Tolerance", convtol_);
684 if (convTest_ != Teuchos::null)
688 if (params->isParameter(
"Show Maximum Residual Norm Only")) {
689 showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,
"Show Maximum Residual Norm Only");
692 params_->set(
"Show Maximum Residual Norm Only", showMaxResNormOnly_);
693 if (convTest_ != Teuchos::null)
694 convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
700 if (maxIterTest_ == Teuchos::null)
704 if (convTest_ == Teuchos::null)
705 convTest_ = Teuchos::rcp(
new StatusTestResNorm_t( convtol_, 1 ) );
707 if (sTest_ == Teuchos::null)
708 sTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
710 if (outputTest_ == Teuchos::null) {
718 std::string solverDesc =
" Recycling CG ";
719 outputTest_->setSolverDesc( solverDesc );
723 if (timerSolve_ == Teuchos::null) {
724 std::string solveLabel = label_ +
": RCGSolMgr total solve time";
725 #ifdef BELOS_TEUCHOS_TIME_MONITOR 726 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
735 template<
class ScalarType,
class MV,
class OP>
736 Teuchos::RCP<const Teuchos::ParameterList>
739 static Teuchos::RCP<const Teuchos::ParameterList> validPL;
742 if(is_null(validPL)) {
743 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
745 "The relative residual tolerance that needs to be achieved by the\n" 746 "iterative solver in order for the linear system to be declared converged.");
747 pl->set(
"Maximum Iterations", static_cast<int>(maxIters_default_),
748 "The maximum number of block iterations allowed for each\n" 749 "set of RHS solved.");
750 pl->set(
"Block Size", static_cast<int>(blockSize_default_),
751 "Block Size Parameter -- currently must be 1 for RCG");
752 pl->set(
"Num Blocks", static_cast<int>(numBlocks_default_),
753 "The length of a cycle (and this max number of search vectors kept)\n");
754 pl->set(
"Num Recycled Blocks", static_cast<int>(recycleBlocks_default_),
755 "The number of vectors in the recycle subspace.");
756 pl->set(
"Verbosity", static_cast<int>(verbosity_default_),
757 "What type(s) of solver information should be outputted\n" 758 "to the output stream.");
759 pl->set(
"Output Style", static_cast<int>(outputStyle_default_),
760 "What style is used for the solver information outputted\n" 761 "to the output stream.");
762 pl->set(
"Output Frequency", static_cast<int>(outputFreq_default_),
763 "How often convergence information should be outputted\n" 764 "to the output stream.");
765 pl->set(
"Output Stream", Teuchos::rcp(outputStream_default_,
false),
766 "A reference-counted pointer to the output stream where all\n" 767 "solver output is sent.");
768 pl->set(
"Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
769 "When convergence information is printed, only show the maximum\n" 770 "relative residual norm when the block size is greater than one.");
771 pl->set(
"Timer Label", static_cast<const char *>(label_default_),
772 "The string to use as a prefix for the timer labels.");
779 template<
class ScalarType,
class MV,
class OP>
783 Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
784 if (rhsMV == Teuchos::null) {
791 TEUCHOS_TEST_FOR_EXCEPTION(static_cast<ptrdiff_t>(numBlocks_) > MVT::GetGlobalLength(*rhsMV),std::invalid_argument,
792 "Belos::RCGSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!");
795 if (P_ == Teuchos::null) {
796 P_ = MVT::Clone( *rhsMV, numBlocks_+2 );
800 if (MVT::GetNumberVecs(*P_) < numBlocks_+2) {
801 Teuchos::RCP<const MV> tmp = P_;
802 P_ = MVT::Clone( *tmp, numBlocks_+2 );
807 if (Ap_ == Teuchos::null)
808 Ap_ = MVT::Clone( *rhsMV, 1 );
811 if (r_ == Teuchos::null)
812 r_ = MVT::Clone( *rhsMV, 1 );
815 if (z_ == Teuchos::null)
816 z_ = MVT::Clone( *rhsMV, 1 );
819 if (U_ == Teuchos::null) {
820 U_ = MVT::Clone( *rhsMV, recycleBlocks_ );
824 if (MVT::GetNumberVecs(*U_) < recycleBlocks_) {
825 Teuchos::RCP<const MV> tmp = U_;
826 U_ = MVT::Clone( *tmp, recycleBlocks_ );
831 if (AU_ == Teuchos::null) {
832 AU_ = MVT::Clone( *rhsMV, recycleBlocks_ );
836 if (MVT::GetNumberVecs(*AU_) < recycleBlocks_) {
837 Teuchos::RCP<const MV> tmp = AU_;
838 AU_ = MVT::Clone( *tmp, recycleBlocks_ );
843 if (U1_ == Teuchos::null) {
844 U1_ = MVT::Clone( *rhsMV, recycleBlocks_ );
848 if (MVT::GetNumberVecs(*U1_) < recycleBlocks_) {
849 Teuchos::RCP<const MV> tmp = U1_;
850 U1_ = MVT::Clone( *tmp, recycleBlocks_ );
855 if (Alpha_ == Teuchos::null)
856 Alpha_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_, 1 ) );
858 if ( (Alpha_->numRows() != numBlocks_) || (Alpha_->numCols() != 1) )
859 Alpha_->reshape( numBlocks_, 1 );
863 if (Beta_ == Teuchos::null)
864 Beta_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + 1, 1 ) );
866 if ( (Beta_->numRows() != (numBlocks_+1)) || (Beta_->numCols() != 1) )
867 Beta_->reshape( numBlocks_ + 1, 1 );
871 if (D_ == Teuchos::null)
872 D_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ , 1 ) );
874 if ( (D_->numRows() != numBlocks_) || (D_->numCols() != 1) )
875 D_->reshape( numBlocks_, 1 );
879 if (Delta_ == Teuchos::null)
880 Delta_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ + 1 ) );
882 if ( (Delta_->numRows() != recycleBlocks_) || (Delta_->numCols() != (numBlocks_ + 1)) )
883 Delta_->reshape( recycleBlocks_, numBlocks_ + 1 );
887 if (UTAU_ == Teuchos::null)
888 UTAU_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
890 if ( (UTAU_->numRows() != recycleBlocks_) || (UTAU_->numCols() != recycleBlocks_) )
891 UTAU_->reshape( recycleBlocks_, recycleBlocks_ );
895 if (LUUTAU_ == Teuchos::null)
896 LUUTAU_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
898 if ( (LUUTAU_->numRows() != recycleBlocks_) || (LUUTAU_->numCols() != recycleBlocks_) )
899 LUUTAU_->reshape( recycleBlocks_, recycleBlocks_ );
903 if (ipiv_ == Teuchos::null)
904 ipiv_ = Teuchos::rcp(
new std::vector<int>(recycleBlocks_) );
906 if ( (
int)ipiv_->size() != recycleBlocks_ )
907 ipiv_->resize(recycleBlocks_);
911 if (AUTAU_ == Teuchos::null)
912 AUTAU_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
914 if ( (AUTAU_->numRows() != recycleBlocks_) || (AUTAU_->numCols() != recycleBlocks_) )
915 AUTAU_->reshape( recycleBlocks_, recycleBlocks_ );
919 if (rTz_old_ == Teuchos::null)
920 rTz_old_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( 1, 1 ) );
922 if ( (rTz_old_->numRows() != 1) || (rTz_old_->numCols() != 1) )
923 rTz_old_->reshape( 1, 1 );
927 if (F_ == Teuchos::null)
928 F_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ) );
930 if ( (F_->numRows() != (numBlocks_+recycleBlocks_)) || (F_->numCols() != numBlocks_+recycleBlocks_) )
931 F_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
935 if (G_ == Teuchos::null)
936 G_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ) );
938 if ( (G_->numRows() != (numBlocks_+recycleBlocks_)) || (G_->numCols() != numBlocks_+recycleBlocks_) )
939 G_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
943 if (Y_ == Teuchos::null)
944 Y_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, recycleBlocks_ ) );
946 if ( (Y_->numRows() != (numBlocks_+recycleBlocks_)) || (Y_->numCols() != numBlocks_+recycleBlocks_) )
947 Y_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
951 if (L2_ == Teuchos::null)
952 L2_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+1, numBlocks_ ) );
954 if ( (L2_->numRows() != (numBlocks_+1)) || (L2_->numCols() != numBlocks_) )
955 L2_->reshape( numBlocks_+1, numBlocks_ );
959 if (DeltaL2_ == Teuchos::null)
960 DeltaL2_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
962 if ( (DeltaL2_->numRows() != (recycleBlocks_)) || (DeltaL2_->numCols() != (numBlocks_) ) )
963 DeltaL2_->reshape( recycleBlocks_, numBlocks_ );
967 if (AU1TUDeltaL2_ == Teuchos::null)
968 AU1TUDeltaL2_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
970 if ( (AU1TUDeltaL2_->numRows() != (recycleBlocks_)) || (AU1TUDeltaL2_->numCols() != (numBlocks_) ) )
971 AU1TUDeltaL2_->reshape( recycleBlocks_, numBlocks_ );
975 if (AU1TAU1_ == Teuchos::null)
976 AU1TAU1_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
978 if ( (AU1TAU1_->numRows() != (recycleBlocks_)) || (AU1TAU1_->numCols() != (recycleBlocks_) ) )
979 AU1TAU1_->reshape( recycleBlocks_, recycleBlocks_ );
983 if (GY_ == Teuchos::null)
984 GY_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + recycleBlocks_, recycleBlocks_ ) );
986 if ( (GY_->numRows() != (numBlocks_ + recycleBlocks_)) || (GY_->numCols() != (recycleBlocks_) ) )
987 GY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ );
991 if (AU1TU1_ == Teuchos::null)
992 AU1TU1_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
994 if ( (AU1TU1_->numRows() != (recycleBlocks_)) || (AU1TU1_->numCols() != (recycleBlocks_) ) )
995 AU1TU1_->reshape( recycleBlocks_, recycleBlocks_ );
999 if (FY_ == Teuchos::null)
1000 FY_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + recycleBlocks_, recycleBlocks_ ) );
1002 if ( (FY_->numRows() != (numBlocks_ + recycleBlocks_)) || (FY_->numCols() != (recycleBlocks_) ) )
1003 FY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ );
1007 if (AU1TAP_ == Teuchos::null)
1008 AU1TAP_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
1010 if ( (AU1TAP_->numRows() != (recycleBlocks_)) || (AU1TAP_->numCols() != (numBlocks_) ) )
1011 AU1TAP_->reshape( recycleBlocks_, numBlocks_ );
1015 if (APTAP_ == Teuchos::null)
1016 APTAP_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_, numBlocks_ ) );
1018 if ( (APTAP_->numRows() != (numBlocks_)) || (APTAP_->numCols() != (numBlocks_) ) )
1019 APTAP_->reshape( numBlocks_, numBlocks_ );
1023 if (U1Y1_ == Teuchos::null) {
1024 U1Y1_ = MVT::Clone( *rhsMV, recycleBlocks_ );
1028 if (MVT::GetNumberVecs(*U1Y1_) < recycleBlocks_) {
1029 Teuchos::RCP<const MV> tmp = U1Y1_;
1030 U1Y1_ = MVT::Clone( *tmp, recycleBlocks_ );
1035 if (PY2_ == Teuchos::null) {
1036 PY2_ = MVT::Clone( *rhsMV, recycleBlocks_ );
1040 if (MVT::GetNumberVecs(*PY2_) < recycleBlocks_) {
1041 Teuchos::RCP<const MV> tmp = PY2_;
1042 PY2_ = MVT::Clone( *tmp, recycleBlocks_ );
1047 if (AUTAP_ == Teuchos::null)
1048 AUTAP_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
1050 if ( (AUTAP_->numRows() != (recycleBlocks_)) || (AUTAP_->numCols() != (numBlocks_) ) )
1051 AUTAP_->reshape( recycleBlocks_, numBlocks_ );
1055 if (AU1TU_ == Teuchos::null)
1056 AU1TU_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
1058 if ( (AU1TU_->numRows() != (recycleBlocks_)) || (AU1TU_->numCols() != (recycleBlocks_) ) )
1059 AU1TU_->reshape( recycleBlocks_, recycleBlocks_ );
1066 template<
class ScalarType,
class MV,
class OP>
1069 Teuchos::LAPACK<int,ScalarType> lapack;
1070 std::vector<int> index(recycleBlocks_);
1071 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1072 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1081 setParameters(Teuchos::parameterList(*getValidParameters()));
1085 "Belos::RCGSolMgr::solve(): Linear problem is not a valid object.");
1087 "Belos::RCGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
1088 TEUCHOS_TEST_FOR_EXCEPTION((problem_->getLeftPrec() != Teuchos::null)&&(problem_->getRightPrec() != Teuchos::null),
1090 "Belos::RCGSolMgr::solve(): RCG does not support split preconditioning, only set left or right preconditioner.");
1093 Teuchos::RCP<OP> precObj;
1094 if (problem_->getLeftPrec() != Teuchos::null) {
1095 precObj = Teuchos::rcp_const_cast<OP>(problem_->getLeftPrec());
1097 else if (problem_->getRightPrec() != Teuchos::null) {
1098 precObj = Teuchos::rcp_const_cast<OP>(problem_->getRightPrec());
1102 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1103 std::vector<int> currIdx(1);
1107 problem_->setLSIndex( currIdx );
1110 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1111 if (numBlocks_ > dim) {
1112 numBlocks_ = Teuchos::asSafe<int>(dim);
1113 params_->set(
"Num Blocks", numBlocks_);
1115 "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl <<
1116 " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl;
1120 initializeStateStorage();
1123 Teuchos::ParameterList plist;
1124 plist.set(
"Num Blocks",numBlocks_);
1125 plist.set(
"Recycled Blocks",recycleBlocks_);
1128 outputTest_->reset();
1131 bool isConverged =
true;
1135 index.resize(recycleBlocks_);
1136 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1137 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1138 index.resize(recycleBlocks_);
1139 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1140 Teuchos::RCP<MV> AUtmp = MVT::CloneViewNonConst( *AU_, index );
1142 problem_->applyOp( *Utmp, *AUtmp );
1144 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1145 MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp );
1147 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
1148 if ( precObj != Teuchos::null ) {
1149 index.resize(recycleBlocks_);
1150 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1151 index.resize(recycleBlocks_);
1152 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1153 Teuchos::RCP<MV> PCAU = MVT::CloneViewNonConst( *U1_, index );
1154 OPT::Apply( *precObj, *AUtmp, *PCAU );
1155 MVT::MvTransMv( one, *AUtmp, *PCAU, AUTAUtmp );
1157 MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp );
1164 Teuchos::RCP<RCGIter<ScalarType,MV,OP> > rcg_iter;
1169 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1170 Teuchos::TimeMonitor slvtimer(*timerSolve_);
1173 while ( numRHS2Solve > 0 ) {
1176 if (printer_->isVerbosity(
Debug ) ) {
1177 if (existU_) printer_->print(
Debug,
"Using recycle space generated from previous call to solve()." );
1178 else printer_->print(
Debug,
"No recycle space exists." );
1182 rcg_iter->resetNumIters();
1185 rcg_iter->setSize( recycleBlocks_, numBlocks_ );
1188 outputTest_->resetNumCalls();
1197 problem_->computeCurrResVec( &*r_ );
1202 Teuchos::SerialDenseMatrix<int,ScalarType> Utr(recycleBlocks_,1);
1203 index.resize(recycleBlocks_);
1204 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1205 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1206 MVT::MvTransMv( one, *Utmp, *r_, Utr );
1207 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1208 Teuchos::SerialDenseMatrix<int,ScalarType> LUUTAUtmp( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ );
1209 LUUTAUtmp.assign(UTAUtmp);
1211 lapack.GESV(recycleBlocks_, 1, LUUTAUtmp.values(), LUUTAUtmp.stride(), &(*ipiv_)[0], Utr.values(), Utr.stride(), &info);
1213 "Belos::RCGSolMgr::solve(): LAPACK GESV failed to compute a solution.");
1216 MVT::MvTimesMatAddMv( one, *Utmp, Utr, one, *problem_->getCurrLHSVec() );
1219 index.resize(recycleBlocks_);
1220 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1221 Teuchos::RCP<const MV> AUtmp = MVT::CloneView( *AU_, index );
1222 MVT::MvTimesMatAddMv( -one, *AUtmp, Utr, one, *r_ );
1225 if ( precObj != Teuchos::null ) {
1226 OPT::Apply( *precObj, *r_, *z_ );
1232 MVT::MvTransMv( one, *r_, *z_, *rTz_old_ );
1236 Teuchos::SerialDenseMatrix<int,ScalarType> mu( Teuchos::View, *Delta_, recycleBlocks_, 1);
1237 index.resize(recycleBlocks_);
1238 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1239 Teuchos::RCP<const MV> AUtmp = MVT::CloneView( *AU_, index );
1240 MVT::MvTransMv( one, *AUtmp, *z_, mu );
1241 Teuchos::SerialDenseMatrix<int,ScalarType> LUUTAUtmp( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ );
1244 lapack.GETRS(
TRANS, recycleBlocks_, 1, LUUTAUtmp.values(), LUUTAUtmp.stride(), &(*ipiv_)[0], mu.values(), mu.stride(), &info );
1246 "Belos::RCGSolMgr::solve(): LAPACK GETRS failed to compute a solution.");
1250 Teuchos::RCP<MV> Ptmp = MVT::CloneViewNonConst( *P_, index );
1251 MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp);
1252 MVT::MvTimesMatAddMv( -one, *U_, mu, one, *Ptmp );
1257 Teuchos::RCP<MV> Ptmp = MVT::CloneViewNonConst( *P_, index );
1258 MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp);
1265 index.resize( numBlocks_+1 );
1266 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1267 newstate.
P = MVT::CloneViewNonConst( *P_, index );
1268 index.resize( recycleBlocks_ );
1269 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1270 newstate.
U = MVT::CloneViewNonConst( *U_, index );
1271 index.resize( recycleBlocks_ );
1272 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1273 newstate.
AU = MVT::CloneViewNonConst( *AU_, index );
1274 newstate.
Alpha = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Alpha_, numBlocks_, 1 ) );
1275 newstate.
Beta = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Beta_, numBlocks_, 1 ) );
1276 newstate.
D = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *D_, numBlocks_, 1 ) );
1277 newstate.
Delta = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_, 0, 1 ) );
1278 newstate.
LUUTAU = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ ) );
1284 newstate.
existU = existU_;
1285 newstate.
ipiv = ipiv_;
1288 rcg_iter->initialize(newstate);
1294 rcg_iter->iterate();
1301 if ( convTest_->getStatus() ==
Passed ) {
1310 else if ( maxIterTest_->getStatus() ==
Passed ) {
1312 isConverged =
false;
1320 else if ( rcg_iter->getCurSubspaceDim() == rcg_iter->getMaxSubspaceDim() ) {
1325 if (recycleBlocks_ > 0) {
1329 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, numBlocks_, numBlocks_ );
1330 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, numBlocks_, numBlocks_ );
1331 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1332 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1333 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_, 1 );
1334 Ftmp.putScalar(zero);
1335 Gtmp.putScalar(zero);
1336 for (
int ii=0;ii<numBlocks_;ii++) {
1337 Gtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0));
1339 Gtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1340 Gtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1342 Ftmp(ii,ii) = Dtmp(ii,0);
1346 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, numBlocks_, recycleBlocks_ );
1347 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1350 index.resize( numBlocks_ );
1351 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; }
1352 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1353 index.resize( recycleBlocks_ );
1354 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1355 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1356 MVT::MvTimesMatAddMv( one, *Ptmp, Ytmp, zero, *U1tmp );
1361 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, numBlocks_, recycleBlocks_ );
1362 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1363 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1364 AU1TAU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1367 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, numBlocks_, recycleBlocks_ );
1368 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1369 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1370 AU1TU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1372 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, 1 );
1374 AU1TAPtmp.putScalar(zero);
1376 ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0);
1377 for (
int ii=0; ii<recycleBlocks_; ++ii) {
1378 AU1TAPtmp(ii,0) = Ytmp(numBlocks_-1,ii) * alphatmp;
1386 lastBeta = numBlocks_-1;
1393 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1394 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, numBlocks_ );
1395 AU1TAPtmp.scale(Dtmp(0,0));
1397 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1398 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_+1, 1 );
1399 Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
1400 APTAPtmp.putScalar(zero);
1401 for (
int ii=0; ii<numBlocks_; ii++) {
1402 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0));
1404 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1405 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1410 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1411 Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
1412 Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1413 Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1414 Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1415 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1416 F11.assign(AU1TU1tmp);
1417 F12.putScalar(zero);
1418 F21.putScalar(zero);
1419 F22.putScalar(zero);
1420 for(
int ii=0;ii<numBlocks_;ii++) {
1421 F22(ii,ii) = Dtmp(ii,0);
1425 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1426 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1427 Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
1428 Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1429 Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1430 Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1431 G11.assign(AU1TAU1tmp);
1432 G12.assign(AU1TAPtmp);
1434 for (
int ii=0;ii<recycleBlocks_;++ii)
1435 for (
int jj=0;jj<numBlocks_;++jj)
1436 G21(jj,ii) = G12(ii,jj);
1437 G22.assign(APTAPtmp);
1440 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
1441 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1444 index.resize( numBlocks_ );
1445 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; }
1446 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1447 index.resize( recycleBlocks_ );
1448 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1449 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index );
1450 Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1451 index.resize( recycleBlocks_ );
1452 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1453 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1454 index.resize( recycleBlocks_ );
1455 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1456 Teuchos::RCP<MV> U1Y1tmp = MVT::CloneViewNonConst( *U1Y1_, index );
1457 Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
1458 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1459 MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp );
1460 MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp);
1465 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1466 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1467 AU1TAU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1471 AU1TAPtmp.putScalar(zero);
1472 ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0);
1473 for (
int ii=0; ii<recycleBlocks_; ++ii) {
1474 AU1TAPtmp(ii,0) = Ytmp(numBlocks_+recycleBlocks_-1,ii) * alphatmp;
1478 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1479 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1481 AU1TU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1484 lastp = numBlocks_+1;
1485 lastBeta = numBlocks_;
1491 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1492 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_, 1 );
1493 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1494 Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
1495 APTAPtmp.putScalar(zero);
1496 for (
int ii=0; ii<numBlocks_; ii++) {
1497 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0));
1499 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1500 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1504 Teuchos::SerialDenseMatrix<int,ScalarType> L2tmp( Teuchos::View, *L2_, numBlocks_+1, numBlocks_ );
1505 L2tmp.putScalar(zero);
1506 for(
int ii=0;ii<numBlocks_;ii++) {
1507 L2tmp(ii,ii) = 1./Alphatmp(ii,0);
1508 L2tmp(ii+1,ii) = -1./Alphatmp(ii,0);
1512 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAPtmp( Teuchos::View, *AUTAP_, recycleBlocks_, numBlocks_ );
1513 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1514 Teuchos::SerialDenseMatrix<int,ScalarType> Deltatmp( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_+1 );
1515 Teuchos::SerialDenseMatrix<int,ScalarType> DeltaL2tmp( Teuchos::View, *DeltaL2_, recycleBlocks_, numBlocks_ );
1516 DeltaL2tmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Deltatmp,L2tmp,zero);
1517 AUTAPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,UTAUtmp,DeltaL2tmp,zero);
1520 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1521 Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
1522 Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1523 Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1524 Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1525 F11.assign(UTAUtmp);
1526 F12.putScalar(zero);
1527 F21.putScalar(zero);
1528 F22.putScalar(zero);
1529 for(
int ii=0;ii<numBlocks_;ii++) {
1530 F22(ii,ii) = Dtmp(ii,0);
1534 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1535 Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
1536 Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1537 Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1538 Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1539 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
1540 G11.assign(AUTAUtmp);
1541 G12.assign(AUTAPtmp);
1543 for (
int ii=0;ii<recycleBlocks_;++ii)
1544 for (
int jj=0;jj<numBlocks_;++jj)
1545 G21(jj,ii) = G12(ii,jj);
1546 G22.assign(APTAPtmp);
1549 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
1550 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1553 index.resize( recycleBlocks_ );
1554 for (
int ii=0; ii<(recycleBlocks_); ++ii) { index[ii] = ii; }
1555 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1556 index.resize( numBlocks_ );
1557 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; }
1558 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1559 index.resize( recycleBlocks_ );
1560 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1561 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index );
1562 Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1563 index.resize( recycleBlocks_ );
1564 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1565 Teuchos::RCP<MV> UY1tmp = MVT::CloneViewNonConst( *U1Y1_, index );
1566 Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
1567 index.resize( recycleBlocks_ );
1568 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1569 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1570 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1571 MVT::MvTimesMatAddMv( one, *Utmp, Y1, zero, *UY1tmp );
1572 MVT::MvAddMv(one,*UY1tmp, one, *PY2tmp, *U1tmp);
1577 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1578 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1579 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1580 AU1TAU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1583 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1584 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1585 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1586 AU1TU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1589 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUtmp( Teuchos::View, *AU1TU_, recycleBlocks_, recycleBlocks_ );
1590 AU1TUtmp.assign(UTAUtmp);
1593 dold = Dtmp(numBlocks_-1,0);
1600 lastBeta = numBlocks_-1;
1603 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1604 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_+1, 1 );
1605 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1606 Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
1607 for (
int ii=0; ii<numBlocks_; ii++) {
1608 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0));
1610 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1611 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1615 Teuchos::SerialDenseMatrix<int,ScalarType> L2tmp( Teuchos::View, *L2_, numBlocks_+1, numBlocks_ );
1616 for(
int ii=0;ii<numBlocks_;ii++) {
1617 L2tmp(ii,ii) = 1./Alphatmp(ii,0);
1618 L2tmp(ii+1,ii) = -1./Alphatmp(ii,0);
1623 Teuchos::SerialDenseMatrix<int,ScalarType> DeltaL2( Teuchos::View, *DeltaL2_, recycleBlocks_, numBlocks_ );
1624 Teuchos::SerialDenseMatrix<int,ScalarType> Deltatmp( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_+1 );
1625 DeltaL2.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Deltatmp,L2tmp,zero);
1626 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUDeltaL2( Teuchos::View, *AU1TUDeltaL2_, recycleBlocks_, numBlocks_ );
1627 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUtmp( Teuchos::View, *AU1TU_, recycleBlocks_, recycleBlocks_ );
1628 AU1TUDeltaL2.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,AU1TUtmp,DeltaL2,zero);
1629 Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
1630 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, numBlocks_ );
1631 AU1TAPtmp.putScalar(zero);
1632 AU1TAPtmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Y1,AU1TUDeltaL2,zero);
1633 Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1634 ScalarType val = dold * (-Betatmp(0,0)/Alphatmp(0,0));
1635 for(
int ii=0;ii<recycleBlocks_;ii++) {
1636 AU1TAPtmp(ii,0) += Y2(numBlocks_-1,ii)*val;
1640 Teuchos::SerialDenseMatrix<int,ScalarType> Y1TAU1TU( Teuchos::View, *GY_, recycleBlocks_, recycleBlocks_ );
1641 Y1TAU1TU.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Y1,AU1TUtmp,zero);
1642 AU1TUtmp.assign(Y1TAU1TU);
1645 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1646 Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
1647 Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1648 Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1649 Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1650 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1651 F11.assign(AU1TU1tmp);
1652 for(
int ii=0;ii<numBlocks_;ii++) {
1653 F22(ii,ii) = Dtmp(ii,0);
1657 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1658 Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
1659 Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1660 Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1661 Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1662 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1663 G11.assign(AU1TAU1tmp);
1664 G12.assign(AU1TAPtmp);
1666 for (
int ii=0;ii<recycleBlocks_;++ii)
1667 for (
int jj=0;jj<numBlocks_;++jj)
1668 G21(jj,ii) = G12(ii,jj);
1669 G22.assign(APTAPtmp);
1672 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
1673 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1676 index.resize( numBlocks_ );
1677 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; }
1678 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1679 index.resize( recycleBlocks_ );
1680 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1681 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index );
1682 index.resize( recycleBlocks_ );
1683 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1684 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1685 index.resize( recycleBlocks_ );
1686 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1687 Teuchos::RCP<MV> U1Y1tmp = MVT::CloneViewNonConst( *U1Y1_, index );
1688 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1689 MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp );
1690 MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp);
1695 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1696 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1697 AU1TAU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1700 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1701 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1702 AU1TU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1705 dold = Dtmp(numBlocks_-1,0);
1708 lastp = numBlocks_+1;
1709 lastBeta = numBlocks_;
1719 index[0] = lastp-1; index[1] = lastp;
1720 Teuchos::RCP<const MV> Ptmp2 = MVT::CloneView( *P_, index );
1721 index[0] = 0; index[1] = 1;
1722 MVT::SetBlock(*Ptmp2,index,*P_);
1725 (*Beta_)(0,0) = (*Beta_)(lastBeta,0);
1729 Teuchos::SerialDenseMatrix<int,ScalarType> mu1( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, 0 );
1730 Teuchos::SerialDenseMatrix<int,ScalarType> mu2( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, numBlocks_ );
1735 newstate.
P = Teuchos::null;
1736 index.resize( numBlocks_+1 );
1737 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii+1; }
1738 newstate.
P = MVT::CloneViewNonConst( *P_, index );
1740 newstate.
Beta = Teuchos::null;
1741 newstate.
Beta = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Beta_, numBlocks_, 1, 1, 0 ) );
1743 newstate.
Delta = Teuchos::null;
1744 newstate.
Delta = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_, 0, 1 ) );
1749 rcg_iter->initialize(newstate);
1762 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
1763 "Belos::RCGSolMgr::solve(): Invalid return from RCGIter::iterate().");
1766 catch (
const std::exception &e) {
1767 printer_->stream(
Errors) <<
"Error! Caught std::exception in RCGIter::iterate() at iteration " 1768 << rcg_iter->getNumIters() << std::endl
1769 << e.what() << std::endl;
1775 problem_->setCurrLS();
1779 if ( numRHS2Solve > 0 ) {
1782 problem_->setLSIndex( currIdx );
1785 currIdx.resize( numRHS2Solve );
1791 index.resize(recycleBlocks_);
1792 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1793 MVT::SetBlock(*U1_,index,*U_);
1796 if (numRHS2Solve > 0) {
1798 newstate.
P = Teuchos::null;
1799 newstate.
Ap = Teuchos::null;
1800 newstate.
r = Teuchos::null;
1801 newstate.
z = Teuchos::null;
1802 newstate.
U = Teuchos::null;
1803 newstate.
AU = Teuchos::null;
1804 newstate.
Alpha = Teuchos::null;
1805 newstate.
Beta = Teuchos::null;
1806 newstate.
D = Teuchos::null;
1807 newstate.
Delta = Teuchos::null;
1808 newstate.
LUUTAU = Teuchos::null;
1809 newstate.
ipiv = Teuchos::null;
1810 newstate.
rTz_old = Teuchos::null;
1813 index.resize(recycleBlocks_);
1814 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1815 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1816 index.resize(recycleBlocks_);
1817 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1818 Teuchos::RCP<MV> AUtmp = MVT::CloneViewNonConst( *AU_, index );
1820 problem_->applyOp( *Utmp, *AUtmp );
1822 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1823 MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp );
1825 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
1826 if ( precObj != Teuchos::null ) {
1827 index.resize(recycleBlocks_);
1828 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1829 index.resize(recycleBlocks_);
1830 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1831 Teuchos::RCP<MV> LeftPCAU = MVT::CloneViewNonConst( *U1_, index );
1832 OPT::Apply( *precObj, *AUtmp, *LeftPCAU );
1833 MVT::MvTransMv( one, *AUtmp, *LeftPCAU, AUTAUtmp );
1835 MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp );
1848 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1853 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails) );
1857 numIters_ = maxIterTest_->getNumIters();
1861 using Teuchos::rcp_dynamic_cast;
1864 const std::vector<MagnitudeType>* pTestValues =
1865 rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
1867 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1868 "Belos::RCGSolMgr::solve(): The convergence test's getTestValue() " 1869 "method returned NULL. Please report this bug to the Belos developers.");
1871 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1872 "Belos::RCGSolMgr::solve(): The convergence test's getTestValue() " 1873 "method returned a vector of length zero. Please report this bug to the " 1874 "Belos developers.");
1879 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1889 template<
class ScalarType,
class MV,
class OP>
1891 const Teuchos::SerialDenseMatrix<int,ScalarType>& G,
1892 Teuchos::SerialDenseMatrix<int,ScalarType>& Y) {
1894 int n = F.numCols();
1897 Teuchos::LAPACK<int,ScalarType> lapack;
1900 std::vector<MagnitudeType> w(n);
1903 std::vector<int> iperm(n);
1910 std::vector<ScalarType> work(1);
1914 Teuchos::SerialDenseMatrix<int,ScalarType> F2( Teuchos::Copy, *F_ );
1915 Teuchos::SerialDenseMatrix<int,ScalarType> G2( Teuchos::Copy, *G_ );
1918 lapack.SYGV(itype, jobz, uplo, n, G2.values(), G2.stride(), F2.values(), F2.stride(), &w[0], &work[0], lwork, &info);
1920 "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to query optimal work size.");
1921 lwork = (int)work[0];
1923 lapack.SYGV(itype, jobz, uplo, n, G2.values(), G2.stride(), F2.values(), F2.stride(), &w[0], &work[0], lwork, &info);
1925 "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to compute eigensolutions.");
1929 this->sort(w,n,iperm);
1932 for(
int i=0; i<recycleBlocks_; i++ ) {
1933 for(
int j=0; j<n; j++ ) {
1934 Y(j,i) = G2(j,iperm[i]);
1941 template<
class ScalarType,
class MV,
class OP>
1942 void RCGSolMgr<ScalarType,MV,OP,true>::sort(std::vector<ScalarType>& dlist,
int n, std::vector<int>& iperm)
1944 int l, r, j, i, flag;
1973 if (dlist[j] > dlist[j - 1]) j = j + 1;
1975 if (dlist[j - 1] > dK) {
1976 dlist[i - 1] = dlist[j - 1];
1977 iperm[i - 1] = iperm[j - 1];
1990 dlist[r] = dlist[0];
1991 iperm[r] = iperm[0];
2006 template<
class ScalarType,
class MV,
class OP>
2009 std::ostringstream oss;
2010 oss <<
"Belos::RCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<
">";
Teuchos::RCP< MV > r
The current residual.
Collection of types and exceptions used within the Belos solvers.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > rTz_old
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
Teuchos::RCP< MV > P
The current Krylov basis.
Teuchos::RCP< std::vector< int > > ipiv
Data from LU factorization of U^T A U.
Class which manages the output and verbosity of the Belos solvers.
RCGSolMgrRecyclingFailure is thrown when any problem occurs in using/creating the recycling subspace...
int getNumIters() const override
Get the iteration count for the most recent call to solve().
virtual ~RCGSolMgr()
Destructor.
Teuchos::RCP< MV > U
The recycled subspace and its image.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > LUUTAU
The LU factorization of the matrix U^T A U.
Base class for Belos::SolverManager subclasses which normally can only compile with real ScalarType t...
A factory class for generating StatusTestOutput objects.
Implementation of the RCG (Recycling Conjugate Gradient) iterative linear solver. ...
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Beta
This class implements the RCG iteration, where a single-std::vector Krylov subspace is constructed...
An implementation of StatusTestResNorm using a family of residual norms.
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
Belos::StatusTest class for specifying a maximum number of iterations.
bool isLOADetected() const override
Return whether a loss of accuracy was detected by this solver during the most current solve...
A factory class for generating StatusTestOutput objects.
Iterated Modified Gram-Schmidt (IMGS) implementation of the Belos::OrthoManager class.
RCGSolMgrLinearProblemFailure(const std::string &what_arg)
RCGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i.e.
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
Structure to contain pointers to RCGIter state variables.
Belos concrete class for performing the RCG iteration.
Classical Gram-Schmidt (with DGKS correction) implementation of the Belos::OrthoManager class...
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
A Belos::StatusTest class for specifying a maximum number of iterations.
ResetType
How to reset the solver.
bool existU
Flag to indicate the recycle space should be used.
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Return a reference to the linear problem being solved by this solver manager.
Pure virtual base class which describes the basic interface for a solver manager. ...
static constexpr double convTol
Default convergence tolerance.
Teuchos::RCP< MV > z
The current preconditioned residual.
Teuchos::RCP< MV > Ap
A times current search vector.
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
Type traits class that says whether Teuchos::LAPACK has a valid implementation for the given ScalarTy...
ReturnType
Whether the Belos solve converged for all linear systems.
Iterated Classical Gram-Schmidt (ICGS) implementation of the Belos::OrthoManager class.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Alpha
Coefficients arising in RCG iteration.
int setTolerance(MagnitudeType tolerance)
Set the value of the tolerance.
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > create(const Teuchos::RCP< OutputManager< ScalarType > > &printer, Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test, int mod, int printStates)
Create the StatusTestOutput object specified by the outputStyle.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
int curDim
The current dimension of the reduction.
RCGSolMgrLAPACKFailure(const std::string &what_arg)
RCGSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
A class for extending the status testing capabilities of Belos via logical combinations.
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.
Class which defines basic traits for the operator type.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > D
Parent class to all Belos exceptions.
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
RCGSolMgrRecyclingFailure(const std::string &what_arg)
Belos header file which uses auto-configuration information to include necessary C++ headers...
RCGSolMgrLAPACKFailure is thrown when a nonzero value is retuned from an LAPACK call.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Delta
Solutions to local least-squares problems.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
void reset(const ResetType type) override
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...