45 #ifndef BELOS_BLOCK_GCRODR_SOLMGR_HPP 46 #define BELOS_BLOCK_GCRODR_SOLMGR_HPP 62 #include "Teuchos_BLAS.hpp" 63 #include "Teuchos_LAPACK.hpp" 64 #include "Teuchos_as.hpp" 65 #ifdef BELOS_TEUCHOS_TIME_MONITOR 66 #include "Teuchos_TimeMonitor.hpp" 67 #endif // BELOS_TEUCHOS_TIME_MONITOR 126 template<
class ScalarType,
class MV,
class OP>
132 typedef Teuchos::ScalarTraits<ScalarType> SCT;
133 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
134 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
135 typedef Teuchos::ScalarTraits<MagnitudeType> SMT;
137 typedef Teuchos::SerialDenseMatrix<int,ScalarType> SDM;
138 typedef Teuchos::SerialDenseVector<int,ScalarType> SDV;
177 const Teuchos::RCP<Teuchos::ParameterList> &pl);
224 TEUCHOS_TEST_FOR_EXCEPTION(problem.is_null(), std::invalid_argument,
225 "Belos::BlockGCRODRSolMgr::setProblem: The input LinearProblem cannot be null.");
228 if (! problem->isProblemSet()) {
229 const bool success = problem->setProblem();
230 TEUCHOS_TEST_FOR_EXCEPTION(success, std::runtime_error,
231 "Belos::BlockGCRODRSolMgr::setProblem: Calling the input LinearProblem's setProblem() method failed. This likely means that the " 232 "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.");
239 void setParameters(
const Teuchos::RCP<Teuchos::ParameterList> ¶ms );
254 problem_->setProblem();
292 void initializeStateStorage();
310 int getHarmonicVecsKryl (
int m,
const SDM& HH, SDM& PP);
315 int getHarmonicVecsAugKryl (
int keff,
318 const Teuchos::RCP<const MV>& VV,
322 void sort (std::vector<MagnitudeType>& dlist,
int n, std::vector<int>& iperm);
325 Teuchos::LAPACK<int,ScalarType> lapack;
328 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
331 Teuchos::RCP<OutputManager<ScalarType> > printer_;
332 Teuchos::RCP<std::ostream> outputStream_;
335 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
336 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
337 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
338 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > expConvTest_, impConvTest_;
339 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
342 ortho_factory_type orthoFactory_;
349 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
352 Teuchos::RCP<Teuchos::ParameterList> params_;
364 mutable Teuchos::RCP<const Teuchos::ParameterList> defaultParams_;
367 static const bool adaptiveBlockSize_default_;
368 static const std::string recycleMethod_default_;
371 MagnitudeType convTol_, orthoKappa_, achievedTol_;
372 int blockSize_, maxRestarts_, maxIters_, numIters_;
373 int verbosity_, outputStyle_, outputFreq_;
374 bool adaptiveBlockSize_;
375 std::string orthoType_, recycleMethod_;
376 std::string impResScale_, expResScale_;
384 int numBlocks_, recycledBlocks_;
395 Teuchos::RCP<MV> U_, C_;
398 Teuchos::RCP<MV> U1_, C1_;
401 Teuchos::RCP<SDM > G_;
402 Teuchos::RCP<SDM > H_;
403 Teuchos::RCP<SDM > B_;
404 Teuchos::RCP<SDM > PP_;
405 Teuchos::RCP<SDM > HP_;
406 std::vector<ScalarType> tau_;
407 std::vector<ScalarType> work_;
408 Teuchos::RCP<SDM > F_;
409 std::vector<int> ipiv_;
412 Teuchos::RCP<Teuchos::Time> timerSolve_;
421 bool builtRecycleSpace_;
427 template<
class ScalarType,
class MV,
class OP>
428 const bool BlockGCRODRSolMgr<ScalarType,MV,OP>::adaptiveBlockSize_default_ =
true;
430 template<
class ScalarType,
class MV,
class OP>
431 const std::string BlockGCRODRSolMgr<ScalarType,MV,OP>::recycleMethod_default_ =
"harmvecs";
437 template<
class ScalarType,
class MV,
class OP>
443 template<
class ScalarType,
class MV,
class OP>
446 const Teuchos::RCP<Teuchos::ParameterList> &pl ) {
451 TEUCHOS_TEST_FOR_EXCEPTION(problem.is_null(), std::invalid_argument,
452 "Belos::BlockGCRODR constructor: The solver manager's constructor needs " 453 "the linear problem argument 'problem' to be nonnull.");
463 template<
class ScalarType,
class MV,
class OP>
465 adaptiveBlockSize_ = adaptiveBlockSize_default_;
466 recycleMethod_ = recycleMethod_default_;
468 loaDetected_ =
false;
469 builtRecycleSpace_ =
false;
540 template<
class ScalarType,
class MV,
class OP>
542 std::ostringstream oss;
543 oss <<
"Belos::BlockGCRODRSolMgr<" << SCT::name() <<
", ...>";
545 oss <<
"Ortho Type='"<<orthoType_ ;
546 oss <<
", Num Blocks=" <<numBlocks_;
547 oss <<
", Block Size=" <<blockSize_;
548 oss <<
", Num Recycle Blocks=" << recycledBlocks_;
549 oss <<
", Max Restarts=" << maxRestarts_;
554 template<
class ScalarType,
class MV,
class OP>
555 Teuchos::RCP<const Teuchos::ParameterList>
557 using Teuchos::ParameterList;
558 using Teuchos::parameterList;
560 using Teuchos::rcpFromRef;
562 if (defaultParams_.is_null()) {
563 RCP<ParameterList> pl = parameterList ();
565 const MagnitudeType convTol = SMT::squareroot (SCT::magnitude (SCT::eps()));
566 const int maxRestarts = 1000;
567 const int maxIters = 5000;
568 const int blockSize = 2;
569 const int numBlocks = 100;
570 const int numRecycledBlocks = 25;
573 const int outputFreq = 1;
574 RCP<std::ostream> outputStream = rcpFromRef (std::cout);
575 const std::string impResScale (
"Norm of Preconditioned Initial Residual");
576 const std::string expResScale (
"Norm of Initial Residual");
577 const std::string timerLabel (
"Belos");
578 const std::string orthoType (
"DGKS");
579 RCP<const ParameterList> orthoParams = orthoFactory_.getDefaultParameters (orthoType);
586 const MagnitudeType orthoKappa = -SMT::one();
589 pl->set (
"Convergence Tolerance", convTol,
590 "The tolerance that the solver needs to achieve in order for " 591 "the linear system(s) to be declared converged. The meaning " 592 "of this tolerance depends on the convergence test details.");
593 pl->set(
"Maximum Restarts", maxRestarts,
594 "The maximum number of restart cycles (not counting the first) " 595 "allowed for each set of right-hand sides solved.");
596 pl->set(
"Maximum Iterations", maxIters,
597 "The maximum number of iterations allowed for each set of " 598 "right-hand sides solved.");
599 pl->set(
"Block Size", blockSize,
600 "How many linear systems to solve at once.");
601 pl->set(
"Num Blocks", numBlocks,
602 "The maximum number of blocks allowed in the Krylov subspace " 603 "for each set of right-hand sides solved. This includes the " 604 "initial block. Total Krylov basis storage, not counting the " 605 "recycled subspace, is \"Num Blocks\" times \"Block Size\".");
606 pl->set(
"Num Recycled Blocks", numRecycledBlocks,
607 "The maximum number of vectors in the recycled subspace." );
608 pl->set(
"Verbosity", verbosity,
609 "What type(s) of solver information should be written " 610 "to the output stream.");
611 pl->set(
"Output Style", outputStyle,
612 "What style is used for the solver information to write " 613 "to the output stream.");
614 pl->set(
"Output Frequency", outputFreq,
615 "How often convergence information should be written " 616 "to the output stream.");
617 pl->set(
"Output Stream", outputStream,
618 "A reference-counted pointer to the output stream where all " 619 "solver output is sent.");
620 pl->set(
"Implicit Residual Scaling", impResScale,
621 "The type of scaling used in the implicit residual convergence test.");
622 pl->set(
"Explicit Residual Scaling", expResScale,
623 "The type of scaling used in the explicit residual convergence test.");
624 pl->set(
"Timer Label", timerLabel,
625 "The string to use as a prefix for the timer labels.");
627 pl->set(
"Orthogonalization", orthoType,
628 "The orthogonalization method to use. Valid options: " +
629 orthoFactory_.validNamesString());
630 pl->set (
"Orthogonalization Parameters", *orthoParams,
631 "Sublist of parameters specific to the orthogonalization method to use.");
632 pl->set(
"Orthogonalization Constant", orthoKappa,
633 "When using DGKS orthogonalization: the \"depTol\" constant, used " 634 "to determine whether another step of classical Gram-Schmidt is " 635 "necessary. Otherwise ignored. Nonpositive values are ignored.");
638 return defaultParams_;
641 template<
class ScalarType,
class MV,
class OP>
645 using Teuchos::isParameterType;
646 using Teuchos::getParameter;
648 using Teuchos::ParameterList;
649 using Teuchos::parameterList;
652 using Teuchos::rcp_dynamic_cast;
653 using Teuchos::rcpFromRef;
654 using Teuchos::Exceptions::InvalidParameter;
655 using Teuchos::Exceptions::InvalidParameterName;
656 using Teuchos::Exceptions::InvalidParameterType;
659 RCP<const ParameterList> defaultParams = getValidParameters();
661 if (params.is_null()) {
666 params_ = parameterList (*defaultParams);
678 params->validateParametersAndSetDefaults (*defaultParams);
694 const int maxRestarts = params_->get<
int> (
"Maximum Restarts");
695 TEUCHOS_TEST_FOR_EXCEPTION(maxRestarts <= 0, std::invalid_argument,
696 "Belos::BlockGCRODRSolMgr: The \"Maximum Restarts\" parameter " 697 "must be nonnegative, but you specified a negative value of " 698 << maxRestarts <<
".");
700 const int maxIters = params_->get<
int> (
"Maximum Iterations");
701 TEUCHOS_TEST_FOR_EXCEPTION(maxIters <= 0, std::invalid_argument,
702 "Belos::BlockGCRODRSolMgr: The \"Maximum Iterations\" parameter " 703 "must be positive, but you specified a nonpositive value of " 706 const int numBlocks = params_->get<
int> (
"Num Blocks");
707 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument,
708 "Belos::BlockGCRODRSolMgr: The \"Num Blocks\" parameter must be " 709 "positive, but you specified a nonpositive value of " << numBlocks
712 const int blockSize = params_->get<
int> (
"Block Size");
713 TEUCHOS_TEST_FOR_EXCEPTION(blockSize <= 0, std::invalid_argument,
714 "Belos::BlockGCRODRSolMgr: The \"Block Size\" parameter must be " 715 "positive, but you specified a nonpositive value of " << blockSize
718 const int recycledBlocks = params_->get<
int> (
"Num Recycled Blocks");
719 TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks <= 0, std::invalid_argument,
720 "Belos::BlockGCRODRSolMgr: The \"Num Recycled Blocks\" parameter must " 721 "be positive, but you specified a nonpositive value of " 722 << recycledBlocks <<
".");
723 TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks >= numBlocks,
724 std::invalid_argument,
"Belos::BlockGCRODRSolMgr: The \"Num Recycled " 725 "Blocks\" parameter must be less than the \"Num Blocks\" parameter, " 726 "but you specified \"Num Recycled Blocks\" = " << recycledBlocks
727 <<
" and \"Num Blocks\" = " << numBlocks <<
".");
731 maxRestarts_ = maxRestarts;
732 maxIters_ = maxIters;
733 numBlocks_ = numBlocks;
734 blockSize_ = blockSize;
735 recycledBlocks_ = recycledBlocks;
742 std::string tempLabel = params_->get<std::string> (
"Timer Label");
743 const bool labelChanged = (tempLabel != label_);
745 #ifdef BELOS_TEUCHOS_TIME_MONITOR 746 std::string solveLabel = label_ +
": BlockGCRODRSolMgr total solve time";
747 if (timerSolve_.is_null()) {
749 timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
750 }
else if (labelChanged) {
756 Teuchos::TimeMonitor::clearCounter (solveLabel);
757 timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
759 #endif // BELOS_TEUCHOS_TIME_MONITOR 765 if (params_->isParameter (
"Verbosity")) {
766 if (isParameterType<int> (*params_,
"Verbosity")) {
767 verbosity_ = params_->get<
int> (
"Verbosity");
770 verbosity_ = (int) getParameter<MsgType> (*params_,
"Verbosity");
777 if (params_->isParameter (
"Output Style")) {
778 if (isParameterType<int> (*params_,
"Output Style")) {
779 outputStyle_ = params_->get<
int> (
"Output Style");
782 outputStyle_ = (int) getParameter<OutputType> (*params_,
"Output Style");
810 if (params_->isParameter (
"Output Stream")) {
812 outputStream_ = getParameter<RCP<std::ostream> > (*params_,
"Output Stream");
814 catch (InvalidParameter&) {
815 outputStream_ = rcpFromRef (std::cout);
822 if (outputStream_.is_null()) {
823 outputStream_ = rcp (
new Teuchos::oblackholestream);
827 outputFreq_ = params_->get<
int> (
"Output Frequency");
830 if (! outputTest_.is_null()) {
831 outputTest_->setOutputFrequency (outputFreq_);
839 if (printer_.is_null()) {
842 printer_->setVerbosity (verbosity_);
843 printer_->setOStream (outputStream_);
854 if (params_->isParameter (
"Orthogonalization")) {
855 const std::string& tempOrthoType =
856 params_->get<std::string> (
"Orthogonalization");
858 if (! orthoFactory_.isValidName (tempOrthoType)) {
859 std::ostringstream os;
860 os <<
"Belos::BlockGCRODRSolMgr: Invalid orthogonalization name \"" 861 << tempOrthoType <<
"\". The following are valid options " 862 <<
"for the \"Orthogonalization\" name parameter: ";
863 orthoFactory_.printValidNames (os);
864 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument, os.str());
866 if (tempOrthoType != orthoType_) {
868 orthoType_ = tempOrthoType;
885 RCP<ParameterList> orthoParams = sublist (params,
"Orthogonalization Parameters",
true);
886 TEUCHOS_TEST_FOR_EXCEPTION(orthoParams.is_null(), std::logic_error,
887 "Failed to get orthogonalization parameters. " 888 "Please report this bug to the Belos developers.");
914 ortho_ = orthoFactory_.makeMatOrthoManager (orthoType_, null, printer_,
915 label_, orthoParams);
927 if (params_->isParameter (
"Orthogonalization Constant")) {
928 const MagnitudeType orthoKappa =
929 params_->get<MagnitudeType> (
"Orthogonalization Constant");
930 if (orthoKappa > 0) {
931 orthoKappa_ = orthoKappa;
934 if (orthoType_ ==
"DGKS" && ! ortho_.is_null()) {
940 rcp_dynamic_cast<ortho_man_type>(ortho_)->setDepTol (orthoKappa_);
950 convTol_ = params_->get<MagnitudeType> (
"Convergence Tolerance");
951 if (! impConvTest_.is_null()) {
954 if (! expConvTest_.is_null()) {
955 expConvTest_->setTolerance (convTol_);
959 if (params_->isParameter (
"Implicit Residual Scaling")) {
960 std::string tempImpResScale =
961 getParameter<std::string> (*params_,
"Implicit Residual Scaling");
964 if (impResScale_ != tempImpResScale) {
966 impResScale_ = tempImpResScale;
968 if (! impConvTest_.is_null()) {
981 if (params_->isParameter(
"Explicit Residual Scaling")) {
982 std::string tempExpResScale =
983 getParameter<std::string> (*params_,
"Explicit Residual Scaling");
986 if (expResScale_ != tempExpResScale) {
988 expResScale_ = tempExpResScale;
990 if (! expConvTest_.is_null()) {
1008 if (maxIterTest_.is_null()) {
1011 maxIterTest_->setMaxIters (maxIters_);
1016 if (impConvTest_.is_null()) {
1017 impConvTest_ = rcp (
new StatusTestResNorm_t (convTol_));
1023 if (expConvTest_.is_null()) {
1024 expConvTest_ = rcp (
new StatusTestResNorm_t (convTol_));
1025 expConvTest_->defineResForm (StatusTestResNorm_t::Explicit,
Belos::TwoNorm);
1031 if (convTest_.is_null()) {
1032 convTest_ = rcp (
new StatusTestCombo_t (StatusTestCombo_t::SEQ,
1040 sTest_ = rcp (
new StatusTestCombo_t (StatusTestCombo_t::OR,
1041 maxIterTest_, convTest_));
1045 outputTest_ = stoFactory.
create (printer_, sTest_, outputFreq_,
1049 std::string solverDesc =
"Block GCRODR ";
1050 outputTest_->setSolverDesc (solverDesc);
1057 template<
class ScalarType,
class MV,
class OP>
1062 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1065 Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
1068 int KrylSpaDim = (numBlocks_ - 1) * blockSize_;
1071 int augSpaDim = KrylSpaDim + recycledBlocks_ + 1;
1074 int augSpaImgDim = KrylSpaDim + blockSize_ + recycledBlocks_+1;
1098 if (U_ == Teuchos::null) {
1099 U_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1103 if (MVT::GetNumberVecs(*U_) < recycledBlocks_+1) {
1104 Teuchos::RCP<const MV> tmp = U_;
1105 U_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1110 if (C_ == Teuchos::null) {
1111 C_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1115 if (MVT::GetNumberVecs(*C_) < recycledBlocks_+1) {
1116 Teuchos::RCP<const MV> tmp = C_;
1117 C_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1122 if (U1_ == Teuchos::null) {
1123 U1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1127 if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1128 Teuchos::RCP<const MV> tmp = U1_;
1129 U1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1134 if (C1_ == Teuchos::null) {
1135 C1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1139 if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1140 Teuchos::RCP<const MV> tmp = C1_;
1141 C1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1146 if (R_ == Teuchos::null){
1147 R_ = MVT::Clone( *rhsMV, blockSize_ );
1153 if (G_ == Teuchos::null){
1154 G_ = Teuchos::rcp(
new SDM( augSpaImgDim, augSpaDim ) );
1157 if ( (G_->numRows() != augSpaImgDim) || (G_->numCols() != augSpaDim) )
1159 G_->reshape( augSpaImgDim, augSpaDim );
1161 G_->putScalar(zero);
1165 if (H_ == Teuchos::null){
1166 H_ = Teuchos::rcp (
new SDM ( Teuchos::View, *G_, KrylSpaDim + blockSize_, KrylSpaDim, recycledBlocks_+1 ,recycledBlocks_+1 ) );
1170 if (F_ == Teuchos::null){
1171 F_ = Teuchos::rcp(
new SDM( recycledBlocks_+1, recycledBlocks_+1 ) );
1174 if ( (F_->numRows() != recycledBlocks_+1) || (F_->numCols() != recycledBlocks_+1) ){
1175 F_->reshape( recycledBlocks_+1, recycledBlocks_+1 );
1178 F_->putScalar(zero);
1181 if (PP_ == Teuchos::null){
1182 PP_ = Teuchos::rcp(
new SDM( augSpaImgDim, recycledBlocks_+1 ) );
1185 if ( (PP_->numRows() != augSpaImgDim) || (PP_->numCols() != recycledBlocks_+1) ){
1186 PP_->reshape( augSpaImgDim, recycledBlocks_+1 );
1191 if (HP_ == Teuchos::null)
1192 HP_ = Teuchos::rcp(
new SDM( augSpaImgDim, augSpaDim ) );
1194 if ( (HP_->numRows() != augSpaImgDim) || (HP_->numCols() != augSpaDim) ){
1195 HP_->reshape( augSpaImgDim, augSpaDim );
1200 tau_.resize(recycledBlocks_+1);
1203 work_.resize(recycledBlocks_+1);
1206 ipiv_.resize(recycledBlocks_+1);
1210 template<
class ScalarType,
class MV,
class OP>
1211 void BlockGCRODRSolMgr<ScalarType,MV,OP>::buildRecycleSpaceKryl(
int& keff, Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > block_gmres_iter){
1213 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1214 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1216 int p = block_gmres_iter->getState().curDim;
1217 std::vector<int> index(keff);
1221 SDM HH(Teuchos::Copy, *H_, p+blockSize_, p);
1222 if(recycledBlocks_ >= p + blockSize_){
1226 for (
int ii=0; ii < p; ++ii) index[ii] = ii;
1227 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1228 MVT::SetBlock(*V_, index, *Utmp);
1233 Teuchos::RCP<SDM > PPtmp = Teuchos::rcp (
new SDM ( Teuchos::View, *PP_, p, recycledBlocks_+1 ) );
1234 if(recycleMethod_ ==
"harmvecs"){
1235 keff = getHarmonicVecsKryl(p, HH, *PPtmp);
1236 printer_->stream(
Debug) <<
"keff = " << keff << std::endl;
1239 PPtmp = Teuchos::rcp (
new SDM ( Teuchos::View, *PP_, p, keff ) );
1242 for (
int ii=0; ii<keff; ++ii) index[ii] = ii;
1243 Teuchos::RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1244 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1245 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1247 for (
int ii=0; ii < p; ++ii) index[ii] = ii;
1248 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1252 MVT::MvTimesMatAddMv( one, *Vtmp, *PPtmp, zero, *U1tmp );
1255 SDM HPtmp( Teuchos::View, *HP_, p+blockSize_, keff );
1256 HPtmp.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, *H_, *PPtmp, zero );
1260 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1261 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
1264 lwork =
static_cast<int> (Teuchos::ScalarTraits<ScalarType>::magnitude (work_[0]));
1265 work_.resize(lwork);
1266 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1267 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
1271 SDM Rtmp( Teuchos::View, *F_, keff, keff );
1272 for(
int ii=0;ii<keff;ii++) {
for(
int jj=ii;jj<keff;jj++) Rtmp(ii,jj) = HPtmp(ii,jj); }
1274 lapack.UNGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1275 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _UNGQR failed to construct the Q factor.");
1279 index.resize( p+blockSize_ );
1280 for (
int ii=0; ii < (p+blockSize_); ++ii) { index[ii] = ii; }
1281 Vtmp = MVT::CloneView( *V_, index );
1282 MVT::MvTimesMatAddMv( one, *Vtmp, HPtmp, zero, *Ctmp );
1289 ipiv_.resize(Rtmp.numRows());
1290 lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
1291 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1293 lwork = Rtmp.numRows();
1294 work_.resize(lwork);
1295 lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1298 MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
1304 template<
class ScalarType,
class MV,
class OP>
1305 void BlockGCRODRSolMgr<ScalarType,MV,OP>::buildRecycleSpaceAugKryl(Teuchos::RCP<BlockGCRODRIter<ScalarType,MV,OP> > block_gcrodr_iter){
1306 const MagnitudeType one = Teuchos::ScalarTraits<ScalarType>::one();
1307 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1309 std::vector<MagnitudeType> d(keff);
1310 std::vector<ScalarType> dscalar(keff);
1311 std::vector<int> index(numBlocks_+1);
1314 BlockGCRODRIterState<ScalarType,MV> oldState = block_gcrodr_iter->getState();
1315 int p = oldState.curDim;
1320 if(recycledBlocks_ >= keff + p){
1323 for (
int ii=0; ii < p; ++ii) { index[ii] = keff+ii; }
1324 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1325 for (
int ii=0; ii < p; ++ii) { index[ii] =ii; }
1326 MVT::SetBlock(*V_, index, *Utmp);
1333 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1334 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1336 dscalar.resize(keff);
1337 MVT::MvNorm( *Utmp, d );
1338 for (
int i=0; i<keff; ++i) {
1340 dscalar[i] = (ScalarType)d[i];
1342 MVT::MvScale( *Utmp, dscalar );
1347 Teuchos::RCP<SDM> Gtmp = Teuchos::rcp(
new SDM( Teuchos::View, *G_, p+keff, p+keff-blockSize_ ) );
1350 for (
int i=0; i<keff; ++i)
1351 (*Gtmp)(i,i) = d[i];
1358 SDM PPtmp( Teuchos::View, *PP_, p+keff-blockSize_, recycledBlocks_+1 );
1359 keff_new = getHarmonicVecsAugKryl( keff, p-blockSize_, *Gtmp, oldState.V, PPtmp );
1366 Teuchos::RCP<MV> U1tmp;
1368 index.resize( keff );
1369 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1370 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1371 index.resize( keff_new );
1372 for (
int ii=0; ii<keff_new; ++ii) { index[ii] = ii; }
1373 U1tmp = MVT::CloneViewNonConst( *U1_, index );
1374 SDM PPtmp( Teuchos::View, *PP_, keff, keff_new );
1375 MVT::MvTimesMatAddMv( one, *Utmp, PPtmp, zero, *U1tmp );
1380 index.resize(p-blockSize_);
1381 for (
int ii=0; ii < p-blockSize_; ii++) { index[ii] = ii; }
1382 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1383 SDM PPtmp( Teuchos::View, *PP_, p-blockSize_, keff_new, keff );
1384 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *U1tmp );
1388 SDM HPtmp( Teuchos::View, *HP_, p+keff, keff_new );
1390 SDM PPtmp( Teuchos::View, *PP_, p-blockSize_+keff, keff_new );
1391 HPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*Gtmp,PPtmp,zero);
1395 int info = 0, lwork = -1;
1396 tau_.resize(keff_new);
1397 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1398 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
1400 lwork =
static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work_[0]));
1401 work_.resize(lwork);
1402 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1403 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
1407 SDM Ftmp( Teuchos::View, *F_, keff_new, keff_new );
1408 for(
int i=0;i<keff_new;i++) {
for(
int j=i;j<keff_new;j++) Ftmp(i,j) = HPtmp(i,j); }
1410 lapack.UNGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1411 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _UNGQR failed to construct the Q factor.");
1418 Teuchos::RCP<MV> C1tmp;
1421 for (
int i=0; i < keff; i++) { index[i] = i; }
1422 Teuchos::RCP<const MV> Ctmp = MVT::CloneView( *C_, index );
1423 index.resize(keff_new);
1424 for (
int i=0; i < keff_new; i++) { index[i] = i; }
1425 C1tmp = MVT::CloneViewNonConst( *C1_, index );
1426 SDM PPtmp( Teuchos::View, *HP_, keff, keff_new );
1427 MVT::MvTimesMatAddMv( one, *Ctmp, PPtmp, zero, *C1tmp );
1432 for (
int i=0; i < p; ++i) { index[i] = i; }
1433 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1434 SDM PPtmp( Teuchos::View, *HP_, p, keff_new, keff, 0 );
1435 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *C1tmp );
1444 ipiv_.resize(Ftmp.numRows());
1445 lapack.GETRF(Ftmp.numRows(),Ftmp.numCols(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&info);
1446 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1449 lwork = Ftmp.numRows();
1450 work_.resize(lwork);
1451 lapack.GETRI(Ftmp.numRows(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1452 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRI failed to compute an LU factorization.");
1455 index.resize(keff_new);
1456 for (
int i=0; i < keff_new; i++) { index[i] = i; }
1457 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1458 MVT::MvTimesMatAddMv( one, *U1tmp, Ftmp, zero, *Utmp );
1463 if (keff != keff_new) {
1465 block_gcrodr_iter->setSize( keff, numBlocks_ );
1467 SDM b1( Teuchos::View, *G_, recycledBlocks_+2, 1, 0, recycledBlocks_ );
1473 template<
class ScalarType,
class MV,
class OP>
1474 int BlockGCRODRSolMgr<ScalarType,MV,OP>::getHarmonicVecsAugKryl(
int keff,
int m,
const SDM& GG,
const Teuchos::RCP<const MV>& VV, SDM& PP){
1476 int m2 = GG.numCols();
1477 bool xtraVec =
false;
1478 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1479 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1480 std::vector<int> index;
1483 std::vector<MagnitudeType> wr(m2), wi(m2);
1486 std::vector<MagnitudeType> w(m2);
1489 SDM vr(m2,m2,
false);
1492 std::vector<int> iperm(m2);
1495 builtRecycleSpace_ =
true;
1505 SDM A_tmp( keff+m+blockSize_, keff+m );
1510 for (
int i=0; i<keff; ++i) { index[i] = i; }
1511 Teuchos::RCP<const MV> Ctmp = MVT::CloneView( *C_, index );
1512 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1513 SDM A11( Teuchos::View, A_tmp, keff, keff );
1514 MVT::MvTransMv( one, *Ctmp, *Utmp, A11 );
1517 SDM A21( Teuchos::View, A_tmp, m+blockSize_, keff, keff );
1518 index.resize(m+blockSize_);
1519 for (i=0; i < m+blockSize_; i++) { index[i] = i; }
1520 Teuchos::RCP<const MV> Vp = MVT::CloneView( *VV, index );
1521 MVT::MvTransMv( one, *Vp, *Utmp, A21 );
1524 for( i=keff; i<keff+m; i++)
1528 SDM A( m2, A_tmp.numCols() );
1529 A.multiply(
Teuchos::TRANS, Teuchos::NO_TRANS, one, GG, A_tmp, zero );
1537 char balanc=
'P', jobvl=
'N', jobvr=
'V', sense=
'N';
1538 int ld = A.numRows();
1540 int ldvl = ld, ldvr = ld;
1541 int info = 0,ilo = 0,ihi = 0;
1542 MagnitudeType abnrm = 0.0, bbnrm = 0.0;
1544 std::vector<ScalarType> beta(ld);
1545 std::vector<ScalarType> work(lwork);
1546 std::vector<MagnitudeType> rwork(lwork);
1547 std::vector<MagnitudeType> lscale(ld), rscale(ld);
1548 std::vector<MagnitudeType> rconde(ld), rcondv(ld);
1549 std::vector<int> iwork(ld+6);
1554 lapack.GGEVX(balanc, jobvl, jobvr, sense, ld, A.values(), ld, B.values(), ld, &wr[0], &wi[0],
1555 &beta[0], vl, ldvl, vr.values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0],
1556 &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &rwork[0],
1557 &iwork[0], bwork, &info);
1558 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK GGEVX failed to compute eigensolutions.");
1562 for( i=0; i<ld; i++ )
1563 w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] ) / std::abs( beta[i] );
1565 this->sort(w,ld,iperm);
1567 bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
1570 for( i=0; i<recycledBlocks_; i++ )
1571 for( j=0; j<ld; j++ )
1572 PP(j,i) = vr(j,iperm[ld-recycledBlocks_+i]);
1574 if(scalarTypeIsComplex==
false) {
1577 if (wi[iperm[ld-recycledBlocks_]] != 0.0) {
1579 for ( i=ld-recycledBlocks_; i<ld; i++ )
1580 if (wi[iperm[i]] != 0.0) countImag++;
1582 if (countImag % 2) xtraVec =
true;
1586 if (wi[iperm[ld-recycledBlocks_]] > 0.0) {
1587 for( j=0; j<ld; j++ )
1588 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]+1);
1591 for( j=0; j<ld; j++ )
1592 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]-1);
1600 return recycledBlocks_+1;
1602 return recycledBlocks_;
1606 template<
class ScalarType,
class MV,
class OP>
1607 int BlockGCRODRSolMgr<ScalarType,MV,OP>::getHarmonicVecsKryl(
int m,
const SDM& HH, SDM& PP){
1608 bool xtraVec =
false;
1609 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1610 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1613 std::vector<MagnitudeType> wr(m), wi(m);
1619 std::vector<MagnitudeType> w(m);
1622 std::vector<int> iperm(m);
1626 std::vector<ScalarType> work(lwork);
1627 std::vector<MagnitudeType> rwork(lwork);
1633 builtRecycleSpace_ =
true;
1637 Teuchos::RCP<SDM> harmRitzMatrix = Teuchos::rcp(
new SDM( m, blockSize_));
1640 for(
int i=0; i<=blockSize_-1; i++) (*harmRitzMatrix)[blockSize_-1-i][harmRitzMatrix->numRows()-1-i] = 1;
1643 lapack.GESV(m, blockSize_, HHt.values(), HHt.stride(), &iperm[0], harmRitzMatrix->values(), harmRitzMatrix->stride(), &info);
1645 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK GESV failed to compute a solution.");
1649 Teuchos::SerialDenseMatrix<int, ScalarType> H_lbl(Teuchos::View, HH, blockSize_, blockSize_, (HH).numRows()-blockSize_, (HH).numCols()-blockSize_ );
1650 Teuchos::SerialDenseMatrix<int, ScalarType> H_lbl_t( H_lbl,
Teuchos::TRANS );
1655 Teuchos::RCP<SDM> Htemp = Teuchos::null;
1656 Htemp = Teuchos::rcp(
new SDM(H_lbl_t.numRows(), H_lbl_t.numCols()));
1657 Htemp -> multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, H_lbl_t, H_lbl, zero);
1658 H_lbl_t.assign(*Htemp);
1660 Htemp = Teuchos::rcp(
new SDM(harmRitzMatrix -> numRows(), harmRitzMatrix -> numCols()));
1661 Htemp -> multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, *harmRitzMatrix, H_lbl_t, zero);
1662 harmRitzMatrix -> assign(*Htemp);
1665 int harmColIndex, HHColIndex;
1666 Htemp = Teuchos::rcp(
new SDM(Teuchos::Copy,HH,HH.numRows()-blockSize_,HH.numCols()));
1667 for(
int i = 0; i<blockSize_; i++){
1668 harmColIndex = harmRitzMatrix -> numCols() - i -1;
1670 for(
int j=0; j<m; j++) (*Htemp)[HHColIndex][j] += (*harmRitzMatrix)[harmColIndex][j];
1672 harmRitzMatrix = Htemp;
1682 lapack.GEEV(
'N',
'V', m, harmRitzMatrix -> values(), harmRitzMatrix -> stride(), &wr[0], &wi[0],
1683 vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &rwork[0], &info);
1685 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK GEEV failed to compute eigensolutions.");
1688 for(
int i=0; i<m; ++i ) w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] );
1690 this->sort(w, m, iperm);
1692 bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
1695 for(
int i=0; i<recycledBlocks_; ++i )
1696 for(
int j=0; j<m; j++ )
1697 PP(j,i) = vr(j,iperm[i]);
1699 if(scalarTypeIsComplex==
false) {
1702 if (wi[iperm[recycledBlocks_-1]] != 0.0) {
1704 for (
int i=0; i<recycledBlocks_; ++i )
1705 if (wi[iperm[i]] != 0.0) countImag++;
1707 if (countImag % 2) xtraVec =
true;
1711 if (wi[iperm[recycledBlocks_-1]] > 0.0) {
1712 for(
int j=0; j<m; ++j )
1713 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]+1);
1716 for(
int j=0; j<m; ++j )
1717 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]-1);
1725 printer_->stream(
Debug) <<
"Recycled " << recycledBlocks_+1 <<
" vectors" << std::endl;
1726 return recycledBlocks_+1;
1729 printer_->stream(
Debug) <<
"Recycled " << recycledBlocks_ <<
" vectors" << std::endl;
1730 return recycledBlocks_;
1735 template<
class ScalarType,
class MV,
class OP>
1736 void BlockGCRODRSolMgr<ScalarType,MV,OP>::sort(std::vector<MagnitudeType>& dlist,
int n, std::vector<int>& iperm) {
1737 int l, r, j, i, flag;
1739 MagnitudeType dRR, dK;
1764 if (dlist[j] > dlist[j - 1]) j = j + 1;
1765 if (dlist[j - 1] > dK) {
1766 dlist[i - 1] = dlist[j - 1];
1767 iperm[i - 1] = iperm[j - 1];
1780 dlist[r] = dlist[0];
1781 iperm[r] = iperm[0];
1795 template<
class ScalarType,
class MV,
class OP>
1799 using Teuchos::rcp_const_cast;
1803 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1804 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1805 std::vector<int> index(numBlocks_+1);
1821 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1822 int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1826 std::vector<int> currIdx;
1828 if ( adaptiveBlockSize_ ) {
1829 blockSize_ = numCurrRHS;
1830 currIdx.resize( numCurrRHS );
1831 for (
int i=0; i<numCurrRHS; ++i)
1832 currIdx[i] = startPtr+i;
1835 currIdx.resize( blockSize_ );
1836 for (
int i=0; i<numCurrRHS; ++i)
1837 currIdx[i] = startPtr+i;
1838 for (
int i=numCurrRHS; i<blockSize_; ++i)
1843 problem_->setLSIndex( currIdx );
1849 loaDetected_ =
false;
1852 bool isConverged =
true;
1855 initializeStateStorage();
1858 Teuchos::ParameterList plist;
1860 while (numRHS2Solve > 0){
1870 printer_->stream(
Debug) <<
" Now solving RHS index " << currIdx[0] <<
" using recycled subspace of dimension " << keff << std::endl << std::endl;
1874 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1875 RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1876 RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1877 problem_->apply( *Utmp, *Ctmp );
1879 RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1883 SDM Ftmp( Teuchos::View, *F_, keff, keff );
1884 int rank = ortho_->normalize(*Ctmp, rcp(&Ftmp,
false));
1886 TEUCHOS_TEST_FOR_EXCEPTION(rank != keff,
BlockGCRODRSolMgrOrthoFailure,
"Belos::BlockGCRODRSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
1891 ipiv_.resize(Ftmp.numRows());
1892 lapack.GETRF(Ftmp.numRows(),Ftmp.numCols(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&info);
1893 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1896 int lwork = Ftmp.numRows();
1897 work_.resize(lwork);
1898 lapack.GETRI(Ftmp.numRows(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1899 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
BlockGCRODRSolMgrLAPACKFailure,
"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
1902 MVT::MvTimesMatAddMv( one, *Utmp, Ftmp, zero, *U1tmp );
1907 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1908 Ctmp = MVT::CloneViewNonConst( *C_, index );
1909 Utmp = MVT::CloneView( *U_, index );
1912 SDM Ctr(keff,blockSize_);
1913 problem_->computeCurrPrecResVec( &*R_ );
1914 MVT::MvTransMv( one, *Ctmp, *R_, Ctr );
1917 RCP<MV> update = MVT::Clone( *problem_->getCurrLHSVec(), blockSize_ );
1918 MVT::MvInit( *update, 0.0 );
1919 MVT::MvTimesMatAddMv( one, *Utmp, Ctr, one, *update );
1920 problem_->updateSolution( update,
true );
1923 MVT::MvTimesMatAddMv( -one, *Ctmp, Ctr, one, *R_ );
1931 if (V_ == Teuchos::null) {
1933 Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
1934 V_ = MVT::Clone( *rhsMV, (numBlocks_+1)*blockSize_ );
1938 if (MVT::GetNumberVecs(*V_) < (numBlocks_+1)*blockSize_ ) {
1939 Teuchos::RCP<const MV> tmp = V_;
1940 V_ = MVT::Clone( *tmp, (numBlocks_+1)*blockSize_ );
1945 printer_->stream(
Debug) <<
" No recycled subspace available for RHS index " << std::endl << std::endl;
1947 Teuchos::ParameterList primeList;
1950 primeList.set(
"Num Blocks",numBlocks_-1);
1951 primeList.set(
"Block Size",blockSize_);
1952 primeList.set(
"Recycled Blocks",0);
1953 primeList.set(
"Keep Hessenberg",
true);
1954 primeList.set(
"Initialize Hessenberg",
true);
1956 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1957 if (blockSize_*static_cast<ptrdiff_t>(numBlocks_) > dim) {
1958 ptrdiff_t tmpNumBlocks = 0;
1959 if (blockSize_ == 1)
1960 tmpNumBlocks = dim / blockSize_;
1962 tmpNumBlocks = ( dim - blockSize_) / blockSize_;
1963 printer_->stream(
Warnings) <<
"Belos::BlockGmresSolMgr::solve(): Warning! Requested Krylov subspace dimension is larger than operator dimension!" 1964 << std::endl <<
"The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << tmpNumBlocks << std::endl;
1965 primeList.set(
"Num Blocks",Teuchos::as<int>(tmpNumBlocks));
1969 primeList.set(
"Num Blocks",numBlocks_-1);
1972 Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > block_gmres_iter;
1976 block_gmres_iter->setSize( blockSize_, numBlocks_-1 );
1979 Teuchos::RCP<MV> V_0;
1980 if (currIdx[blockSize_-1] == -1) {
1981 V_0 = MVT::Clone( *(problem_->getInitPrecResVec()), blockSize_ );
1982 problem_->computeCurrPrecResVec( &*V_0 );
1985 V_0 = MVT::CloneCopy( *(problem_->getInitPrecResVec()), currIdx );
1989 Teuchos::RCP<SDM > z_0 = Teuchos::rcp(
new SDM( blockSize_, blockSize_ ) );
1992 int rank = ortho_->normalize( *V_0, z_0 );
2001 block_gmres_iter->initializeGmres(newstate);
2003 bool primeConverged =
false;
2006 printer_->stream(
Debug) <<
" Preparing to Iterate!!!!" << std::endl << std::endl;
2007 block_gmres_iter->iterate();
2012 if ( convTest_->getStatus() ==
Passed ) {
2013 printer_->stream(
Debug) <<
"We converged during the prime the pump stage" << std::endl << std::endl;
2014 primeConverged = !(expConvTest_->getLOADetected());
2015 if ( expConvTest_->getLOADetected() ) {
2017 loaDetected_ =
true;
2018 printer_->stream(
Warnings) <<
"Belos::BlockGmresSolMgr::solve(): Warning! Solver has experienced a loss of accuracy!" << std::endl;
2024 else if( maxIterTest_->getStatus() ==
Passed ) {
2026 primeConverged =
false;
2032 printer_->stream(
Debug) <<
" We did not converge on priming cycle of Block GMRES. Therefore we recycle and restart. " << std::endl << std::endl;
2037 if (blockSize_ != 1) {
2038 printer_->stream(
Errors) <<
"Error! Caught std::exception in BlockGmresIter::iterate() at iteration " 2039 << block_gmres_iter->getNumIters() << std::endl
2040 << e.what() << std::endl;
2041 if (convTest_->getStatus() !=
Passed)
2042 primeConverged =
false;
2046 block_gmres_iter->updateLSQR( block_gmres_iter->getCurSubspaceDim() );
2048 sTest_->checkStatus( &*block_gmres_iter );
2049 if (convTest_->getStatus() !=
Passed)
2050 isConverged =
false;
2053 catch (
const std::exception &e) {
2054 printer_->stream(
Errors) <<
"Error! Caught std::exception in BlockGmresIter::iterate() at iteration " 2055 << block_gmres_iter->getNumIters() << std::endl
2056 << e.what() << std::endl;
2064 RCP<MV> update = block_gmres_iter->getCurrentUpdate();
2065 problem_->updateSolution( update,
true );
2069 problem_->computeCurrPrecResVec( &*R_ );
2072 newstate = block_gmres_iter->getState();
2075 H_->assign(*(newstate.
H));
2084 V_ = rcp_const_cast<MV>(newstate.
V);
2085 newstate.
V.release();
2087 buildRecycleSpaceKryl(keff, block_gmres_iter);
2088 printer_->stream(
Debug) <<
"Generated recycled subspace using RHS index " << currIdx[0] <<
" of dimension " << keff << std::endl << std::endl;
2092 if (primeConverged) {
2113 problem_->setCurrLS();
2116 startPtr += numCurrRHS;
2117 numRHS2Solve -= numCurrRHS;
2118 if ( numRHS2Solve > 0 ) {
2119 numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
2120 if ( adaptiveBlockSize_ ) {
2121 blockSize_ = numCurrRHS;
2122 currIdx.resize( numCurrRHS );
2123 for (
int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2126 currIdx.resize( blockSize_ );
2127 for (
int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2128 for (
int i=numCurrRHS; i<blockSize_; ++i) currIdx[i] = -1;
2131 problem_->setLSIndex( currIdx );
2134 currIdx.resize( numRHS2Solve );
2144 Teuchos::ParameterList blockgcrodrList;
2145 blockgcrodrList.set(
"Num Blocks",numBlocks_);
2146 blockgcrodrList.set(
"Block Size",blockSize_);
2147 blockgcrodrList.set(
"Recycled Blocks",keff);
2149 Teuchos::RCP<BlockGCRODRIter<ScalarType,MV,OP> > block_gcrodr_iter;
2153 index.resize( blockSize_ );
2154 for(
int ii = 0; ii < blockSize_; ii++) index[ii] = ii;
2155 Teuchos::RCP<MV> V0 = MVT::CloneViewNonConst( *V_, index );
2158 MVT::Assign(*R_,*V0);
2161 for(
int i=0; i < keff; i++){ index[i] = i;};
2162 B_ = rcp(
new SDM(Teuchos::View, *G_, keff, numBlocks_*blockSize_, 0, keff));
2163 H_ = rcp(
new SDM(Teuchos::View, *G_, (numBlocks_-1)*blockSize_ + blockSize_, (numBlocks_-1)*blockSize_, keff ,keff ));
2167 newstate.
U = MVT::CloneViewNonConst(*U_, index);
2168 newstate.
C = MVT::CloneViewNonConst(*C_, index);
2170 newstate.
curDim = blockSize_;
2171 block_gcrodr_iter -> initialize(newstate);
2173 int numRestarts = 0;
2177 block_gcrodr_iter -> iterate();
2182 if( convTest_->getStatus() ==
Passed ) {
2190 else if(maxIterTest_->getStatus() ==
Passed ){
2192 isConverged =
false;
2199 else if (block_gcrodr_iter->getCurSubspaceDim() == block_gcrodr_iter->getMaxSubspaceDim()){
2204 Teuchos::RCP<MV> update = block_gcrodr_iter->getCurrentUpdate();
2205 problem_->updateSolution(update,
true);
2206 buildRecycleSpaceAugKryl(block_gcrodr_iter);
2208 printer_->stream(
Debug) <<
" Generated new recycled subspace using RHS index " << currIdx[0] <<
" of dimension " << keff << std::endl << std::endl;
2210 if(numRestarts >= maxRestarts_) {
2211 isConverged =
false;
2217 printer_ -> stream(
Debug) <<
" Performing restart number " << numRestarts <<
" of " << maxRestarts_ << std::endl << std::endl;
2220 problem_->computeCurrPrecResVec( &*R_ );
2221 index.resize( blockSize_ );
2222 for (
int ii=0; ii<blockSize_; ++ii) index[ii] = ii;
2223 Teuchos::RCP<MV> V0 = MVT::CloneViewNonConst( *V_, index );
2224 MVT::SetBlock(*R_,index,*V0);
2228 index.resize( numBlocks_*blockSize_ );
2229 for (
int ii=0; ii<(numBlocks_*blockSize_); ++ii) index[ii] = ii;
2230 restartState.
V = MVT::CloneViewNonConst( *V_, index );
2231 index.resize( keff );
2232 for (
int ii=0; ii<keff; ++ii) index[ii] = ii;
2233 restartState.
U = MVT::CloneViewNonConst( *U_, index );
2234 restartState.
C = MVT::CloneViewNonConst( *C_, index );
2235 B_ = rcp(
new SDM(Teuchos::View, *G_, keff, (numBlocks_-1)*blockSize_, 0, keff));
2236 H_ = rcp(
new SDM(Teuchos::View, *G_, numBlocks_*blockSize_, (numBlocks_-1)*blockSize_, keff ,keff ));
2237 restartState.
B = B_;
2238 restartState.
H = H_;
2239 restartState.
curDim = blockSize_;
2240 block_gcrodr_iter->initialize(restartState);
2249 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Belos::BlockGCRODRSolMgr::solve(): Invalid return from BlockGCRODRIter::iterate().");
2255 block_gcrodr_iter->updateLSQR( block_gcrodr_iter->getCurSubspaceDim() );
2257 sTest_->checkStatus( &*block_gcrodr_iter );
2258 if (convTest_->getStatus() !=
Passed) isConverged =
false;
2261 catch(
const std::exception &e){
2262 printer_->stream(
Errors) <<
"Error! Caught exception in BlockGCRODRIter::iterate() at iteration " 2263 << block_gcrodr_iter->getNumIters() << std::endl
2264 << e.what() << std::endl;
2271 Teuchos::RCP<MV> update = block_gcrodr_iter->getCurrentUpdate();
2272 problem_->updateSolution( update,
true );
2291 problem_->setCurrLS();
2294 startPtr += numCurrRHS;
2295 numRHS2Solve -= numCurrRHS;
2296 if ( numRHS2Solve > 0 ) {
2297 numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
2298 if ( adaptiveBlockSize_ ) {
2299 blockSize_ = numCurrRHS;
2300 currIdx.resize( numCurrRHS );
2301 for (
int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2304 currIdx.resize( blockSize_ );
2305 for (
int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2306 for (
int i=numCurrRHS; i<blockSize_; ++i) currIdx[i] = -1;
2309 problem_->setLSIndex( currIdx );
2312 currIdx.resize( numRHS2Solve );
2316 if (!builtRecycleSpace_) {
2317 buildRecycleSpaceAugKryl(block_gcrodr_iter);
2318 printer_->stream(
Debug) <<
" Generated new recycled subspace using RHS index " << currIdx[0] <<
" of dimension " << keff << std::endl << std::endl;
2326 #ifdef BELOS_TEUCHOS_TIME_MONITOR 2332 numIters_ = maxIterTest_->getNumIters();
2335 const std::vector<MagnitudeType>* pTestValues = NULL;
2336 pTestValues = impConvTest_->getTestValue();
2337 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
2338 "Belos::BlockGCRODRSolMgr::solve(): The implicit convergence test's " 2339 "getTestValue() method returned NULL. Please report this bug to the " 2340 "Belos developers.");
2341 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
2342 "Belos::BlockGCRODRSolMgr::solve(): The implicit convergence test's " 2343 "getTestValue() method returned a vector of length zero. Please report " 2344 "this bug to the Belos developers.");
2348 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
Collection of types and exceptions used within the Belos solvers.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
Class which manages the output and verbosity of the Belos solvers.
Structure to contain pointers to BlockGCRODRIter state variables.
virtual ~BlockGCRODRSolMgr()
Destructor.
ScaleType
The type of scaling to use on the residual norm value.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B
The projection of the Krylov subspace against the recycled subspace *
Exception thrown to signal error in a status test during Belos::StatusTest::checkStatus().
This class implements the block GMRES iteration, where a block Krylov subspace is constructed...
Teuchos::RCP< const MV > V
The current Krylov basis.
BlockGCRODRIterOrthoFailure is thrown when the BlockGCRODRIter object is unable to compute independen...
Thrown when the linear problem was not set up correctly.
A factory class for generating StatusTestOutput objects.
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > ¶ms)
Set the parameters the solver should use to solve the linear problem.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > z
The current right-hand side of the least squares system RY = Z.
An implementation of StatusTestResNorm using a family of residual norms.
Structure to contain pointers to GmresIteration state variables.
int getNumIters() const
Get the iteration count for the most recent call to solve().
Belos::StatusTest class for specifying a maximum number of iterations.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Get a parameter list containing the valid parameters for this object.
BlockGCRODRSolMgrRecyclingFailure(const std::string &what_arg)
void reset(const ResetType type)
Performs a reset of the solver manager specified by the ResetType.
Thrown if any problem occurs in using or creating the recycle subspace.
A factory class for generating StatusTestOutput objects.
Teuchos::RCP< MV > V
The current Krylov basis.
BlockGCRODRSolMgrOrthoFailure(const std::string &what_arg)
int curDim
The current dimension of the reduction.
Traits class which defines basic operations on multivectors.
std::string description() const
A description of the Block GCRODR solver manager.
Belos::StatusTest for logically combining several status tests.
A Belos::StatusTest class for specifying a maximum number of iterations.
ResetType
How to reset the solver.
BlockGCRODRSolMgrLAPACKFailure(const std::string &what_arg)
A solver manager for the Block GCRO-DR (Block Recycling GMRES) linear solver.
Pure virtual base class which describes the basic interface for a solver manager. ...
int curDim
The current dimension of the reduction.
Belos concrete class for performing the block GCRO-DR (block GMRES with recycling) iteration...
bool isLOADetected() const
Whether a loss of accuracy was detected during the most recent solve.
BlockGCRODRSolMgr()
Default constructor.
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem)
Set the linear problem to solve on the next call to solve().
MagnitudeType achievedTol() const
Get the residual for the most recent call to solve().
ReturnType
Whether the Belos solve converged for all linear systems.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get current linear problem being solved for in this object.
The Belos::SolverManager is a templated virtual base class that defines the basic interface that any ...
Thrown when the solution manager was unable to orthogonalize the basis vectors.
Teuchos::RCP< MV > U
The recycled subspace and its projection.
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const
Get a parameter list containing the current parameters for this object.
Belos concrete class for performing the block GMRES iteration.
Implementation of the Block GCRO-DR (Block Recycling GMRES) 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.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
A class for extending the status testing capabilities of Belos via logical combinations.
Class which defines basic traits for the operator type.
Parent class to all Belos exceptions.
Pure virtual base class which augments the basic interface for a Gmres linear solver iteration...
GmresIterationOrthoFailure is thrown when the GmresIteration object is unable to compute independent ...
Belos header file which uses auto-configuration information to include necessary C++ headers...
BlockGCRODRSolMgrLinearProblemFailure(const std::string &what_arg)
ReturnType solve()
Solve the current linear problem.
Thrown when an LAPACK call fails.
Belos concrete class for performing the block, flexible GMRES iteration.
OutputType
Style of output used to display status test information.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.