42 #ifndef BELOS_RCG_SOLMGR_HPP 43 #define BELOS_RCG_SOLMGR_HPP 64 #include "Teuchos_BLAS.hpp" 65 #include "Teuchos_LAPACK.hpp" 66 #include "Teuchos_as.hpp" 67 #ifdef BELOS_TEUCHOS_TIME_MONITOR 68 #include "Teuchos_TimeMonitor.hpp" 150 template<
class ScalarType,
class MV,
class OP,
151 const bool supportsScalarType =
153 ! Teuchos::ScalarTraits<ScalarType>::isComplex>
156 Belos::Details::LapackSupportsScalar<ScalarType>::value &&
157 ! Teuchos::ScalarTraits<ScalarType>::isComplex>
159 static const bool scalarTypeIsSupported =
161 ! Teuchos::ScalarTraits<ScalarType>::isComplex;
170 const Teuchos::RCP<Teuchos::ParameterList> &pl) :
179 template<
class ScalarType,
class MV,
class OP>
185 typedef Teuchos::ScalarTraits<ScalarType> SCT;
186 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
187 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
223 const Teuchos::RCP<Teuchos::ParameterList> &pl );
237 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters()
const;
247 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
248 return Teuchos::tuple(timerSolve_);
276 void setParameters(
const Teuchos::RCP<Teuchos::ParameterList> ¶ms );
289 if ((type &
Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem();
322 std::string description()
const;
333 void getHarmonicVecs(
const Teuchos::SerialDenseMatrix<int,ScalarType> &F,
334 const Teuchos::SerialDenseMatrix<int,ScalarType> &G,
335 Teuchos::SerialDenseMatrix<int,ScalarType>& Y);
338 void sort(std::vector<ScalarType>& dlist,
int n, std::vector<int>& iperm);
341 void initializeStateStorage();
344 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
347 Teuchos::RCP<OutputManager<ScalarType> > printer_;
348 Teuchos::RCP<std::ostream> outputStream_;
351 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
352 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
353 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
354 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
357 Teuchos::RCP<Teuchos::ParameterList> params_;
360 static const MagnitudeType convtol_default_;
361 static const int maxIters_default_;
362 static const int blockSize_default_;
363 static const int numBlocks_default_;
364 static const int recycleBlocks_default_;
365 static const bool showMaxResNormOnly_default_;
366 static const int verbosity_default_;
367 static const int outputStyle_default_;
368 static const int outputFreq_default_;
369 static const std::string label_default_;
370 static const Teuchos::RCP<std::ostream> outputStream_default_;
377 MagnitudeType convtol_;
383 MagnitudeType achievedTol_;
391 int numBlocks_, recycleBlocks_;
392 bool showMaxResNormOnly_;
393 int verbosity_, outputStyle_, outputFreq_;
402 Teuchos::RCP<MV> Ap_;
417 Teuchos::RCP<MV> U_, AU_;
420 Teuchos::RCP<MV> U1_;
423 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Alpha_;
424 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Beta_;
425 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > D_;
428 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Delta_;
431 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > UTAU_;
434 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > LUUTAU_;
437 Teuchos::RCP<std::vector<int> > ipiv_;
440 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AUTAU_;
443 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > rTz_old_;
446 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > F_,G_,Y_;
449 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > L2_,DeltaL2_,AU1TUDeltaL2_;
450 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AU1TAU1_, AU1TU1_, AU1TAP_;
451 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > FY_,GY_;
452 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > APTAP_;
453 Teuchos::RCP<MV> U1Y1_, PY2_;
454 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AUTAP_, AU1TU_;
460 Teuchos::RCP<Teuchos::Time> timerSolve_;
468 template<
class ScalarType,
class MV,
class OP>
469 const typename RCGSolMgr<ScalarType,MV,OP,true>::MagnitudeType
470 RCGSolMgr<ScalarType,MV,OP,true>::convtol_default_ = 1e-8;
472 template<
class ScalarType,
class MV,
class OP>
473 const int RCGSolMgr<ScalarType,MV,OP,true>::maxIters_default_ = 1000;
475 template<
class ScalarType,
class MV,
class OP>
476 const int RCGSolMgr<ScalarType,MV,OP,true>::numBlocks_default_ = 25;
478 template<
class ScalarType,
class MV,
class OP>
479 const int RCGSolMgr<ScalarType,MV,OP,true>::blockSize_default_ = 1;
481 template<
class ScalarType,
class MV,
class OP>
482 const int RCGSolMgr<ScalarType,MV,OP,true>::recycleBlocks_default_ = 3;
484 template<
class ScalarType,
class MV,
class OP>
485 const bool RCGSolMgr<ScalarType,MV,OP,true>::showMaxResNormOnly_default_ =
false;
487 template<
class ScalarType,
class MV,
class OP>
488 const int RCGSolMgr<ScalarType,MV,OP,true>::verbosity_default_ =
Belos::Errors;
490 template<
class ScalarType,
class MV,
class OP>
491 const int RCGSolMgr<ScalarType,MV,OP,true>::outputStyle_default_ =
Belos::General;
493 template<
class ScalarType,
class MV,
class OP>
494 const int RCGSolMgr<ScalarType,MV,OP,true>::outputFreq_default_ = -1;
496 template<
class ScalarType,
class MV,
class OP>
497 const std::string RCGSolMgr<ScalarType,MV,OP,true>::label_default_ =
"Belos";
499 template<
class ScalarType,
class MV,
class OP>
500 const Teuchos::RCP<std::ostream> RCGSolMgr<ScalarType,MV,OP,true>::outputStream_default_ = Teuchos::rcp(&std::cout,
false);
503 template<
class ScalarType,
class MV,
class OP>
512 template<
class ScalarType,
class MV,
class OP>
515 const Teuchos::RCP<Teuchos::ParameterList> &pl ) :
521 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument,
"Problem not given to solver manager.");
524 if ( !is_null(pl) ) {
530 template<
class ScalarType,
class MV,
class OP>
533 outputStream_ = outputStream_default_;
534 convtol_ = convtol_default_;
535 maxIters_ = maxIters_default_;
536 numBlocks_ = numBlocks_default_;
537 recycleBlocks_ = recycleBlocks_default_;
538 verbosity_ = verbosity_default_;
539 outputStyle_= outputStyle_default_;
540 outputFreq_= outputFreq_default_;
541 showMaxResNormOnly_ = showMaxResNormOnly_default_;
542 label_ = label_default_;
553 Alpha_ = Teuchos::null;
554 Beta_ = Teuchos::null;
556 Delta_ = Teuchos::null;
557 UTAU_ = Teuchos::null;
558 LUUTAU_ = Teuchos::null;
559 ipiv_ = Teuchos::null;
560 AUTAU_ = Teuchos::null;
561 rTz_old_ = Teuchos::null;
566 DeltaL2_ = Teuchos::null;
567 AU1TUDeltaL2_ = Teuchos::null;
568 AU1TAU1_ = Teuchos::null;
569 AU1TU1_ = Teuchos::null;
570 AU1TAP_ = Teuchos::null;
573 APTAP_ = Teuchos::null;
574 U1Y1_ = Teuchos::null;
575 PY2_ = Teuchos::null;
576 AUTAP_ = Teuchos::null;
577 AU1TU_ = Teuchos::null;
581 template<
class ScalarType,
class MV,
class OP>
585 if (params_ == Teuchos::null) {
586 params_ = Teuchos::rcp(
new Teuchos::ParameterList(*getValidParameters()) );
589 params->validateParameters(*getValidParameters());
593 if (params->isParameter(
"Maximum Iterations")) {
594 maxIters_ = params->get(
"Maximum Iterations",maxIters_default_);
597 params_->set(
"Maximum Iterations", maxIters_);
598 if (maxIterTest_!=Teuchos::null)
599 maxIterTest_->setMaxIters( maxIters_ );
603 if (params->isParameter(
"Num Blocks")) {
604 numBlocks_ = params->get(
"Num Blocks",numBlocks_default_);
605 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
606 "Belos::RCGSolMgr: \"Num Blocks\" must be strictly positive.");
609 params_->set(
"Num Blocks", numBlocks_);
613 if (params->isParameter(
"Num Recycled Blocks")) {
614 recycleBlocks_ = params->get(
"Num Recycled Blocks",recycleBlocks_default_);
615 TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks_ <= 0, std::invalid_argument,
616 "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be strictly positive.");
618 TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks_ >= numBlocks_, std::invalid_argument,
619 "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be less than \"Num Blocks\".");
622 params_->set(
"Num Recycled Blocks", recycleBlocks_);
626 if (params->isParameter(
"Timer Label")) {
627 std::string tempLabel = params->get(
"Timer Label", label_default_);
630 if (tempLabel != label_) {
632 params_->set(
"Timer Label", label_);
633 std::string solveLabel = label_ +
": RCGSolMgr total solve time";
634 #ifdef BELOS_TEUCHOS_TIME_MONITOR 635 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
641 if (params->isParameter(
"Verbosity")) {
642 if (Teuchos::isParameterType<int>(*params,
"Verbosity")) {
643 verbosity_ = params->get(
"Verbosity", verbosity_default_);
645 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,
"Verbosity");
649 params_->set(
"Verbosity", verbosity_);
650 if (printer_ != Teuchos::null)
651 printer_->setVerbosity(verbosity_);
655 if (params->isParameter(
"Output Style")) {
656 if (Teuchos::isParameterType<int>(*params,
"Output Style")) {
657 outputStyle_ = params->get(
"Output Style", outputStyle_default_);
659 outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,
"Output Style");
663 params_->set(
"Output Style", outputStyle_);
664 outputTest_ = Teuchos::null;
668 if (params->isParameter(
"Output Stream")) {
669 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,
"Output Stream");
672 params_->set(
"Output Stream", outputStream_);
673 if (printer_ != Teuchos::null)
674 printer_->setOStream( outputStream_ );
679 if (params->isParameter(
"Output Frequency")) {
680 outputFreq_ = params->get(
"Output Frequency", outputFreq_default_);
684 params_->set(
"Output Frequency", outputFreq_);
685 if (outputTest_ != Teuchos::null)
686 outputTest_->setOutputFrequency( outputFreq_ );
690 if (printer_ == Teuchos::null) {
699 if (params->isParameter(
"Convergence Tolerance")) {
700 convtol_ = params->get(
"Convergence Tolerance",convtol_default_);
703 params_->set(
"Convergence Tolerance", convtol_);
704 if (convTest_ != Teuchos::null)
708 if (params->isParameter(
"Show Maximum Residual Norm Only")) {
709 showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,
"Show Maximum Residual Norm Only");
712 params_->set(
"Show Maximum Residual Norm Only", showMaxResNormOnly_);
713 if (convTest_ != Teuchos::null)
714 convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
720 if (maxIterTest_ == Teuchos::null)
724 if (convTest_ == Teuchos::null)
725 convTest_ = Teuchos::rcp(
new StatusTestResNorm_t( convtol_, 1 ) );
727 if (sTest_ == Teuchos::null)
728 sTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
730 if (outputTest_ == Teuchos::null) {
738 std::string solverDesc =
" Recycling CG ";
739 outputTest_->setSolverDesc( solverDesc );
743 if (timerSolve_ == Teuchos::null) {
744 std::string solveLabel = label_ +
": RCGSolMgr total solve time";
745 #ifdef BELOS_TEUCHOS_TIME_MONITOR 746 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
755 template<
class ScalarType,
class MV,
class OP>
756 Teuchos::RCP<const Teuchos::ParameterList>
759 static Teuchos::RCP<const Teuchos::ParameterList> validPL;
762 if(is_null(validPL)) {
763 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
764 pl->set(
"Convergence Tolerance", convtol_default_,
765 "The relative residual tolerance that needs to be achieved by the\n" 766 "iterative solver in order for the linear system to be declared converged.");
767 pl->set(
"Maximum Iterations", maxIters_default_,
768 "The maximum number of block iterations allowed for each\n" 769 "set of RHS solved.");
770 pl->set(
"Block Size", blockSize_default_,
771 "Block Size Parameter -- currently must be 1 for RCG");
772 pl->set(
"Num Blocks", numBlocks_default_,
773 "The length of a cycle (and this max number of search vectors kept)\n");
774 pl->set(
"Num Recycled Blocks", recycleBlocks_default_,
775 "The number of vectors in the recycle subspace.");
776 pl->set(
"Verbosity", verbosity_default_,
777 "What type(s) of solver information should be outputted\n" 778 "to the output stream.");
779 pl->set(
"Output Style", outputStyle_default_,
780 "What style is used for the solver information outputted\n" 781 "to the output stream.");
782 pl->set(
"Output Frequency", outputFreq_default_,
783 "How often convergence information should be outputted\n" 784 "to the output stream.");
785 pl->set(
"Output Stream", outputStream_default_,
786 "A reference-counted pointer to the output stream where all\n" 787 "solver output is sent.");
788 pl->set(
"Show Maximum Residual Norm Only", showMaxResNormOnly_default_,
789 "When convergence information is printed, only show the maximum\n" 790 "relative residual norm when the block size is greater than one.");
791 pl->set(
"Timer Label", label_default_,
792 "The string to use as a prefix for the timer labels.");
800 template<
class ScalarType,
class MV,
class OP>
804 Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
805 if (rhsMV == Teuchos::null) {
812 TEUCHOS_TEST_FOR_EXCEPTION(static_cast<ptrdiff_t>(numBlocks_) > MVT::GetGlobalLength(*rhsMV),std::invalid_argument,
813 "Belos::RCGSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!");
816 if (P_ == Teuchos::null) {
817 P_ = MVT::Clone( *rhsMV, numBlocks_+2 );
821 if (MVT::GetNumberVecs(*P_) < numBlocks_+2) {
822 Teuchos::RCP<const MV> tmp = P_;
823 P_ = MVT::Clone( *tmp, numBlocks_+2 );
828 if (Ap_ == Teuchos::null)
829 Ap_ = MVT::Clone( *rhsMV, 1 );
832 if (r_ == Teuchos::null)
833 r_ = MVT::Clone( *rhsMV, 1 );
836 if (z_ == Teuchos::null)
837 z_ = MVT::Clone( *rhsMV, 1 );
840 if (U_ == Teuchos::null) {
841 U_ = MVT::Clone( *rhsMV, recycleBlocks_ );
845 if (MVT::GetNumberVecs(*U_) < recycleBlocks_) {
846 Teuchos::RCP<const MV> tmp = U_;
847 U_ = MVT::Clone( *tmp, recycleBlocks_ );
852 if (AU_ == Teuchos::null) {
853 AU_ = MVT::Clone( *rhsMV, recycleBlocks_ );
857 if (MVT::GetNumberVecs(*AU_) < recycleBlocks_) {
858 Teuchos::RCP<const MV> tmp = AU_;
859 AU_ = MVT::Clone( *tmp, recycleBlocks_ );
864 if (U1_ == Teuchos::null) {
865 U1_ = MVT::Clone( *rhsMV, recycleBlocks_ );
869 if (MVT::GetNumberVecs(*U1_) < recycleBlocks_) {
870 Teuchos::RCP<const MV> tmp = U1_;
871 U1_ = MVT::Clone( *tmp, recycleBlocks_ );
876 if (Alpha_ == Teuchos::null)
877 Alpha_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_, 1 ) );
879 if ( (Alpha_->numRows() != numBlocks_) || (Alpha_->numCols() != 1) )
880 Alpha_->reshape( numBlocks_, 1 );
884 if (Beta_ == Teuchos::null)
885 Beta_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + 1, 1 ) );
887 if ( (Beta_->numRows() != (numBlocks_+1)) || (Beta_->numCols() != 1) )
888 Beta_->reshape( numBlocks_ + 1, 1 );
892 if (D_ == Teuchos::null)
893 D_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ , 1 ) );
895 if ( (D_->numRows() != numBlocks_) || (D_->numCols() != 1) )
896 D_->reshape( numBlocks_, 1 );
900 if (Delta_ == Teuchos::null)
901 Delta_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ + 1 ) );
903 if ( (Delta_->numRows() != recycleBlocks_) || (Delta_->numCols() != (numBlocks_ + 1)) )
904 Delta_->reshape( recycleBlocks_, numBlocks_ + 1 );
908 if (UTAU_ == Teuchos::null)
909 UTAU_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
911 if ( (UTAU_->numRows() != recycleBlocks_) || (UTAU_->numCols() != recycleBlocks_) )
912 UTAU_->reshape( recycleBlocks_, recycleBlocks_ );
916 if (LUUTAU_ == Teuchos::null)
917 LUUTAU_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
919 if ( (LUUTAU_->numRows() != recycleBlocks_) || (LUUTAU_->numCols() != recycleBlocks_) )
920 LUUTAU_->reshape( recycleBlocks_, recycleBlocks_ );
924 if (ipiv_ == Teuchos::null)
925 ipiv_ = Teuchos::rcp(
new std::vector<int>(recycleBlocks_) );
927 if ( (
int)ipiv_->size() != recycleBlocks_ )
928 ipiv_->resize(recycleBlocks_);
932 if (AUTAU_ == Teuchos::null)
933 AUTAU_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
935 if ( (AUTAU_->numRows() != recycleBlocks_) || (AUTAU_->numCols() != recycleBlocks_) )
936 AUTAU_->reshape( recycleBlocks_, recycleBlocks_ );
940 if (rTz_old_ == Teuchos::null)
941 rTz_old_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( 1, 1 ) );
943 if ( (rTz_old_->numRows() != 1) || (rTz_old_->numCols() != 1) )
944 rTz_old_->reshape( 1, 1 );
948 if (F_ == Teuchos::null)
949 F_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ) );
951 if ( (F_->numRows() != (numBlocks_+recycleBlocks_)) || (F_->numCols() != numBlocks_+recycleBlocks_) )
952 F_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
956 if (G_ == Teuchos::null)
957 G_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ) );
959 if ( (G_->numRows() != (numBlocks_+recycleBlocks_)) || (G_->numCols() != numBlocks_+recycleBlocks_) )
960 G_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
964 if (Y_ == Teuchos::null)
965 Y_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, recycleBlocks_ ) );
967 if ( (Y_->numRows() != (numBlocks_+recycleBlocks_)) || (Y_->numCols() != numBlocks_+recycleBlocks_) )
968 Y_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
972 if (L2_ == Teuchos::null)
973 L2_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+1, numBlocks_ ) );
975 if ( (L2_->numRows() != (numBlocks_+1)) || (L2_->numCols() != numBlocks_) )
976 L2_->reshape( numBlocks_+1, numBlocks_ );
980 if (DeltaL2_ == Teuchos::null)
981 DeltaL2_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
983 if ( (DeltaL2_->numRows() != (recycleBlocks_)) || (DeltaL2_->numCols() != (numBlocks_) ) )
984 DeltaL2_->reshape( recycleBlocks_, numBlocks_ );
988 if (AU1TUDeltaL2_ == Teuchos::null)
989 AU1TUDeltaL2_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
991 if ( (AU1TUDeltaL2_->numRows() != (recycleBlocks_)) || (AU1TUDeltaL2_->numCols() != (numBlocks_) ) )
992 AU1TUDeltaL2_->reshape( recycleBlocks_, numBlocks_ );
996 if (AU1TAU1_ == Teuchos::null)
997 AU1TAU1_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
999 if ( (AU1TAU1_->numRows() != (recycleBlocks_)) || (AU1TAU1_->numCols() != (recycleBlocks_) ) )
1000 AU1TAU1_->reshape( recycleBlocks_, recycleBlocks_ );
1004 if (GY_ == Teuchos::null)
1005 GY_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + recycleBlocks_, recycleBlocks_ ) );
1007 if ( (GY_->numRows() != (numBlocks_ + recycleBlocks_)) || (GY_->numCols() != (recycleBlocks_) ) )
1008 GY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ );
1012 if (AU1TU1_ == Teuchos::null)
1013 AU1TU1_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
1015 if ( (AU1TU1_->numRows() != (recycleBlocks_)) || (AU1TU1_->numCols() != (recycleBlocks_) ) )
1016 AU1TU1_->reshape( recycleBlocks_, recycleBlocks_ );
1020 if (FY_ == Teuchos::null)
1021 FY_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + recycleBlocks_, recycleBlocks_ ) );
1023 if ( (FY_->numRows() != (numBlocks_ + recycleBlocks_)) || (FY_->numCols() != (recycleBlocks_) ) )
1024 FY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ );
1028 if (AU1TAP_ == Teuchos::null)
1029 AU1TAP_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
1031 if ( (AU1TAP_->numRows() != (recycleBlocks_)) || (AU1TAP_->numCols() != (numBlocks_) ) )
1032 AU1TAP_->reshape( recycleBlocks_, numBlocks_ );
1036 if (APTAP_ == Teuchos::null)
1037 APTAP_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_, numBlocks_ ) );
1039 if ( (APTAP_->numRows() != (numBlocks_)) || (APTAP_->numCols() != (numBlocks_) ) )
1040 APTAP_->reshape( numBlocks_, numBlocks_ );
1044 if (U1Y1_ == Teuchos::null) {
1045 U1Y1_ = MVT::Clone( *rhsMV, recycleBlocks_ );
1049 if (MVT::GetNumberVecs(*U1Y1_) < recycleBlocks_) {
1050 Teuchos::RCP<const MV> tmp = U1Y1_;
1051 U1Y1_ = MVT::Clone( *tmp, recycleBlocks_ );
1056 if (PY2_ == Teuchos::null) {
1057 PY2_ = MVT::Clone( *rhsMV, recycleBlocks_ );
1061 if (MVT::GetNumberVecs(*PY2_) < recycleBlocks_) {
1062 Teuchos::RCP<const MV> tmp = PY2_;
1063 PY2_ = MVT::Clone( *tmp, recycleBlocks_ );
1068 if (AUTAP_ == Teuchos::null)
1069 AUTAP_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
1071 if ( (AUTAP_->numRows() != (recycleBlocks_)) || (AUTAP_->numCols() != (numBlocks_) ) )
1072 AUTAP_->reshape( recycleBlocks_, numBlocks_ );
1076 if (AU1TU_ == Teuchos::null)
1077 AU1TU_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
1079 if ( (AU1TU_->numRows() != (recycleBlocks_)) || (AU1TU_->numCols() != (recycleBlocks_) ) )
1080 AU1TU_->reshape( recycleBlocks_, recycleBlocks_ );
1087 template<
class ScalarType,
class MV,
class OP>
1090 Teuchos::LAPACK<int,ScalarType> lapack;
1091 std::vector<int> index(recycleBlocks_);
1092 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1093 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1102 setParameters(Teuchos::parameterList(*getValidParameters()));
1106 "Belos::RCGSolMgr::solve(): Linear problem is not a valid object.");
1108 "Belos::RCGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
1109 TEUCHOS_TEST_FOR_EXCEPTION((problem_->getLeftPrec() != Teuchos::null)&&(problem_->getRightPrec() != Teuchos::null),
1111 "Belos::RCGSolMgr::solve(): RCG does not support split preconditioning, only set left or right preconditioner.");
1114 Teuchos::RCP<OP> precObj;
1115 if (problem_->getLeftPrec() != Teuchos::null) {
1116 precObj = Teuchos::rcp_const_cast<OP>(problem_->getLeftPrec());
1118 else if (problem_->getRightPrec() != Teuchos::null) {
1119 precObj = Teuchos::rcp_const_cast<OP>(problem_->getRightPrec());
1123 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1124 std::vector<int> currIdx(1);
1128 problem_->setLSIndex( currIdx );
1131 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1132 if (numBlocks_ > dim) {
1133 numBlocks_ = Teuchos::asSafe<int>(dim);
1134 params_->set(
"Num Blocks", numBlocks_);
1136 "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl <<
1137 " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl;
1141 initializeStateStorage();
1144 Teuchos::ParameterList plist;
1145 plist.set(
"Num Blocks",numBlocks_);
1146 plist.set(
"Recycled Blocks",recycleBlocks_);
1149 outputTest_->reset();
1152 bool isConverged =
true;
1156 index.resize(recycleBlocks_);
1157 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1158 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1159 index.resize(recycleBlocks_);
1160 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1161 Teuchos::RCP<MV> AUtmp = MVT::CloneViewNonConst( *AU_, index );
1163 problem_->applyOp( *Utmp, *AUtmp );
1165 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1166 MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp );
1168 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
1169 if ( precObj != Teuchos::null ) {
1170 index.resize(recycleBlocks_);
1171 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1172 index.resize(recycleBlocks_);
1173 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1174 Teuchos::RCP<MV> PCAU = MVT::CloneViewNonConst( *U1_, index );
1175 OPT::Apply( *precObj, *AUtmp, *PCAU );
1176 MVT::MvTransMv( one, *AUtmp, *PCAU, AUTAUtmp );
1178 MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp );
1185 Teuchos::RCP<RCGIter<ScalarType,MV,OP> > rcg_iter;
1190 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1191 Teuchos::TimeMonitor slvtimer(*timerSolve_);
1194 while ( numRHS2Solve > 0 ) {
1197 if (printer_->isVerbosity(
Debug ) ) {
1198 if (existU_) printer_->print(
Debug,
"Using recycle space generated from previous call to solve()." );
1199 else printer_->print(
Debug,
"No recycle space exists." );
1203 rcg_iter->resetNumIters();
1206 rcg_iter->setSize( recycleBlocks_, numBlocks_ );
1209 outputTest_->resetNumCalls();
1218 problem_->computeCurrResVec( &*r_ );
1223 Teuchos::SerialDenseMatrix<int,ScalarType> Utr(recycleBlocks_,1);
1224 index.resize(recycleBlocks_);
1225 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1226 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1227 MVT::MvTransMv( one, *Utmp, *r_, Utr );
1228 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1229 Teuchos::SerialDenseMatrix<int,ScalarType> LUUTAUtmp( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ );
1230 LUUTAUtmp.assign(UTAUtmp);
1232 lapack.GESV(recycleBlocks_, 1, LUUTAUtmp.values(), LUUTAUtmp.stride(), &(*ipiv_)[0], Utr.values(), Utr.stride(), &info);
1234 "Belos::RCGSolMgr::solve(): LAPACK GESV failed to compute a solution.");
1237 MVT::MvTimesMatAddMv( one, *Utmp, Utr, one, *problem_->getCurrLHSVec() );
1240 index.resize(recycleBlocks_);
1241 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1242 Teuchos::RCP<const MV> AUtmp = MVT::CloneView( *AU_, index );
1243 MVT::MvTimesMatAddMv( -one, *AUtmp, Utr, one, *r_ );
1246 if ( precObj != Teuchos::null ) {
1247 OPT::Apply( *precObj, *r_, *z_ );
1253 MVT::MvTransMv( one, *r_, *z_, *rTz_old_ );
1257 Teuchos::SerialDenseMatrix<int,ScalarType> mu( Teuchos::View, *Delta_, recycleBlocks_, 1);
1258 index.resize(recycleBlocks_);
1259 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1260 Teuchos::RCP<const MV> AUtmp = MVT::CloneView( *AU_, index );
1261 MVT::MvTransMv( one, *AUtmp, *z_, mu );
1262 Teuchos::SerialDenseMatrix<int,ScalarType> LUUTAUtmp( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ );
1265 lapack.GETRS(
TRANS, recycleBlocks_, 1, LUUTAUtmp.values(), LUUTAUtmp.stride(), &(*ipiv_)[0], mu.values(), mu.stride(), &info );
1267 "Belos::RCGSolMgr::solve(): LAPACK GETRS failed to compute a solution.");
1271 Teuchos::RCP<MV> Ptmp = MVT::CloneViewNonConst( *P_, index );
1272 MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp);
1273 MVT::MvTimesMatAddMv( -one, *U_, mu, one, *Ptmp );
1278 Teuchos::RCP<MV> Ptmp = MVT::CloneViewNonConst( *P_, index );
1279 MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp);
1286 index.resize( numBlocks_+1 );
1287 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1288 newstate.
P = MVT::CloneViewNonConst( *P_, index );
1289 index.resize( recycleBlocks_ );
1290 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1291 newstate.
U = MVT::CloneViewNonConst( *U_, index );
1292 index.resize( recycleBlocks_ );
1293 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1294 newstate.
AU = MVT::CloneViewNonConst( *AU_, index );
1295 newstate.
Alpha = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Alpha_, numBlocks_, 1 ) );
1296 newstate.
Beta = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Beta_, numBlocks_, 1 ) );
1297 newstate.
D = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *D_, numBlocks_, 1 ) );
1298 newstate.
Delta = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_, 0, 1 ) );
1299 newstate.
LUUTAU = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ ) );
1305 newstate.
existU = existU_;
1306 newstate.
ipiv = ipiv_;
1309 rcg_iter->initialize(newstate);
1315 rcg_iter->iterate();
1322 if ( convTest_->getStatus() ==
Passed ) {
1331 else if ( maxIterTest_->getStatus() ==
Passed ) {
1333 isConverged =
false;
1341 else if ( rcg_iter->getCurSubspaceDim() == rcg_iter->getMaxSubspaceDim() ) {
1346 if (recycleBlocks_ > 0) {
1350 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, numBlocks_, numBlocks_ );
1351 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, numBlocks_, numBlocks_ );
1352 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1353 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1354 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_, 1 );
1355 Ftmp.putScalar(zero);
1356 Gtmp.putScalar(zero);
1357 for (
int ii=0;ii<numBlocks_;ii++) {
1358 Gtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0));
1360 Gtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1361 Gtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1363 Ftmp(ii,ii) = Dtmp(ii,0);
1367 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, numBlocks_, recycleBlocks_ );
1368 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1371 index.resize( numBlocks_ );
1372 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; }
1373 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1374 index.resize( recycleBlocks_ );
1375 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1376 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1377 MVT::MvTimesMatAddMv( one, *Ptmp, Ytmp, zero, *U1tmp );
1382 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, numBlocks_, recycleBlocks_ );
1383 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1384 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1385 AU1TAU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1388 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, numBlocks_, recycleBlocks_ );
1389 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1390 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1391 AU1TU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1393 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, 1 );
1395 AU1TAPtmp.putScalar(zero);
1397 ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0);
1398 for (
int ii=0; ii<recycleBlocks_; ++ii) {
1399 AU1TAPtmp(ii,0) = Ytmp(numBlocks_-1,ii) * alphatmp;
1407 lastBeta = numBlocks_-1;
1414 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1415 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, numBlocks_ );
1416 AU1TAPtmp.scale(Dtmp(0,0));
1418 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1419 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_+1, 1 );
1420 Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
1421 APTAPtmp.putScalar(zero);
1422 for (
int ii=0; ii<numBlocks_; ii++) {
1423 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0));
1425 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1426 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1431 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1432 Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
1433 Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1434 Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1435 Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1436 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1437 F11.assign(AU1TU1tmp);
1438 F12.putScalar(zero);
1439 F21.putScalar(zero);
1440 F22.putScalar(zero);
1441 for(
int ii=0;ii<numBlocks_;ii++) {
1442 F22(ii,ii) = Dtmp(ii,0);
1446 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1447 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1448 Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
1449 Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1450 Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1451 Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1452 G11.assign(AU1TAU1tmp);
1453 G12.assign(AU1TAPtmp);
1455 for (
int ii=0;ii<recycleBlocks_;++ii)
1456 for (
int jj=0;jj<numBlocks_;++jj)
1457 G21(jj,ii) = G12(ii,jj);
1458 G22.assign(APTAPtmp);
1461 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
1462 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1465 index.resize( numBlocks_ );
1466 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; }
1467 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1468 index.resize( recycleBlocks_ );
1469 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1470 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index );
1471 Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1472 index.resize( recycleBlocks_ );
1473 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1474 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1475 index.resize( recycleBlocks_ );
1476 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1477 Teuchos::RCP<MV> U1Y1tmp = MVT::CloneViewNonConst( *U1Y1_, index );
1478 Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
1479 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1480 MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp );
1481 MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp);
1486 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1487 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1488 AU1TAU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1492 AU1TAPtmp.putScalar(zero);
1493 ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0);
1494 for (
int ii=0; ii<recycleBlocks_; ++ii) {
1495 AU1TAPtmp(ii,0) = Ytmp(numBlocks_+recycleBlocks_-1,ii) * alphatmp;
1499 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1500 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1502 AU1TU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1505 lastp = numBlocks_+1;
1506 lastBeta = numBlocks_;
1512 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1513 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_, 1 );
1514 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1515 Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
1516 APTAPtmp.putScalar(zero);
1517 for (
int ii=0; ii<numBlocks_; ii++) {
1518 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0));
1520 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1521 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1525 Teuchos::SerialDenseMatrix<int,ScalarType> L2tmp( Teuchos::View, *L2_, numBlocks_+1, numBlocks_ );
1526 L2tmp.putScalar(zero);
1527 for(
int ii=0;ii<numBlocks_;ii++) {
1528 L2tmp(ii,ii) = 1./Alphatmp(ii,0);
1529 L2tmp(ii+1,ii) = -1./Alphatmp(ii,0);
1533 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAPtmp( Teuchos::View, *AUTAP_, recycleBlocks_, numBlocks_ );
1534 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1535 Teuchos::SerialDenseMatrix<int,ScalarType> Deltatmp( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_+1 );
1536 Teuchos::SerialDenseMatrix<int,ScalarType> DeltaL2tmp( Teuchos::View, *DeltaL2_, recycleBlocks_, numBlocks_ );
1537 DeltaL2tmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Deltatmp,L2tmp,zero);
1538 AUTAPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,UTAUtmp,DeltaL2tmp,zero);
1541 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1542 Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
1543 Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1544 Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1545 Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1546 F11.assign(UTAUtmp);
1547 F12.putScalar(zero);
1548 F21.putScalar(zero);
1549 F22.putScalar(zero);
1550 for(
int ii=0;ii<numBlocks_;ii++) {
1551 F22(ii,ii) = Dtmp(ii,0);
1555 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1556 Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
1557 Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1558 Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1559 Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1560 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
1561 G11.assign(AUTAUtmp);
1562 G12.assign(AUTAPtmp);
1564 for (
int ii=0;ii<recycleBlocks_;++ii)
1565 for (
int jj=0;jj<numBlocks_;++jj)
1566 G21(jj,ii) = G12(ii,jj);
1567 G22.assign(APTAPtmp);
1570 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
1571 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1574 index.resize( recycleBlocks_ );
1575 for (
int ii=0; ii<(recycleBlocks_); ++ii) { index[ii] = ii; }
1576 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1577 index.resize( numBlocks_ );
1578 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; }
1579 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1580 index.resize( recycleBlocks_ );
1581 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1582 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index );
1583 Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1584 index.resize( recycleBlocks_ );
1585 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1586 Teuchos::RCP<MV> UY1tmp = MVT::CloneViewNonConst( *U1Y1_, index );
1587 Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
1588 index.resize( recycleBlocks_ );
1589 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1590 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1591 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1592 MVT::MvTimesMatAddMv( one, *Utmp, Y1, zero, *UY1tmp );
1593 MVT::MvAddMv(one,*UY1tmp, one, *PY2tmp, *U1tmp);
1598 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1599 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1600 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1601 AU1TAU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1604 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1605 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1606 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1607 AU1TU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1610 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUtmp( Teuchos::View, *AU1TU_, recycleBlocks_, recycleBlocks_ );
1611 AU1TUtmp.assign(UTAUtmp);
1614 dold = Dtmp(numBlocks_-1,0);
1621 lastBeta = numBlocks_-1;
1624 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1625 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_+1, 1 );
1626 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1627 Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
1628 for (
int ii=0; ii<numBlocks_; ii++) {
1629 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0));
1631 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1632 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1636 Teuchos::SerialDenseMatrix<int,ScalarType> L2tmp( Teuchos::View, *L2_, numBlocks_+1, numBlocks_ );
1637 for(
int ii=0;ii<numBlocks_;ii++) {
1638 L2tmp(ii,ii) = 1./Alphatmp(ii,0);
1639 L2tmp(ii+1,ii) = -1./Alphatmp(ii,0);
1644 Teuchos::SerialDenseMatrix<int,ScalarType> DeltaL2( Teuchos::View, *DeltaL2_, recycleBlocks_, numBlocks_ );
1645 Teuchos::SerialDenseMatrix<int,ScalarType> Deltatmp( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_+1 );
1646 DeltaL2.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Deltatmp,L2tmp,zero);
1647 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUDeltaL2( Teuchos::View, *AU1TUDeltaL2_, recycleBlocks_, numBlocks_ );
1648 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUtmp( Teuchos::View, *AU1TU_, recycleBlocks_, recycleBlocks_ );
1649 AU1TUDeltaL2.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,AU1TUtmp,DeltaL2,zero);
1650 Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
1651 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, numBlocks_ );
1652 AU1TAPtmp.putScalar(zero);
1653 AU1TAPtmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Y1,AU1TUDeltaL2,zero);
1654 Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1655 ScalarType val = dold * (-Betatmp(0,0)/Alphatmp(0,0));
1656 for(
int ii=0;ii<recycleBlocks_;ii++) {
1657 AU1TAPtmp(ii,0) += Y2(numBlocks_-1,ii)*val;
1661 Teuchos::SerialDenseMatrix<int,ScalarType> Y1TAU1TU( Teuchos::View, *GY_, recycleBlocks_, recycleBlocks_ );
1662 Y1TAU1TU.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Y1,AU1TUtmp,zero);
1663 AU1TUtmp.assign(Y1TAU1TU);
1666 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1667 Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
1668 Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1669 Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1670 Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1671 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1672 F11.assign(AU1TU1tmp);
1673 for(
int ii=0;ii<numBlocks_;ii++) {
1674 F22(ii,ii) = Dtmp(ii,0);
1678 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1679 Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
1680 Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1681 Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1682 Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1683 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1684 G11.assign(AU1TAU1tmp);
1685 G12.assign(AU1TAPtmp);
1687 for (
int ii=0;ii<recycleBlocks_;++ii)
1688 for (
int jj=0;jj<numBlocks_;++jj)
1689 G21(jj,ii) = G12(ii,jj);
1690 G22.assign(APTAPtmp);
1693 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
1694 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1697 index.resize( numBlocks_ );
1698 for (
int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; }
1699 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1700 index.resize( recycleBlocks_ );
1701 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1702 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index );
1703 index.resize( recycleBlocks_ );
1704 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1705 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1706 index.resize( recycleBlocks_ );
1707 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1708 Teuchos::RCP<MV> U1Y1tmp = MVT::CloneViewNonConst( *U1Y1_, index );
1709 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1710 MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp );
1711 MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp);
1716 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1717 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1718 AU1TAU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1721 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1722 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1723 AU1TU1tmp.multiply(
Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1726 dold = Dtmp(numBlocks_-1,0);
1729 lastp = numBlocks_+1;
1730 lastBeta = numBlocks_;
1740 index[0] = lastp-1; index[1] = lastp;
1741 Teuchos::RCP<const MV> Ptmp2 = MVT::CloneView( *P_, index );
1742 index[0] = 0; index[1] = 1;
1743 MVT::SetBlock(*Ptmp2,index,*P_);
1746 (*Beta_)(0,0) = (*Beta_)(lastBeta,0);
1750 Teuchos::SerialDenseMatrix<int,ScalarType> mu1( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, 0 );
1751 Teuchos::SerialDenseMatrix<int,ScalarType> mu2( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, numBlocks_ );
1756 newstate.
P = Teuchos::null;
1757 index.resize( numBlocks_+1 );
1758 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii+1; }
1759 newstate.
P = MVT::CloneViewNonConst( *P_, index );
1761 newstate.
Beta = Teuchos::null;
1762 newstate.
Beta = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Beta_, numBlocks_, 1, 1, 0 ) );
1764 newstate.
Delta = Teuchos::null;
1765 newstate.
Delta = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_, 0, 1 ) );
1770 rcg_iter->initialize(newstate);
1783 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
1784 "Belos::RCGSolMgr::solve(): Invalid return from RCGIter::iterate().");
1787 catch (
const std::exception &e) {
1788 printer_->stream(
Errors) <<
"Error! Caught std::exception in RCGIter::iterate() at iteration " 1789 << rcg_iter->getNumIters() << std::endl
1790 << e.what() << std::endl;
1796 problem_->setCurrLS();
1800 if ( numRHS2Solve > 0 ) {
1803 problem_->setLSIndex( currIdx );
1806 currIdx.resize( numRHS2Solve );
1812 index.resize(recycleBlocks_);
1813 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1814 MVT::SetBlock(*U1_,index,*U_);
1817 if (numRHS2Solve > 0) {
1819 newstate.
P = Teuchos::null;
1820 newstate.
Ap = Teuchos::null;
1821 newstate.
r = Teuchos::null;
1822 newstate.
z = Teuchos::null;
1823 newstate.
U = Teuchos::null;
1824 newstate.
AU = Teuchos::null;
1825 newstate.
Alpha = Teuchos::null;
1826 newstate.
Beta = Teuchos::null;
1827 newstate.
D = Teuchos::null;
1828 newstate.
Delta = Teuchos::null;
1829 newstate.
LUUTAU = Teuchos::null;
1830 newstate.
ipiv = Teuchos::null;
1831 newstate.
rTz_old = Teuchos::null;
1834 index.resize(recycleBlocks_);
1835 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1836 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1837 index.resize(recycleBlocks_);
1838 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1839 Teuchos::RCP<MV> AUtmp = MVT::CloneViewNonConst( *AU_, index );
1841 problem_->applyOp( *Utmp, *AUtmp );
1843 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1844 MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp );
1846 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
1847 if ( precObj != Teuchos::null ) {
1848 index.resize(recycleBlocks_);
1849 for (
int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1850 index.resize(recycleBlocks_);
1851 for (
int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1852 Teuchos::RCP<MV> LeftPCAU = MVT::CloneViewNonConst( *U1_, index );
1853 OPT::Apply( *precObj, *AUtmp, *LeftPCAU );
1854 MVT::MvTransMv( one, *AUtmp, *LeftPCAU, AUTAUtmp );
1856 MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp );
1869 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1874 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails) );
1878 numIters_ = maxIterTest_->getNumIters();
1882 using Teuchos::rcp_dynamic_cast;
1885 const std::vector<MagnitudeType>* pTestValues =
1886 rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
1888 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1889 "Belos::RCGSolMgr::solve(): The convergence test's getTestValue() " 1890 "method returned NULL. Please report this bug to the Belos developers.");
1892 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1893 "Belos::RCGSolMgr::solve(): The convergence test's getTestValue() " 1894 "method returned a vector of length zero. Please report this bug to the " 1895 "Belos developers.");
1900 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1910 template<
class ScalarType,
class MV,
class OP>
1912 const Teuchos::SerialDenseMatrix<int,ScalarType>& G,
1913 Teuchos::SerialDenseMatrix<int,ScalarType>& Y) {
1915 int n = F.numCols();
1918 Teuchos::LAPACK<int,ScalarType> lapack;
1921 std::vector<MagnitudeType> w(n);
1924 std::vector<int> iperm(n);
1931 std::vector<ScalarType> work(1);
1935 Teuchos::SerialDenseMatrix<int,ScalarType> F2( Teuchos::Copy, *F_ );
1936 Teuchos::SerialDenseMatrix<int,ScalarType> G2( Teuchos::Copy, *G_ );
1939 lapack.SYGV(itype, jobz, uplo, n, G2.values(), G2.stride(), F2.values(), F2.stride(), &w[0], &work[0], lwork, &info);
1941 "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to query optimal work size.");
1942 lwork = (int)work[0];
1944 lapack.SYGV(itype, jobz, uplo, n, G2.values(), G2.stride(), F2.values(), F2.stride(), &w[0], &work[0], lwork, &info);
1946 "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to compute eigensolutions.");
1950 this->sort(w,n,iperm);
1953 for(
int i=0; i<recycleBlocks_; i++ ) {
1954 for(
int j=0; j<n; j++ ) {
1955 Y(j,i) = G2(j,iperm[i]);
1962 template<
class ScalarType,
class MV,
class OP>
1963 void RCGSolMgr<ScalarType,MV,OP,true>::sort(std::vector<ScalarType>& dlist,
int n, std::vector<int>& iperm)
1965 int l, r, j, i, flag;
1994 if (dlist[j] > dlist[j - 1]) j = j + 1;
1996 if (dlist[j - 1] > dK) {
1997 dlist[i - 1] = dlist[j - 1];
1998 iperm[i - 1] = iperm[j - 1];
2011 dlist[r] = dlist[0];
2012 iperm[r] = iperm[0];
2027 template<
class ScalarType,
class MV,
class OP>
2030 std::ostringstream oss;
2031 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
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.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem)
Set the linear problem that needs to be solved.
RCGSolMgrRecyclingFailure is thrown when any problem occurs in using/creating the recycling subspace...
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.
int getNumIters() const
Get the iteration count for the most recent call to solve().
Belos::StatusTest class for specifying a maximum number of iterations.
A factory class for generating StatusTestOutput objects.
Iterated Modified Gram-Schmidt (IMGS) implementation of the Belos::OrthoManager class.
RCGSolMgrLinearProblemFailure(const std::string &what_arg)
RCGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i.e.
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
Structure to contain pointers to RCGIter state variables.
Belos concrete class for performing the RCG iteration.
Classical Gram-Schmidt (with DGKS correction) implementation of the Belos::OrthoManager class...
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
A Belos::StatusTest class for specifying a maximum number of iterations.
ResetType
How to reset the solver.
bool existU
Flag to indicate the recycle space should be used.
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.
MagnitudeType achievedTol() const
Tolerance achieved by the last solve() invocation.
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.
bool isLOADetected() const
Return whether a loss of accuracy was detected by this solver during the most current solve...
Iterated Classical Gram-Schmidt (ICGS) implementation of the Belos::OrthoManager class.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Alpha
Coefficients arising in RCG iteration.
int setTolerance(MagnitudeType tolerance)
Set the value of the tolerance.
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > create(const Teuchos::RCP< OutputManager< ScalarType > > &printer, Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test, int mod, int printStates)
Create the StatusTestOutput object specified by the outputStyle.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
int curDim
The current dimension of the reduction.
RCGSolMgrLAPACKFailure(const std::string &what_arg)
void reset(const ResetType type)
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
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.
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const
Get a parameter list containing the current parameters for this object.
Class which defines basic traits for the operator type.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > D
Parent class to all Belos exceptions.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Return a reference to the linear problem being solved by this solver manager.
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.