184#if defined(HAVE_TEUCHOSCORE_CXX11)
185# if defined(HAVE_TEUCHOS_COMPLEX)
186 #if defined(HAVE_TEUCHOS_LONG_DOUBLE)
187 static_assert (std::is_same<ScalarType, std::complex<float> >::value ||
188 std::is_same<ScalarType, std::complex<double> >::value ||
189 std::is_same<ScalarType, float>::value ||
190 std::is_same<ScalarType, double>::value ||
191 std::is_same<ScalarType, long double>::value,
192 "Belos::GCRODRSolMgr: ScalarType must be one of the four "
193 "types (S,D,C,Z) supported by LAPACK or long double (largely not impl'd).");
195 static_assert (std::is_same<ScalarType, std::complex<float> >::value ||
196 std::is_same<ScalarType, std::complex<double> >::value ||
197 std::is_same<ScalarType, float>::value ||
198 std::is_same<ScalarType, double>::value,
199 "Belos::GCRODRSolMgr: ScalarType must be one of the four "
200 "types (S,D,C,Z) supported by LAPACK.");
203 #if defined(HAVE_TEUCHOS_LONG_DOUBLE)
204 static_assert (std::is_same<ScalarType, float>::value ||
205 std::is_same<ScalarType, double>::value ||
206 std::is_same<ScalarType, long double>::value,
207 "Belos::GCRODRSolMgr: ScalarType must be float, double or long double. "
208 "Complex arithmetic support is currently disabled. To "
209 "enable it, set Teuchos_ENABLE_COMPLEX=ON.");
211 static_assert (std::is_same<ScalarType, float>::value ||
212 std::is_same<ScalarType, double>::value,
213 "Belos::GCRODRSolMgr: ScalarType must be float or double. "
214 "Complex arithmetic support is currently disabled. To "
215 "enable it, set Teuchos_ENABLE_COMPLEX=ON.");
223 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
224 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
MagnitudeType;
225 typedef Teuchos::ScalarTraits<MagnitudeType>
MT;
292 const Teuchos::RCP<Teuchos::ParameterList> &pl);
298 Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const override {
314 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters()
const override;
327 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
328 return Teuchos::tuple(timerSolve_);
360 void setParameters(
const Teuchos::RCP<Teuchos::ParameterList> ¶ms )
override;
372 bool set = problem_->setProblem();
374 throw "Could not set problem.";
418 std::string description()
const override;
428 void initializeStateStorage();
437 int getHarmonicVecs1(
int m,
438 const Teuchos::SerialDenseMatrix<int,ScalarType>& HH,
439 Teuchos::SerialDenseMatrix<int,ScalarType>& PP);
446 int getHarmonicVecs2(
int keff,
int m,
447 const Teuchos::SerialDenseMatrix<int,ScalarType>& HH,
448 const Teuchos::RCP<const MV>& VV,
449 Teuchos::SerialDenseMatrix<int,ScalarType>& PP);
452 void sort(std::vector<MagnitudeType>& dlist,
int n, std::vector<int>& iperm);
458 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >
problem_;
465 Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
sTest_;
468 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> >
expConvTest_, impConvTest_;
474 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> >
ortho_;
480 static constexpr double orthoKappa_default_ = 0.0;
481 static constexpr int maxRestarts_default_ = 100;
482 static constexpr int maxIters_default_ = 1000;
483 static constexpr int numBlocks_default_ = 50;
484 static constexpr int blockSize_default_ = 1;
485 static constexpr int recycledBlocks_default_ = 5;
488 static constexpr int outputFreq_default_ = -1;
489 static constexpr const char * impResScale_default_ =
"Norm of Preconditioned Initial Residual";
490 static constexpr const char * expResScale_default_ =
"Norm of Initial Residual";
491 static constexpr const char * label_default_ =
"Belos";
492 static constexpr const char * orthoType_default_ =
"ICGS";
517 Teuchos::RCP<MV> U_,
C_;
520 Teuchos::RCP<MV> U1_,
C1_;
523 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
H2_;
524 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
H_;
525 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
B_;
526 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
PP_;
527 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
HP_;
530 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
R_;
622setParameters (
const Teuchos::RCP<Teuchos::ParameterList> ¶ms)
624 using Teuchos::isParameterType;
625 using Teuchos::getParameter;
627 using Teuchos::ParameterList;
628 using Teuchos::parameterList;
631 using Teuchos::rcp_dynamic_cast;
632 using Teuchos::rcpFromRef;
633 using Teuchos::Exceptions::InvalidParameter;
634 using Teuchos::Exceptions::InvalidParameterName;
635 using Teuchos::Exceptions::InvalidParameterType;
639 RCP<const ParameterList> defaultParams = getValidParameters();
657 if (params_.is_null()) {
658 params_ = parameterList (*defaultParams);
666 if (params_ != params) {
672 params_ = parameterList (*params);
707 params_->validateParametersAndSetDefaults (*defaultParams);
711 if (params->isParameter (
"Maximum Restarts")) {
712 maxRestarts_ = params->get(
"Maximum Restarts", maxRestarts_default_);
715 params_->set (
"Maximum Restarts", maxRestarts_);
719 if (params->isParameter (
"Maximum Iterations")) {
720 maxIters_ = params->get (
"Maximum Iterations", maxIters_default_);
723 params_->set (
"Maximum Iterations", maxIters_);
724 if (! maxIterTest_.is_null())
725 maxIterTest_->setMaxIters (maxIters_);
729 if (params->isParameter (
"Num Blocks")) {
730 numBlocks_ = params->get (
"Num Blocks", numBlocks_default_);
731 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
732 "Belos::GCRODRSolMgr: The \"Num Blocks\" parameter must "
733 "be strictly positive, but you specified a value of "
734 << numBlocks_ <<
".");
736 params_->set (
"Num Blocks", numBlocks_);
740 if (params->isParameter (
"Num Recycled Blocks")) {
741 recycledBlocks_ = params->get (
"Num Recycled Blocks",
742 recycledBlocks_default_);
743 TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks_ <= 0, std::invalid_argument,
744 "Belos::GCRODRSolMgr: The \"Num Recycled Blocks\" "
745 "parameter must be strictly positive, but you specified "
746 "a value of " << recycledBlocks_ <<
".");
747 TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks_ >= numBlocks_, std::invalid_argument,
748 "Belos::GCRODRSolMgr: The \"Num Recycled Blocks\" "
749 "parameter must be less than the \"Num Blocks\" "
750 "parameter, but you specified \"Num Recycled Blocks\" "
751 "= " << recycledBlocks_ <<
" and \"Num Blocks\" = "
752 << numBlocks_ <<
".");
754 params_->set(
"Num Recycled Blocks", recycledBlocks_);
760 if (params->isParameter (
"Timer Label")) {
761 std::string tempLabel = params->get (
"Timer Label", label_default_);
764 if (tempLabel != label_) {
766 params_->set (
"Timer Label", label_);
767 std::string solveLabel = label_ +
": GCRODRSolMgr total solve time";
768#ifdef BELOS_TEUCHOS_TIME_MONITOR
769 timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
771 if (ortho_ != Teuchos::null) {
772 ortho_->setLabel( label_ );
778 if (params->isParameter (
"Verbosity")) {
779 if (isParameterType<int> (*params,
"Verbosity")) {
780 verbosity_ = params->get (
"Verbosity", verbosity_default_);
782 verbosity_ = (int) getParameter<Belos::MsgType> (*params,
"Verbosity");
785 params_->set (
"Verbosity", verbosity_);
788 if (! printer_.is_null())
789 printer_->setVerbosity (verbosity_);
793 if (params->isParameter (
"Output Style")) {
794 if (isParameterType<int> (*params,
"Output Style")) {
795 outputStyle_ = params->get (
"Output Style", outputStyle_default_);
797 outputStyle_ = (int) getParameter<OutputType> (*params,
"Output Style");
801 params_->set (
"Output Style", outputStyle_);
819 if (params->isParameter (
"Output Stream")) {
821 outputStream_ = getParameter<RCP<std::ostream> > (*params,
"Output Stream");
822 }
catch (InvalidParameter&) {
823 outputStream_ = rcpFromRef (std::cout);
830 if (outputStream_.is_null()) {
831 outputStream_ = rcp (
new Teuchos::oblackholestream);
834 params_->set (
"Output Stream", outputStream_);
837 if (! printer_.is_null()) {
838 printer_->setOStream (outputStream_);
844 if (params->isParameter (
"Output Frequency")) {
845 outputFreq_ = params->get (
"Output Frequency", outputFreq_default_);
849 params_->set(
"Output Frequency", outputFreq_);
850 if (! outputTest_.is_null())
851 outputTest_->setOutputFrequency (outputFreq_);
858 if (printer_.is_null()) {
869 bool changedOrthoType =
false;
870 if (params->isParameter (
"Orthogonalization")) {
871 const std::string& tempOrthoType =
872 params->get (
"Orthogonalization", orthoType_default_);
875 std::ostringstream os;
876 os <<
"Belos::GCRODRSolMgr: Invalid orthogonalization name \""
877 << tempOrthoType <<
"\". The following are valid options "
878 <<
"for the \"Orthogonalization\" name parameter: ";
880 throw std::invalid_argument (os.str());
882 if (tempOrthoType != orthoType_) {
883 changedOrthoType =
true;
884 orthoType_ = tempOrthoType;
886 params_->set (
"Orthogonalization", orthoType_);
902 RCP<ParameterList> orthoParams;
905 using Teuchos::sublist;
907 const std::string paramName (
"Orthogonalization Parameters");
910 orthoParams = sublist (params_, paramName,
true);
911 }
catch (InvalidParameter&) {
918 orthoParams = sublist (params_, paramName,
true);
921 TEUCHOS_TEST_FOR_EXCEPTION(orthoParams.is_null(), std::logic_error,
922 "Failed to get orthogonalization parameters. "
923 "Please report this bug to the Belos developers.");
928 if (ortho_.is_null() || changedOrthoType) {
934 label_, orthoParams);
942 typedef Teuchos::ParameterListAcceptor PLA;
943 RCP<PLA> pla = rcp_dynamic_cast<PLA> (ortho_);
949 label_, orthoParams);
951 pla->setParameterList (orthoParams);
963 if (params->isParameter (
"Orthogonalization Constant")) {
965 if (params->isType<
MagnitudeType> (
"Orthogonalization Constant")) {
966 orthoKappa = params->get (
"Orthogonalization Constant", orthoKappa);
969 orthoKappa = params->get (
"Orthogonalization Constant", orthoKappa_default_);
972 if (orthoKappa > 0) {
973 orthoKappa_ = orthoKappa;
975 params_->set(
"Orthogonalization Constant", orthoKappa_);
977 if (orthoType_ ==
"DGKS" && ! ortho_.is_null()) {
984 rcp_dynamic_cast<ortho_man_type>(ortho_)->
setDepTol (orthoKappa_);
994 if (params->isParameter(
"Convergence Tolerance")) {
995 if (params->isType<
MagnitudeType> (
"Convergence Tolerance")) {
996 convTol_ = params->get (
"Convergence Tolerance",
1004 params_->set (
"Convergence Tolerance", convTol_);
1005 if (! impConvTest_.is_null())
1007 if (! expConvTest_.is_null())
1008 expConvTest_->setTolerance (convTol_);
1012 if (params->isParameter (
"Implicit Residual Scaling")) {
1013 std::string tempImpResScale =
1014 getParameter<std::string> (*params,
"Implicit Residual Scaling");
1017 if (impResScale_ != tempImpResScale) {
1019 impResScale_ = tempImpResScale;
1022 params_->set(
"Implicit Residual Scaling", impResScale_);
1032 if (! impConvTest_.is_null()) {
1038 impConvTest_ = null;
1045 if (params->isParameter(
"Explicit Residual Scaling")) {
1046 std::string tempExpResScale =
1047 getParameter<std::string> (*params,
"Explicit Residual Scaling");
1050 if (expResScale_ != tempExpResScale) {
1052 expResScale_ = tempExpResScale;
1055 params_->set(
"Explicit Residual Scaling", expResScale_);
1058 if (! expConvTest_.is_null()) {
1064 expConvTest_ = null;
1075 if (maxIterTest_.is_null())
1080 if (impConvTest_.is_null()) {
1081 impConvTest_ = rcp (
new StatusTestResNorm_t (convTol_));
1087 if (expConvTest_.is_null()) {
1088 expConvTest_ = rcp (
new StatusTestResNorm_t (convTol_));
1089 expConvTest_->defineResForm (StatusTestResNorm_t::Explicit,
Belos::TwoNorm);
1095 if (convTest_.is_null()) {
1096 convTest_ = rcp (
new StatusTestCombo_t (StatusTestCombo_t::SEQ,
1104 sTest_ = rcp (
new StatusTestCombo_t (StatusTestCombo_t::OR,
1110 outputTest_ = stoFactory.
create (printer_, sTest_, outputFreq_,
1114 std::string solverDesc =
" GCRODR ";
1115 outputTest_->setSolverDesc( solverDesc );
1118 if (timerSolve_.is_null()) {
1119 std::string solveLabel = label_ +
": GCRODRSolMgr total solve time";
1120#ifdef BELOS_TEUCHOS_TIME_MONITOR
1121 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
1203 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1206 Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
1207 if (rhsMV == Teuchos::null) {
1214 TEUCHOS_TEST_FOR_EXCEPTION(
static_cast<ptrdiff_t
>(numBlocks_) > MVT::GetGlobalLength(*rhsMV),std::invalid_argument,
1215 "Belos::GCRODRSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!");
1218 if (U_ == Teuchos::null) {
1219 U_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1223 if (MVT::GetNumberVecs(*U_) < recycledBlocks_+1) {
1224 Teuchos::RCP<const MV> tmp = U_;
1225 U_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1230 if (C_ == Teuchos::null) {
1231 C_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1235 if (MVT::GetNumberVecs(*C_) < recycledBlocks_+1) {
1236 Teuchos::RCP<const MV> tmp = C_;
1237 C_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1242 if (V_ == Teuchos::null) {
1243 V_ = MVT::Clone( *rhsMV, numBlocks_+1 );
1247 if (MVT::GetNumberVecs(*V_) < numBlocks_+1) {
1248 Teuchos::RCP<const MV> tmp = V_;
1249 V_ = MVT::Clone( *tmp, numBlocks_+1 );
1254 if (U1_ == Teuchos::null) {
1255 U1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1259 if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1260 Teuchos::RCP<const MV> tmp = U1_;
1261 U1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1266 if (C1_ == Teuchos::null) {
1267 C1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1271 if (MVT::GetNumberVecs(*C1_) < recycledBlocks_+1) {
1272 Teuchos::RCP<const MV> tmp = C1_;
1273 C1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1278 if (r_ == Teuchos::null)
1279 r_ = MVT::Clone( *rhsMV, 1 );
1282 tau_.resize(recycledBlocks_+1);
1285 work_.resize(recycledBlocks_+1);
1288 ipiv_.resize(recycledBlocks_+1);
1291 if (H2_ == Teuchos::null)
1292 H2_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 ) );
1294 if ( (H2_->numRows() != numBlocks_+recycledBlocks_+2) || (H2_->numCols() != numBlocks_+recycledBlocks_+1) )
1295 H2_->reshape( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 );
1297 H2_->putScalar(zero);
1300 if (R_ == Teuchos::null)
1301 R_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( recycledBlocks_+1, recycledBlocks_+1 ) );
1303 if ( (R_->numRows() != recycledBlocks_+1) || (R_->numCols() != recycledBlocks_+1) )
1304 R_->reshape( recycledBlocks_+1, recycledBlocks_+1 );
1306 R_->putScalar(zero);
1309 if (PP_ == Teuchos::null)
1310 PP_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, recycledBlocks_+1 ) );
1312 if ( (PP_->numRows() != numBlocks_+recycledBlocks_+2) || (PP_->numCols() != recycledBlocks_+1) )
1313 PP_->reshape( numBlocks_+recycledBlocks_+2, recycledBlocks_+1 );
1317 if (HP_ == Teuchos::null)
1318 HP_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 ) );
1320 if ( (HP_->numRows() != numBlocks_+recycledBlocks_+2) || (HP_->numCols() != numBlocks_+recycledBlocks_+1) )
1321 HP_->reshape( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 );
1337 if (!isSet_) { setParameters( params_ ); }
1339 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1340 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1341 std::vector<int> index(numBlocks_+1);
1343 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,
GCRODRSolMgrLinearProblemFailure,
"Belos::GCRODRSolMgr::solve(): Linear problem is not a valid object.");
1345 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),
GCRODRSolMgrLinearProblemFailure,
"Belos::GCRODRSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
1348 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1349 std::vector<int> currIdx(1);
1353 problem_->setLSIndex( currIdx );
1356 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1357 if (
static_cast<ptrdiff_t
>(numBlocks_) > dim) {
1358 numBlocks_ = Teuchos::as<int>(dim);
1360 "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl <<
1361 " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl;
1362 params_->set(
"Num Blocks", numBlocks_);
1366 bool isConverged =
true;
1369 initializeStateStorage();
1373 Teuchos::ParameterList plist;
1375 plist.set(
"Num Blocks",numBlocks_);
1376 plist.set(
"Recycled Blocks",recycledBlocks_);
1381 RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_iter;
1384 int prime_iterations = 0;
1388#ifdef BELOS_TEUCHOS_TIME_MONITOR
1389 Teuchos::TimeMonitor slvtimer(*timerSolve_);
1392 while ( numRHS2Solve > 0 ) {
1395 builtRecycleSpace_ =
false;
1398 outputTest_->reset();
1406 "Belos::GCRODRSolMgr::solve(): Requested size of recycled subspace is not consistent with the current recycle subspace.");
1408 printer_->stream(
Debug) <<
" Now solving RHS index " << currIdx[0] <<
" using recycled subspace of dimension " << keff << std::endl << std::endl;
1411 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1412 RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1413 RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1414 problem_->apply( *Utmp, *Ctmp );
1416 RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1420 Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff, keff );
1421 int rank = ortho_->normalize(*Ctmp, rcp(&Rtmp,
false));
1423 TEUCHOS_TEST_FOR_EXCEPTION(rank != keff,
GCRODRSolMgrOrthoFailure,
"Belos::GCRODRSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
1428 ipiv_.resize(Rtmp.numRows());
1429 lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
1430 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1433 int lwork = Rtmp.numRows();
1434 work_.resize(lwork);
1435 lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1436 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
1439 MVT::MvTimesMatAddMv( one, *Utmp, Rtmp, zero, *U1tmp );
1444 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1445 Ctmp = MVT::CloneViewNonConst( *C_, index );
1446 Utmp = MVT::CloneView( *U_, index );
1449 Teuchos::SerialDenseMatrix<int,ScalarType> Ctr(keff,1);
1450 problem_->computeCurrPrecResVec( &*r_ );
1451 MVT::MvTransMv( one, *Ctmp, *r_, Ctr );
1454 RCP<MV> update = MVT::Clone( *problem_->getCurrLHSVec(), 1 );
1455 MVT::MvInit( *update, 0.0 );
1456 MVT::MvTimesMatAddMv( one, *Utmp, Ctr, one, *update );
1457 problem_->updateSolution( update,
true );
1460 MVT::MvTimesMatAddMv( -one, *Ctmp, Ctr, one, *r_ );
1463 prime_iterations = 0;
1469 printer_->stream(
Debug) <<
" No recycled subspace available for RHS index " << currIdx[0] << std::endl << std::endl;
1471 Teuchos::ParameterList primeList;
1474 primeList.set(
"Num Blocks",numBlocks_);
1475 primeList.set(
"Recycled Blocks",0);
1478 RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_prime_iter;
1482 problem_->computeCurrPrecResVec( &*r_ );
1483 index.resize( 1 ); index[0] = 0;
1484 RCP<MV> v0 = MVT::CloneViewNonConst( *V_, index );
1485 MVT::SetBlock(*r_,index,*v0);
1489 index.resize( numBlocks_+1 );
1490 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1491 newstate.
V = MVT::CloneViewNonConst( *V_, index );
1492 newstate.
U = Teuchos::null;
1493 newstate.
C = Teuchos::null;
1494 newstate.
H = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, recycledBlocks_+1, recycledBlocks_+1 ) );
1495 newstate.
B = Teuchos::null;
1497 gcrodr_prime_iter->initialize(newstate);
1500 bool primeConverged =
false;
1502 gcrodr_prime_iter->iterate();
1505 if ( convTest_->getStatus() ==
Passed ) {
1507 primeConverged =
true;
1512 gcrodr_prime_iter->updateLSQR( gcrodr_prime_iter->getCurSubspaceDim() );
1515 sTest_->checkStatus( &*gcrodr_prime_iter );
1516 if (convTest_->getStatus() ==
Passed)
1517 primeConverged =
true;
1519 catch (
const std::exception &e) {
1520 printer_->stream(
Errors) <<
"Error! Caught exception in GCRODRIter::iterate() at iteration "
1521 << gcrodr_prime_iter->getNumIters() << std::endl
1522 << e.what() << std::endl;
1526 prime_iterations = gcrodr_prime_iter->getNumIters();
1529 RCP<MV> update = gcrodr_prime_iter->getCurrentUpdate();
1530 problem_->updateSolution( update,
true );
1533 newstate = gcrodr_prime_iter->getState();
1541 if (recycledBlocks_ < p+1) {
1543 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PPtmp = rcp (
new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View, *PP_, p, recycledBlocks_+1 ) );
1545 keff = getHarmonicVecs1( p, *newstate.
H, *PPtmp );
1547 PPtmp = rcp (
new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View, *PP_, p, keff ) );
1550 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1551 RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1552 RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1553 RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1555 for (
int ii=0; ii < p; ++ii) { index[ii] = ii; }
1556 RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1560 MVT::MvTimesMatAddMv( one, *Vtmp, *PPtmp, zero, *U1tmp );
1567 Teuchos::SerialDenseMatrix<int,ScalarType> Htmp( Teuchos::View, *H2_, p+1, p, recycledBlocks_+1,recycledBlocks_+1);
1568 Teuchos::SerialDenseMatrix<int,ScalarType> HPtmp( Teuchos::View, *HP_, p+1, keff );
1569 HPtmp.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, Htmp, *PPtmp, zero );
1576 lapack.GEQRF (HPtmp.numRows (), HPtmp.numCols (), HPtmp.values (),
1577 HPtmp.stride (), &tau_[0], &work_[0], lwork, &info);
1578 TEUCHOS_TEST_FOR_EXCEPTION(
1580 " LAPACK's _GEQRF failed to compute a workspace size.");
1588 lwork = std::abs (
static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work_[0])));
1589 work_.resize (lwork);
1590 lapack.GEQRF (HPtmp.numRows (), HPtmp.numCols (), HPtmp.values (),
1591 HPtmp.stride (), &tau_[0], &work_[0], lwork, &info);
1592 TEUCHOS_TEST_FOR_EXCEPTION(
1594 " LAPACK's _GEQRF failed to compute a QR factorization.");
1598 Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff, keff );
1599 for (
int ii = 0; ii < keff; ++ii) {
1600 for (
int jj = ii; jj < keff; ++jj) {
1601 Rtmp(ii,jj) = HPtmp(ii,jj);
1608 lapack.UNGQR (HPtmp.numRows (), HPtmp.numCols (), HPtmp.numCols (),
1609 HPtmp.values (), HPtmp.stride (), &tau_[0], &work_[0],
1611 TEUCHOS_TEST_FOR_EXCEPTION(
1613 "LAPACK's _UNGQR failed to construct the Q factor.");
1618 index.resize (p + 1);
1619 for (
int ii = 0; ii < (p+1); ++ii) {
1622 Vtmp = MVT::CloneView( *V_, index );
1623 MVT::MvTimesMatAddMv( one, *Vtmp, HPtmp, zero, *Ctmp );
1630 ipiv_.resize(Rtmp.numRows());
1631 lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
1632 TEUCHOS_TEST_FOR_EXCEPTION(
1634 "LAPACK's _GETRF failed to compute an LU factorization.");
1643 lwork = Rtmp.numRows();
1644 work_.resize(lwork);
1645 lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1646 TEUCHOS_TEST_FOR_EXCEPTION(
1648 "LAPACK's _GETRI failed to invert triangular matrix.");
1651 MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
1653 printer_->stream(
Debug)
1654 <<
" Generated recycled subspace using RHS index " << currIdx[0]
1655 <<
" of dimension " << keff << std::endl << std::endl;
1660 if (primeConverged) {
1662 problem_->setCurrLS();
1666 if (numRHS2Solve > 0) {
1668 problem_->setLSIndex (currIdx);
1671 currIdx.resize (numRHS2Solve);
1681 gcrodr_iter->setSize( keff, numBlocks_ );
1684 gcrodr_iter->resetNumIters(prime_iterations);
1687 outputTest_->resetNumCalls();
1690 problem_->computeCurrPrecResVec( &*r_ );
1691 index.resize( 1 ); index[0] = 0;
1692 RCP<MV> v0 = MVT::CloneViewNonConst( *V_, index );
1693 MVT::SetBlock(*r_,index,*v0);
1697 index.resize( numBlocks_+1 );
1698 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1699 newstate.
V = MVT::CloneViewNonConst( *V_, index );
1700 index.resize( keff );
1701 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1702 newstate.
C = MVT::CloneViewNonConst( *C_, index );
1703 newstate.
U = MVT::CloneViewNonConst( *U_, index );
1704 newstate.
B = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, keff, numBlocks_, 0, keff ) );
1705 newstate.
H = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, keff, keff ) );
1707 gcrodr_iter->initialize(newstate);
1710 int numRestarts = 0;
1715 gcrodr_iter->iterate();
1722 if ( convTest_->getStatus() ==
Passed ) {
1731 else if ( maxIterTest_->getStatus() ==
Passed ) {
1733 isConverged =
false;
1741 else if ( gcrodr_iter->getCurSubspaceDim() == gcrodr_iter->getMaxSubspaceDim() ) {
1746 RCP<MV> update = gcrodr_iter->getCurrentUpdate();
1747 problem_->updateSolution( update,
true );
1749 buildRecycleSpace2(gcrodr_iter);
1751 printer_->stream(
Debug)
1752 <<
" Generated new recycled subspace using RHS index "
1753 << currIdx[0] <<
" of dimension " << keff << std::endl
1757 if (numRestarts >= maxRestarts_) {
1758 isConverged =
false;
1763 printer_->stream(
Debug)
1764 <<
" Performing restart number " << numRestarts <<
" of "
1765 << maxRestarts_ << std::endl << std::endl;
1768 problem_->computeCurrPrecResVec( &*r_ );
1769 index.resize( 1 ); index[0] = 0;
1770 RCP<MV> v00 = MVT::CloneViewNonConst( *V_, index );
1771 MVT::SetBlock(*r_,index,*v00);
1775 index.resize( numBlocks_+1 );
1776 for (
int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1777 restartState.
V = MVT::CloneViewNonConst( *V_, index );
1778 index.resize( keff );
1779 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1780 restartState.
U = MVT::CloneViewNonConst( *U_, index );
1781 restartState.
C = MVT::CloneViewNonConst( *C_, index );
1782 restartState.
B = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, keff, numBlocks_, 0, keff ) );
1783 restartState.
H = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, keff, keff ) );
1785 gcrodr_iter->initialize(restartState);
1798 TEUCHOS_TEST_FOR_EXCEPTION(
1799 true, std::logic_error,
"Belos::GCRODRSolMgr::solve: "
1800 "Invalid return from GCRODRIter::iterate().");
1805 gcrodr_iter->updateLSQR( gcrodr_iter->getCurSubspaceDim() );
1808 sTest_->checkStatus( &*gcrodr_iter );
1809 if (convTest_->getStatus() !=
Passed)
1810 isConverged =
false;
1813 catch (
const std::exception& e) {
1815 <<
"Error! Caught exception in GCRODRIter::iterate() at iteration "
1816 << gcrodr_iter->getNumIters() << std::endl << e.what() << std::endl;
1823 RCP<MV> update = gcrodr_iter->getCurrentUpdate();
1824 problem_->updateSolution( update,
true );
1827 problem_->setCurrLS();
1832 if (!builtRecycleSpace_) {
1833 buildRecycleSpace2(gcrodr_iter);
1834 printer_->stream(
Debug)
1835 <<
" Generated new recycled subspace using RHS index " << currIdx[0]
1836 <<
" of dimension " << keff << std::endl << std::endl;
1841 if (numRHS2Solve > 0) {
1843 problem_->setLSIndex (currIdx);
1846 currIdx.resize (numRHS2Solve);
1854#ifdef BELOS_TEUCHOS_TIME_MONITOR
1859 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails) );
1863 numIters_ = maxIterTest_->getNumIters ();
1875 const std::vector<MagnitudeType>* pTestValues = expConvTest_->getTestValue();
1876 if (pTestValues == NULL || pTestValues->size() < 1) {
1877 pTestValues = impConvTest_->getTestValue();
1879 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1880 "Belos::GCRODRSolMgr::solve(): The implicit convergence test's getTestValue() "
1881 "method returned NULL. Please report this bug to the Belos developers.");
1882 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1883 "Belos::GCRODRSolMgr::solve(): The implicit convergence test's getTestValue() "
1884 "method returned a vector of length zero. Please report this bug to the "
1885 "Belos developers.");
1890 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1900 MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
1901 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1903 std::vector<MagnitudeType> d(keff);
1904 std::vector<ScalarType> dscalar(keff);
1905 std::vector<int> index(numBlocks_+1);
1917 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1918 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1920 dscalar.resize(keff);
1921 MVT::MvNorm( *Utmp, d );
1922 for (
int i=0; i<keff; ++i) {
1924 dscalar[i] = (ScalarType)d[i];
1926 MVT::MvScale( *Utmp, dscalar );
1930 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H2tmp = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, p+keff+1, p+keff ) );
1933 for (
int i=0; i<keff; ++i) {
1934 (*H2tmp)(i,i) = d[i];
1941 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p+keff, recycledBlocks_+1 );
1942 keff_new = getHarmonicVecs2( keff, p, *H2tmp, oldState.
V, PPtmp );
1949 Teuchos::RCP<MV> U1tmp;
1951 index.resize( keff );
1952 for (
int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1953 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1954 index.resize( keff_new );
1955 for (
int ii=0; ii<keff_new; ++ii) { index[ii] = ii; }
1956 U1tmp = MVT::CloneViewNonConst( *U1_, index );
1957 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, keff, keff_new );
1958 MVT::MvTimesMatAddMv( one, *Utmp, PPtmp, zero, *U1tmp );
1964 for (
int ii=0; ii < p; ii++) { index[ii] = ii; }
1965 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1966 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p, keff_new, keff );
1967 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *U1tmp );
1971 Teuchos::SerialDenseMatrix<int,ScalarType> HPtmp( Teuchos::View, *HP_, p+keff+1, keff_new );
1973 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p+keff, keff_new );
1974 HPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*H2tmp,PPtmp,zero);
1978 int info = 0, lwork = -1;
1979 tau_.resize (keff_new);
1980 lapack.GEQRF (HPtmp.numRows (), HPtmp.numCols (), HPtmp.values (),
1981 HPtmp.stride (), &tau_[0], &work_[0], lwork, &info);
1982 TEUCHOS_TEST_FOR_EXCEPTION(
1984 "LAPACK's _GEQRF failed to compute a workspace size.");
1990 lwork = std::abs (
static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work_[0])));
1991 work_.resize (lwork);
1992 lapack.GEQRF (HPtmp.numRows (), HPtmp.numCols (), HPtmp.values (),
1993 HPtmp.stride (), &tau_[0], &work_[0], lwork, &info);
1994 TEUCHOS_TEST_FOR_EXCEPTION(
1996 "LAPACK's _GEQRF failed to compute a QR factorization.");
2000 Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff_new, keff_new );
2001 for(
int i=0;i<keff_new;i++) {
for(
int j=i;j<keff_new;j++) Rtmp(i,j) = HPtmp(i,j); }
2007 lapack.UNGQR (HPtmp.numRows (), HPtmp.numCols (), HPtmp.numCols (),
2008 HPtmp.values (), HPtmp.stride (), &tau_[0], &work_[0],
2010 TEUCHOS_TEST_FOR_EXCEPTION(
2012 "LAPACK's _UNGQR failed to construct the Q factor.");
2019 Teuchos::RCP<MV> C1tmp;
2022 for (
int i=0; i < keff; i++) { index[i] = i; }
2023 Teuchos::RCP<const MV> Ctmp = MVT::CloneView( *C_, index );
2024 index.resize(keff_new);
2025 for (
int i=0; i < keff_new; i++) { index[i] = i; }
2026 C1tmp = MVT::CloneViewNonConst( *C1_, index );
2027 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *HP_, keff, keff_new );
2028 MVT::MvTimesMatAddMv( one, *Ctmp, PPtmp, zero, *C1tmp );
2032 index.resize( p+1 );
2033 for (
int i=0; i < p+1; ++i) { index[i] = i; }
2034 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
2035 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *HP_, p+1, keff_new, keff, 0 );
2036 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *C1tmp );
2045 ipiv_.resize(Rtmp.numRows());
2046 lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
2047 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
2050 lwork = Rtmp.numRows();
2051 work_.resize(lwork);
2052 lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
2053 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to compute an LU factorization.");
2056 index.resize(keff_new);
2057 for (
int i=0; i < keff_new; i++) { index[i] = i; }
2058 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
2059 MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
2063 if (keff != keff_new) {
2065 gcrodr_iter->setSize( keff, numBlocks_ );
2067 Teuchos::SerialDenseMatrix<int,ScalarType> b1( Teuchos::View, *H2_, recycledBlocks_+2, 1, 0, recycledBlocks_ );
2077 const Teuchos::SerialDenseMatrix<int,ScalarType>& HH,
2078 Teuchos::SerialDenseMatrix<int,ScalarType>& PP) {
2080 bool xtraVec =
false;
2081 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
2084 std::vector<MagnitudeType> wr(m), wi(m);
2087 Teuchos::SerialDenseMatrix<int,ScalarType> vr(m,m,
false);
2090 std::vector<MagnitudeType> w(m);
2093 std::vector<int> iperm(m);
2099 builtRecycleSpace_ =
true;
2102 Teuchos::SerialDenseMatrix<int, ScalarType> HHt( HH, Teuchos::TRANS );
2103 Teuchos::SerialDenseVector<int, ScalarType> e_m( m );
2105 lapack.GESV(m, 1, HHt.values(), HHt.stride(), &iperm[0], e_m.values(), e_m.stride(), &info);
2106 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK GESV failed to compute a solution.");
2109 ScalarType d = HH(m, m-1) * HH(m, m-1);
2110 Teuchos::SerialDenseMatrix<int, ScalarType> harmHH( Teuchos::Copy, HH, m, m );
2111 for( i=0; i<m; ++i )
2112 harmHH(i, m-1) += d * e_m[i];
2121 std::vector<ScalarType> work(1);
2122 std::vector<MagnitudeType> rwork(2*m);
2125 lapack.GEEV(
'N',
'V', m, harmHH.values(), harmHH.stride(), &wr[0], &wi[0],
2126 vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &rwork[0], &info);
2128 lwork = std::abs (
static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work[0])));
2129 work.resize( lwork );
2131 lapack.GEEV(
'N',
'V', m, harmHH.values(), harmHH.stride(), &wr[0], &wi[0],
2132 vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &rwork[0], &info);
2133 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK GEEV failed to compute eigensolutions.");
2136 for( i=0; i<m; ++i )
2137 w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] );
2140 this->sort(w, m, iperm);
2142 const bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
2145 for( i=0; i<recycledBlocks_; ++i ) {
2146 for( j=0; j<m; j++ ) {
2147 PP(j,i) = vr(j,iperm[i]);
2151 if(!scalarTypeIsComplex) {
2154 if (wi[iperm[recycledBlocks_-1]] != 0.0) {
2156 for ( i=0; i<recycledBlocks_; ++i ) {
2157 if (wi[iperm[i]] != 0.0)
2166 if (wi[iperm[recycledBlocks_-1]] > 0.0) {
2167 for( j=0; j<m; ++j ) {
2168 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]+1);
2172 for( j=0; j<m; ++j ) {
2173 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]-1);
2182 return recycledBlocks_+1;
2185 return recycledBlocks_;
2193 const Teuchos::SerialDenseMatrix<int,ScalarType>& HH,
2194 const Teuchos::RCP<const MV>& VV,
2195 Teuchos::SerialDenseMatrix<int,ScalarType>& PP) {
2197 int m2 = HH.numCols();
2198 bool xtraVec =
false;
2199 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
2200 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
2201 std::vector<int> index;
2204 std::vector<MagnitudeType> wr(m2), wi(m2);
2207 std::vector<MagnitudeType> w(m2);
2210 Teuchos::SerialDenseMatrix<int,ScalarType> vr(m2,m2,
false);
2213 std::vector<int> iperm(m2);
2216 builtRecycleSpace_ =
true;
2221 Teuchos::SerialDenseMatrix<int,ScalarType> B(m2,m2,
false);
2222 B.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,HH,HH,zero);
2226 Teuchos::SerialDenseMatrix<int,ScalarType> A_tmp( keffloc+m+1, keffloc+m );
2229 index.resize(keffloc);
2230 for (i=0; i<keffloc; ++i) { index[i] = i; }
2231 Teuchos::RCP<const MV> Ctmp = MVT::CloneView( *C_, index );
2232 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
2233 Teuchos::SerialDenseMatrix<int,ScalarType> A11( Teuchos::View, A_tmp, keffloc, keffloc );
2234 MVT::MvTransMv( one, *Ctmp, *Utmp, A11 );
2237 Teuchos::SerialDenseMatrix<int,ScalarType> A21( Teuchos::View, A_tmp, m+1, keffloc, keffloc );
2239 for (i=0; i < m+1; i++) { index[i] = i; }
2240 Teuchos::RCP<const MV> Vp = MVT::CloneView( *VV, index );
2241 MVT::MvTransMv( one, *Vp, *Utmp, A21 );
2244 for( i=keffloc; i<keffloc+m; i++ ) {
2249 Teuchos::SerialDenseMatrix<int,ScalarType> A( m2, A_tmp.numCols() );
2250 A.multiply( Teuchos::TRANS, Teuchos::NO_TRANS, one, HH, A_tmp, zero );
2258 char balanc=
'P', jobvl=
'N', jobvr=
'V', sense=
'N';
2259 int ld = A.numRows();
2261 int ldvl = ld, ldvr = ld;
2262 int info = 0,ilo = 0,ihi = 0;
2265 std::vector<ScalarType> beta(ld);
2266 std::vector<ScalarType> work(lwork);
2267 std::vector<MagnitudeType> rwork(lwork);
2268 std::vector<MagnitudeType> lscale(ld), rscale(ld);
2269 std::vector<MagnitudeType> rconde(ld), rcondv(ld);
2270 std::vector<int> iwork(ld+6);
2275 lapack.GGEVX(balanc, jobvl, jobvr, sense, ld, A.values(), ld, B.values(), ld, &wr[0], &wi[0],
2276 &beta[0], vl, ldvl, vr.values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0],
2277 &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &rwork[0],
2278 &iwork[0], bwork, &info);
2279 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
GCRODRSolMgrLAPACKFailure,
"Belos::GCRODRSolMgr::solve(): LAPACK GGEVX failed to compute eigensolutions.");
2283 for( i=0; i<ld; i++ ) {
2284 w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot (wr[i]*wr[i] + wi[i]*wi[i]) /
2285 Teuchos::ScalarTraits<ScalarType>::magnitude (beta[i]);
2289 this->sort(w,ld,iperm);
2291 const bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
2294 for( i=0; i<recycledBlocks_; i++ ) {
2295 for( j=0; j<ld; j++ ) {
2296 PP(j,i) = vr(j,iperm[ld-recycledBlocks_+i]);
2300 if(!scalarTypeIsComplex) {
2303 if (wi[iperm[ld-recycledBlocks_]] != 0.0) {
2305 for ( i=ld-recycledBlocks_; i<ld; i++ ) {
2306 if (wi[iperm[i]] != 0.0)
2315 if (wi[iperm[ld-recycledBlocks_]] > 0.0) {
2316 for( j=0; j<ld; j++ ) {
2317 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]+1);
2321 for( j=0; j<ld; j++ ) {
2322 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]-1);
2331 return recycledBlocks_+1;
2334 return recycledBlocks_;