42#ifndef BELOS_RCG_SOLMGR_HPP
43#define BELOS_RCG_SOLMGR_HPP
61#include "Teuchos_LAPACK.hpp"
62#include "Teuchos_as.hpp"
63#ifdef BELOS_TEUCHOS_TIME_MONITOR
64#include "Teuchos_TimeMonitor.hpp"
136 template<
class ScalarType,
class MV,
class OP,
137 const bool supportsScalarType =
139 ! Teuchos::ScalarTraits<ScalarType>::isComplex>
142 Belos::Details::LapackSupportsScalar<ScalarType>::value &&
143 ! Teuchos::ScalarTraits<ScalarType>::isComplex>
147 ! Teuchos::ScalarTraits<ScalarType>::isComplex;
156 const Teuchos::RCP<Teuchos::ParameterList> &pl) :
162 Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const override {
170 template<
class ScalarType,
class MV,
class OP>
176 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
177 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
MagnitudeType;
178 typedef Teuchos::ScalarTraits<MagnitudeType>
MT;
214 const Teuchos::RCP<Teuchos::ParameterList> &pl );
220 Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const override {
233 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters()
const override;
243 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
244 return Teuchos::tuple(timerSolve_);
272 void setParameters(
const Teuchos::RCP<Teuchos::ParameterList> ¶ms )
override;
285 if ((type &
Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem();
318 std::string description()
const override;
329 void getHarmonicVecs(
const Teuchos::SerialDenseMatrix<int,ScalarType> &F,
330 const Teuchos::SerialDenseMatrix<int,ScalarType> &G,
331 Teuchos::SerialDenseMatrix<int,ScalarType>& Y);
334 void sort(std::vector<ScalarType>& dlist,
int n, std::vector<int>& iperm);
337 void initializeStateStorage();
340 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >
problem_;
347 Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
sTest_;
349 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> >
convTest_;
356 static constexpr int maxIters_default_ = 1000;
357 static constexpr int blockSize_default_ = 1;
358 static constexpr int numBlocks_default_ = 25;
359 static constexpr int recycleBlocks_default_ = 3;
360 static constexpr bool showMaxResNormOnly_default_ =
false;
363 static constexpr int outputFreq_default_ = -1;
364 static constexpr const char * label_default_ =
"Belos";
417 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
Alpha_;
418 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
Beta_;
419 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
D_;
422 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
Delta_;
425 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
UTAU_;
428 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
LUUTAU_;
431 Teuchos::RCP<std::vector<int> >
ipiv_;
434 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
AUTAU_;
437 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
rTz_old_;
440 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
F_,G_,Y_;
443 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > L2_,DeltaL2_,
AU1TUDeltaL2_;
444 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AU1TAU1_, AU1TU1_,
AU1TAP_;
445 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
FY_,GY_;
446 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
APTAP_;
448 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AUTAP_,
AU1TU_;
462template<
class ScalarType,
class MV,
class OP>
471template<
class ScalarType,
class MV,
class OP>
474 const Teuchos::RCP<Teuchos::ParameterList> &pl ) :
480 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument,
"Problem not given to solver manager.");
483 if ( !is_null(pl) ) {
489template<
class ScalarType,
class MV,
class OP>
492 outputStream_ = Teuchos::rcpFromRef(std::cout);
494 maxIters_ = maxIters_default_;
495 numBlocks_ = numBlocks_default_;
496 recycleBlocks_ = recycleBlocks_default_;
497 verbosity_ = verbosity_default_;
498 outputStyle_= outputStyle_default_;
499 outputFreq_= outputFreq_default_;
500 showMaxResNormOnly_ = showMaxResNormOnly_default_;
501 label_ = label_default_;
512 Alpha_ = Teuchos::null;
513 Beta_ = Teuchos::null;
515 Delta_ = Teuchos::null;
516 UTAU_ = Teuchos::null;
517 LUUTAU_ = Teuchos::null;
518 ipiv_ = Teuchos::null;
519 AUTAU_ = Teuchos::null;
520 rTz_old_ = Teuchos::null;
525 DeltaL2_ = Teuchos::null;
526 AU1TUDeltaL2_ = Teuchos::null;
527 AU1TAU1_ = Teuchos::null;
528 AU1TU1_ = Teuchos::null;
529 AU1TAP_ = Teuchos::null;
532 APTAP_ = Teuchos::null;
533 U1Y1_ = Teuchos::null;
534 PY2_ = Teuchos::null;
535 AUTAP_ = Teuchos::null;
536 AU1TU_ = Teuchos::null;
540template<
class ScalarType,
class MV,
class OP>
544 if (params_ == Teuchos::null) {
545 params_ = Teuchos::rcp(
new Teuchos::ParameterList(*getValidParameters()) );
548 params->validateParameters(*getValidParameters());
552 if (params->isParameter(
"Maximum Iterations")) {
553 maxIters_ = params->get(
"Maximum Iterations",maxIters_default_);
556 params_->set(
"Maximum Iterations", maxIters_);
557 if (maxIterTest_!=Teuchos::null)
558 maxIterTest_->setMaxIters( maxIters_ );
562 if (params->isParameter(
"Num Blocks")) {
563 numBlocks_ = params->get(
"Num Blocks",numBlocks_default_);
564 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
565 "Belos::RCGSolMgr: \"Num Blocks\" must be strictly positive.");
568 params_->set(
"Num Blocks", numBlocks_);
572 if (params->isParameter(
"Num Recycled Blocks")) {
573 recycleBlocks_ = params->get(
"Num Recycled Blocks",recycleBlocks_default_);
574 TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks_ <= 0, std::invalid_argument,
575 "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be strictly positive.");
577 TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks_ >= numBlocks_, std::invalid_argument,
578 "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be less than \"Num Blocks\".");
581 params_->set(
"Num Recycled Blocks", recycleBlocks_);
585 if (params->isParameter(
"Timer Label")) {
586 std::string tempLabel = params->get(
"Timer Label", label_default_);
589 if (tempLabel != label_) {
591 params_->set(
"Timer Label", label_);
592 std::string solveLabel = label_ +
": RCGSolMgr total solve time";
593#ifdef BELOS_TEUCHOS_TIME_MONITOR
594 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
600 if (params->isParameter(
"Verbosity")) {
601 if (Teuchos::isParameterType<int>(*params,
"Verbosity")) {
602 verbosity_ = params->get(
"Verbosity", verbosity_default_);
604 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,
"Verbosity");
608 params_->set(
"Verbosity", verbosity_);
609 if (printer_ != Teuchos::null)
610 printer_->setVerbosity(verbosity_);
614 if (params->isParameter(
"Output Style")) {
615 if (Teuchos::isParameterType<int>(*params,
"Output Style")) {
616 outputStyle_ = params->get(
"Output Style", outputStyle_default_);
618 outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,
"Output Style");
622 params_->set(
"Output Style", outputStyle_);
623 outputTest_ = Teuchos::null;
627 if (params->isParameter(
"Output Stream")) {
628 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,
"Output Stream");
631 params_->set(
"Output Stream", outputStream_);
632 if (printer_ != Teuchos::null)
633 printer_->setOStream( outputStream_ );
638 if (params->isParameter(
"Output Frequency")) {
639 outputFreq_ = params->get(
"Output Frequency", outputFreq_default_);
643 params_->set(
"Output Frequency", outputFreq_);
644 if (outputTest_ != Teuchos::null)
645 outputTest_->setOutputFrequency( outputFreq_ );
649 if (printer_ == Teuchos::null) {
658 if (params->isParameter(
"Convergence Tolerance")) {
659 if (params->isType<
MagnitudeType> (
"Convergence Tolerance")) {
660 convtol_ = params->get (
"Convergence Tolerance",
668 params_->set(
"Convergence Tolerance", convtol_);
669 if (convTest_ != Teuchos::null)
673 if (params->isParameter(
"Show Maximum Residual Norm Only")) {
674 showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,
"Show Maximum Residual Norm Only");
677 params_->set(
"Show Maximum Residual Norm Only", showMaxResNormOnly_);
678 if (convTest_ != Teuchos::null)
679 convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
685 if (maxIterTest_ == Teuchos::null)
689 if (convTest_ == Teuchos::null)
690 convTest_ = Teuchos::rcp(
new StatusTestResNorm_t( convtol_, 1 ) );
692 if (sTest_ == Teuchos::null)
693 sTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
695 if (outputTest_ == Teuchos::null) {
703 std::string solverDesc =
" Recycling CG ";
704 outputTest_->setSolverDesc( solverDesc );
708 if (timerSolve_ == Teuchos::null) {
709 std::string solveLabel = label_ +
": RCGSolMgr total solve time";
710#ifdef BELOS_TEUCHOS_TIME_MONITOR
711 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
720template<
class ScalarType,
class MV,
class OP>
721Teuchos::RCP<const Teuchos::ParameterList>
724 static Teuchos::RCP<const Teuchos::ParameterList> validPL;
727 if(is_null(validPL)) {
728 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
730 "The relative residual tolerance that needs to be achieved by the\n"
731 "iterative solver in order for the linear system to be declared converged.");
732 pl->set(
"Maximum Iterations",
static_cast<int>(maxIters_default_),
733 "The maximum number of block iterations allowed for each\n"
734 "set of RHS solved.");
735 pl->set(
"Block Size",
static_cast<int>(blockSize_default_),
736 "Block Size Parameter -- currently must be 1 for RCG");
737 pl->set(
"Num Blocks",
static_cast<int>(numBlocks_default_),
738 "The length of a cycle (and this max number of search vectors kept)\n");
739 pl->set(
"Num Recycled Blocks",
static_cast<int>(recycleBlocks_default_),
740 "The number of vectors in the recycle subspace.");
741 pl->set(
"Verbosity",
static_cast<int>(verbosity_default_),
742 "What type(s) of solver information should be outputted\n"
743 "to the output stream.");
744 pl->set(
"Output Style",
static_cast<int>(outputStyle_default_),
745 "What style is used for the solver information outputted\n"
746 "to the output stream.");
747 pl->set(
"Output Frequency",
static_cast<int>(outputFreq_default_),
748 "How often convergence information should be outputted\n"
749 "to the output stream.");
750 pl->set(
"Output Stream", Teuchos::rcpFromRef(std::cout),
751 "A reference-counted pointer to the output stream where all\n"
752 "solver output is sent.");
753 pl->set(
"Show Maximum Residual Norm Only",
static_cast<bool>(showMaxResNormOnly_default_),
754 "When convergence information is printed, only show the maximum\n"
755 "relative residual norm when the block size is greater than one.");
756 pl->set(
"Timer Label",
static_cast<const char *
>(label_default_),
757 "The string to use as a prefix for the timer labels.");
764template<
class ScalarType,
class MV,
class OP>
768 Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
769 if (rhsMV == Teuchos::null) {
776 TEUCHOS_TEST_FOR_EXCEPTION(
static_cast<ptrdiff_t
>(numBlocks_) > MVT::GetGlobalLength(*rhsMV),std::invalid_argument,
777 "Belos::RCGSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!");
780 if (P_ == Teuchos::null) {
781 P_ = MVT::Clone( *rhsMV, numBlocks_+2 );
785 if (MVT::GetNumberVecs(*P_) < numBlocks_+2) {
786 Teuchos::RCP<const MV> tmp = P_;
787 P_ = MVT::Clone( *tmp, numBlocks_+2 );
792 if (Ap_ == Teuchos::null)
793 Ap_ = MVT::Clone( *rhsMV, 1 );
796 if (r_ == Teuchos::null)
797 r_ = MVT::Clone( *rhsMV, 1 );
800 if (z_ == Teuchos::null)
801 z_ = MVT::Clone( *rhsMV, 1 );
804 if (U_ == Teuchos::null) {
805 U_ = MVT::Clone( *rhsMV, recycleBlocks_ );
809 if (MVT::GetNumberVecs(*U_) < recycleBlocks_) {
810 Teuchos::RCP<const MV> tmp = U_;
811 U_ = MVT::Clone( *tmp, recycleBlocks_ );
816 if (AU_ == Teuchos::null) {
817 AU_ = MVT::Clone( *rhsMV, recycleBlocks_ );
821 if (MVT::GetNumberVecs(*AU_) < recycleBlocks_) {
822 Teuchos::RCP<const MV> tmp = AU_;
823 AU_ = MVT::Clone( *tmp, recycleBlocks_ );
828 if (U1_ == Teuchos::null) {
829 U1_ = MVT::Clone( *rhsMV, recycleBlocks_ );
833 if (MVT::GetNumberVecs(*U1_) < recycleBlocks_) {
834 Teuchos::RCP<const MV> tmp = U1_;
835 U1_ = MVT::Clone( *tmp, recycleBlocks_ );
840 if (Alpha_ == Teuchos::null)
841 Alpha_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_, 1 ) );
843 if ( (Alpha_->numRows() != numBlocks_) || (Alpha_->numCols() != 1) )
844 Alpha_->reshape( numBlocks_, 1 );
848 if (Beta_ == Teuchos::null)
849 Beta_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + 1, 1 ) );
851 if ( (Beta_->numRows() != (numBlocks_+1)) || (Beta_->numCols() != 1) )
852 Beta_->reshape( numBlocks_ + 1, 1 );
856 if (D_ == Teuchos::null)
857 D_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ , 1 ) );
859 if ( (D_->numRows() != numBlocks_) || (D_->numCols() != 1) )
860 D_->reshape( numBlocks_, 1 );
864 if (Delta_ == Teuchos::null)
865 Delta_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ + 1 ) );
867 if ( (Delta_->numRows() != recycleBlocks_) || (Delta_->numCols() != (numBlocks_ + 1)) )
868 Delta_->reshape( recycleBlocks_, numBlocks_ + 1 );
872 if (UTAU_ == Teuchos::null)
873 UTAU_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
875 if ( (UTAU_->numRows() != recycleBlocks_) || (UTAU_->numCols() != recycleBlocks_) )
876 UTAU_->reshape( recycleBlocks_, recycleBlocks_ );
880 if (LUUTAU_ == Teuchos::null)
881 LUUTAU_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
883 if ( (LUUTAU_->numRows() != recycleBlocks_) || (LUUTAU_->numCols() != recycleBlocks_) )
884 LUUTAU_->reshape( recycleBlocks_, recycleBlocks_ );
888 if (ipiv_ == Teuchos::null)
889 ipiv_ = Teuchos::rcp(
new std::vector<int>(recycleBlocks_) );
891 if ( (
int)ipiv_->size() != recycleBlocks_ )
892 ipiv_->resize(recycleBlocks_);
896 if (AUTAU_ == Teuchos::null)
897 AUTAU_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
899 if ( (AUTAU_->numRows() != recycleBlocks_) || (AUTAU_->numCols() != recycleBlocks_) )
900 AUTAU_->reshape( recycleBlocks_, recycleBlocks_ );
904 if (rTz_old_ == Teuchos::null)
905 rTz_old_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( 1, 1 ) );
907 if ( (rTz_old_->numRows() != 1) || (rTz_old_->numCols() != 1) )
908 rTz_old_->reshape( 1, 1 );
912 if (F_ == Teuchos::null)
913 F_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ) );
915 if ( (F_->numRows() != (numBlocks_+recycleBlocks_)) || (F_->numCols() != numBlocks_+recycleBlocks_) )
916 F_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
920 if (G_ == Teuchos::null)
921 G_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ) );
923 if ( (G_->numRows() != (numBlocks_+recycleBlocks_)) || (G_->numCols() != numBlocks_+recycleBlocks_) )
924 G_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
928 if (Y_ == Teuchos::null)
929 Y_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, recycleBlocks_ ) );
931 if ( (Y_->numRows() != (numBlocks_+recycleBlocks_)) || (Y_->numCols() != numBlocks_+recycleBlocks_) )
932 Y_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
936 if (L2_ == Teuchos::null)
937 L2_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+1, numBlocks_ ) );
939 if ( (L2_->numRows() != (numBlocks_+1)) || (L2_->numCols() != numBlocks_) )
940 L2_->reshape( numBlocks_+1, numBlocks_ );
944 if (DeltaL2_ == Teuchos::null)
945 DeltaL2_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
947 if ( (DeltaL2_->numRows() != (recycleBlocks_)) || (DeltaL2_->numCols() != (numBlocks_) ) )
948 DeltaL2_->reshape( recycleBlocks_, numBlocks_ );
952 if (AU1TUDeltaL2_ == Teuchos::null)
953 AU1TUDeltaL2_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
955 if ( (AU1TUDeltaL2_->numRows() != (recycleBlocks_)) || (AU1TUDeltaL2_->numCols() != (numBlocks_) ) )
956 AU1TUDeltaL2_->reshape( recycleBlocks_, numBlocks_ );
960 if (AU1TAU1_ == Teuchos::null)
961 AU1TAU1_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
963 if ( (AU1TAU1_->numRows() != (recycleBlocks_)) || (AU1TAU1_->numCols() != (recycleBlocks_) ) )
964 AU1TAU1_->reshape( recycleBlocks_, recycleBlocks_ );
968 if (GY_ == Teuchos::null)
969 GY_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + recycleBlocks_, recycleBlocks_ ) );
971 if ( (GY_->numRows() != (numBlocks_ + recycleBlocks_)) || (GY_->numCols() != (recycleBlocks_) ) )
972 GY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ );
976 if (AU1TU1_ == Teuchos::null)
977 AU1TU1_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
979 if ( (AU1TU1_->numRows() != (recycleBlocks_)) || (AU1TU1_->numCols() != (recycleBlocks_) ) )
980 AU1TU1_->reshape( recycleBlocks_, recycleBlocks_ );
984 if (FY_ == Teuchos::null)
985 FY_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + recycleBlocks_, recycleBlocks_ ) );
987 if ( (FY_->numRows() != (numBlocks_ + recycleBlocks_)) || (FY_->numCols() != (recycleBlocks_) ) )
988 FY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ );
992 if (AU1TAP_ == Teuchos::null)
993 AU1TAP_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
995 if ( (AU1TAP_->numRows() != (recycleBlocks_)) || (AU1TAP_->numCols() != (numBlocks_) ) )
996 AU1TAP_->reshape( recycleBlocks_, numBlocks_ );
1000 if (APTAP_ == Teuchos::null)
1001 APTAP_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_, numBlocks_ ) );
1003 if ( (APTAP_->numRows() != (numBlocks_)) || (APTAP_->numCols() != (numBlocks_) ) )
1004 APTAP_->reshape( numBlocks_, numBlocks_ );
1008 if (U1Y1_ == Teuchos::null) {
1009 U1Y1_ = MVT::Clone( *rhsMV, recycleBlocks_ );
1013 if (MVT::GetNumberVecs(*U1Y1_) < recycleBlocks_) {
1014 Teuchos::RCP<const MV> tmp = U1Y1_;
1015 U1Y1_ = MVT::Clone( *tmp, recycleBlocks_ );
1020 if (PY2_ == Teuchos::null) {
1021 PY2_ = MVT::Clone( *rhsMV, recycleBlocks_ );
1025 if (MVT::GetNumberVecs(*PY2_) < recycleBlocks_) {
1026 Teuchos::RCP<const MV> tmp = PY2_;
1027 PY2_ = MVT::Clone( *tmp, recycleBlocks_ );
1032 if (AUTAP_ == Teuchos::null)
1033 AUTAP_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
1035 if ( (AUTAP_->numRows() != (recycleBlocks_)) || (AUTAP_->numCols() != (numBlocks_) ) )
1036 AUTAP_->reshape( recycleBlocks_, numBlocks_ );
1040 if (AU1TU_ == Teuchos::null)
1041 AU1TU_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
1043 if ( (AU1TU_->numRows() != (recycleBlocks_)) || (AU1TU_->numCols() != (recycleBlocks_) ) )
1044 AU1TU_->reshape( recycleBlocks_, recycleBlocks_ );
1051template<
class ScalarType,
class MV,
class OP>
1054 Teuchos::LAPACK<int,ScalarType> lapack;
1055 std::vector<int> index(recycleBlocks_);
1056 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1057 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1066 setParameters(Teuchos::parameterList(*getValidParameters()));
1070 "Belos::RCGSolMgr::solve(): Linear problem is not a valid object.");
1072 "Belos::RCGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
1073 TEUCHOS_TEST_FOR_EXCEPTION((problem_->getLeftPrec() != Teuchos::null)&&(problem_->getRightPrec() != Teuchos::null),
1075 "Belos::RCGSolMgr::solve(): RCG does not support split preconditioning, only set left or right preconditioner.");
1078 Teuchos::RCP<OP> precObj;
1079 if (problem_->getLeftPrec() != Teuchos::null) {
1080 precObj = Teuchos::rcp_const_cast<OP>(problem_->getLeftPrec());
1082 else if (problem_->getRightPrec() != Teuchos::null) {
1083 precObj = Teuchos::rcp_const_cast<OP>(problem_->getRightPrec());
1087 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1088 std::vector<int> currIdx(1);
1092 problem_->setLSIndex( currIdx );
1095 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1096 if (numBlocks_ > dim) {
1097 numBlocks_ = Teuchos::asSafe<int>(dim);
1098 params_->set(
"Num Blocks", numBlocks_);
1100 "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl <<
1101 " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl;
1105 initializeStateStorage();
1108 Teuchos::ParameterList plist;
1109 plist.set(
"Num Blocks",numBlocks_);
1110 plist.set(
"Recycled Blocks",recycleBlocks_);
1113 outputTest_->reset();
1116 bool isConverged =
true;
1120 index.resize(recycleBlocks_);
1121 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1122 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1123 index.resize(recycleBlocks_);
1124 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1125 Teuchos::RCP<MV> AUtmp = MVT::CloneViewNonConst( *AU_, index );
1127 problem_->applyOp( *Utmp, *AUtmp );
1129 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1130 MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp );
1132 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
1133 if ( precObj != Teuchos::null ) {
1134 index.resize(recycleBlocks_);
1135 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1136 index.resize(recycleBlocks_);
1137 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1138 Teuchos::RCP<MV> PCAU = MVT::CloneViewNonConst( *U1_, index );
1139 OPT::Apply( *precObj, *AUtmp, *PCAU );
1140 MVT::MvTransMv( one, *AUtmp, *PCAU, AUTAUtmp );
1142 MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp );
1149 Teuchos::RCP<RCGIter<ScalarType,MV,OP> > rcg_iter;
1154#ifdef BELOS_TEUCHOS_TIME_MONITOR
1155 Teuchos::TimeMonitor slvtimer(*timerSolve_);
1158 while ( numRHS2Solve > 0 ) {
1161 if (printer_->isVerbosity(
Debug ) ) {
1162 if (existU_) printer_->print(
Debug,
"Using recycle space generated from previous call to solve()." );
1163 else printer_->print(
Debug,
"No recycle space exists." );
1167 rcg_iter->resetNumIters();
1170 rcg_iter->setSize( recycleBlocks_, numBlocks_ );
1173 outputTest_->resetNumCalls();
1182 problem_->computeCurrResVec( &*r_ );
1187 Teuchos::SerialDenseMatrix<int,ScalarType> Utr(recycleBlocks_,1);
1188 index.resize(recycleBlocks_);
1189 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1190 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1191 MVT::MvTransMv( one, *Utmp, *r_, Utr );
1192 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1193 Teuchos::SerialDenseMatrix<int,ScalarType> LUUTAUtmp( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ );
1194 LUUTAUtmp.assign(UTAUtmp);
1196 lapack.GESV(recycleBlocks_, 1, LUUTAUtmp.values(), LUUTAUtmp.stride(), &(*ipiv_)[0], Utr.values(), Utr.stride(), &info);
1198 "Belos::RCGSolMgr::solve(): LAPACK GESV failed to compute a solution.");
1201 MVT::MvTimesMatAddMv( one, *Utmp, Utr, one, *problem_->getCurrLHSVec() );
1204 index.resize(recycleBlocks_);
1205 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1206 Teuchos::RCP<const MV> AUtmp = MVT::CloneView( *AU_, index );
1207 MVT::MvTimesMatAddMv( -one, *AUtmp, Utr, one, *r_ );
1210 if ( precObj != Teuchos::null ) {
1211 OPT::Apply( *precObj, *r_, *z_ );
1217 MVT::MvTransMv( one, *r_, *z_, *rTz_old_ );
1221 Teuchos::SerialDenseMatrix<int,ScalarType> mu( Teuchos::View, *Delta_, recycleBlocks_, 1);
1222 index.resize(recycleBlocks_);
1223 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1224 Teuchos::RCP<const MV> AUtmp = MVT::CloneView( *AU_, index );
1225 MVT::MvTransMv( one, *AUtmp, *z_, mu );
1226 Teuchos::SerialDenseMatrix<int,ScalarType> LUUTAUtmp( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ );
1229 lapack.GETRS(
TRANS, recycleBlocks_, 1, LUUTAUtmp.values(), LUUTAUtmp.stride(), &(*ipiv_)[0], mu.values(), mu.stride(), &info );
1231 "Belos::RCGSolMgr::solve(): LAPACK GETRS failed to compute a solution.");
1235 Teuchos::RCP<MV> Ptmp = MVT::CloneViewNonConst( *P_, index );
1236 MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp);
1237 MVT::MvTimesMatAddMv( -one, *U_, mu, one, *Ptmp );
1242 Teuchos::RCP<MV> Ptmp = MVT::CloneViewNonConst( *P_, index );
1243 MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp);
1250 index.resize( numBlocks_+1 );
1251 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1252 newstate.
P = MVT::CloneViewNonConst( *P_, index );
1253 index.resize( recycleBlocks_ );
1254 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1255 newstate.
U = MVT::CloneViewNonConst( *U_, index );
1256 index.resize( recycleBlocks_ );
1257 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1258 newstate.
AU = MVT::CloneViewNonConst( *AU_, index );
1259 newstate.
Alpha = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Alpha_, numBlocks_, 1 ) );
1260 newstate.
Beta = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Beta_, numBlocks_, 1 ) );
1261 newstate.
D = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *D_, numBlocks_, 1 ) );
1262 newstate.
Delta = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_, 0, 1 ) );
1263 newstate.
LUUTAU = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ ) );
1269 newstate.
existU = existU_;
1270 newstate.
ipiv = ipiv_;
1273 rcg_iter->initialize(newstate);
1279 rcg_iter->iterate();
1286 if ( convTest_->getStatus() ==
Passed ) {
1295 else if ( maxIterTest_->getStatus() ==
Passed ) {
1297 isConverged =
false;
1305 else if ( rcg_iter->getCurSubspaceDim() == rcg_iter->getMaxSubspaceDim() ) {
1310 if (recycleBlocks_ > 0) {
1314 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, numBlocks_, numBlocks_ );
1315 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, numBlocks_, numBlocks_ );
1316 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1317 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1318 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_, 1 );
1319 Ftmp.putScalar(zero);
1320 Gtmp.putScalar(zero);
1321 for (
int ii=0;ii<numBlocks_;ii++) {
1322 Gtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0));
1324 Gtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1325 Gtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1327 Ftmp(ii,ii) = Dtmp(ii,0);
1331 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, numBlocks_, recycleBlocks_ );
1332 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1335 index.resize( numBlocks_ );
1336 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; }
1337 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1338 index.resize( recycleBlocks_ );
1339 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1340 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1341 MVT::MvTimesMatAddMv( one, *Ptmp, Ytmp, zero, *U1tmp );
1346 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, numBlocks_, recycleBlocks_ );
1347 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1348 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1349 AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1352 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, numBlocks_, recycleBlocks_ );
1353 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1354 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1355 AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1357 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, 1 );
1359 AU1TAPtmp.putScalar(zero);
1361 ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0);
1362 for (
int ii=0; ii<recycleBlocks_; ++ii) {
1363 AU1TAPtmp(ii,0) = Ytmp(numBlocks_-1,ii) * alphatmp;
1371 lastBeta = numBlocks_-1;
1378 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1379 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, numBlocks_ );
1380 AU1TAPtmp.scale(Dtmp(0,0));
1382 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1383 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_+1, 1 );
1384 Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
1385 APTAPtmp.putScalar(zero);
1386 for (
int ii=0; ii<numBlocks_; ii++) {
1387 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0));
1389 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1390 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1395 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1396 Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
1397 Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1398 Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1399 Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1400 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1401 F11.assign(AU1TU1tmp);
1402 F12.putScalar(zero);
1403 F21.putScalar(zero);
1404 F22.putScalar(zero);
1405 for(
int ii=0;ii<numBlocks_;ii++) {
1406 F22(ii,ii) = Dtmp(ii,0);
1410 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1411 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1412 Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
1413 Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1414 Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1415 Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1416 G11.assign(AU1TAU1tmp);
1417 G12.assign(AU1TAPtmp);
1419 for (
int ii=0;ii<recycleBlocks_;++ii)
1420 for (
int jj=0;jj<numBlocks_;++jj)
1421 G21(jj,ii) = G12(ii,jj);
1422 G22.assign(APTAPtmp);
1425 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
1426 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1429 index.resize( numBlocks_ );
1430 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; }
1431 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1432 index.resize( recycleBlocks_ );
1433 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1434 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index );
1435 Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1436 index.resize( recycleBlocks_ );
1437 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1438 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1439 index.resize( recycleBlocks_ );
1440 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1441 Teuchos::RCP<MV> U1Y1tmp = MVT::CloneViewNonConst( *U1Y1_, index );
1442 Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
1443 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1444 MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp );
1445 MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp);
1450 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1451 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1452 AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1456 AU1TAPtmp.putScalar(zero);
1457 ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0);
1458 for (
int ii=0; ii<recycleBlocks_; ++ii) {
1459 AU1TAPtmp(ii,0) = Ytmp(numBlocks_+recycleBlocks_-1,ii) * alphatmp;
1463 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1464 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1466 AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1469 lastp = numBlocks_+1;
1470 lastBeta = numBlocks_;
1476 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1477 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_, 1 );
1478 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1479 Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
1480 APTAPtmp.putScalar(zero);
1481 for (
int ii=0; ii<numBlocks_; ii++) {
1482 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0));
1484 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1485 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1489 Teuchos::SerialDenseMatrix<int,ScalarType> L2tmp( Teuchos::View, *L2_, numBlocks_+1, numBlocks_ );
1490 L2tmp.putScalar(zero);
1491 for(
int ii=0;ii<numBlocks_;ii++) {
1492 L2tmp(ii,ii) = 1./Alphatmp(ii,0);
1493 L2tmp(ii+1,ii) = -1./Alphatmp(ii,0);
1497 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAPtmp( Teuchos::View, *AUTAP_, recycleBlocks_, numBlocks_ );
1498 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1499 Teuchos::SerialDenseMatrix<int,ScalarType> Deltatmp( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_+1 );
1500 Teuchos::SerialDenseMatrix<int,ScalarType> DeltaL2tmp( Teuchos::View, *DeltaL2_, recycleBlocks_, numBlocks_ );
1501 DeltaL2tmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Deltatmp,L2tmp,zero);
1502 AUTAPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,UTAUtmp,DeltaL2tmp,zero);
1505 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1506 Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
1507 Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1508 Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1509 Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1510 F11.assign(UTAUtmp);
1511 F12.putScalar(zero);
1512 F21.putScalar(zero);
1513 F22.putScalar(zero);
1514 for(
int ii=0;ii<numBlocks_;ii++) {
1515 F22(ii,ii) = Dtmp(ii,0);
1519 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1520 Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
1521 Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1522 Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1523 Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1524 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
1525 G11.assign(AUTAUtmp);
1526 G12.assign(AUTAPtmp);
1528 for (
int ii=0;ii<recycleBlocks_;++ii)
1529 for (
int jj=0;jj<numBlocks_;++jj)
1530 G21(jj,ii) = G12(ii,jj);
1531 G22.assign(APTAPtmp);
1534 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
1535 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1538 index.resize( recycleBlocks_ );
1539 for (
int ii=0; ii<(recycleBlocks_); ++ii) { index[ii] = ii; }
1540 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1541 index.resize( numBlocks_ );
1542 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; }
1543 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1544 index.resize( recycleBlocks_ );
1545 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1546 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index );
1547 Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1548 index.resize( recycleBlocks_ );
1549 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1550 Teuchos::RCP<MV> UY1tmp = MVT::CloneViewNonConst( *U1Y1_, index );
1551 Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
1552 index.resize( recycleBlocks_ );
1553 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1554 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1555 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1556 MVT::MvTimesMatAddMv( one, *Utmp, Y1, zero, *UY1tmp );
1557 MVT::MvAddMv(one,*UY1tmp, one, *PY2tmp, *U1tmp);
1562 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1563 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1564 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1565 AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1568 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1569 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1570 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1571 AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1574 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUtmp( Teuchos::View, *AU1TU_, recycleBlocks_, recycleBlocks_ );
1575 AU1TUtmp.assign(UTAUtmp);
1578 dold = Dtmp(numBlocks_-1,0);
1585 lastBeta = numBlocks_-1;
1588 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1589 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_+1, 1 );
1590 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1591 Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
1592 for (
int ii=0; ii<numBlocks_; ii++) {
1593 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0));
1595 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1596 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1600 Teuchos::SerialDenseMatrix<int,ScalarType> L2tmp( Teuchos::View, *L2_, numBlocks_+1, numBlocks_ );
1601 for(
int ii=0;ii<numBlocks_;ii++) {
1602 L2tmp(ii,ii) = 1./Alphatmp(ii,0);
1603 L2tmp(ii+1,ii) = -1./Alphatmp(ii,0);
1608 Teuchos::SerialDenseMatrix<int,ScalarType> DeltaL2( Teuchos::View, *DeltaL2_, recycleBlocks_, numBlocks_ );
1609 Teuchos::SerialDenseMatrix<int,ScalarType> Deltatmp( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_+1 );
1610 DeltaL2.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Deltatmp,L2tmp,zero);
1611 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUDeltaL2( Teuchos::View, *AU1TUDeltaL2_, recycleBlocks_, numBlocks_ );
1612 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUtmp( Teuchos::View, *AU1TU_, recycleBlocks_, recycleBlocks_ );
1613 AU1TUDeltaL2.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,AU1TUtmp,DeltaL2,zero);
1614 Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
1615 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, numBlocks_ );
1616 AU1TAPtmp.putScalar(zero);
1617 AU1TAPtmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Y1,AU1TUDeltaL2,zero);
1618 Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1619 ScalarType val = dold * (-Betatmp(0,0)/Alphatmp(0,0));
1620 for(
int ii=0;ii<recycleBlocks_;ii++) {
1621 AU1TAPtmp(ii,0) += Y2(numBlocks_-1,ii)*val;
1625 Teuchos::SerialDenseMatrix<int,ScalarType> Y1TAU1TU( Teuchos::View, *GY_, recycleBlocks_, recycleBlocks_ );
1626 Y1TAU1TU.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Y1,AU1TUtmp,zero);
1627 AU1TUtmp.assign(Y1TAU1TU);
1630 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1631 Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
1632 Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1633 Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1634 Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1635 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1636 F11.assign(AU1TU1tmp);
1637 for(
int ii=0;ii<numBlocks_;ii++) {
1638 F22(ii,ii) = Dtmp(ii,0);
1642 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1643 Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
1644 Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1645 Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1646 Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1647 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1648 G11.assign(AU1TAU1tmp);
1649 G12.assign(AU1TAPtmp);
1651 for (
int ii=0;ii<recycleBlocks_;++ii)
1652 for (
int jj=0;jj<numBlocks_;++jj)
1653 G21(jj,ii) = G12(ii,jj);
1654 G22.assign(APTAPtmp);
1657 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
1658 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1661 index.resize( numBlocks_ );
1662 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; }
1663 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1664 index.resize( recycleBlocks_ );
1665 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1666 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index );
1667 index.resize( recycleBlocks_ );
1668 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1669 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1670 index.resize( recycleBlocks_ );
1671 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1672 Teuchos::RCP<MV> U1Y1tmp = MVT::CloneViewNonConst( *U1Y1_, index );
1673 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1674 MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp );
1675 MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp);
1680 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1681 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1682 AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1685 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1686 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1687 AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1690 dold = Dtmp(numBlocks_-1,0);
1693 lastp = numBlocks_+1;
1694 lastBeta = numBlocks_;
1704 index[0] = lastp-1; index[1] = lastp;
1705 Teuchos::RCP<const MV> Ptmp2 = MVT::CloneView( *P_, index );
1706 index[0] = 0; index[1] = 1;
1707 MVT::SetBlock(*Ptmp2,index,*P_);
1710 (*Beta_)(0,0) = (*Beta_)(lastBeta,0);
1714 Teuchos::SerialDenseMatrix<int,ScalarType> mu1( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, 0 );
1715 Teuchos::SerialDenseMatrix<int,ScalarType> mu2( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, numBlocks_ );
1720 newstate.
P = Teuchos::null;
1721 index.resize( numBlocks_+1 );
1722 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii+1; }
1723 newstate.
P = MVT::CloneViewNonConst( *P_, index );
1725 newstate.
Beta = Teuchos::null;
1726 newstate.
Beta = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Beta_, numBlocks_, 1, 1, 0 ) );
1728 newstate.
Delta = Teuchos::null;
1729 newstate.
Delta = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_, 0, 1 ) );
1734 rcg_iter->initialize(newstate);
1747 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
1748 "Belos::RCGSolMgr::solve(): Invalid return from RCGIter::iterate().");
1751 catch (
const std::exception &e) {
1752 printer_->stream(
Errors) <<
"Error! Caught std::exception in RCGIter::iterate() at iteration "
1753 << rcg_iter->getNumIters() << std::endl
1754 << e.what() << std::endl;
1760 problem_->setCurrLS();
1764 if ( numRHS2Solve > 0 ) {
1767 problem_->setLSIndex( currIdx );
1770 currIdx.resize( numRHS2Solve );
1776 index.resize(recycleBlocks_);
1777 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1778 MVT::SetBlock(*U1_,index,*U_);
1781 if (numRHS2Solve > 0) {
1783 newstate.
P = Teuchos::null;
1784 newstate.
Ap = Teuchos::null;
1785 newstate.
r = Teuchos::null;
1786 newstate.
z = Teuchos::null;
1787 newstate.
U = Teuchos::null;
1788 newstate.
AU = Teuchos::null;
1789 newstate.
Alpha = Teuchos::null;
1790 newstate.
Beta = Teuchos::null;
1791 newstate.
D = Teuchos::null;
1792 newstate.
Delta = Teuchos::null;
1793 newstate.
LUUTAU = Teuchos::null;
1794 newstate.
ipiv = Teuchos::null;
1795 newstate.
rTz_old = Teuchos::null;
1798 index.resize(recycleBlocks_);
1799 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1800 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1801 index.resize(recycleBlocks_);
1802 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1803 Teuchos::RCP<MV> AUtmp = MVT::CloneViewNonConst( *AU_, index );
1805 problem_->applyOp( *Utmp, *AUtmp );
1807 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1808 MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp );
1810 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
1811 if ( precObj != Teuchos::null ) {
1812 index.resize(recycleBlocks_);
1813 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1814 index.resize(recycleBlocks_);
1815 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1816 Teuchos::RCP<MV> LeftPCAU = MVT::CloneViewNonConst( *U1_, index );
1817 OPT::Apply( *precObj, *AUtmp, *LeftPCAU );
1818 MVT::MvTransMv( one, *AUtmp, *LeftPCAU, AUTAUtmp );
1820 MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp );
1833#ifdef BELOS_TEUCHOS_TIME_MONITOR
1838 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails) );
1842 numIters_ = maxIterTest_->getNumIters();
1846 using Teuchos::rcp_dynamic_cast;
1849 const std::vector<MagnitudeType>* pTestValues =
1850 rcp_dynamic_cast<conv_test_type>(convTest_)->
getTestValue();
1852 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1853 "Belos::RCGSolMgr::solve(): The convergence test's getTestValue() "
1854 "method returned NULL. Please report this bug to the Belos developers.");
1856 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1857 "Belos::RCGSolMgr::solve(): The convergence test's getTestValue() "
1858 "method returned a vector of length zero. Please report this bug to the "
1859 "Belos developers.");
1864 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1874template<
class ScalarType,
class MV,
class OP>
1876 const Teuchos::SerialDenseMatrix<int,ScalarType>& ,
1877 Teuchos::SerialDenseMatrix<int,ScalarType>& Y) {
1879 int n = F.numCols();
1882 Teuchos::LAPACK<int,ScalarType> lapack;
1885 std::vector<MagnitudeType> w(n);
1888 std::vector<int> iperm(n);
1895 std::vector<ScalarType> work(1);
1899 Teuchos::SerialDenseMatrix<int,ScalarType> F2( Teuchos::Copy, *F_ );
1900 Teuchos::SerialDenseMatrix<int,ScalarType> G2( Teuchos::Copy, *G_ );
1903 lapack.SYGV(itype, jobz, uplo, n, G2.values(), G2.stride(), F2.values(), F2.stride(), &w[0], &work[0], lwork, &info);
1905 "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to query optimal work size.");
1906 lwork = (int)work[0];
1908 lapack.SYGV(itype, jobz, uplo, n, G2.values(), G2.stride(), F2.values(), F2.stride(), &w[0], &work[0], lwork, &info);
1910 "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to compute eigensolutions.");
1914 this->sort(w,n,iperm);
1917 for(
int i=0; i<recycleBlocks_; i++ ) {
1918 for(
int j=0; j<n; j++ ) {
1919 Y(j,i) = G2(j,iperm[i]);
1926template<
class ScalarType,
class MV,
class OP>
1929 int l, r, j, i, flag;
1958 if (dlist[j] > dlist[j - 1]) j = j + 1;
1960 if (dlist[j - 1] > dK) {
1961 dlist[i - 1] = dlist[j - 1];
1962 iperm[i - 1] = iperm[j - 1];
1975 dlist[r] = dlist[0];
1976 iperm[r] = iperm[0];
1991template<
class ScalarType,
class MV,
class OP>
1994 std::ostringstream oss;
1995 oss <<
"Belos::RCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<
">";
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.
Class which manages the output and verbosity of the Belos solvers.
Belos concrete class for performing the RCG iteration.
Pure virtual base class which describes the basic interface for a solver manager.
Belos::StatusTest for logically combining several status tests.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Belos::StatusTest class for specifying a maximum number of iterations.
A factory class for generating StatusTestOutput objects.
Collection of types and exceptions used within the Belos solvers.
Parent class to all Belos exceptions.
Type traits class that says whether Teuchos::LAPACK has a valid implementation for the given ScalarTy...
Base class for Belos::SolverManager subclasses which normally can only compile with real ScalarType t...
A linear system to solve, and its associated information.
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 RCG iteration, where a single-std::vector Krylov subspace is constructed.
RCGSolMgrLAPACKFailure is thrown when a nonzero value is retuned from an LAPACK call.
RCGSolMgrLAPACKFailure(const std::string &what_arg)
RCGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i.e.
RCGSolMgrLinearProblemFailure(const std::string &what_arg)
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.
void reset(const ResetType type) override
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > convTest_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Alpha_
MagnitudeType convtol_
Convergence tolerance (read from parameter list).
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > rTz_old_
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
Teuchos::ScalarTraits< ScalarType > SCT
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Return a reference to the linear problem being solved by this solver manager.
OperatorTraits< ScalarType, MV, OP > OPT
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
MagnitudeType achievedTol_
Tolerance achieved by the last solve() invocation.
Teuchos::RCP< std::vector< int > > ipiv_
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
int maxIters_
Maximum iteration count (read from parameter list).
int getNumIters() const override
Get the iteration count for the most recent call to solve().
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > APTAP_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > AU1TU_
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
Teuchos::ScalarTraits< MagnitudeType > MT
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > AU1TAP_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > D_
virtual ~RCGSolMgr()
Destructor.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > F_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Beta_
Teuchos::RCP< std::ostream > outputStream_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > UTAU_
int numIters_
Number of iterations taken by the last solve() invocation.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > LUUTAU_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > AU1TUDeltaL2_
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > AUTAU_
bool isLOADetected() const override
Return whether a loss of accuracy was detected by this solver during the most current solve.
MultiVecTraits< ScalarType, MV > MVT
Teuchos::RCP< Teuchos::ParameterList > params_
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
Teuchos::RCP< OutputManager< ScalarType > > printer_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Delta_
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > FY_
Teuchos::RCP< Teuchos::Time > timerSolve_
Implementation of the RCG (Recycling Conjugate Gradient) iterative linear solver.
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
static const bool scalarTypeIsSupported
RCGSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
Details::SolverManagerRequiresRealLapack< ScalarType, MV, OP, scalarTypeIsSupported > base_type
A class for extending the status testing capabilities of Belos via logical combinations.
An implementation of StatusTestResNorm using a family of residual norms.
const std::vector< MagnitudeType > * getTestValue() const
Returns the test value, , computed in most recent call to CheckStatus.
int setTolerance(MagnitudeType tolerance)
Set the value of the tolerance.
A Belos::StatusTest class for specifying a maximum number of iterations.
A factory class for generating StatusTestOutput objects.
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.
ReturnType
Whether the Belos solve converged for all linear systems.
ResetType
How to reset the solver.
static const double convTol
Default convergence tolerance.
Structure to contain pointers to RCGIter state variables.
bool existU
Flag to indicate the recycle space should be used.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Beta
Teuchos::RCP< MV > U
The recycled subspace and its image.
Teuchos::RCP< MV > Ap
A times current search vector.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Delta
Solutions to local least-squares problems.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > D
int curDim
The current dimension of the reduction.
Teuchos::RCP< MV > r
The current residual.
Teuchos::RCP< std::vector< int > > ipiv
Data from LU factorization of U^T A U.
Teuchos::RCP< MV > z
The current preconditioned residual.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > LUUTAU
The LU factorization of the matrix U^T A U
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > rTz_old
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Alpha
Coefficients arising in RCG iteration.
Teuchos::RCP< MV > P
The current Krylov basis.