35 #ifndef ANASAZI_TRACEMIN_BASE_HPP 36 #define ANASAZI_TRACEMIN_BASE_HPP 80 template <
class ScalarType,
class MV>
168 template <
class ScalarType,
class MV,
class OP>
242 void harmonicIterate();
347 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
getResNorms();
355 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
getRes2Norms();
365 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
getRitzRes2Norms();
439 void setSize(
int blockSize,
int numBlocks);
458 typedef TraceMinRitzOp<ScalarType,MV,OP> tracemin_ritz_op_type;
459 typedef SaddleContainer<ScalarType,MV> saddle_container_type;
460 typedef SaddleOperator<ScalarType,MV,tracemin_ritz_op_type> saddle_op_type;
461 const MagnitudeType ONE;
462 const MagnitudeType ZERO;
463 const MagnitudeType NANVAL;
484 timerLocal_, timerCompRes_, timerOrtho_, timerInit_;
490 bool checkV, checkX, checkMX,
491 checkKX, checkQ, checkKK;
492 CheckList() : checkV(
false),checkX(
false),
493 checkMX(
false),checkKX(
false),
494 checkQ(
false),checkKK(
false) {};
499 std::string accuracyCheck(
const CheckList &chk,
const std::string &where)
const;
504 int count_ApplyOp_, count_ApplyM_;
531 RCP<MV> X_, KX_, MX_, KV_, MV_, R_, V_;
545 std::vector<MagnitudeType> theta_, Rnorms_, R2norms_;
548 bool Rnorms_current_, R2norms_current_;
552 enum SaddleSolType saddleSolType_;
553 bool previouslyLeveled_;
554 MagnitudeType previousTrace_;
555 bool posSafeToShift_, negSafeToShift_;
556 MagnitudeType largestSafeShift_;
558 std::vector<ScalarType> ritzShifts_;
564 enum WhenToShiftType whenToShift_;
565 enum HowToShiftType howToShift_;
566 bool useMultipleShifts_;
567 bool considerClusters_;
568 bool projectAllVecs_;
569 bool projectLockedVecs_;
573 MagnitudeType traceThresh_;
574 MagnitudeType alpha_;
578 std::string shiftNorm_;
579 MagnitudeType shiftThresh_;
580 bool useRelShiftThresh_;
584 ScalarType getTrace()
const;
590 std::vector<ScalarType> getClusterResids();
593 void computeRitzShifts(
const std::vector<ScalarType>& clusterResids);
596 std::vector<ScalarType> computeTol();
598 void solveSaddlePointProblem(
RCP<MV> Delta);
600 void solveSaddleProj(
RCP<MV> Delta)
const;
603 void solveSaddleProjPrec(
RCP<MV> Delta)
const;
605 void solveSaddleSchur (
RCP<MV> Delta)
const;
607 void solveSaddleBDPrec (
RCP<MV> Delta)
const;
609 void solveSaddleHSSPrec (
RCP<MV> Delta)
const;
613 void computeRitzPairs();
619 void updateResidual();
623 virtual void harmonicAddToBasis(
const RCP<const MV> Delta) =0;
639 template <
class ScalarType,
class MV,
class OP>
648 ONE(
Teuchos::ScalarTraits<MagnitudeType>::one()),
649 ZERO(
Teuchos::ScalarTraits<MagnitudeType>::zero()),
650 NANVAL(
Teuchos::ScalarTraits<MagnitudeType>::nan()),
658 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
659 timerOp_(
Teuchos::TimeMonitor::getNewTimer(
"Anasazi: TraceMinBase::Operation Op*x")),
660 timerMOp_(
Teuchos::TimeMonitor::getNewTimer(
"Anasazi: TraceMinBase::Operation M*x")),
661 timerSaddle_(
Teuchos::TimeMonitor::getNewTimer(
"Anasazi: TraceMinBase::Solving saddle point problem")),
662 timerSortEval_(
Teuchos::TimeMonitor::getNewTimer(
"Anasazi: TraceMinBase::Sorting eigenvalues")),
663 timerDS_(
Teuchos::TimeMonitor::getNewTimer(
"Anasazi: TraceMinBase::Direct solve")),
664 timerLocal_(
Teuchos::TimeMonitor::getNewTimer(
"Anasazi: TraceMinBase::Local update")),
665 timerCompRes_(
Teuchos::TimeMonitor::getNewTimer(
"Anasazi: TraceMinBase::Computing residuals")),
666 timerOrtho_(
Teuchos::TimeMonitor::getNewTimer(
"Anasazi: TraceMinBase::Orthogonalization")),
667 timerInit_(
Teuchos::TimeMonitor::getNewTimer(
"Anasazi: TraceMinBase::Initialization")),
676 auxVecs_(
Teuchos::Array<
RCP<const MV> >(0) ),
677 MauxVecs_(
Teuchos::Array<
RCP<const MV> >(0) ),
680 Rnorms_current_(false),
681 R2norms_current_(false),
683 previouslyLeveled_(false),
684 previousTrace_(ZERO),
685 posSafeToShift_(false),
686 negSafeToShift_(false),
687 largestSafeShift_(ZERO),
691 "Anasazi::TraceMinBase::constructor: user passed null problem pointer.");
693 "Anasazi::TraceMinBase::constructor: user passed null sort manager pointer.");
695 "Anasazi::TraceMinBase::constructor: user passed null output manager pointer.");
697 "Anasazi::TraceMinBase::constructor: user passed null status test pointer.");
699 "Anasazi::TraceMinBase::constructor: user passed null orthogonalization manager pointer.");
701 "Anasazi::TraceMinBase::constructor: problem is not hermitian.");
704 Op_ = problem_->getOperator();
705 MOp_ = problem_->getM();
706 Prec_ = problem_->getPrec();
707 hasM_ = (MOp_ != Teuchos::null);
710 saddleSolType_ = params.
get(
"Saddle Solver Type", PROJECTED_KRYLOV_SOLVER);
711 TEUCHOS_TEST_FOR_EXCEPTION(saddleSolType_ != PROJECTED_KRYLOV_SOLVER && saddleSolType_ != SCHUR_COMPLEMENT_SOLVER && saddleSolType_ != BD_PREC_MINRES && saddleSolType_ != HSS_PREC_GMRES, std::invalid_argument,
712 "Anasazi::TraceMin::constructor: Invalid value for \"Saddle Solver Type\"; valid options are PROJECTED_KRYLOV_SOLVER, SCHUR_COMPLEMENT_SOLVER, and BD_PREC_MINRES.");
715 whenToShift_ = params.
get(
"When To Shift", ALWAYS_SHIFT);
716 TEUCHOS_TEST_FOR_EXCEPTION(whenToShift_ != NEVER_SHIFT && whenToShift_ != SHIFT_WHEN_TRACE_LEVELS && whenToShift_ != SHIFT_WHEN_RESID_SMALL && whenToShift_ != ALWAYS_SHIFT, std::invalid_argument,
717 "Anasazi::TraceMin::constructor: Invalid value for \"When To Shift\"; valid options are \"NEVER_SHIFT\", \"SHIFT_WHEN_TRACE_LEVELS\", \"SHIFT_WHEN_RESID_SMALL\", and \"ALWAYS_SHIFT\".");
719 traceThresh_ = params.
get(
"Trace Threshold", 2e-2);
721 "Anasazi::TraceMin::constructor: Invalid value for \"Trace Threshold\"; Must be positive.");
723 howToShift_ = params.
get(
"How To Choose Shift", ADJUSTED_RITZ_SHIFT);
724 TEUCHOS_TEST_FOR_EXCEPTION(howToShift_ != LARGEST_CONVERGED_SHIFT && howToShift_ != ADJUSTED_RITZ_SHIFT && howToShift_ != RITZ_VALUES_SHIFT && howToShift_ != EXPERIMENTAL_SHIFT, std::invalid_argument,
725 "Anasazi::TraceMin::constructor: Invalid value for \"How To Choose Shift\"; valid options are \"LARGEST_CONVERGED_SHIFT\", \"ADJUSTED_RITZ_SHIFT\", \"RITZ_VALUES_SHIFT\".");
727 useMultipleShifts_ = params.
get(
"Use Multiple Shifts",
true);
729 shiftThresh_ = params.
get(
"Shift Threshold", 1e-2);
730 useRelShiftThresh_ = params.
get(
"Relative Shift Threshold",
true);
731 shiftNorm_ = params.
get(
"Shift Norm",
"2");
733 "Anasazi::TraceMin::constructor: Invalid value for \"Shift Norm\"; valid options are \"2\", \"M\".");
735 considerClusters_ = params.
get(
"Consider Clusters",
true);
737 projectAllVecs_ = params.
get(
"Project All Vectors",
true);
738 projectLockedVecs_ = params.
get(
"Project Locked Vectors",
true);
739 useRHSR_ = params.
get(
"Use Residual as RHS",
true);
740 useHarmonic_ = params.
get(
"Use Harmonic Ritz Values",
false);
741 computeAllRes_ = params.
get(
"Compute All Residuals",
true);
744 int bs = params.
get(
"Block Size", problem_->getNEV());
745 int nb = params.
get(
"Num Blocks", 1);
748 NEV_ = problem_->getNEV();
751 ritzOp_ =
rcp (
new tracemin_ritz_op_type (Op_, MOp_, Prec_));
754 const int innerMaxIts = params.
get (
"Maximum Krylov Iterations", 200);
755 ritzOp_->setMaxIts (innerMaxIts);
757 alpha_ = params.
get (
"HSS: alpha", ONE);
763 template <
class ScalarType,
class MV,
class OP>
770 template <
class ScalarType,
class MV,
class OP>
773 TEUCHOS_TEST_FOR_EXCEPTION(blockSize < 1, std::invalid_argument,
"Anasazi::TraceMinBase::setSize(blocksize,numblocks): blocksize must be strictly positive.");
774 setSize(blockSize,numBlocks_);
780 template <
class ScalarType,
class MV,
class OP>
788 template <
class ScalarType,
class MV,
class OP>
796 template <
class ScalarType,
class MV,
class OP>
804 template <
class ScalarType,
class MV,
class OP>
806 return blockSize_*numBlocks_;
811 template <
class ScalarType,
class MV,
class OP>
813 if (!initialized_)
return 0;
820 template <
class ScalarType,
class MV,
class OP>
821 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
823 return getRes2Norms();
829 template <
class ScalarType,
class MV,
class OP>
831 std::vector<int> ret(curDim_,0);
838 template <
class ScalarType,
class MV,
class OP>
840 std::vector<Value<ScalarType> > ret(curDim_);
841 for (
int i=0; i<curDim_; ++i) {
842 ret[i].realpart = theta_[i];
843 ret[i].imagpart = ZERO;
851 template <
class ScalarType,
class MV,
class OP>
859 template <
class ScalarType,
class MV,
class OP>
867 template <
class ScalarType,
class MV,
class OP>
875 template <
class ScalarType,
class MV,
class OP>
881 if (Op_ != Teuchos::null) {
886 state.
KV = Teuchos::null;
887 state.
KX = Teuchos::null;
894 state.
MopV = Teuchos::null;
895 state.
MX = Teuchos::null;
899 state.
RV = ritzVecs_;
901 state.
T =
rcp(
new std::vector<MagnitudeType>(theta_.begin(),theta_.begin()+curDim_) );
903 state.
T =
rcp(
new std::vector<MagnitudeType>(0));
905 state.
ritzShifts =
rcp(
new std::vector<MagnitudeType>(ritzShifts_.begin(),ritzShifts_.begin()+blockSize_) );
915 template <
class ScalarType,
class MV,
class OP>
926 if (initialized_ ==
false) {
932 const int searchDim = blockSize_*numBlocks_;
935 RCP<MV> Delta = MVT::Clone(*X_,blockSize_);
940 while (tester_->checkStatus(
this) !=
Passed && (numBlocks_ == 1 || curDim_ < searchDim)) {
943 if (om_->isVerbosity(
Debug)) {
944 currentStatus( om_->stream(
Debug) );
953 solveSaddlePointProblem(Delta);
958 if (om_->isVerbosity(
Debug ) ) {
962 om_->print(
Debug, accuracyCheck(chk,
": after addToBasis(Delta)") );
968 if (om_->isVerbosity(
Debug ) ) {
972 om_->print(
Debug, accuracyCheck(chk,
": after computeKK()") );
981 if (om_->isVerbosity(
Debug ) ) {
985 om_->print(
Debug, accuracyCheck(chk,
": after computeX()") );
991 if (om_->isVerbosity(
Debug ) ) {
996 om_->print(
Debug, accuracyCheck(chk,
": after updateKXMX()") );
1009 template <
class ScalarType,
class MV,
class OP>
1014 if (initialized_ ==
false) {
1020 const int searchDim = blockSize_*numBlocks_;
1023 RCP<MV> Delta = MVT::Clone(*X_,blockSize_);
1028 while (tester_->checkStatus(
this) !=
Passed && (numBlocks_ == 1 || curDim_ < searchDim)) {
1031 if (om_->isVerbosity(
Debug)) {
1032 currentStatus( om_->stream(
Debug) );
1041 solveSaddlePointProblem(Delta);
1044 harmonicAddToBasis(Delta);
1046 if (om_->isVerbosity(
Debug ) ) {
1050 om_->print(
Debug, accuracyCheck(chk,
": after addToBasis(Delta)") );
1056 if (om_->isVerbosity(
Debug ) ) {
1060 om_->print(
Debug, accuracyCheck(chk,
": after computeKK()") );
1075 std::vector<int> dimind(nvecs);
1076 for(
int i=0; i<nvecs; i++)
1078 RCP<MV> lclX = MVT::CloneViewNonConst(*X_,dimind);
1079 std::vector<ScalarType> normvec(nvecs);
1080 orthman_->normMat(*lclX,normvec);
1083 for(
int i=0; i<nvecs; i++)
1084 normvec[i] = ONE/normvec[i];
1085 MVT::MvScale(*lclX,normvec);
1088 for(
int i=0; i<nvecs; i++)
1090 theta_[i] = theta_[i] * normvec[i] * normvec[i];
1093 if (om_->isVerbosity(
Debug ) ) {
1097 om_->print(
Debug, accuracyCheck(chk,
": after computeX()") );
1104 if(Op_ != Teuchos::null)
1106 RCP<MV> lclKX = MVT::CloneViewNonConst(*KX_,dimind);
1107 MVT::MvScale(*lclKX,normvec);
1111 RCP<MV> lclMX = MVT::CloneViewNonConst(*MX_,dimind);
1112 MVT::MvScale(*lclMX,normvec);
1115 if (om_->isVerbosity(
Debug ) ) {
1120 om_->print(
Debug, accuracyCheck(chk,
": after updateKXMX()") );
1132 template <
class ScalarType,
class MV,
class OP>
1153 template <
class ScalarType,
class MV,
class OP>
1160 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1164 previouslyLeveled_ =
false;
1168 harmonicInitialize(newstate);
1172 std::vector<int> bsind(blockSize_);
1173 for (
int i=0; i<blockSize_; ++i) bsind[i] = i;
1201 if (newstate.
V != Teuchos::null) {
1202 om_->stream(
Debug) <<
"Copying V from the user\n";
1205 "Anasazi::TraceMinBase::initialize(newstate): Vector length of V not correct." );
1207 "Anasazi::TraceMinBase::initialize(newstate): Rank of new state must be at least blockSize().");
1209 "Anasazi::TraceMinBase::initialize(newstate): Rank of new state must be less than getMaxSubspaceDim().");
1212 "Anasazi::TraceMinBase::initialize(newstate): Multivector for basis in new state must be as large as specified state rank.");
1214 curDim_ = newstate.
curDim;
1216 curDim_ = (int)(curDim_ / blockSize_)*blockSize_;
1219 "Anasazi::TraceMinBase::initialize(newstate): Rank of new state must be a multiple of getBlockSize().");
1222 std::vector<int> nevind(curDim_);
1223 for (
int i=0; i<curDim_; ++i) nevind[i] = i;
1224 if (newstate.
V != V_) {
1225 if (curDim_ < MVT::GetNumberVecs(*newstate.
V)) {
1226 newstate.
V = MVT::CloneView(*newstate.
V,nevind);
1228 MVT::SetBlock(*newstate.
V,nevind,*V_);
1230 lclV = MVT::CloneViewNonConst(*V_,nevind);
1237 "Anasazi::TraceMinBase::initialize(newstate): Eigenproblem did not specify initial vectors to clone from.");
1239 newstate.
X = Teuchos::null;
1240 newstate.
MX = Teuchos::null;
1241 newstate.
KX = Teuchos::null;
1242 newstate.
R = Teuchos::null;
1243 newstate.
T = Teuchos::null;
1244 newstate.
RV = Teuchos::null;
1245 newstate.
KK = Teuchos::null;
1246 newstate.
KV = Teuchos::null;
1247 newstate.
MopV = Teuchos::null;
1252 curDim_ = newstate.
curDim;
1254 curDim_ = MVT::GetNumberVecs(*ivec);
1257 curDim_ = (int)(curDim_ / blockSize_)*blockSize_;
1258 if (curDim_ > blockSize_*numBlocks_) {
1261 curDim_ = blockSize_*numBlocks_;
1263 bool userand =
false;
1268 curDim_ = blockSize_;
1271 std::vector<int> nevind(curDim_);
1272 for (
int i=0; i<curDim_; ++i) nevind[i] = i;
1278 lclV = MVT::CloneViewNonConst(*V_,nevind);
1282 MVT::MvRandom(*lclV);
1288 if(MVT::GetNumberVecs(*newstate.
V) > curDim_) {
1290 MVT::SetBlock(*helperMV,nevind,*lclV);
1293 MVT::SetBlock(*newstate.
V,nevind,*lclV);
1297 if(MVT::GetNumberVecs(*ivec) > curDim_) {
1298 ivec = MVT::CloneView(*ivec,nevind);
1301 MVT::SetBlock(*ivec,nevind,*lclV);
1307 std::vector<int> dimind(curDim_);
1308 for (
int i=0; i<curDim_; ++i) dimind[i] = i;
1311 if(hasM_ && newstate.
MopV == Teuchos::null)
1313 om_->stream(
Debug) <<
"Computing MV\n";
1315 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1318 count_ApplyM_+= curDim_;
1322 lclMV = MVT::CloneViewNonConst(*MV_,dimind);
1323 OPT::Apply(*MOp_,*lclV,*lclMV);
1328 om_->stream(
Debug) <<
"Copying MV\n";
1331 "Anasazi::TraceMinBase::initialize(newstate): Vector length of MV not correct." );
1334 "Anasazi::TraceMinBase::initialize(newstate): Number of vectors in MV not correct.");
1336 if(newstate.
MopV != MV_) {
1337 if(curDim_ < MVT::GetNumberVecs(*newstate.
MopV)) {
1338 newstate.
MopV = MVT::CloneView(*newstate.
MopV,dimind);
1340 MVT::SetBlock(*newstate.
MopV,dimind,*MV_);
1342 lclMV = MVT::CloneViewNonConst(*MV_,dimind);
1347 om_->stream(
Debug) <<
"There is no MV\n";
1356 om_->stream(
Debug) <<
"Project and normalize\n";
1358 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1363 newstate.
X = Teuchos::null;
1364 newstate.
MX = Teuchos::null;
1365 newstate.
KX = Teuchos::null;
1366 newstate.
R = Teuchos::null;
1367 newstate.
T = Teuchos::null;
1368 newstate.
RV = Teuchos::null;
1369 newstate.
KK = Teuchos::null;
1370 newstate.
KV = Teuchos::null;
1373 if(auxVecs_.size() > 0)
1375 rank = orthman_->projectAndNormalizeMat(*lclV, auxVecs_,
1377 Teuchos::null, lclMV, MauxVecs_);
1381 rank = orthman_->normalizeMat(*lclV,Teuchos::null,lclMV);
1385 "Anasazi::TraceMinBase::initialize(): Couldn't generate initial basis of full rank.");
1389 if(Op_ != Teuchos::null && newstate.
KV == Teuchos::null)
1391 om_->stream(
Debug) <<
"Computing KV\n";
1393 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1396 count_ApplyOp_+= curDim_;
1399 newstate.
X = Teuchos::null;
1400 newstate.
MX = Teuchos::null;
1401 newstate.
KX = Teuchos::null;
1402 newstate.
R = Teuchos::null;
1403 newstate.
T = Teuchos::null;
1404 newstate.
RV = Teuchos::null;
1405 newstate.
KK = Teuchos::null;
1406 newstate.
KV = Teuchos::null;
1408 lclKV = MVT::CloneViewNonConst(*KV_,dimind);
1409 OPT::Apply(*Op_,*lclV,*lclKV);
1412 else if(Op_ != Teuchos::null)
1414 om_->stream(
Debug) <<
"Copying MV\n";
1417 "Anasazi::TraceMinBase::initialize(newstate): Vector length of MV not correct." );
1420 "Anasazi::TraceMinBase::initialize(newstate): Number of vectors in KV not correct.");
1422 if (newstate.
KV != KV_) {
1423 if (curDim_ < MVT::GetNumberVecs(*newstate.
KV)) {
1424 newstate.
KV = MVT::CloneView(*newstate.
KV,dimind);
1426 MVT::SetBlock(*newstate.
KV,dimind,*KV_);
1428 lclKV = MVT::CloneViewNonConst(*KV_,dimind);
1432 om_->stream(
Debug) <<
"There is no KV\n";
1439 if(newstate.
KK == Teuchos::null)
1441 om_->stream(
Debug) <<
"Computing KK\n";
1444 newstate.
X = Teuchos::null;
1445 newstate.
MX = Teuchos::null;
1446 newstate.
KX = Teuchos::null;
1447 newstate.
R = Teuchos::null;
1448 newstate.
T = Teuchos::null;
1449 newstate.
RV = Teuchos::null;
1458 MVT::MvTransMv(ONE,*lclV,*lclKV,*lclKK);
1463 om_->stream(
Debug) <<
"Copying KK\n";
1467 "Anasazi::TraceMinBase::initialize(newstate): Projected matrix in new state must be as large as specified state rank.");
1471 if (newstate.
KK != KK_) {
1472 if (newstate.
KK->numRows() > curDim_ || newstate.
KK->numCols() > curDim_) {
1475 lclKK->assign(*newstate.
KK);
1480 if(newstate.
T == Teuchos::null || newstate.
RV == Teuchos::null)
1482 om_->stream(
Debug) <<
"Computing Ritz pairs\n";
1485 newstate.
X = Teuchos::null;
1486 newstate.
MX = Teuchos::null;
1487 newstate.
KX = Teuchos::null;
1488 newstate.
R = Teuchos::null;
1489 newstate.
T = Teuchos::null;
1490 newstate.
RV = Teuchos::null;
1497 om_->stream(
Debug) <<
"Copying Ritz pairs\n";
1500 std::invalid_argument,
"Anasazi::TraceMinBase::initialize(newstate): Size of T must be consistent with dimension of V.");
1503 "Anasazi::TraceMinBase::initialize(newstate): Ritz vectors in new state must be as large as specified state rank.");
1505 std::copy(newstate.
T->begin(),newstate.
T->end(),theta_.begin());
1508 if (newstate.
RV != ritzVecs_) {
1509 if (newstate.
RV->numRows() > curDim_ || newstate.
RV->numCols() > curDim_) {
1512 lclRV->assign(*newstate.
RV);
1517 if(newstate.
X == Teuchos::null)
1519 om_->stream(
Debug) <<
"Computing X\n";
1522 newstate.
MX = Teuchos::null;
1523 newstate.
KX = Teuchos::null;
1524 newstate.
R = Teuchos::null;
1531 om_->stream(
Debug) <<
"Copying X\n";
1533 if(computeAllRes_ ==
false)
1535 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.
X) != blockSize_ || MVT::GetGlobalLength(*newstate.
X) != MVT::GetGlobalLength(*X_),
1536 std::invalid_argument,
"Anasazi::TraceMinBase::initialize(newstate): Size of X must be consistent with block size and length of V.");
1538 if(newstate.
X != X_) {
1539 MVT::SetBlock(*newstate.
X,bsind,*X_);
1544 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.
X) != curDim_ || MVT::GetGlobalLength(*newstate.
X) != MVT::GetGlobalLength(*X_),
1545 std::invalid_argument,
"Anasazi::TraceMinBase::initialize(newstate): Size of X must be consistent with current dimension and length of V.");
1547 if(newstate.
X != X_) {
1548 MVT::SetBlock(*newstate.
X,dimind,*X_);
1555 if((Op_ != Teuchos::null && newstate.
KX == Teuchos::null) || (hasM_ && newstate.
MX == Teuchos::null))
1557 om_->stream(
Debug) <<
"Computing KX and MX\n";
1560 newstate.
R = Teuchos::null;
1567 om_->stream(
Debug) <<
"Copying KX and MX\n";
1568 if(Op_ != Teuchos::null)
1570 if(computeAllRes_ ==
false)
1573 std::invalid_argument,
"Anasazi::TraceMinBase::initialize(newstate): Size of KX must be consistent with block size and length of V.");
1575 if(newstate.
KX != KX_) {
1576 MVT::SetBlock(*newstate.
KX,bsind,*KX_);
1582 std::invalid_argument,
"Anasazi::TraceMinBase::initialize(newstate): Size of KX must be consistent with current dimension and length of V.");
1584 if (newstate.
KX != KX_) {
1585 MVT::SetBlock(*newstate.
KX,dimind,*KX_);
1592 if(computeAllRes_ ==
false)
1595 std::invalid_argument,
"Anasazi::TraceMinBase::initialize(newstate): Size of MX must be consistent with block size and length of V.");
1597 if (newstate.
MX != MX_) {
1598 MVT::SetBlock(*newstate.
MX,bsind,*MX_);
1604 std::invalid_argument,
"Anasazi::TraceMinBase::initialize(newstate): Size of MX must be consistent with current dimension and length of V.");
1606 if (newstate.
MX != MX_) {
1607 MVT::SetBlock(*newstate.
MX,dimind,*MX_);
1614 if(newstate.
R == Teuchos::null)
1616 om_->stream(
Debug) <<
"Computing R\n";
1623 om_->stream(
Debug) <<
"Copying R\n";
1625 if(computeAllRes_ ==
false)
1628 std::invalid_argument,
"Anasazi::TraceMinBase::initialize(newstate): vector length of newstate.R not correct." );
1630 std::invalid_argument,
"Anasazi::TraceMinBase::initialize(newstate): newstate.R must have at least block size vectors." );
1631 if (newstate.
R != R_) {
1632 MVT::SetBlock(*newstate.
R,bsind,*R_);
1638 std::invalid_argument,
"Anasazi::TraceMinBase::initialize(newstate): vector length of newstate.R not correct." );
1640 std::invalid_argument,
"Anasazi::TraceMinBase::initialize(newstate): newstate.R must have at least curDim vectors." );
1641 if (newstate.
R != R_) {
1642 MVT::SetBlock(*newstate.
R,dimind,*R_);
1648 Rnorms_current_ =
false;
1649 R2norms_current_ =
false;
1657 om_->stream(
Debug) <<
"Copying Ritz shifts\n";
1662 om_->stream(
Debug) <<
"Setting Ritz shifts to 0\n";
1663 for(
size_t i=0; i<ritzShifts_.size(); i++)
1664 ritzShifts_[i] = ZERO;
1667 for(
size_t i=0; i<ritzShifts_.size(); i++)
1668 om_->stream(
Debug) <<
"Ritz shifts[" << i <<
"] = " << ritzShifts_[i] << std::endl;
1671 initialized_ =
true;
1673 if (om_->isVerbosity(
Debug ) ) {
1682 om_->print(
Debug, accuracyCheck(chk,
": after initialize()") );
1686 if (om_->isVerbosity(
Debug)) {
1687 currentStatus( om_->stream(
Debug) );
1709 template <
class ScalarType,
class MV,
class OP>
1716 std::vector<int> bsind(blockSize_);
1717 for (
int i=0; i<blockSize_; ++i) bsind[i] = i;
1745 if (newstate.
V != Teuchos::null) {
1746 om_->stream(
Debug) <<
"Copying V from the user\n";
1749 "Anasazi::TraceMinBase::initialize(newstate): Vector length of V not correct." );
1751 "Anasazi::TraceMinBase::initialize(newstate): Rank of new state must be at least blockSize().");
1753 "Anasazi::TraceMinBase::initialize(newstate): Rank of new state must be less than getMaxSubspaceDim().");
1756 "Anasazi::TraceMinBase::initialize(newstate): Multivector for basis in new state must be as large as specified state rank.");
1758 curDim_ = newstate.
curDim;
1760 curDim_ = (int)(curDim_ / blockSize_)*blockSize_;
1763 "Anasazi::TraceMinBase::initialize(newstate): Rank of new state must be a multiple of getBlockSize().");
1766 std::vector<int> nevind(curDim_);
1767 for (
int i=0; i<curDim_; ++i) nevind[i] = i;
1768 if (newstate.
V != V_) {
1769 if (curDim_ < MVT::GetNumberVecs(*newstate.
V)) {
1770 newstate.
V = MVT::CloneView(*newstate.
V,nevind);
1772 MVT::SetBlock(*newstate.
V,nevind,*V_);
1774 lclV = MVT::CloneViewNonConst(*V_,nevind);
1781 "Anasazi::TraceMinBase::initialize(newstate): Eigenproblem did not specify initial vectors to clone from.");
1783 newstate.
X = Teuchos::null;
1784 newstate.
MX = Teuchos::null;
1785 newstate.
KX = Teuchos::null;
1786 newstate.
R = Teuchos::null;
1787 newstate.
T = Teuchos::null;
1788 newstate.
RV = Teuchos::null;
1789 newstate.
KK = Teuchos::null;
1790 newstate.
KV = Teuchos::null;
1791 newstate.
MopV = Teuchos::null;
1796 curDim_ = newstate.
curDim;
1798 curDim_ = MVT::GetNumberVecs(*ivec);
1801 curDim_ = (int)(curDim_ / blockSize_)*blockSize_;
1802 if (curDim_ > blockSize_*numBlocks_) {
1805 curDim_ = blockSize_*numBlocks_;
1807 bool userand =
false;
1812 curDim_ = blockSize_;
1815 std::vector<int> nevind(curDim_);
1816 for (
int i=0; i<curDim_; ++i) nevind[i] = i;
1822 lclV = MVT::CloneViewNonConst(*V_,nevind);
1826 MVT::MvRandom(*lclV);
1832 if(MVT::GetNumberVecs(*newstate.
V) > curDim_) {
1834 MVT::SetBlock(*helperMV,nevind,*lclV);
1837 MVT::SetBlock(*newstate.
V,nevind,*lclV);
1841 if(MVT::GetNumberVecs(*ivec) > curDim_) {
1842 ivec = MVT::CloneView(*ivec,nevind);
1845 MVT::SetBlock(*ivec,nevind,*lclV);
1853 newstate.
X = Teuchos::null;
1854 newstate.
MX = Teuchos::null;
1855 newstate.
KX = Teuchos::null;
1856 newstate.
R = Teuchos::null;
1857 newstate.
T = Teuchos::null;
1858 newstate.
RV = Teuchos::null;
1859 newstate.
KK = Teuchos::null;
1860 newstate.
KV = Teuchos::null;
1861 newstate.
MopV = Teuchos::null;
1865 std::vector<int> dimind(curDim_);
1866 for (
int i=0; i<curDim_; ++i) dimind[i] = i;
1869 if(auxVecs_.size() > 0)
1870 orthman_->projectMat(*lclV,auxVecs_);
1873 if(Op_ != Teuchos::null && newstate.
KV == Teuchos::null)
1875 om_->stream(
Debug) <<
"Computing KV\n";
1877 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1880 count_ApplyOp_+= curDim_;
1883 newstate.
X = Teuchos::null;
1884 newstate.
MX = Teuchos::null;
1885 newstate.
KX = Teuchos::null;
1886 newstate.
R = Teuchos::null;
1887 newstate.
T = Teuchos::null;
1888 newstate.
RV = Teuchos::null;
1889 newstate.
KK = Teuchos::null;
1891 lclKV = MVT::CloneViewNonConst(*KV_,dimind);
1892 OPT::Apply(*Op_,*lclV,*lclKV);
1895 else if(Op_ != Teuchos::null)
1897 om_->stream(
Debug) <<
"Copying KV\n";
1900 "Anasazi::TraceMinBase::initialize(newstate): Vector length of KV not correct." );
1903 "Anasazi::TraceMinBase::initialize(newstate): Number of vectors in KV not correct.");
1905 if (newstate.
KV != KV_) {
1906 if (curDim_ < MVT::GetNumberVecs(*newstate.
KV)) {
1907 newstate.
KV = MVT::CloneView(*newstate.
KV,dimind);
1909 MVT::SetBlock(*newstate.
KV,dimind,*KV_);
1911 lclKV = MVT::CloneViewNonConst(*KV_,dimind);
1915 om_->stream(
Debug) <<
"There is no KV\n";
1926 om_->stream(
Debug) <<
"Project and normalize\n";
1928 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1933 newstate.
MopV = Teuchos::null;
1934 newstate.
X = Teuchos::null;
1935 newstate.
MX = Teuchos::null;
1936 newstate.
KX = Teuchos::null;
1937 newstate.
R = Teuchos::null;
1938 newstate.
T = Teuchos::null;
1939 newstate.
RV = Teuchos::null;
1940 newstate.
KK = Teuchos::null;
1944 int rank = orthman_->normalizeMat(*lclKV,gamma);
1950 RCP<MV> tempMV = MVT::CloneCopy(*lclV);
1951 MVT::MvTimesMatAddMv(ONE,*tempMV,*gamma,ZERO,*lclV);
1954 "Anasazi::TraceMinBase::initialize(): Couldn't generate initial basis of full rank.");
1958 if(hasM_ && newstate.
MopV == Teuchos::null)
1960 om_->stream(
Debug) <<
"Computing MV\n";
1962 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1965 count_ApplyM_+= curDim_;
1968 lclMV = MVT::CloneViewNonConst(*MV_,dimind);
1969 OPT::Apply(*MOp_,*lclV,*lclMV);
1974 om_->stream(
Debug) <<
"Copying MV\n";
1977 "Anasazi::TraceMinBase::initialize(newstate): Vector length of MV not correct." );
1980 "Anasazi::TraceMinBase::initialize(newstate): Number of vectors in MV not correct.");
1982 if(newstate.
MopV != MV_) {
1983 if(curDim_ < MVT::GetNumberVecs(*newstate.
MopV)) {
1984 newstate.
MopV = MVT::CloneView(*newstate.
MopV,dimind);
1986 MVT::SetBlock(*newstate.
MopV,dimind,*MV_);
1988 lclMV = MVT::CloneViewNonConst(*MV_,dimind);
1993 om_->stream(
Debug) <<
"There is no MV\n";
2000 if(newstate.
KK == Teuchos::null)
2002 om_->stream(
Debug) <<
"Computing KK\n";
2005 newstate.
X = Teuchos::null;
2006 newstate.
MX = Teuchos::null;
2007 newstate.
KX = Teuchos::null;
2008 newstate.
R = Teuchos::null;
2009 newstate.
T = Teuchos::null;
2010 newstate.
RV = Teuchos::null;
2019 MVT::MvTransMv(ONE,*lclV,*lclKV,*lclKK);
2024 om_->stream(
Debug) <<
"Copying KK\n";
2028 "Anasazi::TraceMinBase::initialize(newstate): Projected matrix in new state must be as large as specified state rank.");
2032 if (newstate.
KK != KK_) {
2033 if (newstate.
KK->numRows() > curDim_ || newstate.
KK->numCols() > curDim_) {
2036 lclKK->assign(*newstate.
KK);
2041 if(newstate.
T == Teuchos::null || newstate.
RV == Teuchos::null)
2043 om_->stream(
Debug) <<
"Computing Ritz pairs\n";
2046 newstate.
X = Teuchos::null;
2047 newstate.
MX = Teuchos::null;
2048 newstate.
KX = Teuchos::null;
2049 newstate.
R = Teuchos::null;
2050 newstate.
T = Teuchos::null;
2051 newstate.
RV = Teuchos::null;
2058 om_->stream(
Debug) <<
"Copying Ritz pairs\n";
2061 std::invalid_argument,
"Anasazi::TraceMinBase::initialize(newstate): Size of T must be consistent with dimension of V.");
2064 "Anasazi::TraceMinBase::initialize(newstate): Ritz vectors in new state must be as large as specified state rank.");
2066 std::copy(newstate.
T->begin(),newstate.
T->end(),theta_.begin());
2069 if (newstate.
RV != ritzVecs_) {
2070 if (newstate.
RV->numRows() > curDim_ || newstate.
RV->numCols() > curDim_) {
2073 lclRV->assign(*newstate.
RV);
2078 if(newstate.
X == Teuchos::null)
2080 om_->stream(
Debug) <<
"Computing X\n";
2083 newstate.
MX = Teuchos::null;
2084 newstate.
KX = Teuchos::null;
2085 newstate.
R = Teuchos::null;
2092 om_->stream(
Debug) <<
"Copying X\n";
2094 if(computeAllRes_ ==
false)
2096 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.
X) != blockSize_ || MVT::GetGlobalLength(*newstate.
X) != MVT::GetGlobalLength(*X_),
2097 std::invalid_argument,
"Anasazi::TraceMinBase::initialize(newstate): Size of X must be consistent with block size and length of V.");
2099 if(newstate.
X != X_) {
2100 MVT::SetBlock(*newstate.
X,bsind,*X_);
2105 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.
X) != curDim_ || MVT::GetGlobalLength(*newstate.
X) != MVT::GetGlobalLength(*X_),
2106 std::invalid_argument,
"Anasazi::TraceMinBase::initialize(newstate): Size of X must be consistent with current dimension and length of V.");
2108 if(newstate.
X != X_) {
2109 MVT::SetBlock(*newstate.
X,dimind,*X_);
2116 if((Op_ != Teuchos::null && newstate.
KX == Teuchos::null) || (hasM_ && newstate.
MX == Teuchos::null))
2118 om_->stream(
Debug) <<
"Computing KX and MX\n";
2121 newstate.
R = Teuchos::null;
2128 om_->stream(
Debug) <<
"Copying KX and MX\n";
2129 if(Op_ != Teuchos::null)
2131 if(computeAllRes_ ==
false)
2134 std::invalid_argument,
"Anasazi::TraceMinBase::initialize(newstate): Size of KX must be consistent with block size and length of V.");
2136 if(newstate.
KX != KX_) {
2137 MVT::SetBlock(*newstate.
KX,bsind,*KX_);
2143 std::invalid_argument,
"Anasazi::TraceMinBase::initialize(newstate): Size of KX must be consistent with current dimension and length of V.");
2145 if (newstate.
KX != KX_) {
2146 MVT::SetBlock(*newstate.
KX,dimind,*KX_);
2153 if(computeAllRes_ ==
false)
2156 std::invalid_argument,
"Anasazi::TraceMinBase::initialize(newstate): Size of MX must be consistent with block size and length of V.");
2158 if (newstate.
MX != MX_) {
2159 MVT::SetBlock(*newstate.
MX,bsind,*MX_);
2165 std::invalid_argument,
"Anasazi::TraceMinBase::initialize(newstate): Size of MX must be consistent with current dimension and length of V.");
2167 if (newstate.
MX != MX_) {
2168 MVT::SetBlock(*newstate.
MX,dimind,*MX_);
2177 const int nvecs = computeAllRes_ ? curDim_ : blockSize_;
2179 RCP<MV> lclX = MVT::CloneViewNonConst(*X_, dimind2);
2180 std::vector<ScalarType> normvec(nvecs);
2181 orthman_->normMat(*lclX,normvec);
2184 for (
int i = 0; i < nvecs; ++i) {
2185 normvec[i] = ONE / normvec[i];
2187 MVT::MvScale (*lclX, normvec);
2188 if (Op_ != Teuchos::null) {
2189 RCP<MV> lclKX = MVT::CloneViewNonConst (*KX_, dimind2);
2190 MVT::MvScale (*lclKX, normvec);
2193 RCP<MV> lclMX = MVT::CloneViewNonConst (*MX_, dimind2);
2194 MVT::MvScale (*lclMX, normvec);
2198 for (
int i = 0; i < nvecs; ++i) {
2199 theta_[i] = theta_[i] * normvec[i] * normvec[i];
2204 if(newstate.
R == Teuchos::null)
2206 om_->stream(
Debug) <<
"Computing R\n";
2213 om_->stream(
Debug) <<
"Copying R\n";
2215 if(computeAllRes_ ==
false)
2218 std::invalid_argument,
"Anasazi::TraceMinBase::initialize(newstate): vector length of newstate.R not correct." );
2220 std::invalid_argument,
"Anasazi::TraceMinBase::initialize(newstate): newstate.R must have at least block size vectors." );
2221 if (newstate.
R != R_) {
2222 MVT::SetBlock(*newstate.
R,bsind,*R_);
2228 std::invalid_argument,
"Anasazi::TraceMinBase::initialize(newstate): vector length of newstate.R not correct." );
2230 std::invalid_argument,
"Anasazi::TraceMinBase::initialize(newstate): newstate.R must have at least curDim vectors." );
2231 if (newstate.
R != R_) {
2232 MVT::SetBlock(*newstate.
R,dimind,*R_);
2238 Rnorms_current_ =
false;
2239 R2norms_current_ =
false;
2247 om_->stream(
Debug) <<
"Copying Ritz shifts\n";
2252 om_->stream(
Debug) <<
"Setting Ritz shifts to 0\n";
2253 for(
size_t i=0; i<ritzShifts_.size(); i++)
2254 ritzShifts_[i] = ZERO;
2257 for(
size_t i=0; i<ritzShifts_.size(); i++)
2258 om_->stream(
Debug) <<
"Ritz shifts[" << i <<
"] = " << ritzShifts_[i] << std::endl;
2261 initialized_ =
true;
2263 if (om_->isVerbosity(
Debug ) ) {
2272 om_->print(
Debug, accuracyCheck(chk,
": after initialize()") );
2276 if (om_->isVerbosity(
Debug)) {
2277 currentStatus( om_->stream(
Debug) );
2287 template <
class ScalarType,
class MV,
class OP>
2293 template <
class ScalarType,
class MV,
class OP>
2299 TEUCHOS_TEST_FOR_EXCEPTION(blockSize < 1, std::invalid_argument,
"Anasazi::TraceMinBase::setSize(blocksize,numblocks): blocksize must be strictly positive.");
2301 if (blockSize == blockSize_ && numBlocks == numBlocks_) {
2306 blockSize_ = blockSize;
2307 numBlocks_ = numBlocks;
2314 if (X_ != Teuchos::null) {
2318 tmp = problem_->getInitVec();
2320 "Anasazi::TraceMinBase::setSize(): eigenproblem did not specify initial vectors to clone from.");
2323 TEUCHOS_TEST_FOR_EXCEPTION(numAuxVecs_+blockSize*static_cast<ptrdiff_t>(numBlocks) > MVT::GetGlobalLength(*tmp),std::invalid_argument,
2324 "Anasazi::TraceMinBase::setSize(): max subspace dimension and auxilliary subspace too large. Potentially impossible orthogonality constraints.");
2327 int newsd = blockSize_*numBlocks_;
2332 ritzShifts_.resize(blockSize_,ZERO);
2333 if(computeAllRes_ ==
false)
2336 Rnorms_.resize(blockSize_,NANVAL);
2337 R2norms_.resize(blockSize_,NANVAL);
2343 KX_ = Teuchos::null;
2344 MX_ = Teuchos::null;
2347 KV_ = Teuchos::null;
2348 MV_ = Teuchos::null;
2350 om_->print(
Debug,
" >> Allocating X_\n");
2351 X_ = MVT::Clone(*tmp,blockSize_);
2352 if(Op_ != Teuchos::null) {
2353 om_->print(
Debug,
" >> Allocating KX_\n");
2354 KX_ = MVT::Clone(*tmp,blockSize_);
2360 om_->print(
Debug,
" >> Allocating MX_\n");
2361 MX_ = MVT::Clone(*tmp,blockSize_);
2366 om_->print(
Debug,
" >> Allocating R_\n");
2367 R_ = MVT::Clone(*tmp,blockSize_);
2372 Rnorms_.resize(newsd,NANVAL);
2373 R2norms_.resize(newsd,NANVAL);
2379 KX_ = Teuchos::null;
2380 MX_ = Teuchos::null;
2383 KV_ = Teuchos::null;
2384 MV_ = Teuchos::null;
2386 om_->print(
Debug,
" >> Allocating X_\n");
2387 X_ = MVT::Clone(*tmp,newsd);
2388 if(Op_ != Teuchos::null) {
2389 om_->print(
Debug,
" >> Allocating KX_\n");
2390 KX_ = MVT::Clone(*tmp,newsd);
2396 om_->print(
Debug,
" >> Allocating MX_\n");
2397 MX_ = MVT::Clone(*tmp,newsd);
2402 om_->print(
Debug,
" >> Allocating R_\n");
2403 R_ = MVT::Clone(*tmp,newsd);
2409 theta_.resize(newsd,NANVAL);
2410 om_->print(
Debug,
" >> Allocating V_\n");
2411 V_ = MVT::Clone(*tmp,newsd);
2415 if(Op_ != Teuchos::null) {
2416 om_->print(
Debug,
" >> Allocating KV_\n");
2417 KV_ = MVT::Clone(*tmp,newsd);
2423 om_->print(
Debug,
" >> Allocating MV_\n");
2424 MV_ = MVT::Clone(*tmp,newsd);
2430 om_->print(
Debug,
" >> done allocating.\n");
2432 initialized_ =
false;
2439 template <
class ScalarType,
class MV,
class OP>
2450 for (tarcpmv i=auxVecs_.begin(); i != auxVecs_.end(); ++i) {
2451 numAuxVecs_ += MVT::GetNumberVecs(**i);
2455 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 2458 count_ApplyM_+= MVT::GetNumberVecs(**i);
2460 RCP<MV> helperMV = MVT::Clone(**i,MVT::GetNumberVecs(**i));
2461 OPT::Apply(*MOp_,**i,*helperMV);
2462 MauxVecs_.push_back(helperMV);
2467 if (numAuxVecs_ > 0 && initialized_) {
2468 initialized_ =
false;
2471 if (om_->isVerbosity(
Debug ) ) {
2475 om_->print(
Debug, accuracyCheck(chk,
": after setAuxVecs()") );
2482 template <
class ScalarType,
class MV,
class OP>
2483 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
2485 if (Rnorms_current_ ==
false) {
2489 std::vector<int> curind(curDim_);
2490 for(
int i=0; i<curDim_; i++)
2494 std::vector<ScalarType> locNorms(curDim_);
2495 orthman_->norm(*locR,locNorms);
2497 for(
int i=0; i<curDim_; i++)
2498 Rnorms_[i] = locNorms[i];
2499 for(
int i=curDim_+1; i<blockSize_*numBlocks_; i++)
2500 Rnorms_[i] = NANVAL;
2502 Rnorms_current_ =
true;
2503 locNorms.resize(blockSize_);
2507 orthman_->norm(*R_,Rnorms_);
2508 Rnorms_current_ =
true;
2510 else if(computeAllRes_)
2512 std::vector<ScalarType> locNorms(blockSize_);
2513 for(
int i=0; i<blockSize_; i++)
2514 locNorms[i] = Rnorms_[i];
2524 template <
class ScalarType,
class MV,
class OP>
2525 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
2527 if (R2norms_current_ ==
false) {
2531 std::vector<int> curind(curDim_);
2532 for(
int i=0; i<curDim_; i++)
2536 std::vector<ScalarType> locNorms(curDim_);
2537 MVT::MvNorm(*locR,locNorms);
2539 for(
int i=0; i<curDim_; i++)
2541 R2norms_[i] = locNorms[i];
2543 for(
int i=curDim_+1; i<blockSize_*numBlocks_; i++)
2544 R2norms_[i] = NANVAL;
2546 R2norms_current_ =
true;
2547 locNorms.resize(blockSize_);
2551 MVT::MvNorm(*R_,R2norms_);
2552 R2norms_current_ =
true;
2554 else if(computeAllRes_)
2556 std::vector<ScalarType> locNorms(blockSize_);
2557 for(
int i=0; i<blockSize_; i++)
2558 locNorms[i] = R2norms_[i];
2567 template <
class ScalarType,
class MV,
class OP>
2570 "Anasazi::TraceMinBase::setStatusTest() was passed a null StatusTest.");
2576 template <
class ScalarType,
class MV,
class OP>
2583 template <
class ScalarType,
class MV,
class OP>
2589 os.setf(std::ios::scientific, std::ios::floatfield);
2592 os <<
"================================================================================" << endl;
2594 os <<
" TraceMinBase Solver Status" << endl;
2596 os <<
"The solver is "<<(initialized_ ?
"initialized." :
"not initialized.") << endl;
2597 os <<
"The number of iterations performed is " <<iter_<<endl;
2598 os <<
"The block size is " << blockSize_<<endl;
2599 os <<
"The number of blocks is " << numBlocks_<<endl;
2600 os <<
"The current basis size is " << curDim_<<endl;
2601 os <<
"The number of auxiliary vectors is "<< numAuxVecs_ << endl;
2602 os <<
"The number of operations Op*x is "<<count_ApplyOp_<<endl;
2603 os <<
"The number of operations M*x is "<<count_ApplyM_<<endl;
2605 os.setf(std::ios_base::right, std::ios_base::adjustfield);
2609 os <<
"CURRENT EIGENVALUE ESTIMATES "<<endl;
2610 os << std::setw(20) <<
"Eigenvalue" 2611 << std::setw(20) <<
"Residual(M)" 2612 << std::setw(20) <<
"Residual(2)" 2614 os <<
"--------------------------------------------------------------------------------"<<endl;
2615 for (
int i=0; i<blockSize_; ++i) {
2616 os << std::setw(20) << theta_[i];
2617 if (Rnorms_current_) os << std::setw(20) << Rnorms_[i];
2618 else os << std::setw(20) <<
"not current";
2619 if (R2norms_current_) os << std::setw(20) << R2norms_[i];
2620 else os << std::setw(20) <<
"not current";
2624 os <<
"================================================================================" << endl;
2629 template <
class ScalarType,
class MV,
class OP>
2632 ScalarType currentTrace = ZERO;
2634 for(
int i=0; i < blockSize_; i++)
2635 currentTrace += theta_[i];
2637 return currentTrace;
2641 template <
class ScalarType,
class MV,
class OP>
2642 bool TraceMinBase<ScalarType,MV,OP>::traceLeveled()
2644 ScalarType ratioOfChange = traceThresh_;
2646 if(previouslyLeveled_)
2648 om_->stream(
Debug) <<
"The trace already leveled, so we're not going to check it again\n";
2652 ScalarType currentTrace = getTrace();
2654 om_->stream(
Debug) <<
"The current trace is " << currentTrace << std::endl;
2659 if(previousTrace_ != ZERO)
2661 om_->stream(
Debug) <<
"The previous trace was " << previousTrace_ << std::endl;
2663 ratioOfChange = std::abs(previousTrace_-currentTrace)/std::abs(previousTrace_);
2664 om_->stream(
Debug) <<
"The ratio of change is " << ratioOfChange << std::endl;
2667 previousTrace_ = currentTrace;
2669 if(ratioOfChange < traceThresh_)
2671 previouslyLeveled_ =
true;
2680 template <
class ScalarType,
class MV,
class OP>
2681 std::vector<ScalarType> TraceMinBase<ScalarType,MV,OP>::getClusterResids()
2692 std::vector<ScalarType> clusterResids(nvecs);
2693 std::vector<int> clusterIndices;
2694 if(considerClusters_)
2696 for(
int i=0; i < nvecs; i++)
2699 if(clusterIndices.empty() || (theta_[i-1] + R2norms_[i-1] >= theta_[i] - R2norms_[i]))
2702 if(!clusterIndices.empty()) om_->stream(
Debug) << theta_[i-1] <<
" is in a cluster with " << theta_[i] <<
" because " << theta_[i-1] + R2norms_[i-1] <<
" >= " << theta_[i] - R2norms_[i] << std::endl;
2703 clusterIndices.push_back(i);
2708 om_->stream(
Debug) << theta_[i-1] <<
" is NOT in a cluster with " << theta_[i] <<
" because " << theta_[i-1] + R2norms_[i-1] <<
" < " << theta_[i] - R2norms_[i] << std::endl;
2709 ScalarType totalRes = ZERO;
2710 for(
size_t j=0; j < clusterIndices.size(); j++)
2711 totalRes += (R2norms_[clusterIndices[j]]*R2norms_[clusterIndices[j]]);
2716 if(theta_[clusterIndices[0]] < 0 && theta_[i] < 0)
2717 negSafeToShift_ =
true;
2718 else if(theta_[clusterIndices[0]] > 0 && theta_[i] > 0)
2719 posSafeToShift_ =
true;
2721 for(
size_t j=0; j < clusterIndices.size(); j++)
2722 clusterResids[clusterIndices[j]] = sqrt(totalRes);
2724 clusterIndices.clear();
2725 clusterIndices.push_back(i);
2730 ScalarType totalRes = ZERO;
2731 for(
size_t j=0; j < clusterIndices.size(); j++)
2732 totalRes += R2norms_[clusterIndices[j]];
2733 for(
size_t j=0; j < clusterIndices.size(); j++)
2734 clusterResids[clusterIndices[j]] = totalRes;
2738 for(
int j=0; j < nvecs; j++)
2739 clusterResids[j] = R2norms_[j];
2742 return clusterResids;
2751 template <
class ScalarType,
class MV,
class OP>
2752 void TraceMinBase<ScalarType,MV,OP>::computeRitzShifts(
const std::vector<ScalarType>& clusterResids)
2754 std::vector<ScalarType> thetaMag(theta_);
2755 bool traceHasLeveled = traceLeveled();
2758 for(
int i=0; i<blockSize_; i++)
2759 thetaMag[i] = std::abs(thetaMag[i]);
2767 if(whenToShift_ == ALWAYS_SHIFT || (whenToShift_ == SHIFT_WHEN_TRACE_LEVELS && traceHasLeveled))
2770 if(howToShift_ == LARGEST_CONVERGED_SHIFT)
2772 for(
int i=0; i<blockSize_; i++)
2773 ritzShifts_[i] = largestSafeShift_;
2776 else if(howToShift_ == RITZ_VALUES_SHIFT)
2778 ritzShifts_[0] = theta_[0];
2782 if(useMultipleShifts_)
2784 for(
int i=1; i<blockSize_; i++)
2785 ritzShifts_[i] = theta_[i];
2789 for(
int i=1; i<blockSize_; i++)
2790 ritzShifts_[i] = ritzShifts_[0];
2793 else if(howToShift_ == EXPERIMENTAL_SHIFT)
2795 ritzShifts_[0] = std::max(largestSafeShift_,theta_[0]-clusterResids[0]);
2796 for(
int i=1; i<blockSize_; i++)
2798 ritzShifts_[i] = std::max(ritzShifts_[i-1],theta_[i]-clusterResids[i]);
2802 else if(howToShift_ == ADJUSTED_RITZ_SHIFT)
2804 om_->stream(
Debug) <<
"\nSeeking a shift for theta[0]=" << thetaMag[0] << std::endl;
2808 if((theta_[0] > 0 && posSafeToShift_) || (theta_[0] < 0 && negSafeToShift_) || considerClusters_ ==
false)
2811 ritzShifts_[0] = std::max(largestSafeShift_,thetaMag[0]-clusterResids[0]);
2813 om_->stream(
Debug) <<
"Initializing with a conservative shift, either the most positive converged eigenvalue (" 2814 << largestSafeShift_ <<
") or the eigenvalue adjusted by the residual (" << thetaMag[0] <<
"-" 2815 << clusterResids[0] <<
").\n";
2818 if(R2norms_[0] == clusterResids[0])
2820 ritzShifts_[0] = thetaMag[0];
2821 om_->stream(
Debug) <<
"Since this eigenvalue is NOT in a cluster, we can use the eigenvalue itself as a shift: ritzShifts[0]=" << ritzShifts_[0] << std::endl;
2824 om_->stream(
Debug) <<
"This eigenvalue is in a cluster, so it would not be safe to use the eigenvalue itself as a shift\n";
2828 if(largestSafeShift_ > std::abs(ritzShifts_[0]))
2830 om_->stream(
Debug) <<
"Initializing with a conservative shift...the most positive converged eigenvalue: " << largestSafeShift_ << std::endl;
2831 ritzShifts_[0] = largestSafeShift_;
2834 om_->stream(
Debug) <<
"Using the previous value of ritzShifts[0]=" << ritzShifts_[0];
2838 om_->stream(
Debug) <<
"ritzShifts[0]=" << ritzShifts_[0] << std::endl;
2840 if(useMultipleShifts_)
2844 for(
int i=1; i < blockSize_; i++)
2846 om_->stream(
Debug) <<
"\nSeeking a shift for theta[" << i <<
"]=" << thetaMag[i] << std::endl;
2849 if(ritzShifts_[i-1] == thetaMag[i-1] && i < blockSize_-1 && thetaMag[i] < thetaMag[i+1] - clusterResids[i+1])
2851 ritzShifts_[i] = thetaMag[i];
2852 om_->stream(
Debug) <<
"Using an aggressive shift: ritzShifts_[" << i <<
"]=" << ritzShifts_[i] << std::endl;
2856 if(ritzShifts_[0] > std::abs(ritzShifts_[i]))
2858 om_->stream(
Debug) <<
"It was unsafe to use the aggressive shift. Choose the shift used by theta[0]=" 2859 << thetaMag[0] <<
": ritzShifts[0]=" << ritzShifts_[0] << std::endl;
2862 ritzShifts_[i] = ritzShifts_[0];
2865 om_->stream(
Debug) <<
"It was unsafe to use the aggressive shift. We will use the shift from the previous iteration: " << ritzShifts_[i] << std::endl;
2867 om_->stream(
Debug) <<
"Check whether any less conservative shifts would work (such as the biggest eigenvalue outside of the cluster, namely theta[ell] < " 2868 << thetaMag[i] <<
"-" << clusterResids[i] <<
" (" << thetaMag[i] - clusterResids[i] <<
")\n";
2871 for(
int ell=0; ell < i; ell++)
2873 if(thetaMag[ell] < thetaMag[i] - clusterResids[i])
2875 ritzShifts_[i] = thetaMag[ell];
2876 om_->stream(
Debug) <<
"ritzShifts_[" << i <<
"]=" << ritzShifts_[i] <<
" is valid\n";
2883 om_->stream(
Debug) <<
"ritzShifts[" << i <<
"]=" << ritzShifts_[i] << std::endl;
2888 for(
int i=1; i<blockSize_; i++)
2889 ritzShifts_[i] = ritzShifts_[0];
2895 for(
int i=0; i<blockSize_; i++)
2898 ritzShifts_[i] = -abs(ritzShifts_[i]);
2900 ritzShifts_[i] = abs(ritzShifts_[i]);
2905 template <
class ScalarType,
class MV,
class OP>
2906 std::vector<ScalarType> TraceMinBase<ScalarType,MV,OP>::computeTol()
2909 std::vector<ScalarType> tolerances(blockSize_);
2911 for(
int i=0; i < blockSize_-1; i++)
2913 if(std::abs(theta_[blockSize_-1]) != std::abs(ritzShifts_[i]))
2914 temp1 = std::abs(theta_[i]-ritzShifts_[i])/std::abs(std::abs(theta_[blockSize_-1])-std::abs(ritzShifts_[i]));
2920 tolerances[i] = std::min(temp1*temp1,0.5);
2924 tolerances[blockSize_-1] = tolerances[blockSize_-2];
2930 template <
class ScalarType,
class MV,
class OP>
2931 void TraceMinBase<ScalarType,MV,OP>::solveSaddlePointProblem(
RCP<MV> Delta)
2933 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 2938 if(Op_ == Teuchos::null)
2949 std::vector<int> curind(blockSize_);
2950 for(
int i=0; i<blockSize_; i++)
2957 MVT::MvTransMv(ONE,*lclMX,*lclMX,*lclS);
2965 MVT::Assign(*lclX,*Delta);
2966 MVT::MvTimesMatAddMv( -ONE, *lclMX, *lclS, ONE, *Delta);
2971 MVT::MvTransMv(ONE,*MX_,*MX_,*lclS);
2978 MVT::Assign(*X_,*Delta);
2979 MVT::MvTimesMatAddMv( -ONE, *MX_, *lclS, ONE, *Delta);
2984 std::vector<int> order(curDim_);
2985 std::vector<ScalarType> tempvec(blockSize_);
2989 std::vector<ScalarType> clusterResids;
3002 clusterResids = getClusterResids();
3017 computeRitzShifts(clusterResids);
3020 std::vector<ScalarType> tolerances = computeTol();
3022 for(
int i=0; i<blockSize_; i++)
3024 om_->stream(
IterationDetails) <<
"Choosing Ritz shifts...theta[" << i <<
"]=" 3025 << theta_[i] <<
", resids[" << i <<
"]=" << R2norms_[i] <<
", clusterResids[" << i <<
"]=" << clusterResids[i]
3026 <<
", ritzShifts[" << i <<
"]=" << ritzShifts_[i] <<
", and tol[" << i <<
"]=" << tolerances[i] << std::endl;
3030 ritzOp_->setRitzShifts(ritzShifts_);
3034 ritzOp_->setInnerTol(tolerances);
3037 if(saddleSolType_ == PROJECTED_KRYLOV_SOLVER)
3039 if(Prec_ != Teuchos::null)
3040 solveSaddleProjPrec(Delta);
3042 solveSaddleProj(Delta);
3044 else if(saddleSolType_ == SCHUR_COMPLEMENT_SOLVER)
3046 if(Z_ == Teuchos::null || MVT::GetNumberVecs(*Z_) != blockSize_)
3050 Z_ = MVT::Clone(*X_,blockSize_);
3053 solveSaddleSchur(Delta);
3055 else if(saddleSolType_ == BD_PREC_MINRES)
3057 solveSaddleBDPrec(Delta);
3060 else if(saddleSolType_ == HSS_PREC_GMRES)
3062 solveSaddleHSSPrec(Delta);
3065 std::cout <<
"Invalid saddle solver type\n";
3072 template <
class ScalarType,
class MV,
class OP>
3073 void TraceMinBase<ScalarType,MV,OP>::solveSaddleProj (
RCP<MV> Delta)
const 3080 std::vector<int> curind(blockSize_);
3081 for(
int i=0; i<blockSize_; i++)
3089 if(projectLockedVecs_ && numAuxVecs_ > 0)
3090 projOp =
rcp(
new TraceMinProjRitzOp<ScalarType,MV,OP>(ritzOp_,locMX,orthman_,auxVecs_) );
3092 projOp =
rcp(
new TraceMinProjRitzOp<ScalarType,MV,OP>(ritzOp_,locMX,orthman_) );
3095 projOp =
rcp(
new TraceMinProjRitzOp<ScalarType,MV,OP>(ritzOp_,locMX) );
3099 MVT::MvInit(*Delta);
3104 projOp->ApplyInverse(*locR, *Delta);
3109 projOp->ApplyInverse(*locKX, *Delta);
3117 if(projectLockedVecs_ && numAuxVecs_ > 0)
3118 projOp =
rcp(
new TraceMinProjRitzOp<ScalarType,MV,OP>(ritzOp_,MX_,orthman_,auxVecs_) );
3120 projOp =
rcp(
new TraceMinProjRitzOp<ScalarType,MV,OP>(ritzOp_,MX_,orthman_) );
3123 projOp =
rcp(
new TraceMinProjRitzOp<ScalarType,MV,OP>(ritzOp_,MX_) );
3127 MVT::MvInit(*Delta);
3130 projOp->ApplyInverse(*R_, *Delta);
3133 projOp->ApplyInverse(*KX_, *Delta);
3140 template <
class ScalarType,
class MV,
class OP>
3141 void TraceMinBase<ScalarType,MV,OP>::solveSaddleProjPrec (
RCP<MV> Delta)
const 3146 #ifdef HAVE_ANASAZI_BELOS 3153 dimension = curDim_;
3155 dimension = blockSize_;
3158 std::vector<int> curind(dimension);
3159 for(
int i=0; i<dimension; i++)
3167 if(projectLockedVecs_ && numAuxVecs_ > 0)
3168 projOp =
rcp(
new TraceMinProjRitzOpWithPrec<ScalarType,MV,OP>(ritzOp_,locMX,orthman_,auxVecs_) );
3170 projOp =
rcp(
new TraceMinProjRitzOpWithPrec<ScalarType,MV,OP>(ritzOp_,locMX,orthman_) );
3173 projOp =
rcp(
new TraceMinProjRitzOpWithPrec<ScalarType,MV,OP>(ritzOp_,locMX) );
3177 MVT::MvInit(*Delta);
3179 std::vector<int> dimind(blockSize_);
3180 for(
int i=0; i<blockSize_; i++)
3186 projOp->ApplyInverse(*locR, *Delta);
3187 MVT::MvScale(*Delta, -ONE);
3192 projOp->ApplyInverse(*locKX, *Delta);
3200 if(projectLockedVecs_ && numAuxVecs_ > 0)
3201 projOp =
rcp(
new TraceMinProjRitzOpWithPrec<ScalarType,MV,OP>(ritzOp_,MX_,orthman_,auxVecs_) );
3203 projOp =
rcp(
new TraceMinProjRitzOpWithPrec<ScalarType,MV,OP>(ritzOp_,MX_,orthman_) );
3206 projOp =
rcp(
new TraceMinProjRitzOpWithPrec<ScalarType,MV,OP>(ritzOp_,MX_) );
3210 MVT::MvInit(*Delta);
3214 projOp->ApplyInverse(*R_, *Delta);
3215 MVT::MvScale(*Delta,-ONE);
3218 projOp->ApplyInverse(*KX_, *Delta);
3226 template <
class ScalarType,
class MV,
class OP>
3227 void TraceMinBase<ScalarType,MV,OP>::solveSaddleSchur (
RCP<MV> Delta)
const 3238 std::vector<int> curind(blockSize_);
3239 for(
int i=0; i<blockSize_; i++)
3246 #ifdef USE_APPLY_INVERSE 3247 Op_->ApplyInverse(*lclMX,*Z_);
3249 ritzOp_->ApplyInverse(*lclMX,*Z_);
3253 MVT::MvTransMv(ONE,*Z_,*lclMX,*lclS);
3262 MVT::Assign(*lclX,*Delta);
3263 MVT::MvTimesMatAddMv( -ONE, *Z_, *lclL, ONE, *Delta);
3268 #ifdef USE_APPLY_INVERSE 3269 Op_->ApplyInverse(*MX_,*Z_);
3271 ritzOp_->ApplyInverse(*MX_,*Z_);
3275 MVT::MvTransMv(ONE,*Z_,*MX_,*lclS);
3283 MVT::Assign(*X_,*Delta);
3284 MVT::MvTimesMatAddMv( -ONE, *Z_, *lclL, ONE, *Delta);
3291 template <
class ScalarType,
class MV,
class OP>
3292 void TraceMinBase<ScalarType,MV,OP>::solveSaddleBDPrec (
RCP<MV> Delta)
const 3297 std::vector<int> curind(blockSize_);
3298 for(
int i=0; i<blockSize_; i++)
3300 locKX = MVT::CloneViewNonConst(*KX_, curind);
3301 locMX = MVT::CloneViewNonConst(*MX_, curind);
3320 MVT::MvInit(*Delta);
3325 if(Prec_ != Teuchos::null)
3328 sadSolver =
rcp(
new PseudoBlockMinres<ScalarType,saddle_container_type,saddle_op_type>(sadOp, sadPrec));
3331 sadSolver =
rcp(
new PseudoBlockMinres<ScalarType,saddle_container_type,saddle_op_type>(sadOp));
3335 std::vector<ScalarType> tol;
3336 ritzOp_->getInnerTol(tol);
3337 sadSolver->setTol(tol);
3340 sadSolver->setMaxIter(ritzOp_->getMaxIts());
3343 sadSolver->setSol(sadSol);
3346 sadSolver->setRHS(sadRHS);
3356 template <
class ScalarType,
class MV,
class OP>
3357 void TraceMinBase<ScalarType,MV,OP>::solveSaddleHSSPrec (
RCP<MV> Delta)
const 3359 #ifdef HAVE_ANASAZI_BELOS 3360 typedef Belos::LinearProblem<ScalarType,saddle_container_type,saddle_op_type> LP;
3361 typedef Belos::PseudoBlockGmresSolMgr<ScalarType,saddle_container_type,saddle_op_type> GmresSolMgr;
3366 std::vector<int> curind(blockSize_);
3367 for(
int i=0; i<blockSize_; i++)
3369 locKX = MVT::CloneViewNonConst(*KX_, curind);
3370 locMX = MVT::CloneViewNonConst(*MX_, curind);
3385 MVT::MvInit(*Delta);
3392 std::vector<ScalarType> tol;
3393 ritzOp_->getInnerTol(tol);
3394 pl->
set(
"Convergence Tolerance",tol[0]);
3398 pl->
set(
"Maximum Iterations", ritzOp_->getMaxIts());
3399 pl->
set(
"Num Blocks", ritzOp_->getMaxIts());
3404 pl->
set(
"Block Size", blockSize_);
3411 RCP<LP> problem =
rcp(
new LP(sadOp,sadSol,sadRHS));
3414 if(Prec_ != Teuchos::null)
3417 problem->setLeftPrec(sadPrec);
3421 problem->setProblem();
3429 std::cout <<
"No Belos. This is bad\n";
3438 template <
class ScalarType,
class MV,
class OP>
3439 void TraceMinBase<ScalarType,MV,OP>::computeKK()
3442 std::vector<int> curind(curDim_);
3443 for(
int i=0; i<curDim_; i++)
3450 curind.resize(blockSize_);
3451 for(
int i=0; i<blockSize_; i++)
3452 curind[i] = curDim_-blockSize_+i;
3460 MVT::MvTransMv(ONE,*lclV,*lclKV,*lclKK);
3463 for(
int r=0; r<curDim_; r++)
3465 for(
int c=0; c<r; c++)
3467 (*KK_)(r,c) = (*KK_)(c,r);
3474 template <
class ScalarType,
class MV,
class OP>
3475 void TraceMinBase<ScalarType,MV,OP>::computeRitzPairs()
3487 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 3491 Utils::directSolver(curDim_, *lclKK, Teuchos::null, *lclRV, theta_, rank, 10);
3495 "Anasazi::TraceMinBase::computeRitzPairs(): Failed to compute all eigenpairs of KK.");
3500 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 3504 std::vector<int> order(curDim_);
3510 sm.
sort(theta_, Teuchos::rcpFromRef(order), curDim_);
3514 sm_->sort(theta_, Teuchos::rcpFromRef(order), curDim_);
3518 Utils::permuteVectors(order,*lclRV);
3524 template <
class ScalarType,
class MV,
class OP>
3525 void TraceMinBase<ScalarType,MV,OP>::computeX()
3527 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 3532 std::vector<int> curind(curDim_);
3533 for(
int i=0; i<curDim_; i++)
3546 RCP<MV> lclX = MVT::CloneViewNonConst(*X_,curind);
3547 MVT::MvTimesMatAddMv( ONE, *lclV, *relevantEvecs, ZERO, *lclX );
3556 MVT::MvTimesMatAddMv( ONE, *lclV, *relevantEvecs, ZERO, *X_ );
3562 template <
class ScalarType,
class MV,
class OP>
3563 void TraceMinBase<ScalarType,MV,OP>::updateKXMX()
3565 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 3570 std::vector<int> curind(curDim_);
3571 for(
int i=0; i<curDim_; i++)
3585 RCP<MV> lclKX = MVT::CloneViewNonConst(*KX_,curind);
3586 MVT::MvTimesMatAddMv( ONE, *lclKV, *relevantEvecs, ZERO, *lclKX );
3590 RCP<MV> lclMX = MVT::CloneViewNonConst(*MX_,curind);
3591 MVT::MvTimesMatAddMv( ONE, *lclMV, *relevantEvecs, ZERO, *lclMX );
3601 MVT::MvTimesMatAddMv( ONE, *lclKV, *relevantEvecs, ZERO, *KX_ );
3605 MVT::MvTimesMatAddMv( ONE, *lclMV, *relevantEvecs, ZERO, *MX_ );
3612 template <
class ScalarType,
class MV,
class OP>
3613 void TraceMinBase<ScalarType,MV,OP>::updateResidual () {
3614 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 3621 std::vector<int> curind(curDim_);
3622 for(
int i=0; i<curDim_; i++)
3626 RCP<MV> MXT = MVT::CloneCopy(*MX_,curind);
3629 std::vector<ScalarType> locTheta(curDim_);
3630 for(
int i=0; i<curDim_; i++)
3631 locTheta[i] = theta_[i];
3634 MVT::MvScale(*MXT,locTheta);
3638 RCP<MV> locR = MVT::CloneViewNonConst(*R_,curind);
3639 MVT::MvAddMv(ONE,*locKX,-ONE,*MXT,*locR);
3644 RCP<MV> MXT = MVT::CloneCopy(*MX_);
3647 std::vector<ScalarType> locTheta(blockSize_);
3648 for(
int i=0; i<blockSize_; i++)
3649 locTheta[i] = theta_[i];
3652 MVT::MvScale(*MXT,locTheta);
3655 MVT::MvAddMv(ONE,*KX_,-ONE,*MXT,*R_);
3659 Rnorms_current_ =
false;
3660 R2norms_current_ =
false;
3690 template <
class ScalarType,
class MV,
class OP>
3691 std::string TraceMinBase<ScalarType,MV,OP>::accuracyCheck(
const CheckList &chk,
const std::string &where )
const 3695 std::stringstream os;
3697 os.setf(std::ios::scientific, std::ios::floatfield);
3699 os <<
" Debugging checks: iteration " << iter_ << where << endl;
3702 std::vector<int> lclind(curDim_);
3703 for (
int i=0; i<curDim_; ++i) lclind[i] = i;
3706 lclV = MVT::CloneView(*V_,lclind);
3708 if (chk.checkV && initialized_) {
3709 MagnitudeType err = orthman_->orthonormError(*lclV);
3710 os <<
" >> Error in V^H M V == I : " << err << endl;
3711 for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
3712 err = orthman_->orthogError(*lclV,*auxVecs_[i]);
3713 os <<
" >> Error in V^H M Q[" << i <<
"] == 0 : " << err << endl;
3722 lclX = MVT::CloneView(*X_,lclind);
3727 if (chk.checkX && initialized_) {
3728 MagnitudeType err = orthman_->orthonormError(*lclX);
3729 os <<
" >> Error in X^H M X == I : " << err << endl;
3730 for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
3731 err = orthman_->orthogError(*lclX,*auxVecs_[i]);
3732 os <<
" >> Error in X^H M Q[" << i <<
"] == 0 : " << err << endl;
3735 if (chk.checkMX && hasM_ && initialized_) {
3738 lclMX = MVT::CloneView(*MX_,lclind);
3742 MagnitudeType err = Utils::errorEquality(*lclX, *lclMX, MOp_);
3743 os <<
" >> Error in MX == M*X : " << err << endl;
3745 if (Op_ != Teuchos::null && chk.checkKX && initialized_) {
3748 lclKX = MVT::CloneView(*KX_,lclind);
3752 MagnitudeType err = Utils::errorEquality(*lclX, *lclKX, Op_);
3753 os <<
" >> Error in KX == K*X : " << err << endl;
3757 if (chk.checkKK && initialized_) {
3759 if(Op_ != Teuchos::null) {
3760 RCP<MV> lclKV = MVT::Clone(*V_,curDim_);
3761 OPT::Apply(*Op_,*lclV,*lclKV);
3762 MVT::MvTransMv(ONE,*lclV,*lclKV,curKK);
3765 MVT::MvTransMv(ONE,*lclV,*lclV,curKK);
3769 os <<
" >> Error in V^H K V == KK : " << curKK.normFrobenius() << endl;
3772 for (
int j=0; j<curDim_; ++j) {
3773 for (
int i=0; i<curDim_; ++i) {
3774 SDMerr(i,j) = subKK(i,j) - SCT::conjugate(subKK(j,i));
3777 os <<
" >> Error in KK - KK^H == 0 : " << SDMerr.normFrobenius() << endl;
3782 for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
3783 MagnitudeType err = orthman_->orthonormError(*auxVecs_[i]);
3784 os <<
" >> Error in Q[" << i <<
"]^H M Q[" << i <<
"] == I : " << err << endl;
3785 for (Array_size_type j=i+1; j<auxVecs_.size(); ++j) {
3786 err = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
3787 os <<
" >> Error in Q[" << i <<
"]^H M Q[" << j <<
"] == 0 : " << err << endl;
TraceMinBaseState< ScalarType, MV > getState() const
Get access to the current state of the eigensolver.
RCP< const MV > getRitzVectors()
Get access to the current Ritz vectors.
Teuchos::Array< RCP< const MV > > getAuxVecs() const
Get the auxiliary vectors for the solver.
int NEV
Number of unconverged eigenvalues.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
RCP< const MV > X
The current eigenvectors.
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace. For the trace minimization methods...
This class defines the interface required by an eigensolver and status test class to compute solution...
Structure to contain pointers to TraceMinBase state variables.
An implementation of the Anasazi::SortManager that performs a collection of common sorting techniques...
void setBlockSize(int blockSize)
Set the blocksize.
void sort(std::vector< MagnitudeType > &evals, Teuchos::RCP< std::vector< int > > perm=Teuchos::null, int n=-1) const
Sort real eigenvalues, optionally returning the permutation vector.
RCP< const MV > V
The current basis.
Declaration of basic traits for the multivector type.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getResNorms()
Get the current residual norms, computing the norms if they are not up-to-date with the current resid...
T & get(const std::string &name, T def_value)
An operator of the form [A Y; Y' 0] where A is a sparse matrix and Y a multivector.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
RCP< const MV > MopV
The image of V under M, or Teuchos::null if M was not specified.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
TraceMinBaseOrthoFailure is thrown when the orthogonalization manager is unable to orthogonalize the ...
Virtual base class which defines basic traits for the operator type.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRes2Norms()
Get the current residual 2-norms, computing the norms if they are not up-to-date with the current res...
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current eigenvectors and eigenvalues...
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRitzRes2Norms()
Get the 2-norms of the residuals.
Basic implementation of the Anasazi::SortManager class.
TraceMinBaseInitFailure is thrown when the TraceMinBase solver is unable to generate an initial itera...
An exception class parent to all Anasazi exceptions.
std::vector< int > getRitzIndex()
Get the index used for extracting individual Ritz vectors from getRitzVectors().
void setStatusTest(RCP< StatusTest< ScalarType, MV, OP > > test)
Set a new StatusTest for the solver.
void initialize()
Initialize the solver with the initial vectors from the eigenproblem or random data.
void currentStatus(std::ostream &os)
This method requests that the solver print out its current status to the given output stream...
bool isOrtho
Whether V has been projected and orthonormalized already.
Pure virtual base class which describes the basic interface to the iterative eigensolver.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
Anasazi's templated, static class providing utilities for the solvers.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the eigenvalue problem.
Anasazi's templated virtual class for providing routines for orthogonalization and orthonormalization...
int getNumIters() const
Get the current iteration count.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual ~TraceMinBase()
Anasazi::TraceMinBase destructor.
Output managers remove the need for the eigensolver to know any information about the required output...
Traits class which defines basic operations on multivectors.
void iterate()
This method performs trace minimization iterations until the status test indicates the need to stop o...
Virtual base class which defines basic traits for the operator type.
int curDim
The current dimension of the solver.
RCP< const MV > KV
The image of V under K.
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
void resize(size_type new_size, const value_type &x=value_type())
RCP< StatusTest< ScalarType, MV, OP > > getStatusTest() const
Get the current StatusTest used by the solver.
RCP< const MV > MX
The image of the current eigenvectors under M, or Teuchos::null if M was not specified.
RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > RV
The current Ritz vectors.
ScalarType largestSafeShift
Largest safe shift.
RCP< const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > T
The current Ritz values. This vector is a copy of the internal data.
std::vector< Value< ScalarType > > getRitzValues()
Get the Ritz values for the previous iteration.
RCP< const std::vector< ScalarType > > ritzShifts
Current Ritz shifts.
RCP< const MV > KX
The image of the current eigenvectors under K.
RCP< const MV > R
The current residual vectors.
RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > KK
The current projected K matrix.
Types and exceptions used within Anasazi solvers and interfaces.
void setSize(int blockSize, int numBlocks)
Set the blocksize and number of blocks to be used by the iterative solver in solving this eigenproble...
bool isInitialized() const
Indicates whether the solver has been initialized or not.
This is an abstract base class for the trace minimization eigensolvers.
Anasazi's templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
Stores a set of vectors of the form [V; L] where V is a distributed multivector and L is a serialdens...
Common interface of stopping criteria for Anasazi's solvers.
void setAuxVecs(const Teuchos::Array< RCP< const MV > > &auxvecs)
Set the auxiliary vectors for the solver.
The Eigensolver is a templated virtual base class that defines the basic interface that any eigensolv...
void resetNumIters()
Reset the iteration count.
int getBlockSize() const
Get the blocksize used by the iterative solver.
TraceMinBase(const RCP< Eigenproblem< ScalarType, MV, OP > > &problem, const RCP< SortManager< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > &sorter, const RCP< OutputManager< ScalarType > > &printer, const RCP< StatusTest< ScalarType, MV, OP > > &tester, const RCP< MatOrthoManager< ScalarType, MV, OP > > &ortho, Teuchos::ParameterList ¶ms)
TraceMinBase constructor with eigenproblem, solver utilities, and parameter list of solver options...
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)
Class which provides internal utilities for the Anasazi solvers.