42#ifndef BELOS_PCPG_SOLMGR_HPP
43#define BELOS_PCPG_SOLMGR_HPP
62#include "Teuchos_LAPACK.hpp"
63#ifdef BELOS_TEUCHOS_TIME_MONITOR
64# include "Teuchos_TimeMonitor.hpp"
66#if defined(HAVE_TEUCHOSCORE_CXX11)
67# include <type_traits>
69#include "Teuchos_TypeTraits.hpp"
138 template<
class ScalarType,
class MV,
class OP,
139 const bool supportsScalarType =
141 ! Teuchos::ScalarTraits<ScalarType>::isComplex>
144 Belos::Details::LapackSupportsScalar<ScalarType>::value &&
145 ! Teuchos::ScalarTraits<ScalarType>::isComplex>
147 static const bool scalarTypeIsSupported =
149 ! Teuchos::ScalarTraits<ScalarType>::isComplex;
158 const Teuchos::RCP<Teuchos::ParameterList> &pl) :
164 Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const override {
169 template<
class ScalarType,
class MV,
class OP>
175 typedef Teuchos::ScalarTraits<ScalarType> SCT;
176 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
177 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
227 const Teuchos::RCP<Teuchos::ParameterList> &pl );
233 virtual Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const {
249 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters()
const;
260 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
261 return Teuchos::tuple(timerSolve_);
291 void setParameters(
const Teuchos::RCP<Teuchos::ParameterList> ¶ms );
332 std::string description()
const;
340 int ARRQR(
int numVecs,
int numOrthVecs,
const Teuchos::SerialDenseMatrix<int,ScalarType>& D);
343 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
346 Teuchos::RCP<OutputManager<ScalarType> > printer_;
347 Teuchos::RCP<std::ostream> outputStream_;
350 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
351 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
352 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
353 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
356 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
359 Teuchos::RCP<Teuchos::ParameterList> params_;
362 static constexpr int maxIters_default_ = 1000;
363 static constexpr int deflatedBlocks_default_ = 2;
364 static constexpr int savedBlocks_default_ = 16;
367 static constexpr int outputFreq_default_ = -1;
368 static constexpr const char * label_default_ =
"Belos";
369 static constexpr const char * orthoType_default_ =
"ICGS";
376 MagnitudeType convtol_;
379 MagnitudeType orthoKappa_;
382 MagnitudeType achievedTol_;
390 int deflatedBlocks_, savedBlocks_, verbosity_, outputStyle_, outputFreq_;
391 std::string orthoType_;
394 Teuchos::RCP<MV> U_, C_, R_;
401 Teuchos::RCP<Teuchos::Time> timerSolve_;
409template<
class ScalarType,
class MV,
class OP>
411 outputStream_(Teuchos::rcpFromRef(std::cout)),
414 achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
416 maxIters_(maxIters_default_),
417 deflatedBlocks_(deflatedBlocks_default_),
418 savedBlocks_(savedBlocks_default_),
419 verbosity_(verbosity_default_),
420 outputStyle_(outputStyle_default_),
421 outputFreq_(outputFreq_default_),
422 orthoType_(orthoType_default_),
424 label_(label_default_),
430template<
class ScalarType,
class MV,
class OP>
433 const Teuchos::RCP<Teuchos::ParameterList> &pl ) :
435 outputStream_(Teuchos::rcpFromRef(std::cout)),
439 achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
441 maxIters_(maxIters_default_),
442 deflatedBlocks_(deflatedBlocks_default_),
443 savedBlocks_(savedBlocks_default_),
444 verbosity_(verbosity_default_),
445 outputStyle_(outputStyle_default_),
446 outputFreq_(outputFreq_default_),
447 orthoType_(orthoType_default_),
449 label_(label_default_),
452 TEUCHOS_TEST_FOR_EXCEPTION(
453 problem_.is_null (), std::invalid_argument,
454 "Belos::PCPGSolMgr two-argument constructor: "
455 "'problem' is null. You must supply a non-null Belos::LinearProblem "
456 "instance when calling this constructor.");
458 if (! pl.is_null ()) {
465template<
class ScalarType,
class MV,
class OP>
469 if (params_ == Teuchos::null) {
470 params_ = Teuchos::rcp(
new Teuchos::ParameterList(*getValidParameters()) );
473 params->validateParameters(*getValidParameters());
477 if (params->isParameter(
"Maximum Iterations")) {
478 maxIters_ = params->get(
"Maximum Iterations",maxIters_default_);
481 params_->set(
"Maximum Iterations", maxIters_);
482 if (maxIterTest_!=Teuchos::null)
483 maxIterTest_->setMaxIters( maxIters_ );
487 if (params->isParameter(
"Num Saved Blocks")) {
488 savedBlocks_ = params->get(
"Num Saved Blocks",savedBlocks_default_);
489 TEUCHOS_TEST_FOR_EXCEPTION(savedBlocks_ <= 0, std::invalid_argument,
490 "Belos::PCPGSolMgr: \"Num Saved Blocks\" must be strictly positive.");
497 params_->set(
"Num Saved Blocks", savedBlocks_);
499 if (params->isParameter(
"Num Deflated Blocks")) {
500 deflatedBlocks_ = params->get(
"Num Deflated Blocks",deflatedBlocks_default_);
501 TEUCHOS_TEST_FOR_EXCEPTION(deflatedBlocks_ < 0, std::invalid_argument,
502 "Belos::PCPGSolMgr: \"Num Deflated Blocks\" must be positive.");
504 TEUCHOS_TEST_FOR_EXCEPTION(deflatedBlocks_ > savedBlocks_, std::invalid_argument,
505 "Belos::PCPGSolMgr: \"Num Deflated Blocks\" must be <= \"Num Saved Blocks\".");
509 params_->set(
"Num Deflated Blocks",
static_cast<int>(deflatedBlocks_));
513 if (params->isParameter(
"Timer Label")) {
514 std::string tempLabel = params->get(
"Timer Label", label_default_);
517 if (tempLabel != label_) {
519 params_->set(
"Timer Label", label_);
520 std::string solveLabel = label_ +
": PCPGSolMgr total solve time";
521#ifdef BELOS_TEUCHOS_TIME_MONITOR
522 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
524 if (ortho_ != Teuchos::null) {
525 ortho_->setLabel( label_ );
531 if (params->isParameter(
"Verbosity")) {
532 if (Teuchos::isParameterType<int>(*params,
"Verbosity")) {
533 verbosity_ = params->get(
"Verbosity", verbosity_default_);
535 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,
"Verbosity");
539 params_->set(
"Verbosity", verbosity_);
540 if (printer_ != Teuchos::null)
541 printer_->setVerbosity(verbosity_);
545 if (params->isParameter(
"Output Style")) {
546 if (Teuchos::isParameterType<int>(*params,
"Output Style")) {
547 outputStyle_ = params->get(
"Output Style", outputStyle_default_);
549 outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,
"Output Style");
553 params_->set(
"Output Style", outputStyle_);
554 outputTest_ = Teuchos::null;
558 if (params->isParameter(
"Output Stream")) {
559 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,
"Output Stream");
562 params_->set(
"Output Stream", outputStream_);
563 if (printer_ != Teuchos::null)
564 printer_->setOStream( outputStream_ );
569 if (params->isParameter(
"Output Frequency")) {
570 outputFreq_ = params->get(
"Output Frequency", outputFreq_default_);
574 params_->set(
"Output Frequency", outputFreq_);
575 if (outputTest_ != Teuchos::null)
576 outputTest_->setOutputFrequency( outputFreq_ );
580 if (printer_ == Teuchos::null) {
585 bool changedOrthoType =
false;
586 if (params->isParameter(
"Orthogonalization")) {
587 std::string tempOrthoType = params->get(
"Orthogonalization",orthoType_default_);
588 if (tempOrthoType != orthoType_) {
589 orthoType_ = tempOrthoType;
590 changedOrthoType =
true;
593 params_->set(
"Orthogonalization", orthoType_);
596 if (params->isParameter(
"Orthogonalization Constant")) {
597 if (params->isType<MagnitudeType> (
"Orthogonalization Constant")) {
598 orthoKappa_ = params->get (
"Orthogonalization Constant",
602 orthoKappa_ = params->get (
"Orthogonalization Constant",
607 params_->set(
"Orthogonalization Constant",orthoKappa_);
608 if (orthoType_==
"DGKS") {
609 if (orthoKappa_ > 0 && ortho_ != Teuchos::null && !changedOrthoType) {
610 Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
616 if (ortho_ == Teuchos::null || changedOrthoType) {
618 Teuchos::RCP<Teuchos::ParameterList> paramsOrtho;
619 if (orthoType_==
"DGKS" && orthoKappa_ > 0) {
620 paramsOrtho->set (
"depTol", orthoKappa_ );
623 ortho_ = factory.
makeMatOrthoManager (orthoType_, Teuchos::null, printer_, label_, paramsOrtho);
631 if (params->isParameter(
"Convergence Tolerance")) {
632 if (params->isType<MagnitudeType> (
"Convergence Tolerance")) {
633 convtol_ = params->get (
"Convergence Tolerance",
641 params_->set(
"Convergence Tolerance", convtol_);
642 if (convTest_ != Teuchos::null)
649 if (maxIterTest_ == Teuchos::null)
652 if (convTest_ == Teuchos::null)
653 convTest_ = Teuchos::rcp(
new StatusTestResNorm_t( convtol_, 1 ) );
655 sTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
663 std::string solverDesc =
" PCPG ";
664 outputTest_->setSolverDesc( solverDesc );
667 if (timerSolve_ == Teuchos::null) {
668 std::string solveLabel = label_ +
": PCPGSolMgr total solve time";
669#ifdef BELOS_TEUCHOS_TIME_MONITOR
670 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
679template<
class ScalarType,
class MV,
class OP>
680Teuchos::RCP<const Teuchos::ParameterList>
683 static Teuchos::RCP<const Teuchos::ParameterList> validPL;
684 if (is_null(validPL)) {
685 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
688 "The relative residual tolerance that needs to be achieved by the\n"
689 "iterative solver in order for the linear system to be declared converged.");
690 pl->set(
"Maximum Iterations",
static_cast<int>(maxIters_default_),
691 "The maximum number of iterations allowed for each\n"
692 "set of RHS solved.");
693 pl->set(
"Num Deflated Blocks",
static_cast<int>(deflatedBlocks_default_),
694 "The maximum number of vectors in the seed subspace." );
695 pl->set(
"Num Saved Blocks",
static_cast<int>(savedBlocks_default_),
696 "The maximum number of vectors saved from old Krylov subspaces." );
697 pl->set(
"Verbosity",
static_cast<int>(verbosity_default_),
698 "What type(s) of solver information should be outputted\n"
699 "to the output stream.");
700 pl->set(
"Output Style",
static_cast<int>(outputStyle_default_),
701 "What style is used for the solver information outputted\n"
702 "to the output stream.");
703 pl->set(
"Output Frequency",
static_cast<int>(outputFreq_default_),
704 "How often convergence information should be outputted\n"
705 "to the output stream.");
706 pl->set(
"Output Stream", Teuchos::rcpFromRef(std::cout),
707 "A reference-counted pointer to the output stream where all\n"
708 "solver output is sent.");
709 pl->set(
"Timer Label",
static_cast<const char *
>(label_default_),
710 "The string to use as a prefix for the timer labels.");
711 pl->set(
"Orthogonalization",
static_cast<const char *
>(orthoType_default_),
712 "The type of orthogonalization to use: DGKS, ICGS, IMGS");
714 "The constant used by DGKS orthogonalization to determine\n"
715 "whether another step of classical Gram-Schmidt is necessary.");
723template<
class ScalarType,
class MV,
class OP>
727 if (!isSet_) { setParameters( params_ ); }
729 Teuchos::LAPACK<int,ScalarType> lapack;
730 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
731 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
734 "Belos::PCPGSolMgr::solve(): Linear problem is not a valid object.");
737 "Belos::PCPGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
740 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
741 std::vector<int> currIdx(1);
747 problem_->setLSIndex( currIdx );
750 bool isConverged =
true;
754 Teuchos::ParameterList plist;
755 plist.set(
"Saved Blocks", savedBlocks_);
756 plist.set(
"Block Size", 1);
757 plist.set(
"Keep Diagonal",
true);
758 plist.set(
"Initialize Diagonal",
true);
763 Teuchos::RCP<PCPGIter<ScalarType,MV,OP> > pcpg_iter;
769#ifdef BELOS_TEUCHOS_TIME_MONITOR
770 Teuchos::TimeMonitor slvtimer(*timerSolve_);
772 while ( numRHS2Solve > 0 ) {
775 outputTest_->reset();
778 if (R_ == Teuchos::null)
779 R_ = MVT::Clone( *(problem_->getRHS()), 1 );
781 problem_->computeCurrResVec( &*R_ );
787 if( U_ != Teuchos::null ){
792 Teuchos::RCP<MV> cur_soln_vec = problem_->getCurrLHSVec();
793 std::vector<MagnitudeType> rnorm0(1);
794 MVT::MvNorm( *R_, rnorm0 );
797 std::cout <<
"Solver Manager: dimU_ = " << dimU_ << std::endl;
798 Teuchos::SerialDenseMatrix<int,ScalarType> Z( dimU_, 1 );
800 Teuchos::RCP<const MV> Uactive, Cactive;
801 std::vector<int> active_columns( dimU_ );
802 for (
int i=0; i < dimU_; ++i) active_columns[i] = i;
803 Uactive = MVT::CloneView(*U_, active_columns);
804 Cactive = MVT::CloneView(*C_, active_columns);
807 std::cout <<
" Solver Manager : check duality of seed basis " << std::endl;
808 Teuchos::SerialDenseMatrix<int,ScalarType> H( dimU_, dimU_ );
809 MVT::MvTransMv( one, *Uactive, *Cactive, H );
810 H.print( std::cout );
813 MVT::MvTransMv( one, *Uactive, *R_, Z );
814 Teuchos::RCP<MV> tempU = MVT::Clone( *R_, 1 );
815 MVT::MvTimesMatAddMv( one, *Uactive, Z, zero, *tempU );
816 MVT::MvAddMv( one, *tempU, one, *cur_soln_vec, *cur_soln_vec );
817 MVT::MvTimesMatAddMv( one, *Cactive, Z, zero, *tempU );
818 MVT::MvAddMv( -one, *tempU, one, *R_, *R_ );
819 std::vector<MagnitudeType> rnorm(1);
820 MVT::MvNorm( *R_, rnorm );
821 if( rnorm[0] < rnorm0[0] * .001 ){
822 MVT::MvTransMv( one, *Uactive, *R_, Z );
823 MVT::MvTimesMatAddMv( one, *Uactive, Z, zero, *tempU );
824 MVT::MvAddMv( one, *tempU, one, *cur_soln_vec, *cur_soln_vec );
825 MVT::MvTimesMatAddMv( one, *Cactive, Z, zero, *tempU );
826 MVT::MvAddMv( -one, *tempU, one, *R_, *R_ );
828 Uactive = Teuchos::null;
829 Cactive = Teuchos::null;
830 tempU = Teuchos::null;
841 if( U_ != Teuchos::null ) pcpgState.
U = U_;
842 if( C_ != Teuchos::null ) pcpgState.
C = C_;
843 if( dimU_ > 0 ) pcpgState.
curDim = dimU_;
844 pcpg_iter->initialize(pcpgState);
850 if( !dimU_ ) printer_->stream(
Debug) <<
" No recycled subspace available for RHS index " << currIdx[0] << std::endl << std::endl;
851 pcpg_iter->resetNumIters();
853 if( dimU_ > savedBlocks_ )
854 std::cout <<
"Error: dimU_ = " << dimU_ <<
" > savedBlocks_ = " << savedBlocks_ << std::endl;
860 if( debug ) printf(
"********** Calling iterate...\n");
861 pcpg_iter->iterate();
868 if ( convTest_->getStatus() ==
Passed ) {
877 else if ( maxIterTest_->getStatus() ==
Passed ) {
891 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
892 "Belos::PCPGSolMgr::solve(): Invalid return from PCPGIter::iterate().");
895 catch (
const std::exception &e) {
896 printer_->stream(
Errors) <<
"Error! Caught exception in PCPGIter::iterate() at iteration "
897 << pcpg_iter->getNumIters() << std::endl
898 << e.what() << std::endl;
904 Teuchos::RCP<MV> update = pcpg_iter->getCurrentUpdate();
905 problem_->updateSolution( update,
true );
908 problem_->setCurrLS();
916 std::cout <<
"SolverManager: dimU_ " << dimU_ <<
" prevUdim= " << q << std::endl;
918 if( q > deflatedBlocks_ )
919 std::cout <<
"SolverManager: Error deflatedBlocks = " << deflatedBlocks_ << std::endl;
930 rank = ARRQR(dimU_,q, *oldState.
D );
932 std::cout <<
" rank decreased in ARRQR, something to do? " << std::endl;
938 if( dimU_ > deflatedBlocks_ ){
940 if( !deflatedBlocks_ ){
943 dimU_ = deflatedBlocks_;
947 bool Harmonic =
false;
949 Teuchos::RCP<MV> Uorth;
951 std::vector<int> active_cols( dimU_ );
952 for (
int i=0; i < dimU_; ++i) active_cols[i] = i;
955 Uorth = MVT::CloneCopy(*C_, active_cols);
958 Uorth = MVT::CloneCopy(*U_, active_cols);
962 Teuchos::SerialDenseMatrix<int,ScalarType> R(dimU_,dimU_);
963 rank = ortho_->normalize(*Uorth, Teuchos::rcp(&R,
false));
964 Uorth = Teuchos::null;
970 "Belos::PCPGSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
974 Teuchos::SerialDenseMatrix<int,ScalarType> VT;
975 Teuchos::SerialDenseMatrix<int,ScalarType> Ur;
979 if( problem_->isHermitian() ) lrwork = dimU_;
980 std::vector<ScalarType> work(lwork);
981 std::vector<ScalarType> Svec(dimU_);
982 std::vector<ScalarType> rwork(lrwork);
983 lapack.GESVD(
'N',
'O',
984 R.numRows(),R.numCols(),R.values(), R.numRows(),
992 "Belos::PCPGSolMgr::solve(): LAPACK _GESVD failed to compute singular values.");
994 if( work[0] != 67. * dimU_ )
995 std::cout <<
" SVD " << dimU_ <<
" lwork " << work[0] << std::endl;
996 for(
int i=0; i< dimU_; i++)
997 std::cout << i <<
" " << Svec[i] << std::endl;
999 Teuchos::SerialDenseMatrix<int,ScalarType> wholeV( R, Teuchos::TRANS);
1001 int startRow = 0, startCol = 0;
1003 startCol = dimU_ - deflatedBlocks_;
1005 Teuchos::SerialDenseMatrix<int,ScalarType> V(Teuchos::Copy,
1011 std::vector<int> active_columns( dimU_ );
1012 std::vector<int> def_cols( deflatedBlocks_ );
1013 for (
int i=0; i < dimU_; ++i) active_columns[i] = i;
1014 for (
int i=0; i < deflatedBlocks_; ++i) def_cols[i] = i;
1016 Teuchos::RCP<MV> Uactive = MVT::CloneViewNonConst(*U_, def_cols);
1017 Teuchos::RCP<MV> Ucopy = MVT::CloneCopy( *U_, active_columns );
1018 MVT::MvTimesMatAddMv( one, *Ucopy, V, zero, *Uactive );
1019 Ucopy = Teuchos::null;
1020 Uactive = Teuchos::null;
1021 Teuchos::RCP<MV> Cactive = MVT::CloneViewNonConst(*C_, def_cols);
1022 Teuchos::RCP<MV> Ccopy = MVT::CloneCopy( *C_, active_columns );
1023 MVT::MvTimesMatAddMv( one, *Ccopy, V, zero, *Cactive );
1024 Ccopy = Teuchos::null;
1025 Cactive = Teuchos::null;
1026 dimU_ = deflatedBlocks_;
1028 printer_->stream(
Debug) <<
" Generated recycled subspace using RHS index " << currIdx[0] <<
" of dimension " << dimU_ << std::endl << std::endl;
1031 problem_->setCurrLS();
1035 if ( numRHS2Solve > 0 ) {
1039 problem_->setLSIndex( currIdx );
1042 currIdx.resize( numRHS2Solve );
1051#ifdef BELOS_TEUCHOS_TIME_MONITOR
1056 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails) );
1061 using Teuchos::rcp_dynamic_cast;
1064 const std::vector<MagnitudeType>* pTestValues =
1065 rcp_dynamic_cast<conv_test_type>(convTest_)->
getTestValue();
1067 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1068 "Belos::PCPGSolMgr::solve(): The convergence test's getTestValue() "
1069 "method returned NULL. Please report this bug to the Belos developers.");
1071 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1072 "Belos::PCPGSolMgr::solve(): The convergence test's getTestValue() "
1073 "method returned a vector of length zero. Please report this bug to the "
1074 "Belos developers.");
1079 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1083 numIters_ = maxIterTest_->getNumIters();
1094template<
class ScalarType,
class MV,
class OP>
1098 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1099 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1102 Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
1103 Teuchos::SerialDenseMatrix<int,ScalarType> gamma( 1, 1 );
1104 Teuchos::SerialDenseMatrix<int,ScalarType> anorm( 1, 1 );
1105 std::vector<int> curind(1);
1106 std::vector<int> ipiv(p - q);
1107 std::vector<ScalarType> Pivots(p);
1109 ScalarType rteps = 1.5e-8;
1112 for( i = q ; i < p ; i++ ){
1115 RCP<MV> P = MVT::CloneViewNonConst(*U_,curind);
1116 RCP<MV> AP = MVT::CloneViewNonConst(*C_,curind);
1117 anorm(0,0) = one / Teuchos::ScalarTraits<ScalarType>::squareroot( D(i-q,i-q) ) ;
1118 MVT::MvAddMv( anorm(0,0), *P, zero, *AP, *P );
1119 MVT::MvAddMv( zero, *P, anorm(0,0), *AP, *AP );
1123 for( i = q ; i < p ; i++ ){
1124 if( q < i && i < p-1 ){
1127 for( j = i+1 ; j < p ; j++ ){
1128 const int k = ipiv[j-q];
1129 if( Pivots[k] > Pivots[l] ){
1136 ipiv[imax-q] = ipiv[i-q];
1142 if( Pivots[k] > 1.5625e-2 ){
1143 anorm(0,0) = Pivots[k];
1147 RCP<const MV> P = MVT::CloneView(*U_,curind);
1148 RCP<const MV> AP = MVT::CloneView(*C_,curind);
1149 MVT::MvTransMv( one, *P, *AP, anorm );
1150 anorm(0,0) = Teuchos::ScalarTraits<ScalarType>::squareroot( anorm(0,0) ) ;
1152 if( rteps <= anorm(0,0) && anorm(0,0) < 9.765625e-4){
1160 std::cout <<
"ARRQR: Bad case not implemented" << std::endl;
1162 if( anorm(0,0) < rteps ){
1163 std::cout <<
"ARRQR : deficient case not implemented " << std::endl;
1171 RCP<MV> P = MVT::CloneViewNonConst(*U_,curind);
1172 RCP<MV> AP = MVT::CloneViewNonConst(*C_,curind);
1173 MVT::MvAddMv( anorm(0,0), *P, zero, *AP, *P );
1174 MVT::MvAddMv( zero, *P, anorm(0,0), *AP, *AP );
1178 P = MVT::CloneViewNonConst(*U_,curind);
1179 AP = MVT::CloneViewNonConst(*C_,curind);
1180 for( j = i+1 ; j < p ; j++ ){
1183 RCP<MV> Q = MVT::CloneViewNonConst(*U_,curind);
1184 MVT::MvTransMv( one, *Q, *AP, alpha);
1185 MVT::MvAddMv( -alpha(0,0), *P, one, *Q, *Q );
1187 RCP<MV> AQ = MVT::CloneViewNonConst(*C_,curind);
1188 MVT::MvAddMv( -alpha(0,0), *AP, one, *AQ, *AQ );
1190 gamma(0,0) = ( Pivots[l] - alpha(0,0))*( Pivots[l] + alpha(0,0));
1191 if( gamma(0,0) > 0){
1192 Pivots[l] = Teuchos::ScalarTraits<ScalarType>::squareroot( gamma(0,0) );
1203template<
class ScalarType,
class MV,
class OP>
1206 std::ostringstream oss;
1207 oss <<
"Belos::PCPGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<
">";
1209 oss <<
"Ortho Type='"<<orthoType_;
Belos header file which uses auto-configuration information to include necessary C++ headers.
Class which describes the linear problem to be solved by the iterative solver.
Class which manages the output and verbosity of the Belos solvers.
Belos concrete class to iterate Preconditioned Conjugate Projected Gradients.
Pure virtual base class which describes the basic interface for a solver manager.
Belos::StatusTest for logically combining several status tests.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Belos::StatusTest class for specifying a maximum number of iterations.
A factory class for generating StatusTestOutput objects.
Collection of types and exceptions used within the Belos solvers.
Parent class to all Belos exceptions.
Type traits class that says whether Teuchos::LAPACK has a valid implementation for the given ScalarTy...
Base class for Belos::SolverManager subclasses which normally can only compile with real ScalarType t...
A linear system to solve, and its associated information.
Traits class which defines basic operations on multivectors.
Class which defines basic traits for the operator type.
Enumeration of all valid Belos (Mat)OrthoManager classes.
Teuchos::RCP< Belos::MatOrthoManager< Scalar, MV, OP > > makeMatOrthoManager(const std::string &ortho, const Teuchos::RCP< const OP > &M, const Teuchos::RCP< OutputManager< Scalar > > &, const std::string &label, const Teuchos::RCP< Teuchos::ParameterList > ¶ms)
Return an instance of the specified MatOrthoManager subclass.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
This class implements the PCPG iteration, where a single-std::vector Krylov subspace is constructed....
PCPGSolMgrLAPACKFailure is thrown when a nonzero value is retuned from an LAPACK call.
PCPGSolMgrLAPACKFailure(const std::string &what_arg)
PCPGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i.e.
PCPGSolMgrLinearProblemFailure(const std::string &what_arg)
PCPGSolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate orthonormal...
PCPGSolMgrOrthoFailure(const std::string &what_arg)
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get current linear problem being solved for in this object.
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const
Get a parameter list containing the current parameters for this object.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem)
Set the linear problem that needs to be solved.
bool isLOADetected() const
Return whether a loss of accuracy was detected by this solver during the most current solve.
int getNumIters() const
Get the iteration count for the most recent call to solve().
virtual ~PCPGSolMgr()
Destructor.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
virtual Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const
clone for Inverted Injection (DII)
MagnitudeType achievedTol() const
Tolerance achieved by the last solve() invocation.
void reset(const ResetType type)
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
PCPG iterative linear solver.
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
PCPGSolMgr(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.
An implementation of StatusTestResNorm using a family of residual norms.
const std::vector< MagnitudeType > * getTestValue() const
Returns the test value, , computed in most recent call to CheckStatus.
int setTolerance(MagnitudeType tolerance)
Set the value of the tolerance.
A Belos::StatusTest class for specifying a maximum number of iterations.
A factory class for generating StatusTestOutput objects.
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > create(const Teuchos::RCP< OutputManager< ScalarType > > &printer, Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test, int mod, int printStates)
Create the StatusTestOutput object specified by the outputStyle.
ReturnType
Whether the Belos solve converged for all linear systems.
ResetType
How to reset the solver.
Default parameters common to most Belos solvers.
static const double convTol
Default convergence tolerance.
static const double orthoKappa
DGKS orthogonalization constant.
Structure to contain pointers to PCPGIter state variables.
int prevUdim
Number of block columns in matrices C and U before current iteration.
Teuchos::RCP< MV > U
The recycled subspace.
Teuchos::RCP< MV > C
C = AU, U spans recycled subspace.
int curDim
The current dimension of the reduction.
Teuchos::RCP< MV > R
The current residual.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > D
The current Hessenberg matrix.