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" 146 template<
class ScalarType,
class MV,
class OP,
147 const bool supportsScalarType =
149 ! Teuchos::ScalarTraits<ScalarType>::isComplex>
152 Belos::Details::LapackSupportsScalar<ScalarType>::value &&
153 ! Teuchos::ScalarTraits<ScalarType>::isComplex>
157 ! Teuchos::ScalarTraits<ScalarType>::isComplex;
166 const Teuchos::RCP<Teuchos::ParameterList> &pl) :
172 Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const override {
180 template<
class ScalarType,
class MV,
class OP>
186 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
187 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
MagnitudeType;
188 typedef Teuchos::ScalarTraits<MagnitudeType>
MT;
224 const Teuchos::RCP<Teuchos::ParameterList> &pl );
230 Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const override {
243 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters()
const override;
253 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
254 return Teuchos::tuple(timerSolve_);
282 void setParameters(
const Teuchos::RCP<Teuchos::ParameterList> ¶ms )
override;
295 if ((type &
Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem();
328 std::string description()
const override;
339 void getHarmonicVecs(
const Teuchos::SerialDenseMatrix<int,ScalarType> &F,
340 const Teuchos::SerialDenseMatrix<int,ScalarType> &G,
341 Teuchos::SerialDenseMatrix<int,ScalarType>& Y);
344 void sort(std::vector<ScalarType>& dlist,
int n, std::vector<int>& iperm);
347 void initializeStateStorage();
350 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >
problem_;
357 Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
sTest_;
359 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> >
convTest_;
366 static constexpr
int maxIters_default_ = 1000;
367 static constexpr
int blockSize_default_ = 1;
368 static constexpr
int numBlocks_default_ = 25;
369 static constexpr
int recycleBlocks_default_ = 3;
370 static constexpr
bool showMaxResNormOnly_default_ =
false;
373 static constexpr
int outputFreq_default_ = -1;
374 static constexpr
const char * label_default_ =
"Belos";
376 #if defined(_WIN32) && defined(__clang__) 377 static constexpr std::ostream * outputStream_default_ =
378 __builtin_constant_p(reinterpret_cast<const std::ostream*>(&std::cout));
380 static constexpr std::ostream * outputStream_default_ = &std::cout;
428 Teuchos::RCP<MV>
U_, AU_;
434 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
Alpha_;
435 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
Beta_;
436 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
D_;
439 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
Delta_;
442 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
UTAU_;
445 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
LUUTAU_;
448 Teuchos::RCP<std::vector<int> >
ipiv_;
451 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
AUTAU_;
454 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
rTz_old_;
457 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > F_,G_,
Y_;
460 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
L2_,DeltaL2_,AU1TUDeltaL2_;
461 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AU1TAU1_,
AU1TU1_, AU1TAP_;
462 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > FY_,
GY_;
463 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
APTAP_;
465 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
AUTAP_, AU1TU_;
479 template<
class ScalarType,
class MV,
class OP>
488 template<
class ScalarType,
class MV,
class OP>
491 const Teuchos::RCP<Teuchos::ParameterList> &pl ) :
497 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument,
"Problem not given to solver manager.");
500 if ( !is_null(pl) ) {
506 template<
class ScalarType,
class MV,
class OP>
509 outputStream_ = Teuchos::rcp(outputStream_default_,
false);
511 maxIters_ = maxIters_default_;
512 numBlocks_ = numBlocks_default_;
513 recycleBlocks_ = recycleBlocks_default_;
514 verbosity_ = verbosity_default_;
515 outputStyle_= outputStyle_default_;
516 outputFreq_= outputFreq_default_;
517 showMaxResNormOnly_ = showMaxResNormOnly_default_;
518 label_ = label_default_;
529 Alpha_ = Teuchos::null;
530 Beta_ = Teuchos::null;
532 Delta_ = Teuchos::null;
533 UTAU_ = Teuchos::null;
534 LUUTAU_ = Teuchos::null;
535 ipiv_ = Teuchos::null;
536 AUTAU_ = Teuchos::null;
537 rTz_old_ = Teuchos::null;
542 DeltaL2_ = Teuchos::null;
543 AU1TUDeltaL2_ = Teuchos::null;
544 AU1TAU1_ = Teuchos::null;
545 AU1TU1_ = Teuchos::null;
546 AU1TAP_ = Teuchos::null;
549 APTAP_ = Teuchos::null;
550 U1Y1_ = Teuchos::null;
551 PY2_ = Teuchos::null;
552 AUTAP_ = Teuchos::null;
553 AU1TU_ = Teuchos::null;
557 template<
class ScalarType,
class MV,
class OP>
561 if (params_ == Teuchos::null) {
562 params_ = Teuchos::rcp(
new Teuchos::ParameterList(*getValidParameters()) );
565 params->validateParameters(*getValidParameters());
569 if (params->isParameter(
"Maximum Iterations")) {
570 maxIters_ = params->get(
"Maximum Iterations",maxIters_default_);
573 params_->set(
"Maximum Iterations", maxIters_);
574 if (maxIterTest_!=Teuchos::null)
575 maxIterTest_->setMaxIters( maxIters_ );
579 if (params->isParameter(
"Num Blocks")) {
580 numBlocks_ = params->get(
"Num Blocks",numBlocks_default_);
581 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
582 "Belos::RCGSolMgr: \"Num Blocks\" must be strictly positive.");
585 params_->set(
"Num Blocks", numBlocks_);
589 if (params->isParameter(
"Num Recycled Blocks")) {
590 recycleBlocks_ = params->get(
"Num Recycled Blocks",recycleBlocks_default_);
591 TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks_ <= 0, std::invalid_argument,
592 "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be strictly positive.");
594 TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks_ >= numBlocks_, std::invalid_argument,
595 "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be less than \"Num Blocks\".");
598 params_->set(
"Num Recycled Blocks", recycleBlocks_);
602 if (params->isParameter(
"Timer Label")) {
603 std::string tempLabel = params->get(
"Timer Label", label_default_);
606 if (tempLabel != label_) {
608 params_->set(
"Timer Label", label_);
609 std::string solveLabel = label_ +
": RCGSolMgr total solve time";
610 #ifdef BELOS_TEUCHOS_TIME_MONITOR 611 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
617 if (params->isParameter(
"Verbosity")) {
618 if (Teuchos::isParameterType<int>(*params,
"Verbosity")) {
619 verbosity_ = params->get(
"Verbosity", verbosity_default_);
621 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,
"Verbosity");
625 params_->set(
"Verbosity", verbosity_);
626 if (printer_ != Teuchos::null)
627 printer_->setVerbosity(verbosity_);
631 if (params->isParameter(
"Output Style")) {
632 if (Teuchos::isParameterType<int>(*params,
"Output Style")) {
633 outputStyle_ = params->get(
"Output Style", outputStyle_default_);
635 outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,
"Output Style");
639 params_->set(
"Output Style", outputStyle_);
640 outputTest_ = Teuchos::null;
644 if (params->isParameter(
"Output Stream")) {
645 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,
"Output Stream");
648 params_->set(
"Output Stream", outputStream_);
649 if (printer_ != Teuchos::null)
650 printer_->setOStream( outputStream_ );
655 if (params->isParameter(
"Output Frequency")) {
656 outputFreq_ = params->get(
"Output Frequency", outputFreq_default_);
660 params_->set(
"Output Frequency", outputFreq_);
661 if (outputTest_ != Teuchos::null)
662 outputTest_->setOutputFrequency( outputFreq_ );
666 if (printer_ == Teuchos::null) {
675 if (params->isParameter(
"Convergence Tolerance")) {
676 if (params->isType<
MagnitudeType> (
"Convergence Tolerance")) {
677 convtol_ = params->get (
"Convergence Tolerance",
685 params_->set(
"Convergence Tolerance", convtol_);
686 if (convTest_ != Teuchos::null)
690 if (params->isParameter(
"Show Maximum Residual Norm Only")) {
691 showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,
"Show Maximum Residual Norm Only");
694 params_->set(
"Show Maximum Residual Norm Only", showMaxResNormOnly_);
695 if (convTest_ != Teuchos::null)
696 convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
702 if (maxIterTest_ == Teuchos::null)
706 if (convTest_ == Teuchos::null)
707 convTest_ = Teuchos::rcp(
new StatusTestResNorm_t( convtol_, 1 ) );
709 if (sTest_ == Teuchos::null)
710 sTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
712 if (outputTest_ == Teuchos::null) {
720 std::string solverDesc =
" Recycling CG ";
721 outputTest_->setSolverDesc( solverDesc );
725 if (timerSolve_ == Teuchos::null) {
726 std::string solveLabel = label_ +
": RCGSolMgr total solve time";
727 #ifdef BELOS_TEUCHOS_TIME_MONITOR 728 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
737 template<
class ScalarType,
class MV,
class OP>
738 Teuchos::RCP<const Teuchos::ParameterList>
741 static Teuchos::RCP<const Teuchos::ParameterList> validPL;
744 if(is_null(validPL)) {
745 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
747 "The relative residual tolerance that needs to be achieved by the\n" 748 "iterative solver in order for the linear system to be declared converged.");
749 pl->set(
"Maximum Iterations", static_cast<int>(maxIters_default_),
750 "The maximum number of block iterations allowed for each\n" 751 "set of RHS solved.");
752 pl->set(
"Block Size", static_cast<int>(blockSize_default_),
753 "Block Size Parameter -- currently must be 1 for RCG");
754 pl->set(
"Num Blocks", static_cast<int>(numBlocks_default_),
755 "The length of a cycle (and this max number of search vectors kept)\n");
756 pl->set(
"Num Recycled Blocks", static_cast<int>(recycleBlocks_default_),
757 "The number of vectors in the recycle subspace.");
758 pl->set(
"Verbosity", static_cast<int>(verbosity_default_),
759 "What type(s) of solver information should be outputted\n" 760 "to the output stream.");
761 pl->set(
"Output Style", static_cast<int>(outputStyle_default_),
762 "What style is used for the solver information outputted\n" 763 "to the output stream.");
764 pl->set(
"Output Frequency", static_cast<int>(outputFreq_default_),
765 "How often convergence information should be outputted\n" 766 "to the output stream.");
767 pl->set(
"Output Stream", Teuchos::rcp(outputStream_default_,
false),
768 "A reference-counted pointer to the output stream where all\n" 769 "solver output is sent.");
770 pl->set(
"Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
771 "When convergence information is printed, only show the maximum\n" 772 "relative residual norm when the block size is greater than one.");
773 pl->set(
"Timer Label", static_cast<const char *>(label_default_),
774 "The string to use as a prefix for the timer labels.");
781 template<
class ScalarType,
class MV,
class OP>
785 Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
786 if (rhsMV == Teuchos::null) {
793 TEUCHOS_TEST_FOR_EXCEPTION(static_cast<ptrdiff_t>(numBlocks_) > MVT::GetGlobalLength(*rhsMV),std::invalid_argument,
794 "Belos::RCGSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!");
797 if (P_ == Teuchos::null) {
798 P_ = MVT::Clone( *rhsMV, numBlocks_+2 );
802 if (MVT::GetNumberVecs(*P_) < numBlocks_+2) {
803 Teuchos::RCP<const MV> tmp = P_;
804 P_ = MVT::Clone( *tmp, numBlocks_+2 );
809 if (Ap_ == Teuchos::null)
810 Ap_ = MVT::Clone( *rhsMV, 1 );
813 if (r_ == Teuchos::null)
814 r_ = MVT::Clone( *rhsMV, 1 );
817 if (z_ == Teuchos::null)
818 z_ = MVT::Clone( *rhsMV, 1 );
821 if (U_ == Teuchos::null) {
822 U_ = MVT::Clone( *rhsMV, recycleBlocks_ );
826 if (MVT::GetNumberVecs(*U_) < recycleBlocks_) {
827 Teuchos::RCP<const MV> tmp = U_;
828 U_ = MVT::Clone( *tmp, recycleBlocks_ );
833 if (AU_ == Teuchos::null) {
834 AU_ = MVT::Clone( *rhsMV, recycleBlocks_ );
838 if (MVT::GetNumberVecs(*AU_) < recycleBlocks_) {
839 Teuchos::RCP<const MV> tmp = AU_;
840 AU_ = MVT::Clone( *tmp, recycleBlocks_ );
845 if (U1_ == Teuchos::null) {
846 U1_ = MVT::Clone( *rhsMV, recycleBlocks_ );
850 if (MVT::GetNumberVecs(*U1_) < recycleBlocks_) {
851 Teuchos::RCP<const MV> tmp = U1_;
852 U1_ = MVT::Clone( *tmp, recycleBlocks_ );
857 if (Alpha_ == Teuchos::null)
858 Alpha_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_, 1 ) );
860 if ( (Alpha_->numRows() != numBlocks_) || (Alpha_->numCols() != 1) )
861 Alpha_->reshape( numBlocks_, 1 );
865 if (Beta_ == Teuchos::null)
866 Beta_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + 1, 1 ) );
868 if ( (Beta_->numRows() != (numBlocks_+1)) || (Beta_->numCols() != 1) )
869 Beta_->reshape( numBlocks_ + 1, 1 );
873 if (D_ == Teuchos::null)
874 D_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ , 1 ) );
876 if ( (D_->numRows() != numBlocks_) || (D_->numCols() != 1) )
877 D_->reshape( numBlocks_, 1 );
881 if (Delta_ == Teuchos::null)
882 Delta_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ + 1 ) );
884 if ( (Delta_->numRows() != recycleBlocks_) || (Delta_->numCols() != (numBlocks_ + 1)) )
885 Delta_->reshape( recycleBlocks_, numBlocks_ + 1 );
889 if (UTAU_ == Teuchos::null)
890 UTAU_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
892 if ( (UTAU_->numRows() != recycleBlocks_) || (UTAU_->numCols() != recycleBlocks_) )
893 UTAU_->reshape( recycleBlocks_, recycleBlocks_ );
897 if (LUUTAU_ == Teuchos::null)
898 LUUTAU_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
900 if ( (LUUTAU_->numRows() != recycleBlocks_) || (LUUTAU_->numCols() != recycleBlocks_) )
901 LUUTAU_->reshape( recycleBlocks_, recycleBlocks_ );
905 if (ipiv_ == Teuchos::null)
906 ipiv_ = Teuchos::rcp(
new std::vector<int>(recycleBlocks_) );
908 if ( (
int)ipiv_->size() != recycleBlocks_ )
909 ipiv_->resize(recycleBlocks_);
913 if (AUTAU_ == Teuchos::null)
914 AUTAU_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
916 if ( (AUTAU_->numRows() != recycleBlocks_) || (AUTAU_->numCols() != recycleBlocks_) )
917 AUTAU_->reshape( recycleBlocks_, recycleBlocks_ );
921 if (rTz_old_ == Teuchos::null)
922 rTz_old_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( 1, 1 ) );
924 if ( (rTz_old_->numRows() != 1) || (rTz_old_->numCols() != 1) )
925 rTz_old_->reshape( 1, 1 );
929 if (F_ == Teuchos::null)
930 F_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ) );
932 if ( (F_->numRows() != (numBlocks_+recycleBlocks_)) || (F_->numCols() != numBlocks_+recycleBlocks_) )
933 F_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
937 if (G_ == Teuchos::null)
938 G_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ) );
940 if ( (G_->numRows() != (numBlocks_+recycleBlocks_)) || (G_->numCols() != numBlocks_+recycleBlocks_) )
941 G_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
945 if (Y_ == Teuchos::null)
946 Y_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, recycleBlocks_ ) );
948 if ( (Y_->numRows() != (numBlocks_+recycleBlocks_)) || (Y_->numCols() != numBlocks_+recycleBlocks_) )
949 Y_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
953 if (L2_ == Teuchos::null)
954 L2_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+1, numBlocks_ ) );
956 if ( (L2_->numRows() != (numBlocks_+1)) || (L2_->numCols() != numBlocks_) )
957 L2_->reshape( numBlocks_+1, numBlocks_ );
961 if (DeltaL2_ == Teuchos::null)
962 DeltaL2_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
964 if ( (DeltaL2_->numRows() != (recycleBlocks_)) || (DeltaL2_->numCols() != (numBlocks_) ) )
965 DeltaL2_->reshape( recycleBlocks_, numBlocks_ );
969 if (AU1TUDeltaL2_ == Teuchos::null)
970 AU1TUDeltaL2_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
972 if ( (AU1TUDeltaL2_->numRows() != (recycleBlocks_)) || (AU1TUDeltaL2_->numCols() != (numBlocks_) ) )
973 AU1TUDeltaL2_->reshape( recycleBlocks_, numBlocks_ );
977 if (AU1TAU1_ == Teuchos::null)
978 AU1TAU1_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
980 if ( (AU1TAU1_->numRows() != (recycleBlocks_)) || (AU1TAU1_->numCols() != (recycleBlocks_) ) )
981 AU1TAU1_->reshape( recycleBlocks_, recycleBlocks_ );
985 if (GY_ == Teuchos::null)
986 GY_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + recycleBlocks_, recycleBlocks_ ) );
988 if ( (GY_->numRows() != (numBlocks_ + recycleBlocks_)) || (GY_->numCols() != (recycleBlocks_) ) )
989 GY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ );
993 if (AU1TU1_ == Teuchos::null)
994 AU1TU1_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
996 if ( (AU1TU1_->numRows() != (recycleBlocks_)) || (AU1TU1_->numCols() != (recycleBlocks_) ) )
997 AU1TU1_->reshape( recycleBlocks_, recycleBlocks_ );
1001 if (FY_ == Teuchos::null)
1002 FY_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + recycleBlocks_, recycleBlocks_ ) );
1004 if ( (FY_->numRows() != (numBlocks_ + recycleBlocks_)) || (FY_->numCols() != (recycleBlocks_) ) )
1005 FY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ );
1009 if (AU1TAP_ == Teuchos::null)
1010 AU1TAP_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
1012 if ( (AU1TAP_->numRows() != (recycleBlocks_)) || (AU1TAP_->numCols() != (numBlocks_) ) )
1013 AU1TAP_->reshape( recycleBlocks_, numBlocks_ );
1017 if (APTAP_ == Teuchos::null)
1018 APTAP_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_, numBlocks_ ) );
1020 if ( (APTAP_->numRows() != (numBlocks_)) || (APTAP_->numCols() != (numBlocks_) ) )
1021 APTAP_->reshape( numBlocks_, numBlocks_ );
1025 if (U1Y1_ == Teuchos::null) {
1026 U1Y1_ = MVT::Clone( *rhsMV, recycleBlocks_ );
1030 if (MVT::GetNumberVecs(*U1Y1_) < recycleBlocks_) {
1031 Teuchos::RCP<const MV> tmp = U1Y1_;
1032 U1Y1_ = MVT::Clone( *tmp, recycleBlocks_ );
1037 if (PY2_ == Teuchos::null) {
1038 PY2_ = MVT::Clone( *rhsMV, recycleBlocks_ );
1042 if (MVT::GetNumberVecs(*PY2_) < recycleBlocks_) {
1043 Teuchos::RCP<const MV> tmp = PY2_;
1044 PY2_ = MVT::Clone( *tmp, recycleBlocks_ );
1049 if (AUTAP_ == Teuchos::null)
1050 AUTAP_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
1052 if ( (AUTAP_->numRows() != (recycleBlocks_)) || (AUTAP_->numCols() != (numBlocks_) ) )
1053 AUTAP_->reshape( recycleBlocks_, numBlocks_ );
1057 if (AU1TU_ == Teuchos::null)
1058 AU1TU_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
1060 if ( (AU1TU_->numRows() != (recycleBlocks_)) || (AU1TU_->numCols() != (recycleBlocks_) ) )
1061 AU1TU_->reshape( recycleBlocks_, recycleBlocks_ );
1068 template<
class ScalarType,
class MV,
class OP>
1071 Teuchos::LAPACK<int,ScalarType> lapack;
1072 std::vector<int> index(recycleBlocks_);
1073 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1074 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1083 setParameters(Teuchos::parameterList(*getValidParameters()));
1087 "Belos::RCGSolMgr::solve(): Linear problem is not a valid object.");
1089 "Belos::RCGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
1090 TEUCHOS_TEST_FOR_EXCEPTION((problem_->getLeftPrec() != Teuchos::null)&&(problem_->getRightPrec() != Teuchos::null),
1092 "Belos::RCGSolMgr::solve(): RCG does not support split preconditioning, only set left or right preconditioner.");
1095 Teuchos::RCP<OP> precObj;
1096 if (problem_->getLeftPrec() != Teuchos::null) {
1097 precObj = Teuchos::rcp_const_cast<OP>(problem_->getLeftPrec());
1099 else if (problem_->getRightPrec() != Teuchos::null) {
1100 precObj = Teuchos::rcp_const_cast<OP>(problem_->getRightPrec());
1104 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1105 std::vector<int> currIdx(1);
1109 problem_->setLSIndex( currIdx );
1112 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1113 if (numBlocks_ > dim) {
1114 numBlocks_ = Teuchos::asSafe<int>(dim);
1115 params_->set(
"Num Blocks", numBlocks_);
1117 "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl <<
1118 " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl;
1122 initializeStateStorage();
1125 Teuchos::ParameterList plist;
1126 plist.set(
"Num Blocks",numBlocks_);
1127 plist.set(
"Recycled Blocks",recycleBlocks_);
1130 outputTest_->reset();
1133 bool isConverged =
true;
1137 index.resize(recycleBlocks_);
1138 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1139 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1140 index.resize(recycleBlocks_);
1141 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1142 Teuchos::RCP<MV> AUtmp = MVT::CloneViewNonConst( *AU_, index );
1144 problem_->applyOp( *Utmp, *AUtmp );
1146 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1147 MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp );
1149 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
1150 if ( precObj != Teuchos::null ) {
1151 index.resize(recycleBlocks_);
1152 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1153 index.resize(recycleBlocks_);
1154 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1155 Teuchos::RCP<MV> PCAU = MVT::CloneViewNonConst( *U1_, index );
1156 OPT::Apply( *precObj, *AUtmp, *PCAU );
1157 MVT::MvTransMv( one, *AUtmp, *PCAU, AUTAUtmp );
1159 MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp );
1166 Teuchos::RCP<RCGIter<ScalarType,MV,OP> > rcg_iter;
1171 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1172 Teuchos::TimeMonitor slvtimer(*timerSolve_);
1175 while ( numRHS2Solve > 0 ) {
1178 if (printer_->isVerbosity(
Debug ) ) {
1179 if (existU_) printer_->print(
Debug,
"Using recycle space generated from previous call to solve()." );
1180 else printer_->print(
Debug,
"No recycle space exists." );
1184 rcg_iter->resetNumIters();
1187 rcg_iter->setSize( recycleBlocks_, numBlocks_ );
1190 outputTest_->resetNumCalls();
1199 problem_->computeCurrResVec( &*r_ );
1204 Teuchos::SerialDenseMatrix<int,ScalarType> Utr(recycleBlocks_,1);
1205 index.resize(recycleBlocks_);
1206 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1207 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1208 MVT::MvTransMv( one, *Utmp, *r_, Utr );
1209 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1210 Teuchos::SerialDenseMatrix<int,ScalarType> LUUTAUtmp( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ );
1211 LUUTAUtmp.assign(UTAUtmp);
1213 lapack.GESV(recycleBlocks_, 1, LUUTAUtmp.values(), LUUTAUtmp.stride(), &(*ipiv_)[0], Utr.values(), Utr.stride(), &info);
1215 "Belos::RCGSolMgr::solve(): LAPACK GESV failed to compute a solution.");
1218 MVT::MvTimesMatAddMv( one, *Utmp, Utr, one, *problem_->getCurrLHSVec() );
1221 index.resize(recycleBlocks_);
1222 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1223 Teuchos::RCP<const MV> AUtmp = MVT::CloneView( *AU_, index );
1224 MVT::MvTimesMatAddMv( -one, *AUtmp, Utr, one, *r_ );
1227 if ( precObj != Teuchos::null ) {
1228 OPT::Apply( *precObj, *r_, *z_ );
1234 MVT::MvTransMv( one, *r_, *z_, *rTz_old_ );
1238 Teuchos::SerialDenseMatrix<int,ScalarType> mu( Teuchos::View, *Delta_, recycleBlocks_, 1);
1239 index.resize(recycleBlocks_);
1240 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1241 Teuchos::RCP<const MV> AUtmp = MVT::CloneView( *AU_, index );
1242 MVT::MvTransMv( one, *AUtmp, *z_, mu );
1243 Teuchos::SerialDenseMatrix<int,ScalarType> LUUTAUtmp( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ );
1246 lapack.GETRS(
TRANS, recycleBlocks_, 1, LUUTAUtmp.values(), LUUTAUtmp.stride(), &(*ipiv_)[0], mu.values(), mu.stride(), &info );
1248 "Belos::RCGSolMgr::solve(): LAPACK GETRS failed to compute a solution.");
1252 Teuchos::RCP<MV> Ptmp = MVT::CloneViewNonConst( *P_, index );
1253 MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp);
1254 MVT::MvTimesMatAddMv( -one, *U_, mu, one, *Ptmp );
1259 Teuchos::RCP<MV> Ptmp = MVT::CloneViewNonConst( *P_, index );
1260 MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp);
1267 index.resize( numBlocks_+1 );
1268 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1269 newstate.
P = MVT::CloneViewNonConst( *P_, index );
1270 index.resize( recycleBlocks_ );
1271 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1272 newstate.
U = MVT::CloneViewNonConst( *U_, index );
1273 index.resize( recycleBlocks_ );
1274 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1275 newstate.
AU = MVT::CloneViewNonConst( *AU_, index );
1276 newstate.
Alpha = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Alpha_, numBlocks_, 1 ) );
1277 newstate.
Beta = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Beta_, numBlocks_, 1 ) );
1278 newstate.
D = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *D_, numBlocks_, 1 ) );
1279 newstate.
Delta = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_, 0, 1 ) );
1280 newstate.
LUUTAU = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ ) );
1286 newstate.
existU = existU_;
1287 newstate.
ipiv = ipiv_;
1290 rcg_iter->initialize(newstate);
1296 rcg_iter->iterate();
1303 if ( convTest_->getStatus() ==
Passed ) {
1312 else if ( maxIterTest_->getStatus() ==
Passed ) {
1314 isConverged =
false;
1322 else if ( rcg_iter->getCurSubspaceDim() == rcg_iter->getMaxSubspaceDim() ) {
1327 if (recycleBlocks_ > 0) {
1331 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, numBlocks_, numBlocks_ );
1332 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, numBlocks_, numBlocks_ );
1333 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1334 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1335 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_, 1 );
1336 Ftmp.putScalar(zero);
1337 Gtmp.putScalar(zero);
1338 for (
int ii=0;ii<numBlocks_;ii++) {
1339 Gtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0));
1341 Gtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1342 Gtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1344 Ftmp(ii,ii) = Dtmp(ii,0);
1348 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, numBlocks_, recycleBlocks_ );
1349 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1352 index.resize( numBlocks_ );
1353 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; }
1354 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1355 index.resize( recycleBlocks_ );
1356 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1357 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1358 MVT::MvTimesMatAddMv( one, *Ptmp, Ytmp, zero, *U1tmp );
1363 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, numBlocks_, recycleBlocks_ );
1364 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1365 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1366 AU1TAU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1369 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, numBlocks_, recycleBlocks_ );
1370 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1371 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1372 AU1TU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1374 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, 1 );
1376 AU1TAPtmp.putScalar(zero);
1378 ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0);
1379 for (
int ii=0; ii<recycleBlocks_; ++ii) {
1380 AU1TAPtmp(ii,0) = Ytmp(numBlocks_-1,ii) * alphatmp;
1388 lastBeta = numBlocks_-1;
1395 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1396 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, numBlocks_ );
1397 AU1TAPtmp.scale(Dtmp(0,0));
1399 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1400 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_+1, 1 );
1401 Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
1402 APTAPtmp.putScalar(zero);
1403 for (
int ii=0; ii<numBlocks_; ii++) {
1404 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0));
1406 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1407 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1412 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1413 Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
1414 Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1415 Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1416 Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1417 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1418 F11.assign(AU1TU1tmp);
1419 F12.putScalar(zero);
1420 F21.putScalar(zero);
1421 F22.putScalar(zero);
1422 for(
int ii=0;ii<numBlocks_;ii++) {
1423 F22(ii,ii) = Dtmp(ii,0);
1427 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1428 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1429 Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
1430 Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1431 Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1432 Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1433 G11.assign(AU1TAU1tmp);
1434 G12.assign(AU1TAPtmp);
1436 for (
int ii=0;ii<recycleBlocks_;++ii)
1437 for (
int jj=0;jj<numBlocks_;++jj)
1438 G21(jj,ii) = G12(ii,jj);
1439 G22.assign(APTAPtmp);
1442 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
1443 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1446 index.resize( numBlocks_ );
1447 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; }
1448 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1449 index.resize( recycleBlocks_ );
1450 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1451 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index );
1452 Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1453 index.resize( recycleBlocks_ );
1454 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1455 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1456 index.resize( recycleBlocks_ );
1457 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1458 Teuchos::RCP<MV> U1Y1tmp = MVT::CloneViewNonConst( *U1Y1_, index );
1459 Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
1460 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1461 MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp );
1462 MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp);
1467 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1468 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1469 AU1TAU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1473 AU1TAPtmp.putScalar(zero);
1474 ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0);
1475 for (
int ii=0; ii<recycleBlocks_; ++ii) {
1476 AU1TAPtmp(ii,0) = Ytmp(numBlocks_+recycleBlocks_-1,ii) * alphatmp;
1480 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1481 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1483 AU1TU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1486 lastp = numBlocks_+1;
1487 lastBeta = numBlocks_;
1493 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1494 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_, 1 );
1495 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1496 Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
1497 APTAPtmp.putScalar(zero);
1498 for (
int ii=0; ii<numBlocks_; ii++) {
1499 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0));
1501 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1502 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1506 Teuchos::SerialDenseMatrix<int,ScalarType> L2tmp( Teuchos::View, *L2_, numBlocks_+1, numBlocks_ );
1507 L2tmp.putScalar(zero);
1508 for(
int ii=0;ii<numBlocks_;ii++) {
1509 L2tmp(ii,ii) = 1./Alphatmp(ii,0);
1510 L2tmp(ii+1,ii) = -1./Alphatmp(ii,0);
1514 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAPtmp( Teuchos::View, *AUTAP_, recycleBlocks_, numBlocks_ );
1515 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1516 Teuchos::SerialDenseMatrix<int,ScalarType> Deltatmp( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_+1 );
1517 Teuchos::SerialDenseMatrix<int,ScalarType> DeltaL2tmp( Teuchos::View, *DeltaL2_, recycleBlocks_, numBlocks_ );
1518 DeltaL2tmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Deltatmp,L2tmp,zero);
1519 AUTAPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,UTAUtmp,DeltaL2tmp,zero);
1522 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1523 Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
1524 Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1525 Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1526 Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1527 F11.assign(UTAUtmp);
1528 F12.putScalar(zero);
1529 F21.putScalar(zero);
1530 F22.putScalar(zero);
1531 for(
int ii=0;ii<numBlocks_;ii++) {
1532 F22(ii,ii) = Dtmp(ii,0);
1536 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1537 Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
1538 Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1539 Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1540 Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1541 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
1542 G11.assign(AUTAUtmp);
1543 G12.assign(AUTAPtmp);
1545 for (
int ii=0;ii<recycleBlocks_;++ii)
1546 for (
int jj=0;jj<numBlocks_;++jj)
1547 G21(jj,ii) = G12(ii,jj);
1548 G22.assign(APTAPtmp);
1551 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
1552 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1555 index.resize( recycleBlocks_ );
1556 for (
int ii=0; ii<(recycleBlocks_); ++ii) { index[ii] = ii; }
1557 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1558 index.resize( numBlocks_ );
1559 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; }
1560 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1561 index.resize( recycleBlocks_ );
1562 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1563 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index );
1564 Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1565 index.resize( recycleBlocks_ );
1566 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1567 Teuchos::RCP<MV> UY1tmp = MVT::CloneViewNonConst( *U1Y1_, index );
1568 Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
1569 index.resize( recycleBlocks_ );
1570 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1571 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1572 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1573 MVT::MvTimesMatAddMv( one, *Utmp, Y1, zero, *UY1tmp );
1574 MVT::MvAddMv(one,*UY1tmp, one, *PY2tmp, *U1tmp);
1579 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1580 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1581 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1582 AU1TAU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1585 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1586 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1587 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1588 AU1TU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1591 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUtmp( Teuchos::View, *AU1TU_, recycleBlocks_, recycleBlocks_ );
1592 AU1TUtmp.assign(UTAUtmp);
1595 dold = Dtmp(numBlocks_-1,0);
1602 lastBeta = numBlocks_-1;
1605 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1606 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_+1, 1 );
1607 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1608 Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
1609 for (
int ii=0; ii<numBlocks_; ii++) {
1610 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0));
1612 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1613 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1617 Teuchos::SerialDenseMatrix<int,ScalarType> L2tmp( Teuchos::View, *L2_, numBlocks_+1, numBlocks_ );
1618 for(
int ii=0;ii<numBlocks_;ii++) {
1619 L2tmp(ii,ii) = 1./Alphatmp(ii,0);
1620 L2tmp(ii+1,ii) = -1./Alphatmp(ii,0);
1625 Teuchos::SerialDenseMatrix<int,ScalarType> DeltaL2( Teuchos::View, *DeltaL2_, recycleBlocks_, numBlocks_ );
1626 Teuchos::SerialDenseMatrix<int,ScalarType> Deltatmp( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_+1 );
1627 DeltaL2.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Deltatmp,L2tmp,zero);
1628 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUDeltaL2( Teuchos::View, *AU1TUDeltaL2_, recycleBlocks_, numBlocks_ );
1629 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUtmp( Teuchos::View, *AU1TU_, recycleBlocks_, recycleBlocks_ );
1630 AU1TUDeltaL2.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,AU1TUtmp,DeltaL2,zero);
1631 Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
1632 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, numBlocks_ );
1633 AU1TAPtmp.putScalar(zero);
1634 AU1TAPtmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Y1,AU1TUDeltaL2,zero);
1635 Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1636 ScalarType val = dold * (-Betatmp(0,0)/Alphatmp(0,0));
1637 for(
int ii=0;ii<recycleBlocks_;ii++) {
1638 AU1TAPtmp(ii,0) += Y2(numBlocks_-1,ii)*val;
1642 Teuchos::SerialDenseMatrix<int,ScalarType> Y1TAU1TU( Teuchos::View, *GY_, recycleBlocks_, recycleBlocks_ );
1643 Y1TAU1TU.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Y1,AU1TUtmp,zero);
1644 AU1TUtmp.assign(Y1TAU1TU);
1647 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1648 Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
1649 Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1650 Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1651 Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1652 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1653 F11.assign(AU1TU1tmp);
1654 for(
int ii=0;ii<numBlocks_;ii++) {
1655 F22(ii,ii) = Dtmp(ii,0);
1659 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1660 Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
1661 Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1662 Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1663 Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1664 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1665 G11.assign(AU1TAU1tmp);
1666 G12.assign(AU1TAPtmp);
1668 for (
int ii=0;ii<recycleBlocks_;++ii)
1669 for (
int jj=0;jj<numBlocks_;++jj)
1670 G21(jj,ii) = G12(ii,jj);
1671 G22.assign(APTAPtmp);
1674 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
1675 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1678 index.resize( numBlocks_ );
1679 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; }
1680 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1681 index.resize( recycleBlocks_ );
1682 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1683 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index );
1684 index.resize( recycleBlocks_ );
1685 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1686 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1687 index.resize( recycleBlocks_ );
1688 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1689 Teuchos::RCP<MV> U1Y1tmp = MVT::CloneViewNonConst( *U1Y1_, index );
1690 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1691 MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp );
1692 MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp);
1697 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1698 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1699 AU1TAU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1702 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1703 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1704 AU1TU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1707 dold = Dtmp(numBlocks_-1,0);
1710 lastp = numBlocks_+1;
1711 lastBeta = numBlocks_;
1721 index[0] = lastp-1; index[1] = lastp;
1722 Teuchos::RCP<const MV> Ptmp2 = MVT::CloneView( *P_, index );
1723 index[0] = 0; index[1] = 1;
1724 MVT::SetBlock(*Ptmp2,index,*P_);
1727 (*Beta_)(0,0) = (*Beta_)(lastBeta,0);
1731 Teuchos::SerialDenseMatrix<int,ScalarType> mu1( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, 0 );
1732 Teuchos::SerialDenseMatrix<int,ScalarType> mu2( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, numBlocks_ );
1737 newstate.
P = Teuchos::null;
1738 index.resize( numBlocks_+1 );
1739 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii+1; }
1740 newstate.
P = MVT::CloneViewNonConst( *P_, index );
1742 newstate.
Beta = Teuchos::null;
1743 newstate.
Beta = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Beta_, numBlocks_, 1, 1, 0 ) );
1745 newstate.
Delta = Teuchos::null;
1746 newstate.
Delta = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_, 0, 1 ) );
1751 rcg_iter->initialize(newstate);
1764 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
1765 "Belos::RCGSolMgr::solve(): Invalid return from RCGIter::iterate().");
1768 catch (
const std::exception &e) {
1769 printer_->stream(
Errors) <<
"Error! Caught std::exception in RCGIter::iterate() at iteration " 1770 << rcg_iter->getNumIters() << std::endl
1771 << e.what() << std::endl;
1777 problem_->setCurrLS();
1781 if ( numRHS2Solve > 0 ) {
1784 problem_->setLSIndex( currIdx );
1787 currIdx.resize( numRHS2Solve );
1793 index.resize(recycleBlocks_);
1794 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1795 MVT::SetBlock(*U1_,index,*U_);
1798 if (numRHS2Solve > 0) {
1800 newstate.
P = Teuchos::null;
1801 newstate.
Ap = Teuchos::null;
1802 newstate.
r = Teuchos::null;
1803 newstate.
z = Teuchos::null;
1804 newstate.
U = Teuchos::null;
1805 newstate.
AU = Teuchos::null;
1806 newstate.
Alpha = Teuchos::null;
1807 newstate.
Beta = Teuchos::null;
1808 newstate.
D = Teuchos::null;
1809 newstate.
Delta = Teuchos::null;
1810 newstate.
LUUTAU = Teuchos::null;
1811 newstate.
ipiv = Teuchos::null;
1812 newstate.
rTz_old = Teuchos::null;
1815 index.resize(recycleBlocks_);
1816 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1817 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1818 index.resize(recycleBlocks_);
1819 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1820 Teuchos::RCP<MV> AUtmp = MVT::CloneViewNonConst( *AU_, index );
1822 problem_->applyOp( *Utmp, *AUtmp );
1824 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1825 MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp );
1827 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
1828 if ( precObj != Teuchos::null ) {
1829 index.resize(recycleBlocks_);
1830 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1831 index.resize(recycleBlocks_);
1832 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1833 Teuchos::RCP<MV> LeftPCAU = MVT::CloneViewNonConst( *U1_, index );
1834 OPT::Apply( *precObj, *AUtmp, *LeftPCAU );
1835 MVT::MvTransMv( one, *AUtmp, *LeftPCAU, AUTAUtmp );
1837 MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp );
1850 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1855 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails) );
1859 numIters_ = maxIterTest_->getNumIters();
1863 using Teuchos::rcp_dynamic_cast;
1866 const std::vector<MagnitudeType>* pTestValues =
1867 rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
1869 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1870 "Belos::RCGSolMgr::solve(): The convergence test's getTestValue() " 1871 "method returned NULL. Please report this bug to the Belos developers.");
1873 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1874 "Belos::RCGSolMgr::solve(): The convergence test's getTestValue() " 1875 "method returned a vector of length zero. Please report this bug to the " 1876 "Belos developers.");
1881 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1891 template<
class ScalarType,
class MV,
class OP>
1893 const Teuchos::SerialDenseMatrix<int,ScalarType>& ,
1894 Teuchos::SerialDenseMatrix<int,ScalarType>& Y) {
1896 int n = F.numCols();
1899 Teuchos::LAPACK<int,ScalarType> lapack;
1902 std::vector<MagnitudeType> w(n);
1905 std::vector<int> iperm(n);
1912 std::vector<ScalarType> work(1);
1916 Teuchos::SerialDenseMatrix<int,ScalarType> F2( Teuchos::Copy, *F_ );
1917 Teuchos::SerialDenseMatrix<int,ScalarType> G2( Teuchos::Copy, *G_ );
1920 lapack.SYGV(itype, jobz, uplo, n, G2.values(), G2.stride(), F2.values(), F2.stride(), &w[0], &work[0], lwork, &info);
1922 "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to query optimal work size.");
1923 lwork = (int)work[0];
1925 lapack.SYGV(itype, jobz, uplo, n, G2.values(), G2.stride(), F2.values(), F2.stride(), &w[0], &work[0], lwork, &info);
1927 "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to compute eigensolutions.");
1931 this->sort(w,n,iperm);
1934 for(
int i=0; i<recycleBlocks_; i++ ) {
1935 for(
int j=0; j<n; j++ ) {
1936 Y(j,i) = G2(j,iperm[i]);
1943 template<
class ScalarType,
class MV,
class OP>
1946 int l, r, j, i, flag;
1975 if (dlist[j] > dlist[j - 1]) j = j + 1;
1977 if (dlist[j - 1] > dK) {
1978 dlist[i - 1] = dlist[j - 1];
1979 iperm[i - 1] = iperm[j - 1];
1992 dlist[r] = dlist[0];
1993 iperm[r] = iperm[0];
2008 template<
class ScalarType,
class MV,
class OP>
2011 std::ostringstream oss;
2012 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< std::ostream > outputStream_
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::ScalarTraits< MagnitudeType > MT
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.
Teuchos::RCP< Teuchos::Time > timerSolve_
RCGSolMgrRecyclingFailure is thrown when any problem occurs in using/creating the recycling subspace...
static const bool scalarTypeIsSupported
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
MagnitudeType achievedTol_
Tolerance achieved by the last solve() invocation.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > D_
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
int getNumIters() const override
Get the iteration count for the most recent call to solve().
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > AUTAU_
virtual ~RCGSolMgr()
Destructor.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > GY_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > L2_
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...
int numIters_
Number of iterations taken by the last solve() invocation.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > AUTAP_
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
An implementation of StatusTestResNorm using a family of residual norms.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > UTAU_
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > rTz_old_
static const double convTol
Default convergence tolerance.
Belos::StatusTest class for specifying a maximum number of iterations.
Teuchos::RCP< Teuchos::ParameterList > params_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > AU1TU1_
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.
RCGSolMgrLinearProblemFailure(const std::string &what_arg)
RCGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i.e.
MagnitudeType convtol_
Convergence tolerance (read from parameter list).
Traits class which defines basic operations on multivectors.
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
Belos::StatusTest for logically combining several status tests.
Structure to contain pointers to RCGIter state variables.
Belos concrete class for performing the RCG iteration.
MultiVecTraits< ScalarType, MV > MVT
int maxIters_
Maximum iteration count (read from parameter list).
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. ...
Teuchos::RCP< MV > z
The current preconditioned residual.
Teuchos::RCP< MV > Ap
A times current search vector.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Alpha_
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...
OperatorTraits< ScalarType, MV, OP > OPT
ReturnType
Whether the Belos solve converged for all linear systems.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > APTAP_
Teuchos::RCP< std::vector< int > > ipiv_
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > convTest_
Teuchos::RCP< OutputManager< ScalarType > > printer_
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.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Beta_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Y_
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
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)
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > LUUTAU_
A class for extending the status testing capabilities of Belos via logical combinations.
Details::SolverManagerRequiresRealLapack< ScalarType, MV, OP, scalarTypeIsSupported > base_type
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
Teuchos::ScalarTraits< ScalarType > SCT
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.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Delta_
void reset(const ResetType type) override
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...