42 #ifndef BELOS_RCG_SOLMGR_HPP 43 #define BELOS_RCG_SOLMGR_HPP 61 #include "Teuchos_BLAS.hpp" 62 #include "Teuchos_LAPACK.hpp" 63 #include "Teuchos_as.hpp" 64 #ifdef BELOS_TEUCHOS_TIME_MONITOR 65 #include "Teuchos_TimeMonitor.hpp" 147 template<
class ScalarType,
class MV,
class OP,
148 const bool supportsScalarType =
150 ! Teuchos::ScalarTraits<ScalarType>::isComplex>
153 Belos::Details::LapackSupportsScalar<ScalarType>::value &&
154 ! Teuchos::ScalarTraits<ScalarType>::isComplex>
156 static const bool scalarTypeIsSupported =
158 ! Teuchos::ScalarTraits<ScalarType>::isComplex;
167 const Teuchos::RCP<Teuchos::ParameterList> &pl) :
173 Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const override {
181 template<
class ScalarType,
class MV,
class OP>
187 typedef Teuchos::ScalarTraits<ScalarType> SCT;
188 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
189 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
225 const Teuchos::RCP<Teuchos::ParameterList> &pl );
231 Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const override {
244 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters()
const override;
254 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
255 return Teuchos::tuple(timerSolve_);
283 void setParameters(
const Teuchos::RCP<Teuchos::ParameterList> ¶ms )
override;
296 if ((type &
Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem();
329 std::string description()
const override;
340 void getHarmonicVecs(
const Teuchos::SerialDenseMatrix<int,ScalarType> &F,
341 const Teuchos::SerialDenseMatrix<int,ScalarType> &G,
342 Teuchos::SerialDenseMatrix<int,ScalarType>& Y);
345 void sort(std::vector<ScalarType>& dlist,
int n, std::vector<int>& iperm);
348 void initializeStateStorage();
351 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
354 Teuchos::RCP<OutputManager<ScalarType> > printer_;
355 Teuchos::RCP<std::ostream> outputStream_;
358 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
359 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
360 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
361 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
364 Teuchos::RCP<Teuchos::ParameterList> params_;
367 static constexpr
int maxIters_default_ = 1000;
368 static constexpr
int blockSize_default_ = 1;
369 static constexpr
int numBlocks_default_ = 25;
370 static constexpr
int recycleBlocks_default_ = 3;
371 static constexpr
bool showMaxResNormOnly_default_ =
false;
374 static constexpr
int outputFreq_default_ = -1;
375 static constexpr
const char * label_default_ =
"Belos";
376 static constexpr std::ostream * outputStream_default_ = &std::cout;
383 MagnitudeType convtol_;
389 MagnitudeType achievedTol_;
397 int numBlocks_, recycleBlocks_;
398 bool showMaxResNormOnly_;
399 int verbosity_, outputStyle_, outputFreq_;
408 Teuchos::RCP<MV> Ap_;
423 Teuchos::RCP<MV> U_, AU_;
426 Teuchos::RCP<MV> U1_;
429 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Alpha_;
430 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Beta_;
431 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > D_;
434 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Delta_;
437 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > UTAU_;
440 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > LUUTAU_;
443 Teuchos::RCP<std::vector<int> > ipiv_;
446 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AUTAU_;
449 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > rTz_old_;
452 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > F_,G_,Y_;
455 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > L2_,DeltaL2_,AU1TUDeltaL2_;
456 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AU1TAU1_, AU1TU1_, AU1TAP_;
457 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > FY_,GY_;
458 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > APTAP_;
459 Teuchos::RCP<MV> U1Y1_, PY2_;
460 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AUTAP_, AU1TU_;
466 Teuchos::RCP<Teuchos::Time> timerSolve_;
474 template<
class ScalarType,
class MV,
class OP>
483 template<
class ScalarType,
class MV,
class OP>
486 const Teuchos::RCP<Teuchos::ParameterList> &pl ) :
492 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument,
"Problem not given to solver manager.");
495 if ( !is_null(pl) ) {
501 template<
class ScalarType,
class MV,
class OP>
504 outputStream_ = Teuchos::rcp(outputStream_default_,
false);
506 maxIters_ = maxIters_default_;
507 numBlocks_ = numBlocks_default_;
508 recycleBlocks_ = recycleBlocks_default_;
509 verbosity_ = verbosity_default_;
510 outputStyle_= outputStyle_default_;
511 outputFreq_= outputFreq_default_;
512 showMaxResNormOnly_ = showMaxResNormOnly_default_;
513 label_ = label_default_;
524 Alpha_ = Teuchos::null;
525 Beta_ = Teuchos::null;
527 Delta_ = Teuchos::null;
528 UTAU_ = Teuchos::null;
529 LUUTAU_ = Teuchos::null;
530 ipiv_ = Teuchos::null;
531 AUTAU_ = Teuchos::null;
532 rTz_old_ = Teuchos::null;
537 DeltaL2_ = Teuchos::null;
538 AU1TUDeltaL2_ = Teuchos::null;
539 AU1TAU1_ = Teuchos::null;
540 AU1TU1_ = Teuchos::null;
541 AU1TAP_ = Teuchos::null;
544 APTAP_ = Teuchos::null;
545 U1Y1_ = Teuchos::null;
546 PY2_ = Teuchos::null;
547 AUTAP_ = Teuchos::null;
548 AU1TU_ = Teuchos::null;
552 template<
class ScalarType,
class MV,
class OP>
556 if (params_ == Teuchos::null) {
557 params_ = Teuchos::rcp(
new Teuchos::ParameterList(*getValidParameters()) );
560 params->validateParameters(*getValidParameters());
564 if (params->isParameter(
"Maximum Iterations")) {
565 maxIters_ = params->get(
"Maximum Iterations",maxIters_default_);
568 params_->set(
"Maximum Iterations", maxIters_);
569 if (maxIterTest_!=Teuchos::null)
570 maxIterTest_->setMaxIters( maxIters_ );
574 if (params->isParameter(
"Num Blocks")) {
575 numBlocks_ = params->get(
"Num Blocks",numBlocks_default_);
576 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
577 "Belos::RCGSolMgr: \"Num Blocks\" must be strictly positive.");
580 params_->set(
"Num Blocks", numBlocks_);
584 if (params->isParameter(
"Num Recycled Blocks")) {
585 recycleBlocks_ = params->get(
"Num Recycled Blocks",recycleBlocks_default_);
586 TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks_ <= 0, std::invalid_argument,
587 "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be strictly positive.");
589 TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks_ >= numBlocks_, std::invalid_argument,
590 "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be less than \"Num Blocks\".");
593 params_->set(
"Num Recycled Blocks", recycleBlocks_);
597 if (params->isParameter(
"Timer Label")) {
598 std::string tempLabel = params->get(
"Timer Label", label_default_);
601 if (tempLabel != label_) {
603 params_->set(
"Timer Label", label_);
604 std::string solveLabel = label_ +
": RCGSolMgr total solve time";
605 #ifdef BELOS_TEUCHOS_TIME_MONITOR 606 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
612 if (params->isParameter(
"Verbosity")) {
613 if (Teuchos::isParameterType<int>(*params,
"Verbosity")) {
614 verbosity_ = params->get(
"Verbosity", verbosity_default_);
616 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,
"Verbosity");
620 params_->set(
"Verbosity", verbosity_);
621 if (printer_ != Teuchos::null)
622 printer_->setVerbosity(verbosity_);
626 if (params->isParameter(
"Output Style")) {
627 if (Teuchos::isParameterType<int>(*params,
"Output Style")) {
628 outputStyle_ = params->get(
"Output Style", outputStyle_default_);
630 outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,
"Output Style");
634 params_->set(
"Output Style", outputStyle_);
635 outputTest_ = Teuchos::null;
639 if (params->isParameter(
"Output Stream")) {
640 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,
"Output Stream");
643 params_->set(
"Output Stream", outputStream_);
644 if (printer_ != Teuchos::null)
645 printer_->setOStream( outputStream_ );
650 if (params->isParameter(
"Output Frequency")) {
651 outputFreq_ = params->get(
"Output Frequency", outputFreq_default_);
655 params_->set(
"Output Frequency", outputFreq_);
656 if (outputTest_ != Teuchos::null)
657 outputTest_->setOutputFrequency( outputFreq_ );
661 if (printer_ == Teuchos::null) {
670 if (params->isParameter(
"Convergence Tolerance")) {
671 if (params->isType<MagnitudeType> (
"Convergence Tolerance")) {
672 convtol_ = params->get (
"Convergence Tolerance",
680 params_->set(
"Convergence Tolerance", convtol_);
681 if (convTest_ != Teuchos::null)
685 if (params->isParameter(
"Show Maximum Residual Norm Only")) {
686 showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,
"Show Maximum Residual Norm Only");
689 params_->set(
"Show Maximum Residual Norm Only", showMaxResNormOnly_);
690 if (convTest_ != Teuchos::null)
691 convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
697 if (maxIterTest_ == Teuchos::null)
701 if (convTest_ == Teuchos::null)
702 convTest_ = Teuchos::rcp(
new StatusTestResNorm_t( convtol_, 1 ) );
704 if (sTest_ == Teuchos::null)
705 sTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
707 if (outputTest_ == Teuchos::null) {
715 std::string solverDesc =
" Recycling CG ";
716 outputTest_->setSolverDesc( solverDesc );
720 if (timerSolve_ == Teuchos::null) {
721 std::string solveLabel = label_ +
": RCGSolMgr total solve time";
722 #ifdef BELOS_TEUCHOS_TIME_MONITOR 723 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
732 template<
class ScalarType,
class MV,
class OP>
733 Teuchos::RCP<const Teuchos::ParameterList>
736 static Teuchos::RCP<const Teuchos::ParameterList> validPL;
739 if(is_null(validPL)) {
740 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
742 "The relative residual tolerance that needs to be achieved by the\n" 743 "iterative solver in order for the linear system to be declared converged.");
744 pl->set(
"Maximum Iterations", static_cast<int>(maxIters_default_),
745 "The maximum number of block iterations allowed for each\n" 746 "set of RHS solved.");
747 pl->set(
"Block Size", static_cast<int>(blockSize_default_),
748 "Block Size Parameter -- currently must be 1 for RCG");
749 pl->set(
"Num Blocks", static_cast<int>(numBlocks_default_),
750 "The length of a cycle (and this max number of search vectors kept)\n");
751 pl->set(
"Num Recycled Blocks", static_cast<int>(recycleBlocks_default_),
752 "The number of vectors in the recycle subspace.");
753 pl->set(
"Verbosity", static_cast<int>(verbosity_default_),
754 "What type(s) of solver information should be outputted\n" 755 "to the output stream.");
756 pl->set(
"Output Style", static_cast<int>(outputStyle_default_),
757 "What style is used for the solver information outputted\n" 758 "to the output stream.");
759 pl->set(
"Output Frequency", static_cast<int>(outputFreq_default_),
760 "How often convergence information should be outputted\n" 761 "to the output stream.");
762 pl->set(
"Output Stream", Teuchos::rcp(outputStream_default_,
false),
763 "A reference-counted pointer to the output stream where all\n" 764 "solver output is sent.");
765 pl->set(
"Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
766 "When convergence information is printed, only show the maximum\n" 767 "relative residual norm when the block size is greater than one.");
768 pl->set(
"Timer Label", static_cast<const char *>(label_default_),
769 "The string to use as a prefix for the timer labels.");
776 template<
class ScalarType,
class MV,
class OP>
780 Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
781 if (rhsMV == Teuchos::null) {
788 TEUCHOS_TEST_FOR_EXCEPTION(static_cast<ptrdiff_t>(numBlocks_) > MVT::GetGlobalLength(*rhsMV),std::invalid_argument,
789 "Belos::RCGSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!");
792 if (P_ == Teuchos::null) {
793 P_ = MVT::Clone( *rhsMV, numBlocks_+2 );
797 if (MVT::GetNumberVecs(*P_) < numBlocks_+2) {
798 Teuchos::RCP<const MV> tmp = P_;
799 P_ = MVT::Clone( *tmp, numBlocks_+2 );
804 if (Ap_ == Teuchos::null)
805 Ap_ = MVT::Clone( *rhsMV, 1 );
808 if (r_ == Teuchos::null)
809 r_ = MVT::Clone( *rhsMV, 1 );
812 if (z_ == Teuchos::null)
813 z_ = MVT::Clone( *rhsMV, 1 );
816 if (U_ == Teuchos::null) {
817 U_ = MVT::Clone( *rhsMV, recycleBlocks_ );
821 if (MVT::GetNumberVecs(*U_) < recycleBlocks_) {
822 Teuchos::RCP<const MV> tmp = U_;
823 U_ = MVT::Clone( *tmp, recycleBlocks_ );
828 if (AU_ == Teuchos::null) {
829 AU_ = MVT::Clone( *rhsMV, recycleBlocks_ );
833 if (MVT::GetNumberVecs(*AU_) < recycleBlocks_) {
834 Teuchos::RCP<const MV> tmp = AU_;
835 AU_ = MVT::Clone( *tmp, recycleBlocks_ );
840 if (U1_ == Teuchos::null) {
841 U1_ = MVT::Clone( *rhsMV, recycleBlocks_ );
845 if (MVT::GetNumberVecs(*U1_) < recycleBlocks_) {
846 Teuchos::RCP<const MV> tmp = U1_;
847 U1_ = MVT::Clone( *tmp, recycleBlocks_ );
852 if (Alpha_ == Teuchos::null)
853 Alpha_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_, 1 ) );
855 if ( (Alpha_->numRows() != numBlocks_) || (Alpha_->numCols() != 1) )
856 Alpha_->reshape( numBlocks_, 1 );
860 if (Beta_ == Teuchos::null)
861 Beta_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + 1, 1 ) );
863 if ( (Beta_->numRows() != (numBlocks_+1)) || (Beta_->numCols() != 1) )
864 Beta_->reshape( numBlocks_ + 1, 1 );
868 if (D_ == Teuchos::null)
869 D_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ , 1 ) );
871 if ( (D_->numRows() != numBlocks_) || (D_->numCols() != 1) )
872 D_->reshape( numBlocks_, 1 );
876 if (Delta_ == Teuchos::null)
877 Delta_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ + 1 ) );
879 if ( (Delta_->numRows() != recycleBlocks_) || (Delta_->numCols() != (numBlocks_ + 1)) )
880 Delta_->reshape( recycleBlocks_, numBlocks_ + 1 );
884 if (UTAU_ == Teuchos::null)
885 UTAU_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
887 if ( (UTAU_->numRows() != recycleBlocks_) || (UTAU_->numCols() != recycleBlocks_) )
888 UTAU_->reshape( recycleBlocks_, recycleBlocks_ );
892 if (LUUTAU_ == Teuchos::null)
893 LUUTAU_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
895 if ( (LUUTAU_->numRows() != recycleBlocks_) || (LUUTAU_->numCols() != recycleBlocks_) )
896 LUUTAU_->reshape( recycleBlocks_, recycleBlocks_ );
900 if (ipiv_ == Teuchos::null)
901 ipiv_ = Teuchos::rcp(
new std::vector<int>(recycleBlocks_) );
903 if ( (
int)ipiv_->size() != recycleBlocks_ )
904 ipiv_->resize(recycleBlocks_);
908 if (AUTAU_ == Teuchos::null)
909 AUTAU_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
911 if ( (AUTAU_->numRows() != recycleBlocks_) || (AUTAU_->numCols() != recycleBlocks_) )
912 AUTAU_->reshape( recycleBlocks_, recycleBlocks_ );
916 if (rTz_old_ == Teuchos::null)
917 rTz_old_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( 1, 1 ) );
919 if ( (rTz_old_->numRows() != 1) || (rTz_old_->numCols() != 1) )
920 rTz_old_->reshape( 1, 1 );
924 if (F_ == Teuchos::null)
925 F_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ) );
927 if ( (F_->numRows() != (numBlocks_+recycleBlocks_)) || (F_->numCols() != numBlocks_+recycleBlocks_) )
928 F_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
932 if (G_ == Teuchos::null)
933 G_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ) );
935 if ( (G_->numRows() != (numBlocks_+recycleBlocks_)) || (G_->numCols() != numBlocks_+recycleBlocks_) )
936 G_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
940 if (Y_ == Teuchos::null)
941 Y_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, recycleBlocks_ ) );
943 if ( (Y_->numRows() != (numBlocks_+recycleBlocks_)) || (Y_->numCols() != numBlocks_+recycleBlocks_) )
944 Y_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
948 if (L2_ == Teuchos::null)
949 L2_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+1, numBlocks_ ) );
951 if ( (L2_->numRows() != (numBlocks_+1)) || (L2_->numCols() != numBlocks_) )
952 L2_->reshape( numBlocks_+1, numBlocks_ );
956 if (DeltaL2_ == Teuchos::null)
957 DeltaL2_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
959 if ( (DeltaL2_->numRows() != (recycleBlocks_)) || (DeltaL2_->numCols() != (numBlocks_) ) )
960 DeltaL2_->reshape( recycleBlocks_, numBlocks_ );
964 if (AU1TUDeltaL2_ == Teuchos::null)
965 AU1TUDeltaL2_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
967 if ( (AU1TUDeltaL2_->numRows() != (recycleBlocks_)) || (AU1TUDeltaL2_->numCols() != (numBlocks_) ) )
968 AU1TUDeltaL2_->reshape( recycleBlocks_, numBlocks_ );
972 if (AU1TAU1_ == Teuchos::null)
973 AU1TAU1_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
975 if ( (AU1TAU1_->numRows() != (recycleBlocks_)) || (AU1TAU1_->numCols() != (recycleBlocks_) ) )
976 AU1TAU1_->reshape( recycleBlocks_, recycleBlocks_ );
980 if (GY_ == Teuchos::null)
981 GY_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + recycleBlocks_, recycleBlocks_ ) );
983 if ( (GY_->numRows() != (numBlocks_ + recycleBlocks_)) || (GY_->numCols() != (recycleBlocks_) ) )
984 GY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ );
988 if (AU1TU1_ == Teuchos::null)
989 AU1TU1_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
991 if ( (AU1TU1_->numRows() != (recycleBlocks_)) || (AU1TU1_->numCols() != (recycleBlocks_) ) )
992 AU1TU1_->reshape( recycleBlocks_, recycleBlocks_ );
996 if (FY_ == Teuchos::null)
997 FY_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + recycleBlocks_, recycleBlocks_ ) );
999 if ( (FY_->numRows() != (numBlocks_ + recycleBlocks_)) || (FY_->numCols() != (recycleBlocks_) ) )
1000 FY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ );
1004 if (AU1TAP_ == Teuchos::null)
1005 AU1TAP_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
1007 if ( (AU1TAP_->numRows() != (recycleBlocks_)) || (AU1TAP_->numCols() != (numBlocks_) ) )
1008 AU1TAP_->reshape( recycleBlocks_, numBlocks_ );
1012 if (APTAP_ == Teuchos::null)
1013 APTAP_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_, numBlocks_ ) );
1015 if ( (APTAP_->numRows() != (numBlocks_)) || (APTAP_->numCols() != (numBlocks_) ) )
1016 APTAP_->reshape( numBlocks_, numBlocks_ );
1020 if (U1Y1_ == Teuchos::null) {
1021 U1Y1_ = MVT::Clone( *rhsMV, recycleBlocks_ );
1025 if (MVT::GetNumberVecs(*U1Y1_) < recycleBlocks_) {
1026 Teuchos::RCP<const MV> tmp = U1Y1_;
1027 U1Y1_ = MVT::Clone( *tmp, recycleBlocks_ );
1032 if (PY2_ == Teuchos::null) {
1033 PY2_ = MVT::Clone( *rhsMV, recycleBlocks_ );
1037 if (MVT::GetNumberVecs(*PY2_) < recycleBlocks_) {
1038 Teuchos::RCP<const MV> tmp = PY2_;
1039 PY2_ = MVT::Clone( *tmp, recycleBlocks_ );
1044 if (AUTAP_ == Teuchos::null)
1045 AUTAP_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
1047 if ( (AUTAP_->numRows() != (recycleBlocks_)) || (AUTAP_->numCols() != (numBlocks_) ) )
1048 AUTAP_->reshape( recycleBlocks_, numBlocks_ );
1052 if (AU1TU_ == Teuchos::null)
1053 AU1TU_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
1055 if ( (AU1TU_->numRows() != (recycleBlocks_)) || (AU1TU_->numCols() != (recycleBlocks_) ) )
1056 AU1TU_->reshape( recycleBlocks_, recycleBlocks_ );
1063 template<
class ScalarType,
class MV,
class OP>
1066 Teuchos::LAPACK<int,ScalarType> lapack;
1067 std::vector<int> index(recycleBlocks_);
1068 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1069 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1078 setParameters(Teuchos::parameterList(*getValidParameters()));
1082 "Belos::RCGSolMgr::solve(): Linear problem is not a valid object.");
1084 "Belos::RCGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
1085 TEUCHOS_TEST_FOR_EXCEPTION((problem_->getLeftPrec() != Teuchos::null)&&(problem_->getRightPrec() != Teuchos::null),
1087 "Belos::RCGSolMgr::solve(): RCG does not support split preconditioning, only set left or right preconditioner.");
1090 Teuchos::RCP<OP> precObj;
1091 if (problem_->getLeftPrec() != Teuchos::null) {
1092 precObj = Teuchos::rcp_const_cast<OP>(problem_->getLeftPrec());
1094 else if (problem_->getRightPrec() != Teuchos::null) {
1095 precObj = Teuchos::rcp_const_cast<OP>(problem_->getRightPrec());
1099 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1100 std::vector<int> currIdx(1);
1104 problem_->setLSIndex( currIdx );
1107 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1108 if (numBlocks_ > dim) {
1109 numBlocks_ = Teuchos::asSafe<int>(dim);
1110 params_->set(
"Num Blocks", numBlocks_);
1112 "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl <<
1113 " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl;
1117 initializeStateStorage();
1120 Teuchos::ParameterList plist;
1121 plist.set(
"Num Blocks",numBlocks_);
1122 plist.set(
"Recycled Blocks",recycleBlocks_);
1125 outputTest_->reset();
1128 bool isConverged =
true;
1132 index.resize(recycleBlocks_);
1133 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1134 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1135 index.resize(recycleBlocks_);
1136 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1137 Teuchos::RCP<MV> AUtmp = MVT::CloneViewNonConst( *AU_, index );
1139 problem_->applyOp( *Utmp, *AUtmp );
1141 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1142 MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp );
1144 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
1145 if ( precObj != Teuchos::null ) {
1146 index.resize(recycleBlocks_);
1147 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1148 index.resize(recycleBlocks_);
1149 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1150 Teuchos::RCP<MV> PCAU = MVT::CloneViewNonConst( *U1_, index );
1151 OPT::Apply( *precObj, *AUtmp, *PCAU );
1152 MVT::MvTransMv( one, *AUtmp, *PCAU, AUTAUtmp );
1154 MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp );
1161 Teuchos::RCP<RCGIter<ScalarType,MV,OP> > rcg_iter;
1166 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1167 Teuchos::TimeMonitor slvtimer(*timerSolve_);
1170 while ( numRHS2Solve > 0 ) {
1173 if (printer_->isVerbosity(
Debug ) ) {
1174 if (existU_) printer_->print(
Debug,
"Using recycle space generated from previous call to solve()." );
1175 else printer_->print(
Debug,
"No recycle space exists." );
1179 rcg_iter->resetNumIters();
1182 rcg_iter->setSize( recycleBlocks_, numBlocks_ );
1185 outputTest_->resetNumCalls();
1194 problem_->computeCurrResVec( &*r_ );
1199 Teuchos::SerialDenseMatrix<int,ScalarType> Utr(recycleBlocks_,1);
1200 index.resize(recycleBlocks_);
1201 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1202 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1203 MVT::MvTransMv( one, *Utmp, *r_, Utr );
1204 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1205 Teuchos::SerialDenseMatrix<int,ScalarType> LUUTAUtmp( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ );
1206 LUUTAUtmp.assign(UTAUtmp);
1208 lapack.GESV(recycleBlocks_, 1, LUUTAUtmp.values(), LUUTAUtmp.stride(), &(*ipiv_)[0], Utr.values(), Utr.stride(), &info);
1210 "Belos::RCGSolMgr::solve(): LAPACK GESV failed to compute a solution.");
1213 MVT::MvTimesMatAddMv( one, *Utmp, Utr, one, *problem_->getCurrLHSVec() );
1216 index.resize(recycleBlocks_);
1217 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1218 Teuchos::RCP<const MV> AUtmp = MVT::CloneView( *AU_, index );
1219 MVT::MvTimesMatAddMv( -one, *AUtmp, Utr, one, *r_ );
1222 if ( precObj != Teuchos::null ) {
1223 OPT::Apply( *precObj, *r_, *z_ );
1229 MVT::MvTransMv( one, *r_, *z_, *rTz_old_ );
1233 Teuchos::SerialDenseMatrix<int,ScalarType> mu( Teuchos::View, *Delta_, recycleBlocks_, 1);
1234 index.resize(recycleBlocks_);
1235 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1236 Teuchos::RCP<const MV> AUtmp = MVT::CloneView( *AU_, index );
1237 MVT::MvTransMv( one, *AUtmp, *z_, mu );
1238 Teuchos::SerialDenseMatrix<int,ScalarType> LUUTAUtmp( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ );
1241 lapack.GETRS(
TRANS, recycleBlocks_, 1, LUUTAUtmp.values(), LUUTAUtmp.stride(), &(*ipiv_)[0], mu.values(), mu.stride(), &info );
1243 "Belos::RCGSolMgr::solve(): LAPACK GETRS failed to compute a solution.");
1247 Teuchos::RCP<MV> Ptmp = MVT::CloneViewNonConst( *P_, index );
1248 MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp);
1249 MVT::MvTimesMatAddMv( -one, *U_, mu, one, *Ptmp );
1254 Teuchos::RCP<MV> Ptmp = MVT::CloneViewNonConst( *P_, index );
1255 MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp);
1262 index.resize( numBlocks_+1 );
1263 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1264 newstate.
P = MVT::CloneViewNonConst( *P_, index );
1265 index.resize( recycleBlocks_ );
1266 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1267 newstate.
U = MVT::CloneViewNonConst( *U_, index );
1268 index.resize( recycleBlocks_ );
1269 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1270 newstate.
AU = MVT::CloneViewNonConst( *AU_, index );
1271 newstate.
Alpha = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Alpha_, numBlocks_, 1 ) );
1272 newstate.
Beta = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Beta_, numBlocks_, 1 ) );
1273 newstate.
D = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *D_, numBlocks_, 1 ) );
1274 newstate.
Delta = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_, 0, 1 ) );
1275 newstate.
LUUTAU = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ ) );
1281 newstate.
existU = existU_;
1282 newstate.
ipiv = ipiv_;
1285 rcg_iter->initialize(newstate);
1291 rcg_iter->iterate();
1298 if ( convTest_->getStatus() ==
Passed ) {
1307 else if ( maxIterTest_->getStatus() ==
Passed ) {
1309 isConverged =
false;
1317 else if ( rcg_iter->getCurSubspaceDim() == rcg_iter->getMaxSubspaceDim() ) {
1322 if (recycleBlocks_ > 0) {
1326 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, numBlocks_, numBlocks_ );
1327 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, numBlocks_, numBlocks_ );
1328 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1329 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1330 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_, 1 );
1331 Ftmp.putScalar(zero);
1332 Gtmp.putScalar(zero);
1333 for (
int ii=0;ii<numBlocks_;ii++) {
1334 Gtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0));
1336 Gtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1337 Gtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1339 Ftmp(ii,ii) = Dtmp(ii,0);
1343 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, numBlocks_, recycleBlocks_ );
1344 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1347 index.resize( numBlocks_ );
1348 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; }
1349 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1350 index.resize( recycleBlocks_ );
1351 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1352 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1353 MVT::MvTimesMatAddMv( one, *Ptmp, Ytmp, zero, *U1tmp );
1358 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, numBlocks_, recycleBlocks_ );
1359 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1360 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1361 AU1TAU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1364 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, numBlocks_, recycleBlocks_ );
1365 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1366 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1367 AU1TU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1369 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, 1 );
1371 AU1TAPtmp.putScalar(zero);
1373 ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0);
1374 for (
int ii=0; ii<recycleBlocks_; ++ii) {
1375 AU1TAPtmp(ii,0) = Ytmp(numBlocks_-1,ii) * alphatmp;
1383 lastBeta = numBlocks_-1;
1390 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1391 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, numBlocks_ );
1392 AU1TAPtmp.scale(Dtmp(0,0));
1394 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1395 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_+1, 1 );
1396 Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
1397 APTAPtmp.putScalar(zero);
1398 for (
int ii=0; ii<numBlocks_; ii++) {
1399 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0));
1401 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1402 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1407 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1408 Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
1409 Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1410 Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1411 Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1412 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1413 F11.assign(AU1TU1tmp);
1414 F12.putScalar(zero);
1415 F21.putScalar(zero);
1416 F22.putScalar(zero);
1417 for(
int ii=0;ii<numBlocks_;ii++) {
1418 F22(ii,ii) = Dtmp(ii,0);
1422 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1423 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1424 Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
1425 Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1426 Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1427 Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1428 G11.assign(AU1TAU1tmp);
1429 G12.assign(AU1TAPtmp);
1431 for (
int ii=0;ii<recycleBlocks_;++ii)
1432 for (
int jj=0;jj<numBlocks_;++jj)
1433 G21(jj,ii) = G12(ii,jj);
1434 G22.assign(APTAPtmp);
1437 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
1438 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1441 index.resize( numBlocks_ );
1442 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; }
1443 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1444 index.resize( recycleBlocks_ );
1445 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1446 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index );
1447 Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1448 index.resize( recycleBlocks_ );
1449 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1450 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1451 index.resize( recycleBlocks_ );
1452 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1453 Teuchos::RCP<MV> U1Y1tmp = MVT::CloneViewNonConst( *U1Y1_, index );
1454 Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
1455 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1456 MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp );
1457 MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp);
1462 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1463 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1464 AU1TAU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1468 AU1TAPtmp.putScalar(zero);
1469 ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0);
1470 for (
int ii=0; ii<recycleBlocks_; ++ii) {
1471 AU1TAPtmp(ii,0) = Ytmp(numBlocks_+recycleBlocks_-1,ii) * alphatmp;
1475 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1476 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1478 AU1TU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1481 lastp = numBlocks_+1;
1482 lastBeta = numBlocks_;
1488 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1489 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_, 1 );
1490 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1491 Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
1492 APTAPtmp.putScalar(zero);
1493 for (
int ii=0; ii<numBlocks_; ii++) {
1494 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0));
1496 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1497 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1501 Teuchos::SerialDenseMatrix<int,ScalarType> L2tmp( Teuchos::View, *L2_, numBlocks_+1, numBlocks_ );
1502 L2tmp.putScalar(zero);
1503 for(
int ii=0;ii<numBlocks_;ii++) {
1504 L2tmp(ii,ii) = 1./Alphatmp(ii,0);
1505 L2tmp(ii+1,ii) = -1./Alphatmp(ii,0);
1509 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAPtmp( Teuchos::View, *AUTAP_, recycleBlocks_, numBlocks_ );
1510 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1511 Teuchos::SerialDenseMatrix<int,ScalarType> Deltatmp( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_+1 );
1512 Teuchos::SerialDenseMatrix<int,ScalarType> DeltaL2tmp( Teuchos::View, *DeltaL2_, recycleBlocks_, numBlocks_ );
1513 DeltaL2tmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Deltatmp,L2tmp,zero);
1514 AUTAPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,UTAUtmp,DeltaL2tmp,zero);
1517 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1518 Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
1519 Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1520 Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1521 Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1522 F11.assign(UTAUtmp);
1523 F12.putScalar(zero);
1524 F21.putScalar(zero);
1525 F22.putScalar(zero);
1526 for(
int ii=0;ii<numBlocks_;ii++) {
1527 F22(ii,ii) = Dtmp(ii,0);
1531 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1532 Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
1533 Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1534 Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1535 Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1536 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
1537 G11.assign(AUTAUtmp);
1538 G12.assign(AUTAPtmp);
1540 for (
int ii=0;ii<recycleBlocks_;++ii)
1541 for (
int jj=0;jj<numBlocks_;++jj)
1542 G21(jj,ii) = G12(ii,jj);
1543 G22.assign(APTAPtmp);
1546 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
1547 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1550 index.resize( recycleBlocks_ );
1551 for (
int ii=0; ii<(recycleBlocks_); ++ii) { index[ii] = ii; }
1552 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1553 index.resize( numBlocks_ );
1554 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; }
1555 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1556 index.resize( recycleBlocks_ );
1557 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1558 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index );
1559 Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1560 index.resize( recycleBlocks_ );
1561 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1562 Teuchos::RCP<MV> UY1tmp = MVT::CloneViewNonConst( *U1Y1_, index );
1563 Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
1564 index.resize( recycleBlocks_ );
1565 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1566 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1567 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1568 MVT::MvTimesMatAddMv( one, *Utmp, Y1, zero, *UY1tmp );
1569 MVT::MvAddMv(one,*UY1tmp, one, *PY2tmp, *U1tmp);
1574 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1575 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1576 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1577 AU1TAU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1580 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1581 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1582 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1583 AU1TU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1586 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUtmp( Teuchos::View, *AU1TU_, recycleBlocks_, recycleBlocks_ );
1587 AU1TUtmp.assign(UTAUtmp);
1590 dold = Dtmp(numBlocks_-1,0);
1597 lastBeta = numBlocks_-1;
1600 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1601 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_+1, 1 );
1602 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1603 Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
1604 for (
int ii=0; ii<numBlocks_; ii++) {
1605 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0));
1607 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1608 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1612 Teuchos::SerialDenseMatrix<int,ScalarType> L2tmp( Teuchos::View, *L2_, numBlocks_+1, numBlocks_ );
1613 for(
int ii=0;ii<numBlocks_;ii++) {
1614 L2tmp(ii,ii) = 1./Alphatmp(ii,0);
1615 L2tmp(ii+1,ii) = -1./Alphatmp(ii,0);
1620 Teuchos::SerialDenseMatrix<int,ScalarType> DeltaL2( Teuchos::View, *DeltaL2_, recycleBlocks_, numBlocks_ );
1621 Teuchos::SerialDenseMatrix<int,ScalarType> Deltatmp( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_+1 );
1622 DeltaL2.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Deltatmp,L2tmp,zero);
1623 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUDeltaL2( Teuchos::View, *AU1TUDeltaL2_, recycleBlocks_, numBlocks_ );
1624 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUtmp( Teuchos::View, *AU1TU_, recycleBlocks_, recycleBlocks_ );
1625 AU1TUDeltaL2.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,AU1TUtmp,DeltaL2,zero);
1626 Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
1627 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, numBlocks_ );
1628 AU1TAPtmp.putScalar(zero);
1629 AU1TAPtmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Y1,AU1TUDeltaL2,zero);
1630 Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1631 ScalarType val = dold * (-Betatmp(0,0)/Alphatmp(0,0));
1632 for(
int ii=0;ii<recycleBlocks_;ii++) {
1633 AU1TAPtmp(ii,0) += Y2(numBlocks_-1,ii)*val;
1637 Teuchos::SerialDenseMatrix<int,ScalarType> Y1TAU1TU( Teuchos::View, *GY_, recycleBlocks_, recycleBlocks_ );
1638 Y1TAU1TU.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Y1,AU1TUtmp,zero);
1639 AU1TUtmp.assign(Y1TAU1TU);
1642 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1643 Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
1644 Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1645 Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1646 Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1647 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1648 F11.assign(AU1TU1tmp);
1649 for(
int ii=0;ii<numBlocks_;ii++) {
1650 F22(ii,ii) = Dtmp(ii,0);
1654 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1655 Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
1656 Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1657 Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1658 Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1659 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1660 G11.assign(AU1TAU1tmp);
1661 G12.assign(AU1TAPtmp);
1663 for (
int ii=0;ii<recycleBlocks_;++ii)
1664 for (
int jj=0;jj<numBlocks_;++jj)
1665 G21(jj,ii) = G12(ii,jj);
1666 G22.assign(APTAPtmp);
1669 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
1670 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1673 index.resize( numBlocks_ );
1674 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; }
1675 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1676 index.resize( recycleBlocks_ );
1677 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1678 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index );
1679 index.resize( recycleBlocks_ );
1680 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1681 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1682 index.resize( recycleBlocks_ );
1683 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1684 Teuchos::RCP<MV> U1Y1tmp = MVT::CloneViewNonConst( *U1Y1_, index );
1685 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1686 MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp );
1687 MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp);
1692 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1693 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1694 AU1TAU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1697 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1698 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1699 AU1TU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1702 dold = Dtmp(numBlocks_-1,0);
1705 lastp = numBlocks_+1;
1706 lastBeta = numBlocks_;
1716 index[0] = lastp-1; index[1] = lastp;
1717 Teuchos::RCP<const MV> Ptmp2 = MVT::CloneView( *P_, index );
1718 index[0] = 0; index[1] = 1;
1719 MVT::SetBlock(*Ptmp2,index,*P_);
1722 (*Beta_)(0,0) = (*Beta_)(lastBeta,0);
1726 Teuchos::SerialDenseMatrix<int,ScalarType> mu1( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, 0 );
1727 Teuchos::SerialDenseMatrix<int,ScalarType> mu2( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, numBlocks_ );
1732 newstate.
P = Teuchos::null;
1733 index.resize( numBlocks_+1 );
1734 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii+1; }
1735 newstate.
P = MVT::CloneViewNonConst( *P_, index );
1737 newstate.
Beta = Teuchos::null;
1738 newstate.
Beta = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Beta_, numBlocks_, 1, 1, 0 ) );
1740 newstate.
Delta = Teuchos::null;
1741 newstate.
Delta = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_, 0, 1 ) );
1746 rcg_iter->initialize(newstate);
1759 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
1760 "Belos::RCGSolMgr::solve(): Invalid return from RCGIter::iterate().");
1763 catch (
const std::exception &e) {
1764 printer_->stream(
Errors) <<
"Error! Caught std::exception in RCGIter::iterate() at iteration " 1765 << rcg_iter->getNumIters() << std::endl
1766 << e.what() << std::endl;
1772 problem_->setCurrLS();
1776 if ( numRHS2Solve > 0 ) {
1779 problem_->setLSIndex( currIdx );
1782 currIdx.resize( numRHS2Solve );
1788 index.resize(recycleBlocks_);
1789 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1790 MVT::SetBlock(*U1_,index,*U_);
1793 if (numRHS2Solve > 0) {
1795 newstate.
P = Teuchos::null;
1796 newstate.
Ap = Teuchos::null;
1797 newstate.
r = Teuchos::null;
1798 newstate.
z = Teuchos::null;
1799 newstate.
U = Teuchos::null;
1800 newstate.
AU = Teuchos::null;
1801 newstate.
Alpha = Teuchos::null;
1802 newstate.
Beta = Teuchos::null;
1803 newstate.
D = Teuchos::null;
1804 newstate.
Delta = Teuchos::null;
1805 newstate.
LUUTAU = Teuchos::null;
1806 newstate.
ipiv = Teuchos::null;
1807 newstate.
rTz_old = Teuchos::null;
1810 index.resize(recycleBlocks_);
1811 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1812 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1813 index.resize(recycleBlocks_);
1814 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1815 Teuchos::RCP<MV> AUtmp = MVT::CloneViewNonConst( *AU_, index );
1817 problem_->applyOp( *Utmp, *AUtmp );
1819 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1820 MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp );
1822 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
1823 if ( precObj != Teuchos::null ) {
1824 index.resize(recycleBlocks_);
1825 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1826 index.resize(recycleBlocks_);
1827 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1828 Teuchos::RCP<MV> LeftPCAU = MVT::CloneViewNonConst( *U1_, index );
1829 OPT::Apply( *precObj, *AUtmp, *LeftPCAU );
1830 MVT::MvTransMv( one, *AUtmp, *LeftPCAU, AUTAUtmp );
1832 MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp );
1845 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1850 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails) );
1854 numIters_ = maxIterTest_->getNumIters();
1858 using Teuchos::rcp_dynamic_cast;
1861 const std::vector<MagnitudeType>* pTestValues =
1862 rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
1864 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1865 "Belos::RCGSolMgr::solve(): The convergence test's getTestValue() " 1866 "method returned NULL. Please report this bug to the Belos developers.");
1868 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1869 "Belos::RCGSolMgr::solve(): The convergence test's getTestValue() " 1870 "method returned a vector of length zero. Please report this bug to the " 1871 "Belos developers.");
1876 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1886 template<
class ScalarType,
class MV,
class OP>
1888 const Teuchos::SerialDenseMatrix<int,ScalarType>& ,
1889 Teuchos::SerialDenseMatrix<int,ScalarType>& Y) {
1891 int n = F.numCols();
1894 Teuchos::LAPACK<int,ScalarType> lapack;
1897 std::vector<MagnitudeType> w(n);
1900 std::vector<int> iperm(n);
1907 std::vector<ScalarType> work(1);
1911 Teuchos::SerialDenseMatrix<int,ScalarType> F2( Teuchos::Copy, *F_ );
1912 Teuchos::SerialDenseMatrix<int,ScalarType> G2( Teuchos::Copy, *G_ );
1915 lapack.SYGV(itype, jobz, uplo, n, G2.values(), G2.stride(), F2.values(), F2.stride(), &w[0], &work[0], lwork, &info);
1917 "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to query optimal work size.");
1918 lwork = (int)work[0];
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 compute eigensolutions.");
1926 this->sort(w,n,iperm);
1929 for(
int i=0; i<recycleBlocks_; i++ ) {
1930 for(
int j=0; j<n; j++ ) {
1931 Y(j,i) = G2(j,iperm[i]);
1938 template<
class ScalarType,
class MV,
class OP>
1939 void RCGSolMgr<ScalarType,MV,OP,true>::sort(std::vector<ScalarType>& dlist,
int n, std::vector<int>& iperm)
1941 int l, r, j, i, flag;
1970 if (dlist[j] > dlist[j - 1]) j = j + 1;
1972 if (dlist[j - 1] > dK) {
1973 dlist[i - 1] = dlist[j - 1];
1974 iperm[i - 1] = iperm[j - 1];
1987 dlist[r] = dlist[0];
1988 iperm[r] = iperm[0];
2003 template<
class ScalarType,
class MV,
class OP>
2006 std::ostringstream oss;
2007 oss <<
"Belos::RCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<
">";
Teuchos::RCP< MV > r
The current residual.
Collection of types and exceptions used within the Belos solvers.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > rTz_old
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
Teuchos::RCP< MV > P
The current Krylov basis.
Teuchos::RCP< std::vector< int > > ipiv
Data from LU factorization of U^T A U.
Class which manages the output and verbosity of the Belos solvers.
RCGSolMgrRecyclingFailure is thrown when any problem occurs in using/creating the recycling subspace...
int getNumIters() const override
Get the iteration count for the most recent call to solve().
virtual ~RCGSolMgr()
Destructor.
Teuchos::RCP< MV > U
The recycled subspace and its image.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > LUUTAU
The LU factorization of the matrix U^T A U.
Base class for Belos::SolverManager subclasses which normally can only compile with real ScalarType t...
A factory class for generating StatusTestOutput objects.
Implementation of the RCG (Recycling Conjugate Gradient) iterative linear solver. ...
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Beta
This class implements the RCG iteration, where a single-std::vector Krylov subspace is constructed...
An implementation of StatusTestResNorm using a family of residual norms.
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
static const double convTol
Default convergence tolerance.
Belos::StatusTest class for specifying a maximum number of iterations.
bool isLOADetected() const override
Return whether a loss of accuracy was detected by this solver during the most current solve...
A factory class for generating StatusTestOutput objects.
RCGSolMgrLinearProblemFailure(const std::string &what_arg)
RCGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i.e.
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
Structure to contain pointers to RCGIter state variables.
Belos concrete class for performing the RCG iteration.
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.
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
Type traits class that says whether Teuchos::LAPACK has a valid implementation for the given ScalarTy...
ReturnType
Whether the Belos solve converged for all linear systems.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Alpha
Coefficients arising in RCG iteration.
int setTolerance(MagnitudeType tolerance)
Set the value of the tolerance.
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > create(const Teuchos::RCP< OutputManager< ScalarType > > &printer, Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test, int mod, int printStates)
Create the StatusTestOutput object specified by the outputStyle.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
int curDim
The current dimension of the reduction.
RCGSolMgrLAPACKFailure(const std::string &what_arg)
RCGSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
A class for extending the status testing capabilities of Belos via logical combinations.
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.
Class which defines basic traits for the operator type.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > D
Parent class to all Belos exceptions.
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
RCGSolMgrRecyclingFailure(const std::string &what_arg)
Belos header file which uses auto-configuration information to include necessary C++ headers...
RCGSolMgrLAPACKFailure is thrown when a nonzero value is retuned from an LAPACK call.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Delta
Solutions to local least-squares problems.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
void reset(const ResetType type) override
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...