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>
155 static const bool scalarTypeIsSupported =
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_;
353 Teuchos::RCP<OutputManager<ScalarType> > printer_;
354 Teuchos::RCP<std::ostream> outputStream_;
357 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
358 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
359 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
360 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
363 Teuchos::RCP<Teuchos::ParameterList> params_;
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;
388 MagnitudeType convtol_;
394 MagnitudeType achievedTol_;
402 int numBlocks_, recycleBlocks_;
403 bool showMaxResNormOnly_;
404 int verbosity_, outputStyle_, outputFreq_;
413 Teuchos::RCP<MV> Ap_;
428 Teuchos::RCP<MV> U_, AU_;
431 Teuchos::RCP<MV> U1_;
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_;
464 Teuchos::RCP<MV> U1Y1_, PY2_;
465 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AUTAP_, AU1TU_;
471 Teuchos::RCP<Teuchos::Time> timerSolve_;
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>
1944 void RCGSolMgr<ScalarType,MV,OP,true>::sort(std::vector<ScalarType>& dlist,
int n, std::vector<int>& iperm)
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< 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...