45#ifndef BELOS_BLOCK_GCRODR_SOLMGR_HPP
46#define BELOS_BLOCK_GCRODR_SOLMGR_HPP
62#include "Teuchos_LAPACK.hpp"
63#include "Teuchos_as.hpp"
64#ifdef BELOS_TEUCHOS_TIME_MONITOR
65#include "Teuchos_TimeMonitor.hpp"
125template<
class ScalarType,
class MV,
class OP>
131 typedef Teuchos::ScalarTraits<ScalarType> SCT;
132 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
133 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
134 typedef Teuchos::ScalarTraits<MagnitudeType> SMT;
136 typedef Teuchos::SerialDenseMatrix<int,ScalarType> SDM;
137 typedef Teuchos::SerialDenseVector<int,ScalarType> SDV;
176 const Teuchos::RCP<Teuchos::ParameterList> &pl);
223 TEUCHOS_TEST_FOR_EXCEPTION(problem.is_null(), std::invalid_argument,
224 "Belos::BlockGCRODRSolMgr::setProblem: The input LinearProblem cannot be null.");
227 if (! problem->isProblemSet()) {
228 const bool success = problem->setProblem();
229 TEUCHOS_TEST_FOR_EXCEPTION(success, std::runtime_error,
230 "Belos::BlockGCRODRSolMgr::setProblem: Calling the input LinearProblem's setProblem() method failed. This likely means that the "
231 "LinearProblem has a missing (null) matrix A, solution vector X, or right-hand side vector B. Please set these items in the LinearProblem and try again.");
238 void setParameters(
const Teuchos::RCP<Teuchos::ParameterList> ¶ms );
253 problem_->setProblem();
291 void initializeStateStorage();
309 int getHarmonicVecsKryl (
int m,
const SDM& HH, SDM& PP);
314 int getHarmonicVecsAugKryl (
int keff,
317 const Teuchos::RCP<const MV>& VV,
321 void sort (std::vector<MagnitudeType>& dlist,
int n, std::vector<int>& iperm);
324 Teuchos::LAPACK<int,ScalarType> lapack;
327 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
330 Teuchos::RCP<OutputManager<ScalarType> > printer_;
331 Teuchos::RCP<std::ostream> outputStream_;
334 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
335 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
336 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
337 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > expConvTest_, impConvTest_;
338 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
341 ortho_factory_type orthoFactory_;
348 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
351 Teuchos::RCP<Teuchos::ParameterList> params_;
363 mutable Teuchos::RCP<const Teuchos::ParameterList> defaultParams_;
366 static const bool adaptiveBlockSize_default_;
367 static const std::string recycleMethod_default_;
370 MagnitudeType convTol_, orthoKappa_, achievedTol_;
371 int blockSize_, maxRestarts_, maxIters_, numIters_;
372 int verbosity_, outputStyle_, outputFreq_;
373 bool adaptiveBlockSize_;
374 std::string orthoType_, recycleMethod_;
375 std::string impResScale_, expResScale_;
383 int numBlocks_, recycledBlocks_;
394 Teuchos::RCP<MV> U_, C_;
397 Teuchos::RCP<MV> U1_, C1_;
400 Teuchos::RCP<SDM > G_;
401 Teuchos::RCP<SDM > H_;
402 Teuchos::RCP<SDM > B_;
403 Teuchos::RCP<SDM > PP_;
404 Teuchos::RCP<SDM > HP_;
405 std::vector<ScalarType> tau_;
406 std::vector<ScalarType> work_;
407 Teuchos::RCP<SDM > F_;
408 std::vector<int> ipiv_;
411 Teuchos::RCP<Teuchos::Time> timerSolve_;
420 bool builtRecycleSpace_;
426 template<
class ScalarType,
class MV,
class OP>
427 const bool BlockGCRODRSolMgr<ScalarType,MV,OP>::adaptiveBlockSize_default_ =
true;
429 template<
class ScalarType,
class MV,
class OP>
430 const std::string BlockGCRODRSolMgr<ScalarType,MV,OP>::recycleMethod_default_ =
"harmvecs";
436 template<
class ScalarType,
class MV,
class OP>
442 template<
class ScalarType,
class MV,
class OP>
445 const Teuchos::RCP<Teuchos::ParameterList> &pl ) {
450 TEUCHOS_TEST_FOR_EXCEPTION(problem.is_null(), std::invalid_argument,
451 "Belos::BlockGCRODR constructor: The solver manager's constructor needs "
452 "the linear problem argument 'problem' to be nonnull.");
462 template<
class ScalarType,
class MV,
class OP>
464 adaptiveBlockSize_ = adaptiveBlockSize_default_;
465 recycleMethod_ = recycleMethod_default_;
467 loaDetected_ =
false;
468 builtRecycleSpace_ =
false;
539 template<
class ScalarType,
class MV,
class OP>
541 std::ostringstream oss;
542 oss <<
"Belos::BlockGCRODRSolMgr<" << SCT::name() <<
", ...>";
544 oss <<
"Ortho Type='"<<orthoType_ ;
545 oss <<
", Num Blocks=" <<numBlocks_;
546 oss <<
", Block Size=" <<blockSize_;
547 oss <<
", Num Recycle Blocks=" << recycledBlocks_;
548 oss <<
", Max Restarts=" << maxRestarts_;
553 template<
class ScalarType,
class MV,
class OP>
554 Teuchos::RCP<const Teuchos::ParameterList>
556 using Teuchos::ParameterList;
557 using Teuchos::parameterList;
559 using Teuchos::rcpFromRef;
561 if (defaultParams_.is_null()) {
562 RCP<ParameterList> pl = parameterList ();
564 const MagnitudeType convTol = SMT::squareroot (SCT::magnitude (SCT::eps()));
565 const int maxRestarts = 1000;
566 const int maxIters = 5000;
567 const int blockSize = 2;
568 const int numBlocks = 100;
569 const int numRecycledBlocks = 25;
572 const int outputFreq = 1;
573 RCP<std::ostream> outputStream = rcpFromRef (std::cout);
574 const std::string impResScale (
"Norm of Preconditioned Initial Residual");
575 const std::string expResScale (
"Norm of Initial Residual");
576 const std::string timerLabel (
"Belos");
577 const std::string orthoType (
"ICGS");
578 RCP<const ParameterList> orthoParams = orthoFactory_.getDefaultParameters (orthoType);
585 const MagnitudeType orthoKappa = -SMT::one();
588 pl->set (
"Convergence Tolerance", convTol,
589 "The tolerance that the solver needs to achieve in order for "
590 "the linear system(s) to be declared converged. The meaning "
591 "of this tolerance depends on the convergence test details.");
592 pl->set(
"Maximum Restarts", maxRestarts,
593 "The maximum number of restart cycles (not counting the first) "
594 "allowed for each set of right-hand sides solved.");
595 pl->set(
"Maximum Iterations", maxIters,
596 "The maximum number of iterations allowed for each set of "
597 "right-hand sides solved.");
598 pl->set(
"Block Size", blockSize,
599 "How many linear systems to solve at once.");
600 pl->set(
"Num Blocks", numBlocks,
601 "The maximum number of blocks allowed in the Krylov subspace "
602 "for each set of right-hand sides solved. This includes the "
603 "initial block. Total Krylov basis storage, not counting the "
604 "recycled subspace, is \"Num Blocks\" times \"Block Size\".");
605 pl->set(
"Num Recycled Blocks", numRecycledBlocks,
606 "The maximum number of vectors in the recycled subspace." );
607 pl->set(
"Verbosity", verbosity,
608 "What type(s) of solver information should be written "
609 "to the output stream.");
610 pl->set(
"Output Style", outputStyle,
611 "What style is used for the solver information to write "
612 "to the output stream.");
613 pl->set(
"Output Frequency", outputFreq,
614 "How often convergence information should be written "
615 "to the output stream.");
616 pl->set(
"Output Stream", outputStream,
617 "A reference-counted pointer to the output stream where all "
618 "solver output is sent.");
619 pl->set(
"Implicit Residual Scaling", impResScale,
620 "The type of scaling used in the implicit residual convergence test.");
621 pl->set(
"Explicit Residual Scaling", expResScale,
622 "The type of scaling used in the explicit residual convergence test.");
623 pl->set(
"Timer Label", timerLabel,
624 "The string to use as a prefix for the timer labels.");
625 pl->set(
"Orthogonalization", orthoType,
626 "The orthogonalization method to use. Valid options: " +
627 orthoFactory_.validNamesString());
628 pl->set (
"Orthogonalization Parameters", *orthoParams,
629 "Sublist of parameters specific to the orthogonalization method to use.");
630 pl->set(
"Orthogonalization Constant", orthoKappa,
631 "When using DGKS orthogonalization: the \"depTol\" constant, used "
632 "to determine whether another step of classical Gram-Schmidt is "
633 "necessary. Otherwise ignored. Nonpositive values are ignored.");
636 return defaultParams_;
639 template<
class ScalarType,
class MV,
class OP>
642 setParameters (
const Teuchos::RCP<Teuchos::ParameterList> ¶ms) {
643 using Teuchos::isParameterType;
644 using Teuchos::getParameter;
646 using Teuchos::ParameterList;
647 using Teuchos::parameterList;
650 using Teuchos::rcp_dynamic_cast;
651 using Teuchos::rcpFromRef;
652 using Teuchos::Exceptions::InvalidParameter;
653 using Teuchos::Exceptions::InvalidParameterName;
654 using Teuchos::Exceptions::InvalidParameterType;
657 RCP<const ParameterList> defaultParams = getValidParameters();
659 if (params.is_null()) {
664 params_ = parameterList (*defaultParams);
676 params->validateParametersAndSetDefaults (*defaultParams);
692 const int maxRestarts = params_->get<
int> (
"Maximum Restarts");
693 TEUCHOS_TEST_FOR_EXCEPTION(maxRestarts <= 0, std::invalid_argument,
694 "Belos::BlockGCRODRSolMgr: The \"Maximum Restarts\" parameter "
695 "must be nonnegative, but you specified a negative value of "
696 << maxRestarts <<
".");
698 const int maxIters = params_->get<
int> (
"Maximum Iterations");
699 TEUCHOS_TEST_FOR_EXCEPTION(maxIters <= 0, std::invalid_argument,
700 "Belos::BlockGCRODRSolMgr: The \"Maximum Iterations\" parameter "
701 "must be positive, but you specified a nonpositive value of "
704 const int numBlocks = params_->get<
int> (
"Num Blocks");
705 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument,
706 "Belos::BlockGCRODRSolMgr: The \"Num Blocks\" parameter must be "
707 "positive, but you specified a nonpositive value of " << numBlocks
710 const int blockSize = params_->get<
int> (
"Block Size");
711 TEUCHOS_TEST_FOR_EXCEPTION(blockSize <= 0, std::invalid_argument,
712 "Belos::BlockGCRODRSolMgr: The \"Block Size\" parameter must be "
713 "positive, but you specified a nonpositive value of " << blockSize
716 const int recycledBlocks = params_->get<
int> (
"Num Recycled Blocks");
717 TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks <= 0, std::invalid_argument,
718 "Belos::BlockGCRODRSolMgr: The \"Num Recycled Blocks\" parameter must "
719 "be positive, but you specified a nonpositive value of "
720 << recycledBlocks <<
".");
721 TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks >= numBlocks,
722 std::invalid_argument,
"Belos::BlockGCRODRSolMgr: The \"Num Recycled "
723 "Blocks\" parameter must be less than the \"Num Blocks\" parameter, "
724 "but you specified \"Num Recycled Blocks\" = " << recycledBlocks
725 <<
" and \"Num Blocks\" = " << numBlocks <<
".");
729 maxRestarts_ = maxRestarts;
730 maxIters_ = maxIters;
731 numBlocks_ = numBlocks;
732 blockSize_ = blockSize;
733 recycledBlocks_ = recycledBlocks;
740 std::string tempLabel = params_->get<std::string> (
"Timer Label");
741 const bool labelChanged = (tempLabel != label_);
743#ifdef BELOS_TEUCHOS_TIME_MONITOR
744 std::string solveLabel = label_ +
": BlockGCRODRSolMgr total solve time";
745 if (timerSolve_.is_null()) {
747 timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
748 }
else if (labelChanged) {
754 Teuchos::TimeMonitor::clearCounter (solveLabel);
755 timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
763 if (params_->isParameter (
"Verbosity")) {
764 if (isParameterType<int> (*params_,
"Verbosity")) {
765 verbosity_ = params_->get<
int> (
"Verbosity");
768 verbosity_ = (int) getParameter<MsgType> (*params_,
"Verbosity");
775 if (params_->isParameter (
"Output Style")) {
776 if (isParameterType<int> (*params_,
"Output Style")) {
777 outputStyle_ = params_->get<
int> (
"Output Style");
780 outputStyle_ = (int) getParameter<OutputType> (*params_,
"Output Style");
808 if (params_->isParameter (
"Output Stream")) {
810 outputStream_ = getParameter<RCP<std::ostream> > (*params_,
"Output Stream");
812 catch (InvalidParameter&) {
813 outputStream_ = rcpFromRef (std::cout);
820 if (outputStream_.is_null()) {
821 outputStream_ = rcp (
new Teuchos::oblackholestream);
825 outputFreq_ = params_->get<
int> (
"Output Frequency");
828 if (! outputTest_.is_null()) {
829 outputTest_->setOutputFrequency (outputFreq_);
837 if (printer_.is_null()) {
840 printer_->setVerbosity (verbosity_);
841 printer_->setOStream (outputStream_);
852 if (params_->isParameter (
"Orthogonalization")) {
853 const std::string& tempOrthoType =
854 params_->get<std::string> (
"Orthogonalization");
856 if (! orthoFactory_.isValidName (tempOrthoType)) {
857 std::ostringstream os;
858 os <<
"Belos::BlockGCRODRSolMgr: Invalid orthogonalization name \""
859 << tempOrthoType <<
"\". The following are valid options "
860 <<
"for the \"Orthogonalization\" name parameter: ";
861 orthoFactory_.printValidNames (os);
862 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument, os.str());
864 if (tempOrthoType != orthoType_) {
866 orthoType_ = tempOrthoType;
883 RCP<ParameterList> orthoParams = sublist (params,
"Orthogonalization Parameters",
true);
884 TEUCHOS_TEST_FOR_EXCEPTION(orthoParams.is_null(), std::logic_error,
885 "Failed to get orthogonalization parameters. "
886 "Please report this bug to the Belos developers.");
912 ortho_ = orthoFactory_.makeMatOrthoManager (orthoType_, null, printer_,
913 label_, orthoParams);
925 if (params_->isParameter (
"Orthogonalization Constant")) {
926 const MagnitudeType orthoKappa =
927 params_->get<MagnitudeType> (
"Orthogonalization Constant");
928 if (orthoKappa > 0) {
929 orthoKappa_ = orthoKappa;
932 if (orthoType_ ==
"DGKS" && ! ortho_.is_null()) {
938 rcp_dynamic_cast<ortho_man_type>(ortho_)->
setDepTol (orthoKappa_);
948 convTol_ = params_->get<MagnitudeType> (
"Convergence Tolerance");
949 if (! impConvTest_.is_null()) {
952 if (! expConvTest_.is_null()) {
953 expConvTest_->setTolerance (convTol_);
957 if (params_->isParameter (
"Implicit Residual Scaling")) {
958 std::string tempImpResScale =
959 getParameter<std::string> (*params_,
"Implicit Residual Scaling");
962 if (impResScale_ != tempImpResScale) {
964 impResScale_ = tempImpResScale;
966 if (! impConvTest_.is_null()) {
979 if (params_->isParameter(
"Explicit Residual Scaling")) {
980 std::string tempExpResScale =
981 getParameter<std::string> (*params_,
"Explicit Residual Scaling");
984 if (expResScale_ != tempExpResScale) {
986 expResScale_ = tempExpResScale;
988 if (! expConvTest_.is_null()) {
1006 if (maxIterTest_.is_null()) {
1009 maxIterTest_->setMaxIters (maxIters_);
1014 if (impConvTest_.is_null()) {
1015 impConvTest_ = rcp (
new StatusTestResNorm_t (convTol_));
1021 if (expConvTest_.is_null()) {
1022 expConvTest_ = rcp (
new StatusTestResNorm_t (convTol_));
1023 expConvTest_->defineResForm (StatusTestResNorm_t::Explicit,
Belos::TwoNorm);
1029 if (convTest_.is_null()) {
1030 convTest_ = rcp (
new StatusTestCombo_t (StatusTestCombo_t::SEQ,
1038 sTest_ = rcp (
new StatusTestCombo_t (StatusTestCombo_t::OR,
1039 maxIterTest_, convTest_));
1043 outputTest_ = stoFactory.
create (printer_, sTest_, outputFreq_,
1047 std::string solverDesc =
"Block GCRODR ";
1048 outputTest_->setSolverDesc (solverDesc);
1055 template<
class ScalarType,
class MV,
class OP>
1060 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1063 Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
1066 int KrylSpaDim = (numBlocks_ - 1) * blockSize_;
1069 int augSpaDim = KrylSpaDim + recycledBlocks_ + 1;
1072 int augSpaImgDim = KrylSpaDim + blockSize_ + recycledBlocks_+1;
1096 if (U_ == Teuchos::null) {
1097 U_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1101 if (MVT::GetNumberVecs(*U_) < recycledBlocks_+1) {
1102 Teuchos::RCP<const MV> tmp = U_;
1103 U_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1108 if (C_ == Teuchos::null) {
1109 C_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1113 if (MVT::GetNumberVecs(*C_) < recycledBlocks_+1) {
1114 Teuchos::RCP<const MV> tmp = C_;
1115 C_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1120 if (U1_ == Teuchos::null) {
1121 U1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1125 if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1126 Teuchos::RCP<const MV> tmp = U1_;
1127 U1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1132 if (C1_ == Teuchos::null) {
1133 C1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1137 if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1138 Teuchos::RCP<const MV> tmp = C1_;
1139 C1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1144 if (R_ == Teuchos::null){
1145 R_ = MVT::Clone( *rhsMV, blockSize_ );
1151 if (G_ == Teuchos::null){
1152 G_ = Teuchos::rcp(
new SDM( augSpaImgDim, augSpaDim ) );
1155 if ( (G_->numRows() != augSpaImgDim) || (G_->numCols() != augSpaDim) )
1157 G_->reshape( augSpaImgDim, augSpaDim );
1159 G_->putScalar(zero);
1163 if (H_ == Teuchos::null){
1164 H_ = Teuchos::rcp (
new SDM ( Teuchos::View, *G_, KrylSpaDim + blockSize_, KrylSpaDim, recycledBlocks_+1 ,recycledBlocks_+1 ) );
1168 if (F_ == Teuchos::null){
1169 F_ = Teuchos::rcp(
new SDM( recycledBlocks_+1, recycledBlocks_+1 ) );
1172 if ( (F_->numRows() != recycledBlocks_+1) || (F_->numCols() != recycledBlocks_+1) ){
1173 F_->reshape( recycledBlocks_+1, recycledBlocks_+1 );
1176 F_->putScalar(zero);
1179 if (PP_ == Teuchos::null){
1180 PP_ = Teuchos::rcp(
new SDM( augSpaImgDim, recycledBlocks_+1 ) );
1183 if ( (PP_->numRows() != augSpaImgDim) || (PP_->numCols() != recycledBlocks_+1) ){
1184 PP_->reshape( augSpaImgDim, recycledBlocks_+1 );
1189 if (HP_ == Teuchos::null)
1190 HP_ = Teuchos::rcp(
new SDM( augSpaImgDim, augSpaDim ) );
1192 if ( (HP_->numRows() != augSpaImgDim) || (HP_->numCols() != augSpaDim) ){
1193 HP_->reshape( augSpaImgDim, augSpaDim );
1198 tau_.resize(recycledBlocks_+1);
1201 work_.resize(recycledBlocks_+1);
1204 ipiv_.resize(recycledBlocks_+1);
1208template<
class ScalarType,
class MV,
class OP>
1209void BlockGCRODRSolMgr<ScalarType,MV,OP>::buildRecycleSpaceKryl(
int& keff, Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > block_gmres_iter){
1211 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1212 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1214 int p = block_gmres_iter->getState().curDim;
1215 std::vector<int> index(keff);
1219 SDM HH(Teuchos::Copy, *H_, p+blockSize_, p);
1220 if(recycledBlocks_ >= p + blockSize_){
1224 for (
int ii=0; ii < p; ++ii) index[ii] = ii;
1225 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1226 MVT::SetBlock(*V_, index, *Utmp);
1231 Teuchos::RCP<SDM > PPtmp = Teuchos::rcp (
new SDM ( Teuchos::View, *PP_, p, recycledBlocks_+1 ) );
1232 if(recycleMethod_ ==
"harmvecs"){
1233 keff = getHarmonicVecsKryl(p, HH, *PPtmp);
1234 printer_->stream(
Debug) <<
"keff = " << keff << std::endl;
1237PPtmp = Teuchos::rcp (
new SDM ( Teuchos::View, *PP_, p, keff ) );
1240for (
int ii=0; ii<keff; ++ii) index[ii] = ii;
1241Teuchos::RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1242Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1243Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1245for (
int ii=0; ii < p; ++ii) index[ii] = ii;
1246Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1250MVT::MvTimesMatAddMv( one, *Vtmp, *PPtmp, zero, *U1tmp );
1253SDM HPtmp( Teuchos::View, *HP_, p+blockSize_, keff );
1254HPtmp.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, *H_, *PPtmp, zero );
1258lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1259TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
1262 lwork =
static_cast<int> (Teuchos::ScalarTraits<ScalarType>::magnitude (work_[0]));
1264lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1265TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
1269SDM Rtmp( Teuchos::View, *F_, keff, keff );
1270for(
int ii=0;ii<keff;ii++) {
for(
int jj=ii;jj<keff;jj++) Rtmp(ii,jj) = HPtmp(ii,jj); }
1272lapack.UNGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1273TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _UNGQR failed to construct the Q factor.");
1277 index.resize( p+blockSize_ );
1278 for (
int ii=0; ii < (p+blockSize_); ++ii) { index[ii] = ii; }
1279 Vtmp = MVT::CloneView( *V_, index );
1280 MVT::MvTimesMatAddMv( one, *Vtmp, HPtmp, zero, *Ctmp );
1287ipiv_.resize(Rtmp.numRows());
1288lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
1289TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1291lwork = Rtmp.numRows();
1293lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1296MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
1302template<
class ScalarType,
class MV,
class OP>
1303void BlockGCRODRSolMgr<ScalarType,MV,OP>::buildRecycleSpaceAugKryl(Teuchos::RCP<BlockGCRODRIter<ScalarType,MV,OP> > block_gcrodr_iter){
1304 const MagnitudeType one = Teuchos::ScalarTraits<ScalarType>::one();
1305 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1307 std::vector<MagnitudeType> d(keff);
1308 std::vector<ScalarType> dscalar(keff);
1309 std::vector<int> index(numBlocks_+1);
1312 BlockGCRODRIterState<ScalarType,MV> oldState = block_gcrodr_iter->getState();
1313 int p = oldState.curDim;
1318 if(recycledBlocks_ >= keff + p){
1321 for (
int ii=0; ii < p; ++ii) { index[ii] = keff+ii; }
1322 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1323 for (
int ii=0; ii < p; ++ii) { index[ii] =ii; }
1324 MVT::SetBlock(*V_, index, *Utmp);
1331 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1332 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1334 dscalar.resize(keff);
1335 MVT::MvNorm( *Utmp, d );
1336 for (
int i=0; i<keff; ++i) {
1338 dscalar[i] = (ScalarType)d[i];
1340 MVT::MvScale( *Utmp, dscalar );
1345 Teuchos::RCP<SDM> Gtmp = Teuchos::rcp(
new SDM( Teuchos::View, *G_, p+keff, p+keff-blockSize_ ) );
1348 for (
int i=0; i<keff; ++i)
1349 (*Gtmp)(i,i) = d[i];
1356 SDM PPtmp( Teuchos::View, *PP_, p+keff-blockSize_, recycledBlocks_+1 );
1357 keff_new = getHarmonicVecsAugKryl( keff, p-blockSize_, *Gtmp, oldState.V, PPtmp );
1364 Teuchos::RCP<MV> U1tmp;
1366 index.resize( keff );
1367 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1368 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1369 index.resize( keff_new );
1370 for (
int ii=0; ii<keff_new; ++ii) { index[ii] = ii; }
1371 U1tmp = MVT::CloneViewNonConst( *U1_, index );
1372 SDM PPtmp( Teuchos::View, *PP_, keff, keff_new );
1373 MVT::MvTimesMatAddMv( one, *Utmp, PPtmp, zero, *U1tmp );
1378 index.resize(p-blockSize_);
1379 for (
int ii=0; ii < p-blockSize_; ii++) { index[ii] = ii; }
1380 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1381 SDM PPtmp( Teuchos::View, *PP_, p-blockSize_, keff_new, keff );
1382 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *U1tmp );
1386 SDM HPtmp( Teuchos::View, *HP_, p+keff, keff_new );
1388 SDM PPtmp( Teuchos::View, *PP_, p-blockSize_+keff, keff_new );
1389 HPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*Gtmp,PPtmp,zero);
1393 int info = 0, lwork = -1;
1394 tau_.resize(keff_new);
1395 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1396 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
1398 lwork =
static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work_[0]));
1399 work_.resize(lwork);
1400 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1401 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
1405 SDM Ftmp( Teuchos::View, *F_, keff_new, keff_new );
1406 for(
int i=0;i<keff_new;i++) {
for(
int j=i;j<keff_new;j++) Ftmp(i,j) = HPtmp(i,j); }
1408 lapack.UNGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1409 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _UNGQR failed to construct the Q factor.");
1416 Teuchos::RCP<MV> C1tmp;
1419 for (
int i=0; i < keff; i++) { index[i] = i; }
1420 Teuchos::RCP<const MV> Ctmp = MVT::CloneView( *C_, index );
1421 index.resize(keff_new);
1422 for (
int i=0; i < keff_new; i++) { index[i] = i; }
1423 C1tmp = MVT::CloneViewNonConst( *C1_, index );
1424 SDM PPtmp( Teuchos::View, *HP_, keff, keff_new );
1425 MVT::MvTimesMatAddMv( one, *Ctmp, PPtmp, zero, *C1tmp );
1430 for (
int i=0; i < p; ++i) { index[i] = i; }
1431 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1432 SDM PPtmp( Teuchos::View, *HP_, p, keff_new, keff, 0 );
1433 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *C1tmp );
1442 ipiv_.resize(Ftmp.numRows());
1443 lapack.GETRF(Ftmp.numRows(),Ftmp.numCols(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&info);
1444 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1447 lwork = Ftmp.numRows();
1448 work_.resize(lwork);
1449 lapack.GETRI(Ftmp.numRows(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1450 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRI failed to compute an LU factorization.");
1453 index.resize(keff_new);
1454 for (
int i=0; i < keff_new; i++) { index[i] = i; }
1455 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1456 MVT::MvTimesMatAddMv( one, *U1tmp, Ftmp, zero, *Utmp );
1461 if (keff != keff_new) {
1463 block_gcrodr_iter->setSize( keff, numBlocks_ );
1465 SDM b1( Teuchos::View, *G_, recycledBlocks_+2, 1, 0, recycledBlocks_ );
1471template<
class ScalarType,
class MV,
class OP>
1472int BlockGCRODRSolMgr<ScalarType,MV,OP>::getHarmonicVecsAugKryl(
int keff,
int m,
const SDM& GG,
const Teuchos::RCP<const MV>& VV, SDM& PP){
1474 int m2 = GG.numCols();
1475 bool xtraVec =
false;
1476 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1477 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1478 std::vector<int> index;
1481 std::vector<MagnitudeType> wr(m2), wi(m2);
1484 std::vector<MagnitudeType> w(m2);
1487 SDM vr(m2,m2,
false);
1490 std::vector<int> iperm(m2);
1493 builtRecycleSpace_ =
true;
1499 B.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,GG,GG,zero);
1503 SDM A_tmp( keff+m+blockSize_, keff+m );
1508 for (
int i=0; i<keff; ++i) { index[i] = i; }
1509 Teuchos::RCP<const MV> Ctmp = MVT::CloneView( *C_, index );
1510 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1511 SDM A11( Teuchos::View, A_tmp, keff, keff );
1512 MVT::MvTransMv( one, *Ctmp, *Utmp, A11 );
1515 SDM A21( Teuchos::View, A_tmp, m+blockSize_, keff, keff );
1516 index.resize(m+blockSize_);
1517 for (i=0; i < m+blockSize_; i++) { index[i] = i; }
1518 Teuchos::RCP<const MV> Vp = MVT::CloneView( *VV, index );
1519 MVT::MvTransMv( one, *Vp, *Utmp, A21 );
1522 for( i=keff; i<keff+m; i++)
1526 SDM A( m2, A_tmp.numCols() );
1527 A.multiply( Teuchos::TRANS, Teuchos::NO_TRANS, one, GG, A_tmp, zero );
1535 char balanc=
'P', jobvl=
'N', jobvr=
'V', sense=
'N';
1536 int ld = A.numRows();
1538 int ldvl = ld, ldvr = ld;
1539 int info = 0,ilo = 0,ihi = 0;
1540 MagnitudeType abnrm = 0.0, bbnrm = 0.0;
1542 std::vector<ScalarType> beta(ld);
1543 std::vector<ScalarType> work(lwork);
1544 std::vector<MagnitudeType> rwork(lwork);
1545 std::vector<MagnitudeType> lscale(ld), rscale(ld);
1546 std::vector<MagnitudeType> rconde(ld), rcondv(ld);
1547 std::vector<int> iwork(ld+6);
1552 lapack.GGEVX(balanc, jobvl, jobvr, sense, ld, A.values(), ld, B.values(), ld, &wr[0], &wi[0],
1553 &beta[0], vl, ldvl, vr.values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0],
1554 &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &rwork[0],
1555 &iwork[0], bwork, &info);
1556 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK GGEVX failed to compute eigensolutions.");
1560 for( i=0; i<ld; i++ )
1561 w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] ) / std::abs( beta[i] );
1563 this->sort(w,ld,iperm);
1565 bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
1568 for( i=0; i<recycledBlocks_; i++ )
1569 for( j=0; j<ld; j++ )
1570 PP(j,i) = vr(j,iperm[ld-recycledBlocks_+i]);
1572 if(scalarTypeIsComplex==
false) {
1575 if (wi[iperm[ld-recycledBlocks_]] != 0.0) {
1577 for ( i=ld-recycledBlocks_; i<ld; i++ )
1578 if (wi[iperm[i]] != 0.0) countImag++;
1580 if (countImag % 2) xtraVec =
true;
1584 if (wi[iperm[ld-recycledBlocks_]] > 0.0) {
1585 for( j=0; j<ld; j++ )
1586 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]+1);
1589 for( j=0; j<ld; j++ )
1590 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]-1);
1598 return recycledBlocks_+1;
1600 return recycledBlocks_;
1604template<
class ScalarType,
class MV,
class OP>
1605int BlockGCRODRSolMgr<ScalarType,MV,OP>::getHarmonicVecsKryl(
int m,
const SDM& HH, SDM& PP){
1606 bool xtraVec =
false;
1607 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1608 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1611 std::vector<MagnitudeType> wr(m), wi(m);
1617 std::vector<MagnitudeType> w(m);
1620 std::vector<int> iperm(m);
1626 builtRecycleSpace_ =
true;
1629 SDM HHt( HH, Teuchos::TRANS );
1630 Teuchos::RCP<SDM> harmRitzMatrix = Teuchos::rcp(
new SDM( m, blockSize_));
1633 for(
int i=0; i<=blockSize_-1; i++) (*harmRitzMatrix)[blockSize_-1-i][harmRitzMatrix->numRows()-1-i] = 1;
1636 lapack.GESV(m, blockSize_, HHt.values(), HHt.stride(), &iperm[0], harmRitzMatrix->values(), harmRitzMatrix->stride(), &info);
1638 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK GESV failed to compute a solution.");
1642 Teuchos::SerialDenseMatrix<int, ScalarType> H_lbl(Teuchos::View, HH, blockSize_, blockSize_, (HH).numRows()-blockSize_, (HH).numCols()-blockSize_ );
1643 Teuchos::SerialDenseMatrix<int, ScalarType> H_lbl_t( H_lbl, Teuchos::TRANS );
1648 Teuchos::RCP<SDM> Htemp = Teuchos::null;
1649 Htemp = Teuchos::rcp(
new SDM(H_lbl_t.numRows(), H_lbl_t.numCols()));
1650 Htemp -> multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, H_lbl_t, H_lbl, zero);
1651 H_lbl_t.assign(*Htemp);
1653 Htemp = Teuchos::rcp(
new SDM(harmRitzMatrix -> numRows(), harmRitzMatrix -> numCols()));
1654 Htemp -> multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, *harmRitzMatrix, H_lbl_t, zero);
1655 harmRitzMatrix -> assign(*Htemp);
1658 int harmColIndex, HHColIndex;
1659 Htemp = Teuchos::rcp(
new SDM(Teuchos::Copy,HH,HH.numRows()-blockSize_,HH.numCols()));
1660 for(
int i = 0; i<blockSize_; i++){
1661 harmColIndex = harmRitzMatrix -> numCols() - i -1;
1663 for(
int j=0; j<m; j++) (*Htemp)[HHColIndex][j] += (*harmRitzMatrix)[harmColIndex][j];
1665 harmRitzMatrix = Htemp;
1673 std::vector<ScalarType> work(1);
1674 std::vector<MagnitudeType> rwork(2*m);
1680 lapack.GEEV(
'N',
'V', m, harmRitzMatrix->values(), harmRitzMatrix->stride(), &wr[0], &wi[0],
1681 vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &rwork[0], &info);
1683 lwork = std::abs (
static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work[0])));
1684 work.resize( lwork );
1686 lapack.GEEV(
'N',
'V', m, harmRitzMatrix->values(), harmRitzMatrix->stride(), &wr[0], &wi[0],
1687 vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &rwork[0], &info);
1689 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK GEEV failed to compute eigensolutions.");
1692 for(
int i=0; i<m; ++i ) w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] );
1694 this->sort(w, m, iperm);
1696 bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
1699 for(
int i=0; i<recycledBlocks_; ++i )
1700 for(
int j=0; j<m; j++ )
1701 PP(j,i) = vr(j,iperm[i]);
1703 if(scalarTypeIsComplex==
false) {
1706 if (wi[iperm[recycledBlocks_-1]] != 0.0) {
1708 for (
int i=0; i<recycledBlocks_; ++i )
1709 if (wi[iperm[i]] != 0.0) countImag++;
1711 if (countImag % 2) xtraVec =
true;
1715 if (wi[iperm[recycledBlocks_-1]] > 0.0) {
1716 for(
int j=0; j<m; ++j )
1717 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]+1);
1720 for(
int j=0; j<m; ++j )
1721 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]-1);
1729 printer_->stream(
Debug) <<
"Recycled " << recycledBlocks_+1 <<
" vectors" << std::endl;
1730 return recycledBlocks_+1;
1733 printer_->stream(
Debug) <<
"Recycled " << recycledBlocks_ <<
" vectors" << std::endl;
1734 return recycledBlocks_;
1739template<
class ScalarType,
class MV,
class OP>
1740void BlockGCRODRSolMgr<ScalarType,MV,OP>::sort(std::vector<MagnitudeType>& dlist,
int n, std::vector<int>& iperm) {
1741 int l, r, j, i, flag;
1743 MagnitudeType dRR, dK;
1768 if (dlist[j] > dlist[j - 1]) j = j + 1;
1769 if (dlist[j - 1] > dK) {
1770 dlist[i - 1] = dlist[j - 1];
1771 iperm[i - 1] = iperm[j - 1];
1784 dlist[r] = dlist[0];
1785 iperm[r] = iperm[0];
1799template<
class ScalarType,
class MV,
class OP>
1803 using Teuchos::rcp_const_cast;
1807 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1808 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1809 std::vector<int> index(numBlocks_+1);
1825 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1826 int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1830 std::vector<int> currIdx;
1832 if ( adaptiveBlockSize_ ) {
1833 blockSize_ = numCurrRHS;
1834 currIdx.resize( numCurrRHS );
1835 for (
int i=0; i<numCurrRHS; ++i)
1836 currIdx[i] = startPtr+i;
1839 currIdx.resize( blockSize_ );
1840 for (
int i=0; i<numCurrRHS; ++i)
1841 currIdx[i] = startPtr+i;
1842 for (
int i=numCurrRHS; i<blockSize_; ++i)
1847 problem_->setLSIndex( currIdx );
1853 loaDetected_ =
false;
1856 bool isConverged =
true;
1859 initializeStateStorage();
1862 Teuchos::ParameterList plist;
1864 while (numRHS2Solve > 0){
1874 printer_->stream(
Debug) <<
" Now solving RHS index " << currIdx[0] <<
" using recycled subspace of dimension " << keff << std::endl << std::endl;
1878 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1879 RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1880 RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1881 problem_->apply( *Utmp, *Ctmp );
1883 RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1887 SDM Ftmp( Teuchos::View, *F_, keff, keff );
1888 int rank = ortho_->normalize(*Ctmp, rcp(&Ftmp,
false));
1890 TEUCHOS_TEST_FOR_EXCEPTION(rank != keff,
BlockGCRODRSolMgrOrthoFailure,
"Belos::BlockGCRODRSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
1895 ipiv_.resize(Ftmp.numRows());
1896 lapack.GETRF(Ftmp.numRows(),Ftmp.numCols(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&info);
1897 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1900 int lwork = Ftmp.numRows();
1901 work_.resize(lwork);
1902 lapack.GETRI(Ftmp.numRows(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1903 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
1906 MVT::MvTimesMatAddMv( one, *Utmp, Ftmp, zero, *U1tmp );
1911 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1912 Ctmp = MVT::CloneViewNonConst( *C_, index );
1913 Utmp = MVT::CloneView( *U_, index );
1916 SDM Ctr(keff,blockSize_);
1917 problem_->computeCurrPrecResVec( &*R_ );
1918 MVT::MvTransMv( one, *Ctmp, *R_, Ctr );
1921 RCP<MV> update = MVT::Clone( *problem_->getCurrLHSVec(), blockSize_ );
1922 MVT::MvInit( *update, 0.0 );
1923 MVT::MvTimesMatAddMv( one, *Utmp, Ctr, one, *update );
1924 problem_->updateSolution( update,
true );
1927 MVT::MvTimesMatAddMv( -one, *Ctmp, Ctr, one, *R_ );
1935 if (V_ == Teuchos::null) {
1937 Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
1938 V_ = MVT::Clone( *rhsMV, (numBlocks_+1)*blockSize_ );
1942 if (MVT::GetNumberVecs(*V_) < (numBlocks_+1)*blockSize_ ) {
1943 Teuchos::RCP<const MV> tmp = V_;
1944 V_ = MVT::Clone( *tmp, (numBlocks_+1)*blockSize_ );
1949 printer_->stream(
Debug) <<
" No recycled subspace available for RHS index " << std::endl << std::endl;
1951 Teuchos::ParameterList primeList;
1954 primeList.set(
"Num Blocks",numBlocks_-1);
1955 primeList.set(
"Block Size",blockSize_);
1956 primeList.set(
"Recycled Blocks",0);
1957 primeList.set(
"Keep Hessenberg",
true);
1958 primeList.set(
"Initialize Hessenberg",
true);
1960 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1961 if (blockSize_*
static_cast<ptrdiff_t
>(numBlocks_) > dim) {
1962 ptrdiff_t tmpNumBlocks = 0;
1963 if (blockSize_ == 1)
1964 tmpNumBlocks = dim / blockSize_;
1966 tmpNumBlocks = ( dim - blockSize_) / blockSize_;
1967 printer_->stream(
Warnings) <<
"Belos::BlockGmresSolMgr::solve(): Warning! Requested Krylov subspace dimension is larger than operator dimension!"
1968 << std::endl <<
"The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << tmpNumBlocks << std::endl;
1969 primeList.set(
"Num Blocks",Teuchos::as<int>(tmpNumBlocks));
1973 primeList.set(
"Num Blocks",numBlocks_-1);
1976 Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > block_gmres_iter;
1980 block_gmres_iter->setSize( blockSize_, numBlocks_-1 );
1983 Teuchos::RCP<MV> V_0;
1984 if (currIdx[blockSize_-1] == -1) {
1985 V_0 = MVT::Clone( *(problem_->getInitPrecResVec()), blockSize_ );
1986 problem_->computeCurrPrecResVec( &*V_0 );
1989 V_0 = MVT::CloneCopy( *(problem_->getInitPrecResVec()), currIdx );
1993 Teuchos::RCP<SDM > z_0 = Teuchos::rcp(
new SDM( blockSize_, blockSize_ ) );
1996 int rank = ortho_->normalize( *V_0, z_0 );
2005 block_gmres_iter->initializeGmres(newstate);
2007 bool primeConverged =
false;
2010 printer_->stream(
Debug) <<
" Preparing to Iterate!!!!" << std::endl << std::endl;
2011 block_gmres_iter->iterate();
2016 if ( convTest_->getStatus() ==
Passed ) {
2017 printer_->stream(
Debug) <<
"We converged during the prime the pump stage" << std::endl << std::endl;
2018 primeConverged = !(expConvTest_->getLOADetected());
2019 if ( expConvTest_->getLOADetected() ) {
2021 loaDetected_ =
true;
2022 printer_->stream(
Warnings) <<
"Belos::BlockGmresSolMgr::solve(): Warning! Solver has experienced a loss of accuracy!" << std::endl;
2028 else if( maxIterTest_->getStatus() ==
Passed ) {
2030 primeConverged =
false;
2036 printer_->stream(
Debug) <<
" We did not converge on priming cycle of Block GMRES. Therefore we recycle and restart. " << std::endl << std::endl;
2041 if (blockSize_ != 1) {
2042 printer_->stream(
Errors) <<
"Error! Caught std::exception in BlockGmresIter::iterate() at iteration "
2043 << block_gmres_iter->getNumIters() << std::endl
2044 << e.what() << std::endl;
2045 if (convTest_->getStatus() !=
Passed)
2046 primeConverged =
false;
2050 block_gmres_iter->updateLSQR( block_gmres_iter->getCurSubspaceDim() );
2052 sTest_->checkStatus( &*block_gmres_iter );
2053 if (convTest_->getStatus() !=
Passed)
2054 isConverged =
false;
2057 catch (
const std::exception &e) {
2058 printer_->stream(
Errors) <<
"Error! Caught std::exception in BlockGmresIter::iterate() at iteration "
2059 << block_gmres_iter->getNumIters() << std::endl
2060 << e.what() << std::endl;
2068 RCP<MV> update = block_gmres_iter->getCurrentUpdate();
2069 problem_->updateSolution( update,
true );
2073 problem_->computeCurrPrecResVec( &*R_ );
2076 newstate = block_gmres_iter->getState();
2079 H_->assign(*(newstate.
H));
2088 V_ = rcp_const_cast<MV>(newstate.
V);
2089 newstate.
V.release();
2091 buildRecycleSpaceKryl(keff, block_gmres_iter);
2092 printer_->stream(
Debug) <<
"Generated recycled subspace using RHS index " << currIdx[0] <<
" of dimension " << keff << std::endl << std::endl;
2096 if (primeConverged) {
2117 problem_->setCurrLS();
2120 startPtr += numCurrRHS;
2121 numRHS2Solve -= numCurrRHS;
2122 if ( numRHS2Solve > 0 ) {
2123 numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
2124 if ( adaptiveBlockSize_ ) {
2125 blockSize_ = numCurrRHS;
2126 currIdx.resize( numCurrRHS );
2127 for (
int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2130 currIdx.resize( blockSize_ );
2131 for (
int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2132 for (
int i=numCurrRHS; i<blockSize_; ++i) currIdx[i] = -1;
2135 problem_->setLSIndex( currIdx );
2138 currIdx.resize( numRHS2Solve );
2148 Teuchos::ParameterList blockgcrodrList;
2149 blockgcrodrList.set(
"Num Blocks",numBlocks_);
2150 blockgcrodrList.set(
"Block Size",blockSize_);
2151 blockgcrodrList.set(
"Recycled Blocks",keff);
2153 Teuchos::RCP<BlockGCRODRIter<ScalarType,MV,OP> > block_gcrodr_iter;
2157 index.resize( blockSize_ );
2158 for(
int ii = 0; ii < blockSize_; ii++) index[ii] = ii;
2159 Teuchos::RCP<MV> V0 = MVT::CloneViewNonConst( *V_, index );
2162 MVT::Assign(*R_,*V0);
2165 for(
int i=0; i < keff; i++){ index[i] = i;};
2166 B_ = rcp(
new SDM(Teuchos::View, *G_, keff, numBlocks_*blockSize_, 0, keff));
2167 H_ = rcp(
new SDM(Teuchos::View, *G_, (numBlocks_-1)*blockSize_ + blockSize_, (numBlocks_-1)*blockSize_, keff ,keff ));
2171 newstate.
U = MVT::CloneViewNonConst(*U_, index);
2172 newstate.
C = MVT::CloneViewNonConst(*C_, index);
2174 newstate.
curDim = blockSize_;
2175 block_gcrodr_iter -> initialize(newstate);
2177 int numRestarts = 0;
2181 block_gcrodr_iter -> iterate();
2186 if( convTest_->getStatus() ==
Passed ) {
2194 else if(maxIterTest_->getStatus() ==
Passed ){
2196 isConverged =
false;
2203 else if (block_gcrodr_iter->getCurSubspaceDim() == block_gcrodr_iter->getMaxSubspaceDim()){
2208 Teuchos::RCP<MV> update = block_gcrodr_iter->getCurrentUpdate();
2209 problem_->updateSolution(update,
true);
2210 buildRecycleSpaceAugKryl(block_gcrodr_iter);
2212 printer_->stream(
Debug) <<
" Generated new recycled subspace using RHS index " << currIdx[0] <<
" of dimension " << keff << std::endl << std::endl;
2214 if(numRestarts >= maxRestarts_) {
2215 isConverged =
false;
2221 printer_ -> stream(
Debug) <<
" Performing restart number " << numRestarts <<
" of " << maxRestarts_ << std::endl << std::endl;
2224 problem_->computeCurrPrecResVec( &*R_ );
2225 index.resize( blockSize_ );
2226 for (
int ii=0; ii<blockSize_; ++ii) index[ii] = ii;
2227 Teuchos::RCP<MV> V0 = MVT::CloneViewNonConst( *V_, index );
2228 MVT::SetBlock(*R_,index,*V0);
2232 index.resize( numBlocks_*blockSize_ );
2233 for (
int ii=0; ii<(numBlocks_*blockSize_); ++ii) index[ii] = ii;
2234 restartState.
V = MVT::CloneViewNonConst( *V_, index );
2235 index.resize( keff );
2236 for (
int ii=0; ii<keff; ++ii) index[ii] = ii;
2237 restartState.
U = MVT::CloneViewNonConst( *U_, index );
2238 restartState.
C = MVT::CloneViewNonConst( *C_, index );
2239 B_ = rcp(
new SDM(Teuchos::View, *G_, keff, (numBlocks_-1)*blockSize_, 0, keff));
2240 H_ = rcp(
new SDM(Teuchos::View, *G_, numBlocks_*blockSize_, (numBlocks_-1)*blockSize_, keff ,keff ));
2241 restartState.
B = B_;
2242 restartState.
H = H_;
2243 restartState.
curDim = blockSize_;
2244 block_gcrodr_iter->initialize(restartState);
2253 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Belos::BlockGCRODRSolMgr::solve(): Invalid return from BlockGCRODRIter::iterate().");
2259 block_gcrodr_iter->updateLSQR( block_gcrodr_iter->getCurSubspaceDim() );
2261 sTest_->checkStatus( &*block_gcrodr_iter );
2262 if (convTest_->getStatus() !=
Passed) isConverged =
false;
2265 catch(
const std::exception &e){
2266 printer_->stream(
Errors) <<
"Error! Caught exception in BlockGCRODRIter::iterate() at iteration "
2267 << block_gcrodr_iter->getNumIters() << std::endl
2268 << e.what() << std::endl;
2275 Teuchos::RCP<MV> update = block_gcrodr_iter->getCurrentUpdate();
2276 problem_->updateSolution( update,
true );
2295 problem_->setCurrLS();
2298 startPtr += numCurrRHS;
2299 numRHS2Solve -= numCurrRHS;
2300 if ( numRHS2Solve > 0 ) {
2301 numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
2302 if ( adaptiveBlockSize_ ) {
2303 blockSize_ = numCurrRHS;
2304 currIdx.resize( numCurrRHS );
2305 for (
int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2308 currIdx.resize( blockSize_ );
2309 for (
int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2310 for (
int i=numCurrRHS; i<blockSize_; ++i) currIdx[i] = -1;
2313 problem_->setLSIndex( currIdx );
2316 currIdx.resize( numRHS2Solve );
2320 if (!builtRecycleSpace_) {
2321 buildRecycleSpaceAugKryl(block_gcrodr_iter);
2322 printer_->stream(
Debug) <<
" Generated new recycled subspace using RHS index " << currIdx[0] <<
" of dimension " << keff << std::endl << std::endl;
2330 #ifdef BELOS_TEUCHOS_TIME_MONITOR
2336 numIters_ = maxIterTest_->getNumIters();
2339 const std::vector<MagnitudeType>* pTestValues = NULL;
2340 pTestValues = impConvTest_->getTestValue();
2341 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
2342 "Belos::BlockGCRODRSolMgr::solve(): The implicit convergence test's "
2343 "getTestValue() method returned NULL. Please report this bug to the "
2344 "Belos developers.");
2345 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
2346 "Belos::BlockGCRODRSolMgr::solve(): The implicit convergence test's "
2347 "getTestValue() method returned a vector of length zero. Please report "
2348 "this bug to the Belos developers.");
2352 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
Belos concrete class for performing the block, flexible GMRES iteration.
Belos concrete class for performing the block GCRO-DR (block GMRES with recycling) iteration.
Belos concrete class for performing the block GMRES iteration.
Belos header file which uses auto-configuration information to include necessary C++ headers.
Pure virtual base class which augments the basic interface for a Gmres linear solver iteration.
Class which describes the linear problem to be solved by the iterative solver.
Class which manages the output and verbosity of the Belos solvers.
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.
BlockGCRODRIterOrthoFailure is thrown when the BlockGCRODRIter object is unable to compute independen...
Implementation of the Block GCRO-DR (Block Recycling GMRES) iteration.
Thrown when an LAPACK call fails.
BlockGCRODRSolMgrLAPACKFailure(const std::string &what_arg)
Thrown when the linear problem was not set up correctly.
BlockGCRODRSolMgrLinearProblemFailure(const std::string &what_arg)
Thrown when the solution manager was unable to orthogonalize the basis vectors.
BlockGCRODRSolMgrOrthoFailure(const std::string &what_arg)
Thrown if any problem occurs in using or creating the recycle subspace.
BlockGCRODRSolMgrRecyclingFailure(const std::string &what_arg)
A solver manager for the Block GCRO-DR (Block Recycling GMRES) linear solver.
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > ¶ms)
Set the parameters the solver should use to solve the linear problem.
BlockGCRODRSolMgr()
Default constructor.
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const
Get a parameter list containing the current parameters for this object.
std::string description() const
A description of the Block GCRODR solver manager.
MagnitudeType achievedTol() const
Get the residual for the most recent call to solve().
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Get a parameter list containing the valid parameters for this object.
void reset(const ResetType type)
Performs a reset of the solver manager specified by the ResetType.
virtual ~BlockGCRODRSolMgr()
Destructor.
bool isLOADetected() const
Whether a loss of accuracy was detected during the most recent solve.
ReturnType solve()
Solve the current linear problem.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem)
Set the linear problem to solve on the next call to solve().
int getNumIters() const
Get the iteration count for the most recent call to solve().
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get current linear problem being solved for in this object.
This class implements the block GMRES iteration, where a block Krylov subspace is constructed....
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
void setDepTol(const MagnitudeType dep_tol)
Set parameter for re-orthogonalization threshhold.
GmresIterationOrthoFailure is thrown when the GmresIteration object is unable to compute independent ...
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.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
The Belos::SolverManager is a templated virtual base class that defines the basic interface that any ...
A class for extending the status testing capabilities of Belos via logical combinations.
Exception thrown to signal error in a status test during Belos::StatusTest::checkStatus().
An implementation of StatusTestResNorm using a family of residual norms.
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.
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
ReturnType
Whether the Belos solve converged for all linear systems.
OutputType
Style of output used to display status test information.
ScaleType
The type of scaling to use on the residual norm value.
ResetType
How to reset the solver.
Structure to contain pointers to BlockGCRODRIter state variables.
Teuchos::RCP< MV > V
The current Krylov basis.
Teuchos::RCP< MV > U
The recycled subspace and its projection.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B
The projection of the Krylov subspace against the recycled subspace *
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
int curDim
The current dimension of the reduction.
Structure to contain pointers to GmresIteration state variables.
Teuchos::RCP< const MV > V
The current Krylov basis.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > z
The current right-hand side of the least squares system RY = Z.
int curDim
The current dimension of the reduction.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.