00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00033
00034
00035 #ifndef ANASAZI_RTRBASE_HPP
00036 #define ANASAZI_RTRBASE_HPP
00037
00038 #include "AnasaziTypes.hpp"
00039
00040 #include "AnasaziEigensolver.hpp"
00041 #include "AnasaziMultiVecTraits.hpp"
00042 #include "AnasaziOperatorTraits.hpp"
00043 #include "Teuchos_ScalarTraits.hpp"
00044
00045 #include "AnasaziGenOrthoManager.hpp"
00046 #include "AnasaziSolverUtils.hpp"
00047
00048 #include "Teuchos_LAPACK.hpp"
00049 #include "Teuchos_BLAS.hpp"
00050 #include "Teuchos_SerialDenseMatrix.hpp"
00051 #include "Teuchos_ParameterList.hpp"
00052 #include "Teuchos_TimeMonitor.hpp"
00053
00103 namespace Anasazi {
00104
00106
00107
00112 template <class ScalarType, class MV>
00113 struct RTRState {
00115 Teuchos::RCP<const MV> X;
00117 Teuchos::RCP<const MV> AX;
00119 Teuchos::RCP<const MV> BX;
00121 Teuchos::RCP<const MV> R;
00123 Teuchos::RCP<const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > T;
00127 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType rho;
00128 RTRState() : X(Teuchos::null),AX(Teuchos::null),BX(Teuchos::null),
00129 R(Teuchos::null),T(Teuchos::null),rho(0) {};
00130 };
00131
00133
00135
00136
00150 class RTRRitzFailure : public AnasaziError {public:
00151 RTRRitzFailure(const std::string& what_arg) : AnasaziError(what_arg)
00152 {}};
00153
00162 class RTRInitFailure : public AnasaziError {public:
00163 RTRInitFailure(const std::string& what_arg) : AnasaziError(what_arg)
00164 {}};
00165
00182 class RTROrthoFailure : public AnasaziError {public:
00183 RTROrthoFailure(const std::string& what_arg) : AnasaziError(what_arg)
00184 {}};
00185
00186
00188
00189
00190 template <class ScalarType, class MV, class OP>
00191 class RTRBase : public Eigensolver<ScalarType,MV,OP> {
00192 public:
00193
00195
00196
00202 RTRBase(const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00203 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
00204 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00205 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00206 const Teuchos::RCP<GenOrthoManager<ScalarType,MV,OP> > &ortho,
00207 Teuchos::ParameterList ¶ms,
00208 const std::string &solverLabel, bool skinnySolver
00209 );
00210
00212 virtual ~RTRBase() {};
00213
00215
00217
00218
00240 virtual void iterate() = 0;
00241
00266 void initialize(RTRState<ScalarType,MV> newstate);
00267
00271 void initialize();
00272
00285 bool isInitialized() const;
00286
00294 RTRState<ScalarType,MV> getState() const;
00295
00297
00299
00300
00302 int getNumIters() const;
00303
00305 void resetNumIters();
00306
00314 Teuchos::RCP<const MV> getRitzVectors();
00315
00321 std::vector<Value<ScalarType> > getRitzValues();
00322
00330 std::vector<int> getRitzIndex();
00331
00337 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getResNorms();
00338
00339
00345 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRes2Norms();
00346
00347
00352 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRitzRes2Norms();
00353
00354
00363 int getCurSubspaceDim() const;
00364
00367 int getMaxSubspaceDim() const;
00368
00370
00372
00373
00375 void setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test);
00376
00378 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > getStatusTest() const;
00379
00381 const Eigenproblem<ScalarType,MV,OP>& getProblem() const;
00382
00383
00392 void setBlockSize(int blockSize);
00393
00394
00396 int getBlockSize() const;
00397
00398
00419 void setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs);
00420
00422 Teuchos::Array<Teuchos::RCP<const MV> > getAuxVecs() const;
00423
00425
00427
00428
00430 virtual void currentStatus(std::ostream &os);
00431
00433
00434 protected:
00435
00436
00437
00438 typedef SolverUtils<ScalarType,MV,OP> Utils;
00439 typedef MultiVecTraits<ScalarType,MV> MVT;
00440 typedef OperatorTraits<ScalarType,MV,OP> OPT;
00441 typedef Teuchos::ScalarTraits<ScalarType> SCT;
00442 typedef typename SCT::magnitudeType MagnitudeType;
00443 typedef Teuchos::ScalarTraits<MagnitudeType> MAT;
00444 const MagnitudeType ONE;
00445 const MagnitudeType ZERO;
00446 const MagnitudeType NANVAL;
00447 typedef typename std::vector<MagnitudeType>::iterator vecMTiter;
00448 typedef typename std::vector<ScalarType>::iterator vecSTiter;
00449
00450
00451
00452 struct CheckList {
00453 bool checkX, checkAX, checkBX;
00454 bool checkEta, checkAEta, checkBEta;
00455 bool checkR, checkQ, checkBR;
00456 bool checkZ, checkPBX;
00457 CheckList() : checkX(false),checkAX(false),checkBX(false),
00458 checkEta(false),checkAEta(false),checkBEta(false),
00459 checkR(false),checkQ(false),checkBR(false),
00460 checkZ(false), checkPBX(false) {};
00461 };
00462
00463
00464
00465 std::string accuracyCheck(const CheckList &chk, const std::string &where) const;
00466
00467 virtual void solveTRSubproblem() = 0;
00468
00469 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType ginner(const MV &xi) const;
00470 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType ginner(const MV &xi, const MV &zeta) const;
00471 void ginnersep(const MV &xi, std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &d) const;
00472 void ginnersep(const MV &xi, const MV &zeta, std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &d) const;
00473
00474
00475
00476 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
00477 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sm_;
00478 const Teuchos::RCP<OutputManager<ScalarType> > om_;
00479 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > tester_;
00480 const Teuchos::RCP<GenOrthoManager<ScalarType,MV,OP> > orthman_;
00481
00482
00483
00484 Teuchos::RCP<const OP> AOp_;
00485 Teuchos::RCP<const OP> BOp_;
00486 Teuchos::RCP<const OP> Prec_;
00487 bool hasBOp_, hasPrec_;
00488
00489
00490
00491 Teuchos::RCP<Teuchos::Time> timerAOp_, timerBOp_, timerPrec_,
00492 timerSort_,
00493 timerLocalProj_, timerDS_,
00494 timerLocalUpdate_, timerCompRes_,
00495 timerOrtho_, timerInit_;
00496
00497
00498
00499
00500 int counterAOp_, counterBOp_, counterPrec_;
00501
00502
00503
00504
00505
00506 int blockSize_;
00507
00508
00509
00510
00511
00512
00513 bool initialized_;
00514
00515
00516
00517 int nevLocal_;
00518
00519
00520 bool isSkinny_;
00521
00522
00523 bool leftMost_;
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563 Teuchos::RCP<MV> V_, BV_, PBV_,
00564 AX_, R_,
00565 eta_, Aeta_, Beta_,
00566 delta_, Adelta_, Bdelta_, Hdelta_,
00567 Z_;
00568 Teuchos::RCP<const MV> X_, BX_;
00569
00570
00571 Teuchos::Array<Teuchos::RCP<const MV> > auxVecs_;
00572 int numAuxVecs_;
00573
00574
00575 int iter_;
00576
00577
00578 std::vector<MagnitudeType> theta_, Rnorms_, R2norms_, ritz2norms_;
00579
00580
00581 bool Rnorms_current_, R2norms_current_;
00582
00583
00584 MagnitudeType conv_kappa_, conv_theta_;
00585 MagnitudeType rho_;
00586
00587
00588 MagnitudeType fx_;
00589 };
00590
00591
00592
00593
00595
00596 template <class ScalarType, class MV, class OP>
00597 RTRBase<ScalarType,MV,OP>::RTRBase(
00598 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00599 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
00600 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00601 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00602 const Teuchos::RCP<GenOrthoManager<ScalarType,MV,OP> > &ortho,
00603 Teuchos::ParameterList ¶ms,
00604 const std::string &solverLabel, bool skinnySolver
00605 ) :
00606 ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
00607 ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
00608 NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan()),
00609
00610 problem_(problem),
00611 sm_(sorter),
00612 om_(printer),
00613 tester_(tester),
00614 orthman_(ortho),
00615
00616 timerAOp_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Operation A*x")),
00617 timerBOp_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Operation B*x")),
00618 timerPrec_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Operation Prec*x")),
00619 timerSort_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Sorting eigenvalues")),
00620 timerLocalProj_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Local projection")),
00621 timerDS_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Direct solve")),
00622 timerLocalUpdate_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Local update")),
00623 timerCompRes_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Computing residuals")),
00624 timerOrtho_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Orthogonalization")),
00625 timerInit_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Initialization")),
00626 counterAOp_(0),
00627 counterBOp_(0),
00628 counterPrec_(0),
00629
00630 blockSize_(0),
00631 initialized_(false),
00632 nevLocal_(0),
00633 isSkinny_(skinnySolver),
00634 leftMost_(true),
00635 auxVecs_( Teuchos::Array<Teuchos::RCP<const MV> >(0) ),
00636 numAuxVecs_(0),
00637 iter_(0),
00638 Rnorms_current_(false),
00639 R2norms_current_(false),
00640 conv_kappa_(.1),
00641 conv_theta_(1),
00642 rho_( MAT::nan() ),
00643 fx_( MAT::nan() )
00644 {
00645 TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument,
00646 "Anasazi::RTRBase::constructor: user passed null problem pointer.");
00647 TEST_FOR_EXCEPTION(sm_ == Teuchos::null,std::invalid_argument,
00648 "Anasazi::RTRBase::constructor: user passed null sort manager pointer.");
00649 TEST_FOR_EXCEPTION(om_ == Teuchos::null,std::invalid_argument,
00650 "Anasazi::RTRBase::constructor: user passed null output manager pointer.");
00651 TEST_FOR_EXCEPTION(tester_ == Teuchos::null,std::invalid_argument,
00652 "Anasazi::RTRBase::constructor: user passed null status test pointer.");
00653 TEST_FOR_EXCEPTION(orthman_ == Teuchos::null,std::invalid_argument,
00654 "Anasazi::RTRBase::constructor: user passed null orthogonalization manager pointer.");
00655 TEST_FOR_EXCEPTION(problem_->isProblemSet() == false, std::invalid_argument,
00656 "Anasazi::RTRBase::constructor: problem is not set.");
00657 TEST_FOR_EXCEPTION(problem_->isHermitian() == false, std::invalid_argument,
00658 "Anasazi::RTRBase::constructor: problem is not Hermitian.");
00659
00660
00661 AOp_ = problem_->getOperator();
00662 TEST_FOR_EXCEPTION(AOp_ == Teuchos::null, std::invalid_argument,
00663 "Anasazi::RTRBase::constructor: problem provides no A matrix.");
00664 BOp_ = problem_->getM();
00665 Prec_ = problem_->getPrec();
00666 hasBOp_ = (BOp_ != Teuchos::null);
00667 hasPrec_ = (Prec_ != Teuchos::null);
00668
00669 TEST_FOR_EXCEPTION(orthman_->getOp() != BOp_,std::invalid_argument,
00670 "Anasazi::RTRBase::constructor: orthogonalization manager must use mass matrix for inner product.");
00671
00672
00673 int bs = params.get("Block Size", problem_->getNEV());
00674 setBlockSize(bs);
00675
00676
00677 leftMost_ = params.get("Leftmost",leftMost_);
00678
00679 conv_kappa_ = params.get("Kappa Convergence",conv_kappa_);
00680 TEST_FOR_EXCEPTION(conv_kappa_ <= 0 || conv_kappa_ >= 1,std::invalid_argument,
00681 "Anasazi::RTRBase::constructor: kappa must be in (0,1).");
00682 conv_theta_ = params.get("Theta Convergence",conv_theta_);
00683 TEST_FOR_EXCEPTION(conv_theta_ <= 0,std::invalid_argument,
00684 "Anasazi::RTRBase::constructor: theta must be strictly postitive.");
00685 }
00686
00687
00689
00690 template <class ScalarType, class MV, class OP>
00691 void RTRBase<ScalarType,MV,OP>::setBlockSize (int blockSize)
00692 {
00693
00694 Teuchos::TimeMonitor lcltimer( *timerInit_ );
00695
00696
00697
00698
00699
00700 Teuchos::RCP<const MV> tmp;
00701
00702
00703
00704
00705
00706 if (blockSize_ > 0) {
00707 tmp = R_;
00708 }
00709 else {
00710 tmp = problem_->getInitVec();
00711 TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::logic_error,
00712 "Anasazi::RTRBase::setBlockSize(): Eigenproblem did not specify initial vectors to clone from");
00713 }
00714
00715 TEST_FOR_EXCEPTION(blockSize <= 0 || blockSize > MVT::GetVecLength(*tmp), std::invalid_argument,
00716 "Anasazi::RTRBase::setBlockSize was passed a non-positive block size");
00717
00718
00719 if (blockSize == blockSize_) {
00720
00721 return;
00722 }
00723
00724
00725 X_ = Teuchos::null;
00726 BX_ = Teuchos::null;
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738 if (blockSize > blockSize_)
00739 {
00740
00741
00742 Teuchos::RCP<const MV> Q;
00743 std::vector<int> indQ(numAuxVecs_);
00744 for (int i=0; i<numAuxVecs_; i++) indQ[i] = i;
00745
00746 TEST_FOR_EXCEPTION(numAuxVecs_ > 0 && blockSize_ == 0, std::logic_error,
00747 "Anasazi::RTRBase::setSize(): logic error. Please contact Anasazi team.");
00748
00749 if (numAuxVecs_ > 0) Q = MVT::CloneView(*V_,indQ);
00750 V_ = MVT::Clone(*tmp,numAuxVecs_ + blockSize);
00751 if (numAuxVecs_ > 0) MVT::SetBlock(*Q,indQ,*V_);
00752
00753 if (hasBOp_) {
00754 if (numAuxVecs_ > 0) Q = MVT::CloneView(*BV_,indQ);
00755 BV_ = MVT::Clone(*tmp,numAuxVecs_ + blockSize);
00756 if (numAuxVecs_ > 0) MVT::SetBlock(*Q,indQ,*BV_);
00757 }
00758 else {
00759 BV_ = V_;
00760 }
00761
00762 if (hasPrec_) {
00763 if (numAuxVecs_ > 0) Q = MVT::CloneView(*PBV_,indQ);
00764 PBV_ = MVT::Clone(*tmp,numAuxVecs_ + blockSize);
00765 if (numAuxVecs_ > 0) MVT::SetBlock(*Q,indQ,*PBV_);
00766 }
00767 }
00768 else
00769 {
00770
00771 std::vector<int> indV(numAuxVecs_+blockSize);
00772 for (int i=0; i<numAuxVecs_+blockSize; i++) indV[i] = i;
00773
00774 V_ = MVT::CloneCopy(*V_,indV);
00775
00776 if (hasBOp_) {
00777 BV_ = MVT::CloneCopy(*BV_,indV);
00778 }
00779 else {
00780 BV_ = V_;
00781 }
00782
00783 if (hasPrec_) {
00784 PBV_ = MVT::CloneCopy(*PBV_,indV);
00785 }
00786 }
00787
00788 if (blockSize < blockSize_) {
00789
00790 blockSize_ = blockSize;
00791
00792 theta_.resize(blockSize_);
00793 ritz2norms_.resize(blockSize_);
00794 Rnorms_.resize(blockSize_);
00795 R2norms_.resize(blockSize_);
00796
00797 if (initialized_) {
00798
00799 std::vector<int> ind(blockSize_);
00800 for (int i=0; i<blockSize_; i++) ind[i] = i;
00801
00802
00803 Z_ = Teuchos::null;
00804
00805
00806 tmp = Teuchos::null;
00807
00808 R_ = MVT::CloneCopy(*R_ ,ind);
00809 eta_ = MVT::CloneCopy(*eta_ ,ind);
00810 delta_ = MVT::CloneCopy(*delta_ ,ind);
00811 Hdelta_ = MVT::CloneCopy(*Hdelta_,ind);
00812 if (!isSkinny_) {
00813 AX_ = MVT::CloneCopy(*AX_ ,ind);
00814 Aeta_ = MVT::CloneCopy(*Aeta_ ,ind);
00815 Adelta_ = MVT::CloneCopy(*Adelta_,ind);
00816 }
00817 else {
00818 AX_ = Teuchos::null;
00819 Aeta_ = Teuchos::null;
00820 Adelta_ = Teuchos::null;
00821 }
00822
00823 if (hasBOp_) {
00824 if (!isSkinny_) {
00825 Beta_ = MVT::CloneCopy(*Beta_,ind);
00826 Bdelta_ = MVT::CloneCopy(*Bdelta_,ind);
00827 }
00828 else {
00829 Beta_ = Teuchos::null;
00830 Bdelta_ = Teuchos::null;
00831 }
00832 }
00833 else {
00834 Beta_ = eta_;
00835 Bdelta_ = delta_;
00836 }
00837
00838
00839
00840 if (hasPrec_ || isSkinny_) {
00841 Z_ = MVT::Clone(*V_,blockSize_);
00842 }
00843 else {
00844 Z_ = R_;
00845 }
00846
00847 }
00848 else {
00849
00850
00851 AX_ = Teuchos::null;
00852 R_ = Teuchos::null;
00853 eta_ = Teuchos::null;
00854 Aeta_ = Teuchos::null;
00855 delta_ = Teuchos::null;
00856 Adelta_ = Teuchos::null;
00857 Hdelta_ = Teuchos::null;
00858 Beta_ = Teuchos::null;
00859 Bdelta_ = Teuchos::null;
00860 Z_ = Teuchos::null;
00861
00862 R_ = MVT::Clone(*tmp,blockSize_);
00863 eta_ = MVT::Clone(*tmp,blockSize_);
00864 delta_ = MVT::Clone(*tmp,blockSize_);
00865 Hdelta_ = MVT::Clone(*tmp,blockSize_);
00866 if (!isSkinny_) {
00867 AX_ = MVT::Clone(*tmp,blockSize_);
00868 Aeta_ = MVT::Clone(*tmp,blockSize_);
00869 Adelta_ = MVT::Clone(*tmp,blockSize_);
00870 }
00871
00872 if (hasBOp_) {
00873 if (!isSkinny_) {
00874 Beta_ = MVT::Clone(*tmp,blockSize_);
00875 Bdelta_ = MVT::Clone(*tmp,blockSize_);
00876 }
00877 }
00878 else {
00879 Beta_ = eta_;
00880 Bdelta_ = delta_;
00881 }
00882
00883
00884
00885 if (hasPrec_ || isSkinny_) {
00886 Z_ = MVT::Clone(*tmp,blockSize_);
00887 }
00888 else {
00889 Z_ = R_;
00890 }
00891 }
00892 }
00893 else {
00894
00895 initialized_ = false;
00896
00897
00898 theta_.resize(blockSize,NANVAL);
00899 ritz2norms_.resize(blockSize,NANVAL);
00900 Rnorms_.resize(blockSize,NANVAL);
00901 R2norms_.resize(blockSize,NANVAL);
00902
00903
00904 AX_ = Teuchos::null;
00905 R_ = Teuchos::null;
00906 eta_ = Teuchos::null;
00907 Aeta_ = Teuchos::null;
00908 delta_ = Teuchos::null;
00909 Adelta_ = Teuchos::null;
00910 Hdelta_ = Teuchos::null;
00911 Beta_ = Teuchos::null;
00912 Bdelta_ = Teuchos::null;
00913 Z_ = Teuchos::null;
00914
00915
00916 R_ = MVT::Clone(*tmp,blockSize);
00917 eta_ = MVT::Clone(*tmp,blockSize);
00918 delta_ = MVT::Clone(*tmp,blockSize);
00919 Hdelta_ = MVT::Clone(*tmp,blockSize);
00920 if (!isSkinny_) {
00921 AX_ = MVT::Clone(*tmp,blockSize);
00922 Aeta_ = MVT::Clone(*tmp,blockSize);
00923 Adelta_ = MVT::Clone(*tmp,blockSize);
00924 }
00925
00926 if (hasBOp_) {
00927 if (!isSkinny_) {
00928 Beta_ = MVT::Clone(*tmp,blockSize);
00929 Bdelta_ = MVT::Clone(*tmp,blockSize);
00930 }
00931 }
00932 else {
00933 Beta_ = eta_;
00934 Bdelta_ = delta_;
00935 }
00936 if (hasPrec_ || isSkinny_) {
00937 Z_ = MVT::Clone(*tmp,blockSize);
00938 }
00939 else {
00940 Z_ = R_;
00941 }
00942 blockSize_ = blockSize;
00943 }
00944
00945
00946
00947 {
00948 std::vector<int> indX(blockSize_);
00949 for (int i=0; i<blockSize_; i++) indX[i] = numAuxVecs_+i;
00950 X_ = MVT::CloneView(*V_,indX);
00951 if (hasBOp_) {
00952 BX_ = MVT::CloneView(*BV_,indX);
00953 }
00954 else {
00955 BX_ = X_;
00956 }
00957 }
00958 }
00959
00960
00962
00963 template <class ScalarType, class MV, class OP>
00964 void RTRBase<ScalarType,MV,OP>::setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test) {
00965 TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument,
00966 "Anasazi::RTRBase::setStatusTest() was passed a null StatusTest.");
00967 tester_ = test;
00968 }
00969
00970
00972
00973 template <class ScalarType, class MV, class OP>
00974 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > RTRBase<ScalarType,MV,OP>::getStatusTest() const {
00975 return tester_;
00976 }
00977
00978
00980
00981 template <class ScalarType, class MV, class OP>
00982 void RTRBase<ScalarType,MV,OP>::setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs) {
00983 typedef typename Teuchos::Array<Teuchos::RCP<const MV> >::const_iterator tarcpmv;
00984
00985
00986 auxVecs_.resize(0);
00987 auxVecs_.reserve(auxvecs.size());
00988
00989 numAuxVecs_ = 0;
00990 for (tarcpmv v=auxvecs.begin(); v != auxvecs.end(); ++v) {
00991 numAuxVecs_ += MVT::GetNumberVecs(**v);
00992 }
00993
00994
00995 if (numAuxVecs_ > 0 && initialized_) {
00996 initialized_ = false;
00997 }
00998
00999
01000 X_ = Teuchos::null;
01001 BX_ = Teuchos::null;
01002
01003 V_ = Teuchos::null;
01004 BV_ = Teuchos::null;
01005 PBV_ = Teuchos::null;
01006
01007
01008 if (numAuxVecs_ > 0) {
01009 V_ = MVT::Clone(*R_,numAuxVecs_ + blockSize_);
01010 int numsofar = 0;
01011 for (tarcpmv v=auxvecs.begin(); v != auxvecs.end(); ++v) {
01012 std::vector<int> ind(MVT::GetNumberVecs(**v));
01013 for (unsigned int j=0; j<ind.size(); j++) ind[j] = numsofar++;
01014 MVT::SetBlock(**v,ind,*V_);
01015 auxVecs_.push_back(MVT::CloneView(*Teuchos::rcp_static_cast<const MV>(V_),ind));
01016 }
01017 TEST_FOR_EXCEPTION(numsofar != numAuxVecs_, std::logic_error,
01018 "Anasazi::RTRBase::setAuxVecs(): Logic error. Please contact Anasazi team.");
01019
01020 if (hasBOp_) {
01021 BV_ = MVT::Clone(*R_,numAuxVecs_ + blockSize_);
01022 OPT::Apply(*BOp_,*V_,*BV_);
01023 }
01024 else {
01025 BV_ = V_;
01026 }
01027 if (hasPrec_) {
01028 PBV_ = MVT::Clone(*R_,numAuxVecs_ + blockSize_);
01029 OPT::Apply(*Prec_,*BV_,*V_);
01030 }
01031 }
01032
01033
01034 if (om_->isVerbosity( Debug ) ) {
01035
01036 CheckList chk;
01037 chk.checkQ = true;
01038 om_->print( Debug, accuracyCheck(chk, "in setAuxVecs()") );
01039 }
01040 }
01041
01042
01044
01045
01046
01047
01048
01049
01050
01051
01052
01053
01054
01055 template <class ScalarType, class MV, class OP>
01056 void RTRBase<ScalarType,MV,OP>::initialize(RTRState<ScalarType,MV> newstate)
01057 {
01058
01059
01060
01061
01062
01063 X_ = Teuchos::null;
01064 BX_ = Teuchos::null;
01065 Teuchos::RCP<MV> X, BX;
01066 {
01067 std::vector<int> ind(blockSize_);
01068 for (int i=0; i<blockSize_; ++i) ind[i] = numAuxVecs_+i;
01069 X = MVT::CloneView(*V_,ind);
01070 if (hasBOp_) {
01071 BX = MVT::CloneView(*BV_,ind);
01072 }
01073 else {
01074 BX = X;
01075 }
01076 }
01077
01078 Teuchos::TimeMonitor inittimer( *timerInit_ );
01079
01080 std::vector<int> bsind(blockSize_);
01081 for (int i=0; i<blockSize_; i++) bsind[i] = i;
01082
01083
01084
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100 if (newstate.X != Teuchos::null) {
01101 TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.X) != MVT::GetVecLength(*X),
01102 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): vector length of newstate.X not correct." );
01103
01104 TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.X) < blockSize_,
01105 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.X must have at least block size vectors.");
01106
01107
01108 MVT::SetBlock(*newstate.X,bsind,*X);
01109
01110
01111
01112
01113
01114 if (isSkinny_) {
01115 AX_ = Z_;
01116 }
01117 if (newstate.AX != Teuchos::null) {
01118 TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.AX) != MVT::GetVecLength(*X),
01119 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): vector length of newstate.AX not correct." );
01120
01121 TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.AX) < blockSize_,
01122 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.AX must have at least block size vectors.");
01123 MVT::SetBlock(*newstate.AX,bsind,*AX_);
01124 }
01125 else {
01126 {
01127 Teuchos::TimeMonitor lcltimer( *timerAOp_ );
01128 OPT::Apply(*AOp_,*X,*AX_);
01129 counterAOp_ += blockSize_;
01130 }
01131
01132 newstate.R = Teuchos::null;
01133 }
01134
01135
01136
01137 if (hasBOp_) {
01138 if (newstate.BX != Teuchos::null) {
01139 TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.BX) != MVT::GetVecLength(*X),
01140 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): vector length of newstate.BX not correct." );
01141
01142 TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.BX) < blockSize_,
01143 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.BX must have at least block size vectors.");
01144 MVT::SetBlock(*newstate.BX,bsind,*BX);
01145 }
01146 else {
01147 {
01148 Teuchos::TimeMonitor lcltimer( *timerBOp_ );
01149 OPT::Apply(*BOp_,*X,*BX);
01150 counterBOp_ += blockSize_;
01151 }
01152
01153 newstate.R = Teuchos::null;
01154 }
01155 }
01156 else {
01157
01158 TEST_FOR_EXCEPTION(BX != X, std::logic_error, "Anasazi::RTRBase::initialize(): solver invariant not satisfied (BX==X).");
01159 }
01160
01161 }
01162 else {
01163
01164
01165
01166 newstate.R = Teuchos::null;
01167 newstate.T = Teuchos::null;
01168
01169
01170 Teuchos::RCP<const MV> ivec = problem_->getInitVec();
01171 TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::logic_error,
01172 "Anasazi::RTRBase::initialize(): Eigenproblem did not specify initial vectors to clone from.");
01173
01174 int initSize = MVT::GetNumberVecs(*ivec);
01175 if (initSize > blockSize_) {
01176
01177 initSize = blockSize_;
01178 std::vector<int> ind(blockSize_);
01179 for (int i=0; i<blockSize_; i++) ind[i] = i;
01180 ivec = MVT::CloneView(*ivec,ind);
01181 }
01182
01183
01184 if (initSize > 0) {
01185 std::vector<int> ind(initSize);
01186 for (int i=0; i<initSize; i++) ind[i] = i;
01187 MVT::SetBlock(*ivec,ind,*X);
01188 }
01189
01190 if (blockSize_ > initSize) {
01191 std::vector<int> ind(blockSize_ - initSize);
01192 for (int i=0; i<blockSize_ - initSize; i++) ind[i] = initSize + i;
01193 Teuchos::RCP<MV> rX = MVT::CloneView(*X,ind);
01194 MVT::MvRandom(*rX);
01195 rX = Teuchos::null;
01196 }
01197
01198
01199 if (hasBOp_) {
01200 Teuchos::TimeMonitor lcltimer( *timerBOp_ );
01201 OPT::Apply(*BOp_,*X,*BX);
01202 counterBOp_ += blockSize_;
01203 }
01204 else {
01205
01206 TEST_FOR_EXCEPTION(BX != X, std::logic_error, "Anasazi::RTRBase::initialize(): solver invariant not satisfied (BX==X).");
01207 }
01208
01209
01210 if (numAuxVecs_ > 0) {
01211 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
01212 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
01213 int rank = orthman_->projectAndNormalizeMat(*X,auxVecs_,dummyC,Teuchos::null,BX);
01214 TEST_FOR_EXCEPTION(rank != blockSize_, RTRInitFailure,
01215 "Anasazi::RTRBase::initialize(): Couldn't generate initial basis of full rank.");
01216 }
01217 else {
01218 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
01219 int rank = orthman_->normalizeMat(*X,Teuchos::null,BX);
01220 TEST_FOR_EXCEPTION(rank != blockSize_, RTRInitFailure,
01221 "Anasazi::RTRBase::initialize(): Couldn't generate initial basis of full rank.");
01222 }
01223
01224
01225 if (isSkinny_) {
01226 AX_ = Z_;
01227 }
01228 {
01229 Teuchos::TimeMonitor lcltimer( *timerAOp_ );
01230 OPT::Apply(*AOp_,*X,*AX_);
01231 counterAOp_ += blockSize_;
01232 }
01233
01234 }
01235
01236
01237
01238 if (newstate.T != Teuchos::null) {
01239 TEST_FOR_EXCEPTION( (signed int)(newstate.T->size()) < blockSize_,
01240 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.T must contain at least block size Ritz values.");
01241 for (int i=0; i<blockSize_; i++) {
01242 theta_[i] = (*newstate.T)[i];
01243 }
01244 }
01245 else {
01246
01247 Teuchos::SerialDenseMatrix<int,ScalarType> AA(blockSize_,blockSize_),
01248 BB(blockSize_,blockSize_),
01249 S(blockSize_,blockSize_);
01250
01251 {
01252 Teuchos::TimeMonitor lcltimer( *timerLocalProj_ );
01253 MVT::MvTransMv(ONE,*X,*AX_,AA);
01254 if (hasBOp_) {
01255 MVT::MvTransMv(ONE,*X,*BX,BB);
01256 }
01257 }
01258 nevLocal_ = blockSize_;
01259
01260
01261
01262 int ret;
01263 if (hasBOp_) {
01264 Teuchos::TimeMonitor lcltimer( *timerDS_ );
01265 ret = Utils::directSolver(blockSize_,AA,Teuchos::rcpFromRef(BB),S,theta_,nevLocal_,1);
01266 }
01267 else {
01268 Teuchos::TimeMonitor lcltimer( *timerDS_ );
01269 ret = Utils::directSolver(blockSize_,AA,Teuchos::null,S,theta_,nevLocal_,10);
01270 }
01271 TEST_FOR_EXCEPTION(ret != 0,RTRInitFailure,
01272 "Anasazi::RTRBase::initialize(): failure solving projected eigenproblem after retraction. LAPACK returns " << ret);
01273 TEST_FOR_EXCEPTION(nevLocal_ != blockSize_,RTRInitFailure,"Anasazi::RTRBase::initialize(): retracted iterate failed in Ritz analysis.");
01274
01275
01276
01277 {
01278 Teuchos::TimeMonitor lcltimer( *timerSort_ );
01279
01280 std::vector<int> order(blockSize_);
01281
01282
01283 sm_->sort(theta_, Teuchos::rcpFromRef(order), blockSize_);
01284
01285
01286 Utils::permuteVectors(order,S);
01287 }
01288
01289
01290 {
01291 Teuchos::BLAS<int,ScalarType> blas;
01292 Teuchos::SerialDenseMatrix<int,ScalarType> RR(blockSize_,blockSize_);
01293
01294 int info;
01295 if (hasBOp_) {
01296 info = RR.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,BB,S,ZERO);
01297 TEST_FOR_EXCEPTION(info != 0, std::logic_error, "Anasazi::RTRBase::initialize(): Logic error calling SerialDenseMatrix::multiply.");
01298 }
01299 else {
01300 RR.assign(S);
01301 }
01302 for (int i=0; i<blockSize_; i++) {
01303 blas.SCAL(blockSize_,theta_[i],RR[i],1);
01304 }
01305 info = RR.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,AA,S,-ONE);
01306 TEST_FOR_EXCEPTION(info != 0, std::logic_error, "Anasazi::RTRBase::initialize(): Logic error calling SerialDenseMatrix::multiply.");
01307 for (int i=0; i<blockSize_; i++) {
01308 ritz2norms_[i] = blas.NRM2(blockSize_,RR[i],1);
01309 }
01310 }
01311
01312
01313
01314
01315 {
01316 Teuchos::TimeMonitor lcltimer( *timerLocalUpdate_ );
01317
01318 MVT::MvAddMv( ONE, *X, ZERO, *X, *R_ );
01319 MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *X );
01320
01321 MVT::MvAddMv( ONE, *AX_, ZERO, *AX_, *R_ );
01322 MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *AX_ );
01323 if (hasBOp_) {
01324
01325 MVT::MvAddMv( ONE, *BX, ZERO, *BX, *R_ );
01326 MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *BX );
01327 }
01328 }
01329 }
01330
01331
01332 X = Teuchos::null;
01333 BX = Teuchos::null;
01334 {
01335 std::vector<int> ind(blockSize_);
01336 for (int i=0; i<blockSize_; ++i) ind[i] = numAuxVecs_+i;
01337 this->X_ = MVT::CloneView(static_cast<const MV&>(*V_),ind);
01338 if (hasBOp_) {
01339 this->BX_ = MVT::CloneView(static_cast<const MV&>(*BV_),ind);
01340 }
01341 else {
01342 this->BX_ = this->X_;
01343 }
01344 }
01345
01346
01347
01348 fx_ = std::accumulate(theta_.begin(),theta_.end(),ZERO);
01349
01350
01351 if (newstate.R != Teuchos::null) {
01352 TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) < blockSize_,
01353 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.R must have blockSize number of vectors." );
01354 TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.R) != MVT::GetVecLength(*R_),
01355 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): vector length of newstate.R not correct." );
01356 MVT::SetBlock(*newstate.R,bsind,*R_);
01357 }
01358 else {
01359 Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
01360
01361 MVT::MvAddMv(ZERO,*AX_,ONE,*AX_,*R_);
01362 Teuchos::SerialDenseMatrix<int,ScalarType> T(blockSize_,blockSize_);
01363 T.putScalar(ZERO);
01364 for (int i=0; i<blockSize_; i++) T(i,i) = theta_[i];
01365 MVT::MvTimesMatAddMv(-ONE,*BX_,T,ONE,*R_);
01366 }
01367
01368
01369 Rnorms_current_ = false;
01370 R2norms_current_ = false;
01371
01372
01373
01374 if (isSkinny_) {
01375 AX_ = Teuchos::null;
01376 }
01377
01378
01379 initialized_ = true;
01380
01381 if (om_->isVerbosity( Debug ) ) {
01382
01383 CheckList chk;
01384 chk.checkX = true;
01385 chk.checkAX = true;
01386 chk.checkBX = true;
01387 chk.checkR = true;
01388 chk.checkQ = true;
01389 om_->print( Debug, accuracyCheck(chk, "after initialize()") );
01390 }
01391 }
01392
01393 template <class ScalarType, class MV, class OP>
01394 void RTRBase<ScalarType,MV,OP>::initialize()
01395 {
01396 RTRState<ScalarType,MV> empty;
01397 initialize(empty);
01398 }
01399
01400
01401
01402
01404
01405 template <class ScalarType, class MV, class OP>
01406 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
01407 RTRBase<ScalarType,MV,OP>::getResNorms() {
01408 if (Rnorms_current_ == false) {
01409
01410 orthman_->norm(*R_,Rnorms_);
01411 Rnorms_current_ = true;
01412 }
01413 return Rnorms_;
01414 }
01415
01416
01418
01419 template <class ScalarType, class MV, class OP>
01420 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
01421 RTRBase<ScalarType,MV,OP>::getRes2Norms() {
01422 if (R2norms_current_ == false) {
01423
01424 MVT::MvNorm(*R_,R2norms_);
01425 R2norms_current_ = true;
01426 }
01427 return R2norms_;
01428 }
01429
01430
01431
01432
01434
01435
01436
01437
01438
01439
01440
01441
01442
01443
01444
01445
01446
01447
01448
01449
01450
01451
01452
01453
01454
01455
01456
01457 template <class ScalarType, class MV, class OP>
01458 std::string RTRBase<ScalarType,MV,OP>::accuracyCheck( const CheckList &chk, const std::string &where ) const
01459 {
01460 using std::setprecision;
01461 using std::scientific;
01462 using std::setw;
01463 using std::endl;
01464 std::stringstream os;
01465 MagnitudeType tmp;
01466
01467 os << " Debugging checks: " << where << endl;
01468
01469
01470 if (chk.checkX && initialized_) {
01471 tmp = orthman_->orthonormError(*X_);
01472 os << " >> Error in X^H B X == I : " << scientific << setprecision(10) << tmp << endl;
01473 for (unsigned int i=0; i<auxVecs_.size(); i++) {
01474 tmp = orthman_->orthogError(*X_,*auxVecs_[i]);
01475 os << " >> Error in X^H B Q[" << i << "] == 0 : " << scientific << setprecision(10) << tmp << endl;
01476 }
01477 }
01478 if (chk.checkAX && !isSkinny_ && initialized_) {
01479 tmp = Utils::errorEquality(*X_, *AX_, AOp_);
01480 os << " >> Error in AX == A*X : " << scientific << setprecision(10) << tmp << endl;
01481 }
01482 if (chk.checkBX && hasBOp_ && initialized_) {
01483 Teuchos::RCP<MV> BX = MVT::Clone(*this->X_,this->blockSize_);
01484 OPT::Apply(*BOp_,*this->X_,*BX);
01485 tmp = Utils::errorEquality(*BX, *BX_);
01486 os << " >> Error in BX == B*X : " << scientific << setprecision(10) << tmp << endl;
01487 }
01488
01489
01490 if (chk.checkEta && initialized_) {
01491 tmp = orthman_->orthogError(*X_,*eta_);
01492 os << " >> Error in X^H B Eta == 0 : " << scientific << setprecision(10) << tmp << endl;
01493 for (unsigned int i=0; i<auxVecs_.size(); i++) {
01494 tmp = orthman_->orthogError(*eta_,*auxVecs_[i]);
01495 os << " >> Error in Eta^H B Q[" << i << "]==0 : " << scientific << setprecision(10) << tmp << endl;
01496 }
01497 }
01498 if (chk.checkAEta && !isSkinny_ && initialized_) {
01499 tmp = Utils::errorEquality(*eta_, *Aeta_, AOp_);
01500 os << " >> Error in AEta == A*Eta : " << scientific << setprecision(10) << tmp << endl;
01501 }
01502 if (chk.checkBEta && !isSkinny_ && hasBOp_ && initialized_) {
01503 tmp = Utils::errorEquality(*eta_, *Beta_, BOp_);
01504 os << " >> Error in BEta == B*Eta : " << scientific << setprecision(10) << tmp << endl;
01505 }
01506
01507
01508 if (chk.checkR && initialized_) {
01509 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
01510 MVT::MvTransMv(ONE,*X_,*R_,xTx);
01511 tmp = xTx.normFrobenius();
01512 os << " >> Error in X^H R == 0 : " << scientific << setprecision(10) << tmp << endl;
01513 }
01514
01515
01516
01517 if (chk.checkBR && hasBOp_ && initialized_) {
01518 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
01519 MVT::MvTransMv(ONE,*BX_,*R_,xTx);
01520 tmp = xTx.normFrobenius();
01521 os << " >> Error in X^H B R == 0 : " << scientific << setprecision(10) << tmp << endl;
01522 }
01523
01524
01525 if (chk.checkZ && initialized_) {
01526 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
01527 MVT::MvTransMv(ONE,*BX_,*Z_,xTx);
01528 tmp = xTx.normFrobenius();
01529 os << " >> Error in X^H B Z == 0 : " << scientific << setprecision(10) << tmp << endl;
01530 }
01531
01532
01533 if (chk.checkQ) {
01534 for (unsigned int i=0; i<auxVecs_.size(); i++) {
01535 tmp = orthman_->orthonormError(*auxVecs_[i]);
01536 os << " >> Error in Q[" << i << "]^H B Q[" << i << "]==I: " << scientific << setprecision(10) << tmp << endl;
01537 for (unsigned int j=i+1; j<auxVecs_.size(); j++) {
01538 tmp = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
01539 os << " >> Error in Q[" << i << "]^H B Q[" << j << "]==0: " << scientific << setprecision(10) << tmp << endl;
01540 }
01541 }
01542 }
01543 os << endl;
01544 return os.str();
01545 }
01546
01547
01549
01550 template <class ScalarType, class MV, class OP>
01551 void
01552 RTRBase<ScalarType,MV,OP>::currentStatus(std::ostream &os)
01553 {
01554 using std::setprecision;
01555 using std::scientific;
01556 using std::setw;
01557 using std::endl;
01558
01559 os <<endl;
01560 os <<"================================================================================" << endl;
01561 os << endl;
01562 os <<" RTRBase Solver Status" << endl;
01563 os << endl;
01564 os <<"The solver is "<<(initialized_ ? "initialized." : "not initialized.") << endl;
01565 os <<"The number of iterations performed is " << iter_ << endl;
01566 os <<"The current block size is " << blockSize_ << endl;
01567 os <<"The number of auxiliary vectors is " << numAuxVecs_ << endl;
01568 os <<"The number of operations A*x is " << counterAOp_ << endl;
01569 os <<"The number of operations B*x is " << counterBOp_ << endl;
01570 os <<"The number of operations Prec*x is " << counterPrec_ << endl;
01571 os <<"The most recent rho was " << scientific << setprecision(10) << rho_ << endl;
01572 os <<"The current value of f(x) is " << scientific << setprecision(10) << fx_ << endl;
01573
01574 if (initialized_) {
01575 os << endl;
01576 os <<"CURRENT EIGENVALUE ESTIMATES "<<endl;
01577 os << setw(20) << "Eigenvalue"
01578 << setw(20) << "Residual(B)"
01579 << setw(20) << "Residual(2)"
01580 << endl;
01581 os <<"--------------------------------------------------------------------------------"<<endl;
01582 for (int i=0; i<blockSize_; i++) {
01583 os << scientific << setprecision(10) << setw(20) << theta_[i];
01584 if (Rnorms_current_) os << scientific << setprecision(10) << setw(20) << Rnorms_[i];
01585 else os << scientific << setprecision(10) << setw(20) << "not current";
01586 if (R2norms_current_) os << scientific << setprecision(10) << setw(20) << R2norms_[i];
01587 else os << scientific << setprecision(10) << setw(20) << "not current";
01588 os << endl;
01589 }
01590 }
01591 os <<"================================================================================" << endl;
01592 os << endl;
01593 }
01594
01595
01597
01598 template <class ScalarType, class MV, class OP>
01599 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
01600 RTRBase<ScalarType,MV,OP>::ginner(const MV &xi) const
01601 {
01602 std::vector<MagnitudeType> d(MVT::GetNumberVecs(xi));
01603 MVT::MvNorm(xi,d);
01604 MagnitudeType ret = 0;
01605 for (vecMTiter i=d.begin(); i != d.end(); ++i) {
01606 ret += (*i)*(*i);
01607 }
01608 return ret;
01609 }
01610
01611
01613
01614 template <class ScalarType, class MV, class OP>
01615 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
01616 RTRBase<ScalarType,MV,OP>::ginner(const MV &xi, const MV &zeta) const
01617 {
01618 std::vector<ScalarType> d(MVT::GetNumberVecs(xi));
01619 MVT::MvDot(xi,zeta,d);
01620 return SCT::real(std::accumulate(d.begin(),d.end(),SCT::zero()));
01621 }
01622
01623
01625
01626 template <class ScalarType, class MV, class OP>
01627 void RTRBase<ScalarType,MV,OP>::ginnersep(
01628 const MV &xi,
01629 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &d) const
01630 {
01631 MVT::MvNorm(xi,d);
01632 for (vecMTiter i=d.begin(); i != d.end(); ++i) {
01633 (*i) = (*i)*(*i);
01634 }
01635 }
01636
01637
01639
01640 template <class ScalarType, class MV, class OP>
01641 void RTRBase<ScalarType,MV,OP>::ginnersep(
01642 const MV &xi, const MV &zeta,
01643 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &d) const
01644 {
01645 std::vector<ScalarType> dC(MVT::GetNumberVecs(xi));
01646 MVT::MvDot(xi,zeta,dC);
01647 vecMTiter iR=d.begin();
01648 vecSTiter iS=dC.begin();
01649 for (; iR != d.end(); iR++, iS++) {
01650 (*iR) = SCT::real(*iS);
01651 }
01652 }
01653
01654 template <class ScalarType, class MV, class OP>
01655 Teuchos::Array<Teuchos::RCP<const MV> > RTRBase<ScalarType,MV,OP>::getAuxVecs() const {
01656 return auxVecs_;
01657 }
01658
01659 template <class ScalarType, class MV, class OP>
01660 int RTRBase<ScalarType,MV,OP>::getBlockSize() const {
01661 return(blockSize_);
01662 }
01663
01664 template <class ScalarType, class MV, class OP>
01665 const Eigenproblem<ScalarType,MV,OP>& RTRBase<ScalarType,MV,OP>::getProblem() const {
01666 return(*problem_);
01667 }
01668
01669 template <class ScalarType, class MV, class OP>
01670 int RTRBase<ScalarType,MV,OP>::getMaxSubspaceDim() const {
01671 return blockSize_;
01672 }
01673
01674 template <class ScalarType, class MV, class OP>
01675 int RTRBase<ScalarType,MV,OP>::getCurSubspaceDim() const
01676 {
01677 if (!initialized_) return 0;
01678 return nevLocal_;
01679 }
01680
01681 template <class ScalarType, class MV, class OP>
01682 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
01683 RTRBase<ScalarType,MV,OP>::getRitzRes2Norms()
01684 {
01685 std::vector<MagnitudeType> ret = ritz2norms_;
01686 ret.resize(nevLocal_);
01687 return ret;
01688 }
01689
01690 template <class ScalarType, class MV, class OP>
01691 std::vector<Value<ScalarType> >
01692 RTRBase<ScalarType,MV,OP>::getRitzValues()
01693 {
01694 std::vector<Value<ScalarType> > ret(nevLocal_);
01695 for (int i=0; i<nevLocal_; i++) {
01696 ret[i].realpart = theta_[i];
01697 ret[i].imagpart = ZERO;
01698 }
01699 return ret;
01700 }
01701
01702 template <class ScalarType, class MV, class OP>
01703 Teuchos::RCP<const MV>
01704 RTRBase<ScalarType,MV,OP>::getRitzVectors()
01705 {
01706 return X_;
01707 }
01708
01709 template <class ScalarType, class MV, class OP>
01710 void RTRBase<ScalarType,MV,OP>::resetNumIters()
01711 {
01712 iter_=0;
01713 }
01714
01715 template <class ScalarType, class MV, class OP>
01716 int RTRBase<ScalarType,MV,OP>::getNumIters() const
01717 {
01718 return iter_;
01719 }
01720
01721 template <class ScalarType, class MV, class OP>
01722 RTRState<ScalarType,MV> RTRBase<ScalarType,MV,OP>::getState() const
01723 {
01724 RTRState<ScalarType,MV> state;
01725 state.X = X_;
01726 state.AX = AX_;
01727 if (hasBOp_) {
01728 state.BX = BX_;
01729 }
01730 else {
01731 state.BX = Teuchos::null;
01732 }
01733 state.rho = rho_;
01734 state.R = R_;
01735 state.T = Teuchos::rcp(new std::vector<MagnitudeType>(theta_));
01736 return state;
01737 }
01738
01739 template <class ScalarType, class MV, class OP>
01740 bool RTRBase<ScalarType,MV,OP>::isInitialized() const
01741 {
01742 return initialized_;
01743 }
01744
01745 template <class ScalarType, class MV, class OP>
01746 std::vector<int> RTRBase<ScalarType,MV,OP>::getRitzIndex()
01747 {
01748 std::vector<int> ret(nevLocal_,0);
01749 return ret;
01750 }
01751
01752
01753 }
01754
01755 #endif // ANASAZI_RTRBASE_HPP