35 #ifndef ANASAZI_TRACEMIN_BASE_HPP 36 #define ANASAZI_TRACEMIN_BASE_HPP 51 #include "Teuchos_ParameterList.hpp" 52 #include "Teuchos_ScalarTraits.hpp" 53 #include "Teuchos_SerialDenseMatrix.hpp" 54 #include "Teuchos_SerialDenseSolver.hpp" 55 #include "Teuchos_TimeMonitor.hpp" 80 template <
class ScalarType,
class MV>
102 RCP<const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> >
T;
108 RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> >
KK;
110 RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> >
RV;
119 TraceMinBaseState() : curDim(0), V(Teuchos::null), KV(Teuchos::null), MopV(Teuchos::null),
120 X(Teuchos::null), KX(Teuchos::null), MX(Teuchos::null), R(Teuchos::null),
121 T(Teuchos::null), KK(Teuchos::null), RV(Teuchos::null), isOrtho(false),
122 NEV(0), largestSafeShift(Teuchos::ScalarTraits<ScalarType>::zero()),
123 ritzShifts(Teuchos::null) {}
168 template <
class ScalarType,
class MV,
class OP>
202 const RCP<
SortManager<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
206 Teuchos::ParameterList ¶ms
242 void harmonicIterate();
290 bool isInitialized()
const;
309 int getNumIters()
const;
312 void resetNumIters();
321 RCP<const MV> getRitzVectors();
328 std::vector<Value<ScalarType> > getRitzValues();
339 std::vector<int> getRitzIndex();
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();
372 int getCurSubspaceDim()
const;
375 int getMaxSubspaceDim()
const;
387 RCP<StatusTest<ScalarType,MV,OP> > getStatusTest()
const;
401 void setBlockSize(
int blockSize);
404 int getBlockSize()
const;
423 void setAuxVecs(
const Teuchos::Array<RCP<const MV> > &auxvecs);
426 Teuchos::Array<RCP<const MV> > getAuxVecs()
const;
439 void setSize(
int blockSize,
int numBlocks);
445 void currentStatus(std::ostream &os);
456 typedef Teuchos::ScalarTraits<ScalarType> SCT;
457 typedef typename SCT::magnitudeType MagnitudeType;
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;
467 const RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
468 const RCP<SortManager<MagnitudeType> > sm_;
469 const RCP<OutputManager<ScalarType> > om_;
470 RCP<StatusTest<ScalarType,MV,OP> > tester_;
471 const RCP<MatOrthoManager<ScalarType,MV,OP> > orthman_;
483 RCP<Teuchos::Time> timerOp_, timerMOp_, timerSaddle_, timerSortEval_, timerDS_,
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_;
535 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > KK_, ritzVecs_;
538 Teuchos::Array<RCP<const MV> > auxVecs_, MauxVecs_;
545 std::vector<MagnitudeType> theta_, Rnorms_, R2norms_;
548 bool Rnorms_current_, R2norms_current_;
551 RCP<tracemin_ritz_op_type> ritzOp_;
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();
622 virtual void addToBasis(
const RCP<const MV> Delta) =0;
623 virtual void harmonicAddToBasis(
const RCP<const MV> Delta) =0;
639 template <
class ScalarType,
class MV,
class OP>
642 const RCP<
SortManager<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
646 Teuchos::ParameterList ¶ms
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),
690 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument,
691 "Anasazi::TraceMinBase::constructor: user passed null problem pointer.");
692 TEUCHOS_TEST_FOR_EXCEPTION(sm_ == Teuchos::null,std::invalid_argument,
693 "Anasazi::TraceMinBase::constructor: user passed null sort manager pointer.");
694 TEUCHOS_TEST_FOR_EXCEPTION(om_ == Teuchos::null,std::invalid_argument,
695 "Anasazi::TraceMinBase::constructor: user passed null output manager pointer.");
696 TEUCHOS_TEST_FOR_EXCEPTION(tester_ == Teuchos::null,std::invalid_argument,
697 "Anasazi::TraceMinBase::constructor: user passed null status test pointer.");
698 TEUCHOS_TEST_FOR_EXCEPTION(orthman_ == Teuchos::null,std::invalid_argument,
699 "Anasazi::TraceMinBase::constructor: user passed null orthogonalization manager pointer.");
700 TEUCHOS_TEST_FOR_EXCEPTION(problem_->isHermitian() ==
false, std::invalid_argument,
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);
720 TEUCHOS_TEST_FOR_EXCEPTION(traceThresh_ < 0, std::invalid_argument,
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, 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");
732 TEUCHOS_TEST_FOR_EXCEPTION(shiftNorm_ !=
"2" && shiftNorm_ !=
"M", std::invalid_argument,
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.");
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>
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_;
940 while (tester_->checkStatus(
this) !=
Passed && (numBlocks_ == 1 || curDim_ < searchDim)) {
943 if (om_->isVerbosity(
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_;
1028 while (tester_->checkStatus(
this) !=
Passed && (numBlocks_ == 1 || curDim_ < searchDim)) {
1031 if (om_->isVerbosity(
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++)
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];
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)
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 1161 Teuchos::TimeMonitor inittimer( *timerInit_ );
1164 previouslyLeveled_ =
false;
1168 harmonicInitialize(newstate);
1172 std::vector<int> bsind(blockSize_);
1173 for (
int i=0; i<blockSize_; ++i) bsind[i] = i;
1197 RCP<MV> lclV, lclKV, lclMV;
1198 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclKK, lclRV;
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." );
1206 TEUCHOS_TEST_FOR_EXCEPTION( newstate.
curDim < blockSize_, std::invalid_argument,
1207 "Anasazi::TraceMinBase::initialize(newstate): Rank of new state must be at least blockSize().");
1208 TEUCHOS_TEST_FOR_EXCEPTION( newstate.
curDim > blockSize_*numBlocks_, std::invalid_argument,
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_;
1218 TEUCHOS_TEST_FOR_EXCEPTION( curDim_ != newstate.
curDim, std::invalid_argument,
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_) {
1235 RCP<const MV> ivec = problem_->getInitVec();
1236 TEUCHOS_TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::invalid_argument,
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;
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;
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 1316 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1318 count_ApplyM_+= curDim_;
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_) {
1347 om_->stream(
Debug) <<
"There is no MV\n";
1356 om_->stream(
Debug) <<
"Project and normalize\n";
1358 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1359 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
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_,
1376 Teuchos::tuple(RCP< Teuchos::SerialDenseMatrix< int, ScalarType > >(Teuchos::null)),
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 1394 Teuchos::TimeMonitor lcltimer( *timerOp_ );
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;
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_) {
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;
1455 lclKK = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
1463 om_->stream(
Debug) <<
"Copying KK\n";
1466 TEUCHOS_TEST_FOR_EXCEPTION( newstate.
KK->numRows() < curDim_ || newstate.
KK->numCols() < curDim_, std::invalid_argument,
1467 "Anasazi::TraceMinBase::initialize(newstate): Projected matrix in new state must be as large as specified state rank.");
1470 lclKK = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
1471 if (newstate.
KK != KK_) {
1472 if (newstate.
KK->numRows() > curDim_ || newstate.
KK->numCols() > curDim_) {
1473 newstate.
KK = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*newstate.
KK,curDim_,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";
1499 TEUCHOS_TEST_FOR_EXCEPTION((
signed int)(newstate.
T->size()) != curDim_,
1500 std::invalid_argument,
"Anasazi::TraceMinBase::initialize(newstate): Size of T must be consistent with dimension of V.");
1502 TEUCHOS_TEST_FOR_EXCEPTION( newstate.
RV->numRows() < curDim_ || newstate.
RV->numCols() < curDim_, std::invalid_argument,
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());
1507 lclRV = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*ritzVecs_,curDim_,curDim_) );
1508 if (newstate.
RV != ritzVecs_) {
1509 if (newstate.
RV->numRows() > curDim_ || newstate.
RV->numCols() > curDim_) {
1510 newstate.
RV = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*newstate.
RV,curDim_,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)
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_) {
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_) {
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_) {
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_) {
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_) {
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_) {
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_) {
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_) {
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)) {
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;
1741 RCP<MV> lclV, lclKV, lclMV;
1742 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclKK, lclRV;
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." );
1750 TEUCHOS_TEST_FOR_EXCEPTION( newstate.
curDim < blockSize_, std::invalid_argument,
1751 "Anasazi::TraceMinBase::initialize(newstate): Rank of new state must be at least blockSize().");
1752 TEUCHOS_TEST_FOR_EXCEPTION( newstate.
curDim > blockSize_*numBlocks_, std::invalid_argument,
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_;
1762 TEUCHOS_TEST_FOR_EXCEPTION( curDim_ != newstate.
curDim, std::invalid_argument,
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_) {
1779 RCP<const MV> ivec = problem_->getInitVec();
1780 TEUCHOS_TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::invalid_argument,
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;
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;
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 1878 Teuchos::TimeMonitor lcltimer( *timerOp_ );
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;
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_) {
1915 om_->stream(
Debug) <<
"There is no KV\n";
1926 om_->stream(
Debug) <<
"Project and normalize\n";
1928 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 1929 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
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;
1943 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > gamma = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(curDim_,curDim_));
1944 int rank = orthman_->normalizeMat(*lclKV,gamma);
1947 Teuchos::SerialDenseSolver<int,ScalarType> SDsolver;
1948 SDsolver.setMatrix(gamma);
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 1963 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1965 count_ApplyM_+= curDim_;
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_) {
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;
2016 lclKK = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
2024 om_->stream(
Debug) <<
"Copying KK\n";
2027 TEUCHOS_TEST_FOR_EXCEPTION( newstate.
KK->numRows() < curDim_ || newstate.
KK->numCols() < curDim_, std::invalid_argument,
2028 "Anasazi::TraceMinBase::initialize(newstate): Projected matrix in new state must be as large as specified state rank.");
2031 lclKK = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
2032 if (newstate.
KK != KK_) {
2033 if (newstate.
KK->numRows() > curDim_ || newstate.
KK->numCols() > curDim_) {
2034 newstate.
KK = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*newstate.
KK,curDim_,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";
2060 TEUCHOS_TEST_FOR_EXCEPTION((
signed int)(newstate.
T->size()) != curDim_,
2061 std::invalid_argument,
"Anasazi::TraceMinBase::initialize(newstate): Size of T must be consistent with dimension of V.");
2063 TEUCHOS_TEST_FOR_EXCEPTION( newstate.
RV->numRows() < curDim_ || newstate.
RV->numCols() < curDim_, std::invalid_argument,
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());
2068 lclRV = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*ritzVecs_,curDim_,curDim_) );
2069 if (newstate.
RV != ritzVecs_) {
2070 if (newstate.
RV->numRows() > curDim_ || newstate.
RV->numCols() > curDim_) {
2071 newstate.
RV = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*newstate.
RV,curDim_,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)
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_) {
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_) {
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_) {
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_) {
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_) {
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_) {
2177 const int nvecs = computeAllRes_ ? curDim_ : blockSize_;
2178 Teuchos::Range1D dimind2 (0, nvecs-1);
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];
2188 if (Op_ != Teuchos::null) {
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_) {
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_) {
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)) {
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();
2319 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
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");
2352 if(Op_ != Teuchos::null) {
2353 om_->print(
Debug,
" >> Allocating KX_\n");
2360 om_->print(
Debug,
" >> Allocating MX_\n");
2366 om_->print(
Debug,
" >> Allocating R_\n");
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");
2388 if(Op_ != Teuchos::null) {
2389 om_->print(
Debug,
" >> Allocating KX_\n");
2396 om_->print(
Debug,
" >> Allocating MX_\n");
2402 om_->print(
Debug,
" >> Allocating R_\n");
2409 theta_.resize(newsd,NANVAL);
2410 om_->print(
Debug,
" >> Allocating V_\n");
2412 KK_ = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(newsd,newsd) );
2413 ritzVecs_ = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(newsd,newsd) );
2415 if(Op_ != Teuchos::null) {
2416 om_->print(
Debug,
" >> Allocating KV_\n");
2423 om_->print(
Debug,
" >> Allocating MV_\n");
2430 om_->print(
Debug,
" >> done allocating.\n");
2432 initialized_ =
false;
2439 template <
class ScalarType,
class MV,
class OP>
2441 typedef typename Teuchos::Array<RCP<const MV> >::iterator tarcpmv;
2447 MauxVecs_.resize(0);
2450 for (tarcpmv i=auxVecs_.begin(); i != auxVecs_.end(); ++i) {
2455 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 2456 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
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_);
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_);
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>
2569 TEUCHOS_TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument,
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>
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>
2692 std::vector<ScalarType> clusterResids(nvecs);
2693 std::vector<int> clusterIndices;
2694 for(
int i=0; i < nvecs; i++)
2697 if(clusterIndices.empty() || (theta_[i-1] + R2norms_[i-1] >= theta_[i] - R2norms_[i]))
2700 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;
2701 clusterIndices.push_back(i);
2706 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;
2707 ScalarType totalRes = ZERO;
2708 for(
size_t j=0; j < clusterIndices.size(); j++)
2709 totalRes += (R2norms_[clusterIndices[j]]*R2norms_[clusterIndices[j]]);
2714 if(theta_[clusterIndices[0]] < 0 && theta_[i] < 0)
2715 negSafeToShift_ =
true;
2716 else if(theta_[clusterIndices[0]] > 0 && theta_[i] > 0)
2717 posSafeToShift_ =
true;
2719 for(
size_t j=0; j < clusterIndices.size(); j++)
2720 clusterResids[clusterIndices[j]] = sqrt(totalRes);
2722 clusterIndices.clear();
2723 clusterIndices.push_back(i);
2728 ScalarType totalRes = ZERO;
2729 for(
size_t j=0; j < clusterIndices.size(); j++)
2730 totalRes += R2norms_[clusterIndices[j]];
2731 for(
size_t j=0; j < clusterIndices.size(); j++)
2732 clusterResids[clusterIndices[j]] = totalRes;
2734 return clusterResids;
2743 template <
class ScalarType,
class MV,
class OP>
2746 std::vector<ScalarType> thetaMag(theta_);
2747 bool traceHasLeveled = traceLeveled();
2750 for(
int i=0; i<blockSize_; i++)
2751 thetaMag[i] = std::abs(thetaMag[i]);
2759 if(whenToShift_ == ALWAYS_SHIFT || (whenToShift_ == SHIFT_WHEN_TRACE_LEVELS && traceHasLeveled))
2762 if(howToShift_ == LARGEST_CONVERGED_SHIFT)
2764 for(
int i=0; i<blockSize_; i++)
2765 ritzShifts_[i] = largestSafeShift_;
2768 else if(howToShift_ == RITZ_VALUES_SHIFT)
2770 ritzShifts_[0] = theta_[0];
2774 if(useMultipleShifts_)
2776 for(
int i=1; i<blockSize_; i++)
2777 ritzShifts_[i] = theta_[i];
2781 for(
int i=1; i<blockSize_; i++)
2782 ritzShifts_[i] = ritzShifts_[0];
2786 else if(howToShift_ == ADJUSTED_RITZ_SHIFT)
2788 om_->stream(
Debug) <<
"\nSeeking a shift for theta[0]=" << thetaMag[0] << std::endl;
2792 if((theta_[0] > 0 && posSafeToShift_) || (theta_[0] < 0 && negSafeToShift_) || considerClusters_ ==
false)
2795 ritzShifts_[0] = std::max(largestSafeShift_,thetaMag[0]-clusterResids[0]);
2797 om_->stream(
Debug) <<
"Initializing with a conservative shift, either the most positive converged eigenvalue (" 2798 << largestSafeShift_ <<
") or the eigenvalue adjusted by the residual (" << thetaMag[0] <<
"-" 2799 << clusterResids[0] <<
").\n";
2802 if(R2norms_[0] == clusterResids[0])
2804 ritzShifts_[0] = thetaMag[0];
2805 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;
2808 om_->stream(
Debug) <<
"This eigenvalue is in a cluster, so it would not be safe to use the eigenvalue itself as a shift\n";
2812 if(largestSafeShift_ > std::abs(ritzShifts_[0]))
2814 om_->stream(
Debug) <<
"Initializing with a conservative shift...the most positive converged eigenvalue: " << largestSafeShift_ << std::endl;
2815 ritzShifts_[0] = largestSafeShift_;
2818 om_->stream(
Debug) <<
"Using the previous value of ritzShifts[0]=" << ritzShifts_[0];
2822 om_->stream(
Debug) <<
"ritzShifts[0]=" << ritzShifts_[0] << std::endl;
2824 if(useMultipleShifts_)
2828 for(
int i=1; i < blockSize_; i++)
2830 om_->stream(
Debug) <<
"\nSeeking a shift for theta[" << i <<
"]=" << thetaMag[i] << std::endl;
2833 if(ritzShifts_[i-1] == thetaMag[i-1] && i < blockSize_-1 && thetaMag[i] < thetaMag[i+1] - clusterResids[i+1])
2835 ritzShifts_[i] = thetaMag[i];
2836 om_->stream(
Debug) <<
"Using an aggressive shift: ritzShifts_[" << i <<
"]=" << ritzShifts_[i] << std::endl;
2840 if(ritzShifts_[0] > std::abs(ritzShifts_[i]))
2842 om_->stream(
Debug) <<
"It was unsafe to use the aggressive shift. Choose the shift used by theta[0]=" 2843 << thetaMag[0] <<
": ritzShifts[0]=" << ritzShifts_[0] << std::endl;
2846 ritzShifts_[i] = ritzShifts_[0];
2849 om_->stream(
Debug) <<
"It was unsafe to use the aggressive shift. We will use the shift from the previous iteration: " << ritzShifts_[i] << std::endl;
2851 om_->stream(
Debug) <<
"Check whether any less conservative shifts would work (such as the biggest eigenvalue outside of the cluster, namely theta[ell] < " 2852 << thetaMag[i] <<
"-" << clusterResids[i] <<
" (" << thetaMag[i] - clusterResids[i] <<
")\n";
2855 for(
int ell=0; ell < i; ell++)
2857 if(thetaMag[ell] < thetaMag[i] - clusterResids[i])
2859 ritzShifts_[i] = thetaMag[ell];
2860 om_->stream(
Debug) <<
"ritzShifts_[" << i <<
"]=" << ritzShifts_[i] <<
" is valid\n";
2867 om_->stream(
Debug) <<
"ritzShifts[" << i <<
"]=" << ritzShifts_[i] << std::endl;
2872 for(
int i=1; i<blockSize_; i++)
2873 ritzShifts_[i] = ritzShifts_[0];
2879 for(
int i=0; i<blockSize_; i++)
2882 ritzShifts_[i] = -abs(ritzShifts_[i]);
2884 ritzShifts_[i] = abs(ritzShifts_[i]);
2889 template <
class ScalarType,
class MV,
class OP>
2892 ScalarType temp1, temp2;
2893 int nvecs = ritzShifts_.size();
2894 std::vector<ScalarType> tolerances;
2896 for(
int i=0; i < nvecs; i++)
2898 if(std::abs(theta_[0]) != std::abs(ritzShifts_[i]))
2899 temp1 = std::abs(theta_[i]-ritzShifts_[i])/std::abs(std::abs(theta_[0])-std::abs(ritzShifts_[i]));
2903 temp2 = pow(2.0,-iter_);
2907 tolerances.push_back(std::max(std::min(temp1*temp1,temp2),1e-8));
2911 tolerances[nvecs-1] = tolerances[nvecs-2];
2917 template <
class ScalarType,
class MV,
class OP>
2920 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 2921 Teuchos::TimeMonitor lcltimer( *timerSaddle_ );
2925 if(Op_ == Teuchos::null)
2928 Teuchos::SerialDenseSolver<int,ScalarType> My_Solver;
2931 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclS = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(blockSize_,blockSize_) );
2936 std::vector<int> curind(blockSize_);
2937 for(
int i=0; i<blockSize_; i++)
2947 My_Solver.setMatrix(lclS);
2961 My_Solver.setMatrix(lclS);
2971 std::vector<int> order(curDim_);
2972 std::vector<ScalarType> tempvec(blockSize_);
2976 std::vector<ScalarType> clusterResids;
2989 clusterResids = getClusterResids();
3004 computeRitzShifts(clusterResids);
3007 std::vector<ScalarType> tolerances = computeTol();
3009 for(
int i=0; i<blockSize_; i++)
3011 om_->stream(
IterationDetails) <<
"Choosing Ritz shifts...theta[" << i <<
"]=" 3012 << theta_[i] <<
", resids[" << i <<
"]=" << R2norms_[i] <<
", clusterResids[" << i <<
"]=" << clusterResids[i]
3013 <<
", ritzShifts[" << i <<
"]=" << ritzShifts_[i] <<
", and tol[" << i <<
"]=" << tolerances[i] << std::endl;
3017 ritzOp_->setRitzShifts(ritzShifts_);
3021 ritzOp_->setInnerTol(tolerances);
3024 if(saddleSolType_ == PROJECTED_KRYLOV_SOLVER)
3026 if(Prec_ != Teuchos::null)
3027 solveSaddleProjPrec(Delta);
3029 solveSaddleProj(Delta);
3031 else if(saddleSolType_ == SCHUR_COMPLEMENT_SOLVER)
3040 solveSaddleSchur(Delta);
3042 else if(saddleSolType_ == BD_PREC_MINRES)
3044 solveSaddleBDPrec(Delta);
3047 else if(saddleSolType_ == HSS_PREC_GMRES)
3049 solveSaddleHSSPrec(Delta);
3052 std::cout <<
"Invalid saddle solver type\n";
3059 template <
class ScalarType,
class MV,
class OP>
3062 RCP<TraceMinProjRitzOp<ScalarType,MV,OP> > projOp;
3067 std::vector<int> curind(blockSize_);
3068 for(
int i=0; i<blockSize_; i++)
3076 if(projectLockedVecs_ && numAuxVecs_ > 0)
3077 projOp = rcp(
new TraceMinProjRitzOp<ScalarType,MV,OP>(ritzOp_,locMX,orthman_,auxVecs_) );
3079 projOp = rcp(
new TraceMinProjRitzOp<ScalarType,MV,OP>(ritzOp_,locMX,orthman_) );
3082 projOp = rcp(
new TraceMinProjRitzOp<ScalarType,MV,OP>(ritzOp_,locMX) );
3091 projOp->ApplyInverse(*locR, *Delta);
3096 projOp->ApplyInverse(*locKX, *Delta);
3104 if(projectLockedVecs_ && numAuxVecs_ > 0)
3105 projOp = rcp(
new TraceMinProjRitzOp<ScalarType,MV,OP>(ritzOp_,MX_,orthman_,auxVecs_) );
3107 projOp = rcp(
new TraceMinProjRitzOp<ScalarType,MV,OP>(ritzOp_,MX_,orthman_) );
3110 projOp = rcp(
new TraceMinProjRitzOp<ScalarType,MV,OP>(ritzOp_,MX_) );
3117 projOp->ApplyInverse(*R_, *Delta);
3120 projOp->ApplyInverse(*KX_, *Delta);
3127 template <
class ScalarType,
class MV,
class OP>
3133 #ifdef HAVE_ANASAZI_BELOS 3134 RCP<TraceMinProjRitzOpWithPrec<ScalarType,MV,OP> > projOp;
3140 dimension = curDim_;
3142 dimension = blockSize_;
3145 std::vector<int> curind(dimension);
3146 for(
int i=0; i<dimension; i++)
3154 if(projectLockedVecs_ && numAuxVecs_ > 0)
3155 projOp = rcp(
new TraceMinProjRitzOpWithPrec<ScalarType,MV,OP>(ritzOp_,locMX,orthman_,auxVecs_) );
3157 projOp = rcp(
new TraceMinProjRitzOpWithPrec<ScalarType,MV,OP>(ritzOp_,locMX,orthman_) );
3160 projOp = rcp(
new TraceMinProjRitzOpWithPrec<ScalarType,MV,OP>(ritzOp_,locMX) );
3166 std::vector<int> dimind(blockSize_);
3167 for(
int i=0; i<blockSize_; i++)
3173 projOp->ApplyInverse(*locR, *Delta);
3179 projOp->ApplyInverse(*locKX, *Delta);
3187 if(projectLockedVecs_ && numAuxVecs_ > 0)
3188 projOp = rcp(
new TraceMinProjRitzOpWithPrec<ScalarType,MV,OP>(ritzOp_,MX_,orthman_,auxVecs_) );
3190 projOp = rcp(
new TraceMinProjRitzOpWithPrec<ScalarType,MV,OP>(ritzOp_,MX_,orthman_) );
3193 projOp = rcp(
new TraceMinProjRitzOpWithPrec<ScalarType,MV,OP>(ritzOp_,MX_) );
3201 projOp->ApplyInverse(*R_, *Delta);
3205 projOp->ApplyInverse(*KX_, *Delta);
3213 template <
class ScalarType,
class MV,
class OP>
3217 Teuchos::SerialDenseSolver<int,ScalarType> My_Solver;
3219 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclL;
3220 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclS = rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(blockSize_,blockSize_) );
3225 std::vector<int> curind(blockSize_);
3226 for(
int i=0; i<blockSize_; i++)
3233 #ifdef USE_APPLY_INVERSE 3234 Op_->ApplyInverse(*lclMX,*Z_);
3236 ritzOp_->ApplyInverse(*lclMX,*Z_);
3243 My_Solver.setMatrix(lclS);
3255 #ifdef USE_APPLY_INVERSE 3256 Op_->ApplyInverse(*MX_,*Z_);
3258 ritzOp_->ApplyInverse(*MX_,*Z_);
3265 My_Solver.setMatrix(lclS);
3278 template <
class ScalarType,
class MV,
class OP>
3281 RCP<MV> locKX, locMX;
3284 std::vector<int> curind(blockSize_);
3285 for(
int i=0; i<blockSize_; i++)
3297 RCP<saddle_op_type> sadOp = rcp(
new saddle_op_type(ritzOp_,locMX));
3300 RCP<saddle_container_type> sadRHS = rcp(
new saddle_container_type(locKX));
3308 RCP<saddle_container_type> sadSol = rcp(
new saddle_container_type(Delta));
3311 RCP<PseudoBlockMinres<ScalarType,saddle_container_type,saddle_op_type > > sadSolver;
3312 if(Prec_ != Teuchos::null)
3314 RCP<saddle_op_type> sadPrec = rcp(
new saddle_op_type(ritzOp_->getPrec(),locMX,BD_PREC));
3315 sadSolver = rcp(
new PseudoBlockMinres<ScalarType,saddle_container_type,saddle_op_type>(sadOp, sadPrec));
3318 sadSolver = rcp(
new PseudoBlockMinres<ScalarType,saddle_container_type,saddle_op_type>(sadOp));
3322 std::vector<ScalarType> tol;
3323 ritzOp_->getInnerTol(tol);
3324 sadSolver->setTol(tol);
3327 sadSolver->setMaxIter(ritzOp_->getMaxIts());
3330 sadSolver->setSol(sadSol);
3333 sadSolver->setRHS(sadRHS);
3343 template <
class ScalarType,
class MV,
class OP>
3346 #ifdef HAVE_ANASAZI_BELOS 3347 typedef Belos::LinearProblem<ScalarType,saddle_container_type,saddle_op_type> LP;
3348 typedef Belos::PseudoBlockGmresSolMgr<ScalarType,saddle_container_type,saddle_op_type> GmresSolMgr;
3350 RCP<MV> locKX, locMX;
3353 std::vector<int> curind(blockSize_);
3354 for(
int i=0; i<blockSize_; i++)
3366 RCP<saddle_op_type> sadOp = rcp(
new saddle_op_type(ritzOp_,locMX,NONSYM));
3369 RCP<saddle_container_type> sadRHS = rcp(
new saddle_container_type(locKX));
3373 RCP<saddle_container_type> sadSol = rcp(
new saddle_container_type(Delta));
3376 RCP<Teuchos::ParameterList> pl = rcp(
new Teuchos::ParameterList());
3379 std::vector<ScalarType> tol;
3380 ritzOp_->getInnerTol(tol);
3381 pl->set(
"Convergence Tolerance",tol[0]);
3385 pl->set(
"Maximum Iterations", ritzOp_->getMaxIts());
3386 pl->set(
"Num Blocks", ritzOp_->getMaxIts());
3391 pl->set(
"Block Size", blockSize_);
3398 RCP<LP> problem = rcp(
new LP(sadOp,sadSol,sadRHS));
3401 if(Prec_ != Teuchos::null)
3403 RCP<saddle_op_type> sadPrec = rcp(
new saddle_op_type(ritzOp_->getPrec(),locMX,HSS_PREC,alpha_));
3404 problem->setLeftPrec(sadPrec);
3408 problem->setProblem();
3411 RCP<GmresSolMgr> sadSolver = rcp(
new GmresSolMgr(problem,pl)) ;
3416 std::cout <<
"No Belos. This is bad\n";
3425 template <
class ScalarType,
class MV,
class OP>
3429 std::vector<int> curind(curDim_);
3430 for(
int i=0; i<curDim_; i++)
3437 curind.resize(blockSize_);
3438 for(
int i=0; i<blockSize_; i++)
3439 curind[i] = curDim_-blockSize_+i;
3443 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > lclKK =
3444 rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,blockSize_,0,curDim_-blockSize_) );
3450 for(
int r=0; r<curDim_; r++)
3452 for(
int c=0; c<r; c++)
3454 (*KK_)(r,c) = (*KK_)(c,r);
3461 template <
class ScalarType,
class MV,
class OP>
3465 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > lclKK =
3466 rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
3469 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > lclRV =
3470 rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*ritzVecs_,curDim_,curDim_) );
3474 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 3475 Teuchos::TimeMonitor lcltimer( *timerDS_ );
3482 "Anasazi::TraceMinBase::computeRitzPairs(): Failed to compute all eigenpairs of KK.");
3487 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 3488 Teuchos::TimeMonitor lcltimer( *timerSortEval_ );
3491 std::vector<int> order(curDim_);
3497 sm.
sort(theta_, Teuchos::rcpFromRef(order), curDim_);
3501 sm_->sort(theta_, Teuchos::rcpFromRef(order), curDim_);
3511 template <
class ScalarType,
class MV,
class OP>
3514 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 3515 Teuchos::TimeMonitor lcltimer( *timerLocal_ );
3519 std::vector<int> curind(curDim_);
3520 for(
int i=0; i<curDim_; i++)
3529 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > relevantEvecs =
3530 rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*ritzVecs_,curDim_,curDim_) );
3539 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > relevantEvecs =
3540 rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*ritzVecs_,curDim_,blockSize_) );
3549 template <
class ScalarType,
class MV,
class OP>
3552 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 3553 Teuchos::TimeMonitor lcltimer( *timerLocal_ );
3557 std::vector<int> curind(curDim_);
3558 for(
int i=0; i<curDim_; i++)
3568 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > relevantEvecs =
3569 rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*ritzVecs_,curDim_,curDim_) );
3584 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > relevantEvecs =
3585 rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*ritzVecs_,curDim_,blockSize_) );
3599 template <
class ScalarType,
class MV,
class OP>
3601 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 3602 Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
3608 std::vector<int> curind(curDim_);
3609 for(
int i=0; i<curDim_; i++)
3616 std::vector<ScalarType> locTheta(curDim_);
3617 for(
int i=0; i<curDim_; i++)
3618 locTheta[i] = theta_[i];
3634 std::vector<ScalarType> locTheta(blockSize_);
3635 for(
int i=0; i<blockSize_; i++)
3636 locTheta[i] = theta_[i];
3646 Rnorms_current_ =
false;
3647 R2norms_current_ =
false;
3677 template <
class ScalarType,
class MV,
class OP>
3682 std::stringstream os;
3684 os.setf(std::ios::scientific, std::ios::floatfield);
3686 os <<
" Debugging checks: iteration " << iter_ << where << endl;
3689 std::vector<int> lclind(curDim_);
3690 for (
int i=0; i<curDim_; ++i) lclind[i] = i;
3695 if (chk.checkV && initialized_) {
3696 MagnitudeType err = orthman_->orthonormError(*lclV);
3697 os <<
" >> Error in V^H M V == I : " << err << endl;
3698 for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
3699 err = orthman_->orthogError(*lclV,*auxVecs_[i]);
3700 os <<
" >> Error in V^H M Q[" << i <<
"] == 0 : " << err << endl;
3714 if (chk.checkX && initialized_) {
3715 MagnitudeType err = orthman_->orthonormError(*lclX);
3716 os <<
" >> Error in X^H M X == I : " << err << endl;
3717 for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
3718 err = orthman_->orthogError(*lclX,*auxVecs_[i]);
3719 os <<
" >> Error in X^H M Q[" << i <<
"] == 0 : " << err << endl;
3722 if (chk.checkMX && hasM_ && initialized_) {
3723 RCP<const MV> lclMX;
3730 os <<
" >> Error in MX == M*X : " << err << endl;
3732 if (Op_ != Teuchos::null && chk.checkKX && initialized_) {
3733 RCP<const MV> lclKX;
3740 os <<
" >> Error in KX == K*X : " << err << endl;
3744 if (chk.checkKK && initialized_) {
3745 Teuchos::SerialDenseMatrix<int,ScalarType> curKK(curDim_,curDim_);
3746 if(Op_ != Teuchos::null) {
3754 Teuchos::SerialDenseMatrix<int,ScalarType> subKK(Teuchos::View,*KK_,curDim_,curDim_);
3756 os <<
" >> Error in V^H K V == KK : " << curKK.normFrobenius() << endl;
3758 Teuchos::SerialDenseMatrix<int,ScalarType> SDMerr(curDim_,curDim_);
3759 for (
int j=0; j<curDim_; ++j) {
3760 for (
int i=0; i<curDim_; ++i) {
3761 SDMerr(i,j) = subKK(i,j) - SCT::conjugate(subKK(j,i));
3764 os <<
" >> Error in KK - KK^H == 0 : " << SDMerr.normFrobenius() << endl;
3769 for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
3770 MagnitudeType err = orthman_->orthonormError(*auxVecs_[i]);
3771 os <<
" >> Error in Q[" << i <<
"]^H M Q[" << i <<
"] == I : " << err << endl;
3772 for (Array_size_type j=i+1; j<auxVecs_.size(); ++j) {
3773 err = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
3774 os <<
" >> Error in Q[" << i <<
"]^H M Q[" << j <<
"] == 0 : " << err << endl;
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
static void MvInit(MV &mv, const ScalarType alpha=Teuchos::ScalarTraits< ScalarType >::zero())
Replace each element of the vectors in mv with alpha.
RCP< const MV > getRitzVectors()
Get access to the current Ritz vectors.
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
int NEV
Number of unconverged eigenvalues.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
RCP< const MV > X
The current eigenvectors.
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.
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...
int getBlockSize() const
Get the blocksize used by the iterative solver.
An operator of the form [A Y; Y' 0] where A is a sparse matrix and Y a multivector.
RCP< const MV > MopV
The image of V under M, or Teuchos::null if M was not specified.
TraceMinBaseState< ScalarType, MV > getState() const
Get access to the current state of the eigensolver.
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 getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace. For the trace minimization methods...
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRitzRes2Norms()
Get the 2-norms of the residuals.
static Teuchos::ScalarTraits< ScalarType >::magnitudeType errorEquality(const MV &X, const MV &MX, Teuchos::RCP< const OP > M=Teuchos::null)
Return the maximum coefficient of the matrix scaled by the maximum coefficient of MX...
static void Apply(const OP &Op, const MV &x, MV &y)
Application method which performs operation y = Op*x. An OperatorError exception is thrown if there i...
Basic implementation of the Anasazi::SortManager class.
TraceMinBaseInitFailure is thrown when the TraceMinBase solver is unable to generate an initial itera...
static void Assign(const MV &A, MV &mv)
mv := A
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.
bool isInitialized() const
Indicates whether the solver has been initialized or not.
Anasazi's templated virtual class for providing routines for orthogonalization and orthonormalization...
virtual ~TraceMinBase()
Anasazi::TraceMinBase destructor.
Output managers remove the need for the eigensolver to know any information about the required output...
static int directSolver(int size, const Teuchos::SerialDenseMatrix< int, ScalarType > &KK, Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > MM, Teuchos::SerialDenseMatrix< int, ScalarType > &EV, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &theta, int &nev, int esType=0)
Routine for computing the first NEV generalized eigenpairs of the Hermitian pencil (KK...
Teuchos::Array< RCP< const MV > > getAuxVecs() const
Get the auxiliary vectors for the solver.
Traits class which defines basic operations on multivectors.
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
void iterate()
This method performs trace minimization iterations until the status test indicates the need to stop o...
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &B, Teuchos::SerialDenseMatrix< int, ScalarType > &C)
Compute C := alpha * A^H B.
Virtual base class which defines basic traits for the operator type.
static void MvScale(MV &mv, const ScalarType alpha)
Scale each element of the vectors in mv with alpha.
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...
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
RCP< const MV > MX
The image of the current eigenvectors under M, or Teuchos::null if M was not specified.
RCP< StatusTest< ScalarType, MV, OP > > getStatusTest() const
Get the current StatusTest used by the solver.
static void MvNorm(const MV &mv, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec)
Compute the 2-norm of each individual vector of mv. Upon return, normvec[i] holds the value of ...
RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > RV
The current Ritz vectors.
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index...
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.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
std::vector< Value< ScalarType > > getRitzValues()
Get the Ritz values for the previous iteration.
RCP< const std::vector< ScalarType > > ritzShifts
Current Ritz shifts.
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.
static void permuteVectors(const int n, const std::vector< int > &perm, MV &Q, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > *resids=0)
Permute the vectors in a multivector according to the permutation vector perm, and optionally the res...
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.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the eigenvalue problem.
void setSize(int blockSize, int numBlocks)
Set the blocksize and number of blocks to be used by the iterative solver in solving this eigenproble...
This is an abstract base class for the trace minimization eigensolvers.
int getNumIters() const
Get the current iteration count.
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.
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current eigenvectors and eigenvalues...
The Eigensolver is a templated virtual base class that defines the basic interface that any eigensolv...
void resetNumIters()
Reset the iteration count.
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...
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
Class which provides internal utilities for the Anasazi solvers.