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 #ifndef ANASAZI_BLOCK_DAVIDSON_HPP
00034 #define ANASAZI_BLOCK_DAVIDSON_HPP
00035
00036 #include "AnasaziTypes.hpp"
00037
00038 #include "AnasaziEigensolver.hpp"
00039 #include "AnasaziMultiVecTraits.hpp"
00040 #include "AnasaziOperatorTraits.hpp"
00041 #include "Teuchos_ScalarTraits.hpp"
00042
00043 #include "AnasaziMatOrthoManager.hpp"
00044 #include "AnasaziSolverUtils.hpp"
00045
00046 #include "Teuchos_LAPACK.hpp"
00047 #include "Teuchos_BLAS.hpp"
00048 #include "Teuchos_SerialDenseMatrix.hpp"
00049 #include "Teuchos_ParameterList.hpp"
00050 #include "Teuchos_TimeMonitor.hpp"
00051
00067 namespace Anasazi {
00068
00070
00071
00076 template <class ScalarType, class MV>
00077 struct BlockDavidsonState {
00082 int curDim;
00087 Teuchos::RCP<const MV> V;
00089 Teuchos::RCP<const MV> X;
00091 Teuchos::RCP<const MV> KX;
00093 Teuchos::RCP<const MV> MX;
00095 Teuchos::RCP<const MV> R;
00100 Teuchos::RCP<const MV> H;
00102 Teuchos::RCP<const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > T;
00108 Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > KK;
00109 BlockDavidsonState() : curDim(0), V(Teuchos::null),
00110 X(Teuchos::null), KX(Teuchos::null), MX(Teuchos::null),
00111 R(Teuchos::null), H(Teuchos::null),
00112 T(Teuchos::null), KK(Teuchos::null) {}
00113 };
00114
00116
00118
00119
00132 class BlockDavidsonInitFailure : public AnasaziError {public:
00133 BlockDavidsonInitFailure(const std::string& what_arg) : AnasaziError(what_arg)
00134 {}};
00135
00143 class BlockDavidsonOrthoFailure : public AnasaziError {public:
00144 BlockDavidsonOrthoFailure(const std::string& what_arg) : AnasaziError(what_arg)
00145 {}};
00146
00148
00149
00150 template <class ScalarType, class MV, class OP>
00151 class BlockDavidson : public Eigensolver<ScalarType,MV,OP> {
00152 public:
00154
00155
00163 BlockDavidson( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00164 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
00165 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00166 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00167 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00168 Teuchos::ParameterList ¶ms
00169 );
00170
00172 virtual ~BlockDavidson();
00174
00175
00177
00178
00202 void iterate();
00203
00237 void initialize(BlockDavidsonState<ScalarType,MV> newstate);
00238
00242 void initialize();
00243
00259 bool isInitialized() const;
00260
00273 BlockDavidsonState<ScalarType,MV> getState() const;
00274
00276
00277
00279
00280
00282 int getNumIters() const;
00283
00285 void resetNumIters();
00286
00294 Teuchos::RCP<const MV> getRitzVectors();
00295
00301 std::vector<Value<ScalarType> > getRitzValues();
00302
00303
00311 std::vector<int> getRitzIndex();
00312
00313
00319 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getResNorms();
00320
00321
00327 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRes2Norms();
00328
00329
00337 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRitzRes2Norms();
00338
00344 int getCurSubspaceDim() const;
00345
00347 int getMaxSubspaceDim() const;
00348
00350
00351
00353
00354
00356 void setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test);
00357
00359 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > getStatusTest() const;
00360
00362 const Eigenproblem<ScalarType,MV,OP>& getProblem() const;
00363
00373 void setBlockSize(int blockSize);
00374
00376 int getBlockSize() const;
00377
00390 void setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs);
00391
00393 Teuchos::Array<Teuchos::RCP<const MV> > getAuxVecs() const;
00394
00396
00398
00399
00409 void setSize(int blockSize, int numBlocks);
00410
00412
00414
00415
00417 void currentStatus(std::ostream &os);
00418
00420
00421 private:
00422
00423
00424
00425 typedef SolverUtils<ScalarType,MV,OP> Utils;
00426 typedef MultiVecTraits<ScalarType,MV> MVT;
00427 typedef OperatorTraits<ScalarType,MV,OP> OPT;
00428 typedef Teuchos::ScalarTraits<ScalarType> SCT;
00429 typedef typename SCT::magnitudeType MagnitudeType;
00430 const MagnitudeType ONE;
00431 const MagnitudeType ZERO;
00432 const MagnitudeType NANVAL;
00433
00434
00435
00436 struct CheckList {
00437 bool checkV;
00438 bool checkX, checkMX, checkKX;
00439 bool checkH, checkMH, checkKH;
00440 bool checkR, checkQ;
00441 bool checkKK;
00442 CheckList() : checkV(false),
00443 checkX(false),checkMX(false),checkKX(false),
00444 checkH(false),checkMH(false),checkKH(false),
00445 checkR(false),checkQ(false),checkKK(false) {};
00446 };
00447
00448
00449
00450 std::string accuracyCheck(const CheckList &chk, const std::string &where) const;
00451
00452
00453
00454 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
00455 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sm_;
00456 const Teuchos::RCP<OutputManager<ScalarType> > om_;
00457 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > tester_;
00458 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > orthman_;
00459
00460
00461
00462 Teuchos::RCP<const OP> Op_;
00463 Teuchos::RCP<const OP> MOp_;
00464 Teuchos::RCP<const OP> Prec_;
00465 bool hasM_;
00466
00467
00468
00469 Teuchos::RCP<Teuchos::Time> timerOp_, timerMOp_, timerPrec_,
00470 timerSortEval_, timerDS_,
00471 timerLocal_, timerCompRes_,
00472 timerOrtho_, timerInit_;
00473
00474
00475
00476 int count_ApplyOp_, count_ApplyM_, count_ApplyPrec_;
00477
00478
00479
00480
00481
00482
00483
00484 int blockSize_;
00485
00486 int numBlocks_;
00487
00488
00489
00490
00491
00492
00493
00494 bool initialized_;
00495
00496
00497
00498
00499 int curDim_;
00500
00501
00502
00503
00504
00505 Teuchos::RCP<MV> X_, KX_, MX_, R_,
00506 H_, KH_, MH_,
00507 V_;
00508
00509
00510
00511 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > KK_;
00512
00513
00514 Teuchos::Array<Teuchos::RCP<const MV> > auxVecs_;
00515 int numAuxVecs_;
00516
00517
00518 int iter_;
00519
00520
00521 std::vector<MagnitudeType> theta_, Rnorms_, R2norms_;
00522
00523
00524 bool Rnorms_current_, R2norms_current_;
00525
00526 };
00527
00530
00531
00532
00535
00536
00538
00539 template <class ScalarType, class MV, class OP>
00540 BlockDavidson<ScalarType,MV,OP>::BlockDavidson(
00541 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00542 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
00543 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00544 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00545 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00546 Teuchos::ParameterList ¶ms
00547 ) :
00548 ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
00549 ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
00550 NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan()),
00551
00552 problem_(problem),
00553 sm_(sorter),
00554 om_(printer),
00555 tester_(tester),
00556 orthman_(ortho),
00557
00558 timerOp_(Teuchos::TimeMonitor::getNewTimer("BlockDavidson::Operation Op*x")),
00559 timerMOp_(Teuchos::TimeMonitor::getNewTimer("BlockDavidson::Operation M*x")),
00560 timerPrec_(Teuchos::TimeMonitor::getNewTimer("BlockDavidson::Operation Prec*x")),
00561 timerSortEval_(Teuchos::TimeMonitor::getNewTimer("BlockDavidson::Sorting eigenvalues")),
00562 timerDS_(Teuchos::TimeMonitor::getNewTimer("BlockDavidson::Direct solve")),
00563 timerLocal_(Teuchos::TimeMonitor::getNewTimer("BlockDavidson::Local update")),
00564 timerCompRes_(Teuchos::TimeMonitor::getNewTimer("BlockDavidson::Computing residuals")),
00565 timerOrtho_(Teuchos::TimeMonitor::getNewTimer("BlockDavidson::Orthogonalization")),
00566 timerInit_(Teuchos::TimeMonitor::getNewTimer("BlockDavidson::Initialization")),
00567 count_ApplyOp_(0),
00568 count_ApplyM_(0),
00569 count_ApplyPrec_(0),
00570
00571 blockSize_(0),
00572 numBlocks_(0),
00573 initialized_(false),
00574 curDim_(0),
00575 auxVecs_( Teuchos::Array<Teuchos::RCP<const MV> >(0) ),
00576 numAuxVecs_(0),
00577 iter_(0),
00578 Rnorms_current_(false),
00579 R2norms_current_(false)
00580 {
00581 TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument,
00582 "Anasazi::BlockDavidson::constructor: user passed null problem pointer.");
00583 TEST_FOR_EXCEPTION(sm_ == Teuchos::null,std::invalid_argument,
00584 "Anasazi::BlockDavidson::constructor: user passed null sort manager pointer.");
00585 TEST_FOR_EXCEPTION(om_ == Teuchos::null,std::invalid_argument,
00586 "Anasazi::BlockDavidson::constructor: user passed null output manager pointer.");
00587 TEST_FOR_EXCEPTION(tester_ == Teuchos::null,std::invalid_argument,
00588 "Anasazi::BlockDavidson::constructor: user passed null status test pointer.");
00589 TEST_FOR_EXCEPTION(orthman_ == Teuchos::null,std::invalid_argument,
00590 "Anasazi::BlockDavidson::constructor: user passed null orthogonalization manager pointer.");
00591 TEST_FOR_EXCEPTION(problem_->isProblemSet() == false, std::invalid_argument,
00592 "Anasazi::BlockDavidson::constructor: problem is not set.");
00593 TEST_FOR_EXCEPTION(problem_->isHermitian() == false, std::invalid_argument,
00594 "Anasazi::BlockDavidson::constructor: problem is not hermitian.");
00595
00596
00597 Op_ = problem_->getOperator();
00598 TEST_FOR_EXCEPTION(Op_ == Teuchos::null, std::invalid_argument,
00599 "Anasazi::BlockDavidson::constructor: problem provides no operator.");
00600 MOp_ = problem_->getM();
00601 Prec_ = problem_->getPrec();
00602 hasM_ = (MOp_ != Teuchos::null);
00603
00604
00605 int bs = params.get("Block Size", problem_->getNEV());
00606 int nb = params.get("Num Blocks", 2);
00607 setSize(bs,nb);
00608 }
00609
00610
00612
00613 template <class ScalarType, class MV, class OP>
00614 BlockDavidson<ScalarType,MV,OP>::~BlockDavidson() {}
00615
00616
00618
00619
00620 template <class ScalarType, class MV, class OP>
00621 void BlockDavidson<ScalarType,MV,OP>::setBlockSize (int blockSize)
00622 {
00623 setSize(blockSize,numBlocks_);
00624 }
00625
00626
00628
00629 template <class ScalarType, class MV, class OP>
00630 Teuchos::Array<Teuchos::RCP<const MV> > BlockDavidson<ScalarType,MV,OP>::getAuxVecs() const {
00631 return auxVecs_;
00632 }
00633
00634
00635
00637
00638 template <class ScalarType, class MV, class OP>
00639 int BlockDavidson<ScalarType,MV,OP>::getBlockSize() const {
00640 return(blockSize_);
00641 }
00642
00643
00645
00646 template <class ScalarType, class MV, class OP>
00647 const Eigenproblem<ScalarType,MV,OP>& BlockDavidson<ScalarType,MV,OP>::getProblem() const {
00648 return(*problem_);
00649 }
00650
00651
00653
00654 template <class ScalarType, class MV, class OP>
00655 int BlockDavidson<ScalarType,MV,OP>::getMaxSubspaceDim() const {
00656 return blockSize_*numBlocks_;
00657 }
00658
00660
00661 template <class ScalarType, class MV, class OP>
00662 int BlockDavidson<ScalarType,MV,OP>::getCurSubspaceDim() const {
00663 if (!initialized_) return 0;
00664 return curDim_;
00665 }
00666
00667
00669
00670 template <class ScalarType, class MV, class OP>
00671 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
00672 BlockDavidson<ScalarType,MV,OP>::getRitzRes2Norms() {
00673 return this->getRes2Norms();
00674 }
00675
00676
00678
00679 template <class ScalarType, class MV, class OP>
00680 std::vector<int> BlockDavidson<ScalarType,MV,OP>::getRitzIndex() {
00681 std::vector<int> ret(curDim_,0);
00682 return ret;
00683 }
00684
00685
00687
00688 template <class ScalarType, class MV, class OP>
00689 std::vector<Value<ScalarType> > BlockDavidson<ScalarType,MV,OP>::getRitzValues() {
00690 std::vector<Value<ScalarType> > ret(curDim_);
00691 for (int i=0; i<curDim_; ++i) {
00692 ret[i].realpart = theta_[i];
00693 ret[i].imagpart = ZERO;
00694 }
00695 return ret;
00696 }
00697
00698
00700
00701 template <class ScalarType, class MV, class OP>
00702 Teuchos::RCP<const MV> BlockDavidson<ScalarType,MV,OP>::getRitzVectors() {
00703 return X_;
00704 }
00705
00706
00708
00709 template <class ScalarType, class MV, class OP>
00710 void BlockDavidson<ScalarType,MV,OP>::resetNumIters() {
00711 iter_=0;
00712 }
00713
00714
00716
00717 template <class ScalarType, class MV, class OP>
00718 int BlockDavidson<ScalarType,MV,OP>::getNumIters() const {
00719 return(iter_);
00720 }
00721
00722
00724
00725 template <class ScalarType, class MV, class OP>
00726 BlockDavidsonState<ScalarType,MV> BlockDavidson<ScalarType,MV,OP>::getState() const {
00727 BlockDavidsonState<ScalarType,MV> state;
00728 state.curDim = curDim_;
00729 state.V = V_;
00730 state.X = X_;
00731 state.KX = KX_;
00732 if (hasM_) {
00733 state.MX = MX_;
00734 }
00735 else {
00736 state.MX = Teuchos::null;
00737 }
00738 state.R = R_;
00739 state.H = H_;
00740 state.KK = KK_;
00741 if (curDim_ > 0) {
00742 state.T = Teuchos::rcp(new std::vector<MagnitudeType>(theta_.begin(),theta_.begin()+curDim_) );
00743 } else {
00744 state.T = Teuchos::rcp(new std::vector<MagnitudeType>(0));
00745 }
00746 return state;
00747 }
00748
00749
00751
00752 template <class ScalarType, class MV, class OP>
00753 bool BlockDavidson<ScalarType,MV,OP>::isInitialized() const { return initialized_; }
00754
00755
00757
00758 template <class ScalarType, class MV, class OP>
00759 void BlockDavidson<ScalarType,MV,OP>::setSize (int blockSize, int numBlocks)
00760 {
00761
00762 Teuchos::TimeMonitor initimer( *timerInit_ );
00763
00764
00765
00766
00767 TEST_FOR_EXCEPTION(blockSize < 1, std::invalid_argument, "Anasazi::BlockDavidson::setSize(blocksize,numblocks): blocksize must be strictly positive.");
00768 TEST_FOR_EXCEPTION(numBlocks < 2, std::invalid_argument, "Anasazi::BlockDavidson::setSize(blocksize,numblocks): numblocks must be greater than one.");
00769 if (blockSize == blockSize_ && numBlocks == numBlocks_) {
00770
00771 return;
00772 }
00773
00774 blockSize_ = blockSize;
00775 numBlocks_ = numBlocks;
00776
00777 Teuchos::RCP<const MV> tmp;
00778
00779
00780
00781
00782 if (X_ != Teuchos::null) {
00783 tmp = X_;
00784 }
00785 else {
00786 tmp = problem_->getInitVec();
00787 TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
00788 "Anasazi::BlockDavidson::setSize(): eigenproblem did not specify initial vectors to clone from.");
00789 }
00790
00791 TEST_FOR_EXCEPTION(numAuxVecs_+blockSize*numBlocks > MVT::GetVecLength(*tmp),std::invalid_argument,
00792 "Anasazi::BlockDavidson::setSize(): max subspace dimension and auxilliary subspace too large.");
00793
00794
00796
00797
00798
00799 Rnorms_.resize(blockSize_,NANVAL);
00800 R2norms_.resize(blockSize_,NANVAL);
00801
00802
00803
00804
00805 X_ = Teuchos::null;
00806 KX_ = Teuchos::null;
00807 MX_ = Teuchos::null;
00808 R_ = Teuchos::null;
00809 V_ = Teuchos::null;
00810
00811 om_->print(Debug," >> Allocating X_\n");
00812 X_ = MVT::Clone(*tmp,blockSize_);
00813 om_->print(Debug," >> Allocating KX_\n");
00814 KX_ = MVT::Clone(*tmp,blockSize_);
00815 if (hasM_) {
00816 om_->print(Debug," >> Allocating MX_\n");
00817 MX_ = MVT::Clone(*tmp,blockSize_);
00818 }
00819 else {
00820 MX_ = X_;
00821 }
00822 om_->print(Debug," >> Allocating R_\n");
00823 R_ = MVT::Clone(*tmp,blockSize_);
00824
00826
00827
00828 int newsd = blockSize_*numBlocks_;
00829 theta_.resize(blockSize_*numBlocks_,NANVAL);
00830 om_->print(Debug," >> Allocating V_\n");
00831 V_ = MVT::Clone(*tmp,newsd);
00832 KK_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(newsd,newsd) );
00833
00834 om_->print(Debug," >> done allocating.\n");
00835
00836 initialized_ = false;
00837 curDim_ = 0;
00838 }
00839
00840
00842
00843 template <class ScalarType, class MV, class OP>
00844 void BlockDavidson<ScalarType,MV,OP>::setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs) {
00845 typedef typename Teuchos::Array<Teuchos::RCP<const MV> >::iterator tarcpmv;
00846
00847
00848 auxVecs_ = auxvecs;
00849 numAuxVecs_ = 0;
00850 for (tarcpmv i=auxVecs_.begin(); i != auxVecs_.end(); ++i) {
00851 numAuxVecs_ += MVT::GetNumberVecs(**i);
00852 }
00853
00854
00855 if (numAuxVecs_ > 0 && initialized_) {
00856 initialized_ = false;
00857 }
00858
00859 if (om_->isVerbosity( Debug ) ) {
00860 CheckList chk;
00861 chk.checkQ = true;
00862 om_->print( Debug, accuracyCheck(chk, ": in setAuxVecs()") );
00863 }
00864 }
00865
00866
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880 template <class ScalarType, class MV, class OP>
00881 void BlockDavidson<ScalarType,MV,OP>::initialize(BlockDavidsonState<ScalarType,MV> newstate)
00882 {
00883
00884
00885
00886 Teuchos::TimeMonitor inittimer( *timerInit_ );
00887
00888 std::vector<int> bsind(blockSize_);
00889 for (int i=0; i<blockSize_; ++i) bsind[i] = i;
00890
00891 Teuchos::BLAS<int,ScalarType> blas;
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915 Teuchos::RCP<MV> lclV;
00916 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclKK;
00917
00918 if (newstate.V != Teuchos::null && newstate.KK != Teuchos::null) {
00919 TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.V) != MVT::GetVecLength(*V_), std::invalid_argument,
00920 "Anasazi::BlockDavidson::initialize(newstate): Vector length of V not correct." );
00921 TEST_FOR_EXCEPTION( newstate.curDim < blockSize_, std::invalid_argument,
00922 "Anasazi::BlockDavidson::initialize(newstate): Rank of new state must be at least blockSize().");
00923 TEST_FOR_EXCEPTION( newstate.curDim > blockSize_*numBlocks_, std::invalid_argument,
00924 "Anasazi::BlockDavidson::initialize(newstate): Rank of new state must be less than getMaxSubspaceDim().");
00925 TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < newstate.curDim, std::invalid_argument,
00926 "Anasazi::BlockDavidson::initialize(newstate): Multivector for basis in new state must be as large as specified state rank.");
00927
00928 curDim_ = newstate.curDim;
00929
00930 curDim_ = (int)(curDim_ / blockSize_)*blockSize_;
00931
00932 TEST_FOR_EXCEPTION( curDim_ != newstate.curDim, std::invalid_argument,
00933 "Anasazi::BlockDavidson::initialize(newstate): Rank of new state must be a multiple of getBlockSize().");
00934
00935
00936 TEST_FOR_EXCEPTION( newstate.KK->numRows() < curDim_ || newstate.KK->numCols() < curDim_, std::invalid_argument,
00937 "Anasazi::BlockDavidson::initialize(newstate): Projected matrix in new state must be as large as specified state rank.");
00938
00939
00940 std::vector<int> nevind(curDim_);
00941 for (int i=0; i<curDim_; ++i) nevind[i] = i;
00942 if (newstate.V != V_) {
00943 if (curDim_ < MVT::GetNumberVecs(*newstate.V)) {
00944 newstate.V = MVT::CloneView(*newstate.V,nevind);
00945 }
00946 MVT::SetBlock(*newstate.V,nevind,*V_);
00947 }
00948 lclV = MVT::CloneView(*V_,nevind);
00949
00950
00951 lclKK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
00952 if (newstate.KK != KK_) {
00953 if (newstate.KK->numRows() > curDim_ || newstate.KK->numCols() > curDim_) {
00954 newstate.KK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*newstate.KK,curDim_,curDim_) );
00955 }
00956 lclKK->assign(*newstate.KK);
00957 }
00958
00959
00960 for (int j=0; j<curDim_-1; ++j) {
00961 for (int i=j+1; i<curDim_; ++i) {
00962 (*lclKK)(i,j) = SCT::conjugate((*lclKK)(j,i));
00963 }
00964 }
00965 }
00966 else {
00967
00968
00969 Teuchos::RCP<const MV> ivec = problem_->getInitVec();
00970 TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::invalid_argument,
00971 "Anasazi::BlockDavdison::initialize(newstate): Eigenproblem did not specify initial vectors to clone from.");
00972
00973 newstate.X = Teuchos::null;
00974 newstate.MX = Teuchos::null;
00975 newstate.KX = Teuchos::null;
00976 newstate.R = Teuchos::null;
00977 newstate.H = Teuchos::null;
00978 newstate.T = Teuchos::null;
00979 newstate.KK = Teuchos::null;
00980 newstate.V = Teuchos::null;
00981 newstate.curDim = 0;
00982
00983 curDim_ = MVT::GetNumberVecs(*ivec);
00984
00985 curDim_ = (int)(curDim_ / blockSize_)*blockSize_;
00986 if (curDim_ > blockSize_*numBlocks_) {
00987
00988
00989 curDim_ = blockSize_*numBlocks_;
00990 }
00991 bool userand = false;
00992 if (curDim_ == 0) {
00993
00994
00995 userand = true;
00996 curDim_ = blockSize_;
00997 }
00998
00999
01000
01001
01002
01003
01004
01005
01006 std::vector<int> dimind(curDim_);
01007 for (int i=0; i<curDim_; ++i) dimind[i] = i;
01008 lclV = MVT::CloneView(*V_,dimind);
01009 if (userand) {
01010
01011 MVT::MvRandom(*lclV);
01012 }
01013 else {
01014 if (MVT::GetNumberVecs(*ivec) > curDim_) {
01015 ivec = MVT::CloneView(*ivec,dimind);
01016 }
01017
01018 MVT::SetBlock(*ivec,dimind,*lclV);
01019 }
01020 Teuchos::RCP<MV> tmpVecs;
01021 if (curDim_*2 <= blockSize_*numBlocks_) {
01022
01023 std::vector<int> block2(curDim_);
01024 for (int i=0; i<curDim_; ++i) block2[i] = curDim_+i;
01025 tmpVecs = MVT::CloneView(*V_,block2);
01026 }
01027 else {
01028
01029 tmpVecs = MVT::Clone(*V_,curDim_);
01030 }
01031
01032
01033 if (hasM_) {
01034 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
01035 OPT::Apply(*MOp_,*lclV,*tmpVecs);
01036 count_ApplyM_ += curDim_;
01037 }
01038
01039
01040 if (auxVecs_.size() > 0) {
01041 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
01042
01043 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
01044 int rank = orthman_->projectAndNormalizeMat(*lclV,auxVecs_,dummyC,Teuchos::null,tmpVecs);
01045 TEST_FOR_EXCEPTION(rank != curDim_,BlockDavidsonInitFailure,
01046 "Anasazi::BlockDavidson::initialize(): Couldn't generate initial basis of full rank.");
01047 }
01048 else {
01049 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
01050
01051 int rank = orthman_->normalizeMat(*lclV,Teuchos::null,tmpVecs);
01052 TEST_FOR_EXCEPTION(rank != curDim_,BlockDavidsonInitFailure,
01053 "Anasazi::BlockDavidson::initialize(): Couldn't generate initial basis of full rank.");
01054 }
01055
01056
01057 {
01058 Teuchos::TimeMonitor lcltimer( *timerOp_ );
01059 OPT::Apply(*Op_,*lclV,*tmpVecs);
01060 count_ApplyOp_ += curDim_;
01061 }
01062
01063
01064 lclKK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
01065 MVT::MvTransMv(ONE,*lclV,*tmpVecs,*lclKK);
01066
01067
01068 tmpVecs = Teuchos::null;
01069 }
01070
01071
01072 if (newstate.X != Teuchos::null && newstate.T != Teuchos::null) {
01073 TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.X) != blockSize_ || MVT::GetVecLength(*newstate.X) != MVT::GetVecLength(*X_),
01074 std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): Size of X must be consistent with block size and length of V.");
01075 TEST_FOR_EXCEPTION((signed int)(newstate.T->size()) != curDim_,
01076 std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): Size of T must be consistent with dimension of V.");
01077
01078 if (newstate.X != X_) {
01079 MVT::SetBlock(*newstate.X,bsind,*X_);
01080 }
01081
01082 std::copy(newstate.T->begin(),newstate.T->end(),theta_.begin());
01083 }
01084 else {
01085
01086 Teuchos::SerialDenseMatrix<int,ScalarType> S(curDim_,curDim_);
01087 {
01088 Teuchos::TimeMonitor lcltimer( *timerDS_ );
01089 int rank = curDim_;
01090 Utils::directSolver(curDim_, *lclKK, Teuchos::null, S, theta_, rank, 10);
01091
01092 TEST_FOR_EXCEPTION(rank != curDim_,BlockDavidsonInitFailure,
01093 "Anasazi::BlockDavidson::initialize(newstate): Not enough Ritz vectors to initialize algorithm.");
01094 }
01095
01096 {
01097 Teuchos::TimeMonitor lcltimer( *timerSortEval_ );
01098
01099 std::vector<int> order(curDim_);
01100
01101
01102 sm_->sort(theta_, Teuchos::rcpFromRef(order), curDim_);
01103
01104
01105 Utils::permuteVectors(order,S);
01106 }
01107
01108
01109 Teuchos::SerialDenseMatrix<int,ScalarType> S1(Teuchos::View,S,curDim_,blockSize_);
01110 {
01111 Teuchos::TimeMonitor lcltimer( *timerLocal_ );
01112
01113
01114 MVT::MvTimesMatAddMv( ONE, *lclV, S1, ZERO, *X_ );
01115 }
01116
01117 newstate.KX = Teuchos::null;
01118 newstate.MX = Teuchos::null;
01119 }
01120
01121
01122 lclV = Teuchos::null;
01123 lclKK = Teuchos::null;
01124
01125
01126 if ( newstate.KX != Teuchos::null ) {
01127 TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.KX) != blockSize_,
01128 std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): vector length of newstate.KX not correct." );
01129 TEST_FOR_EXCEPTION(MVT::GetVecLength(*newstate.KX) != MVT::GetVecLength(*X_),
01130 std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): newstate.KX must have at least block size vectors." );
01131 if (newstate.KX != KX_) {
01132 MVT::SetBlock(*newstate.KX,bsind,*KX_);
01133 }
01134 }
01135 else {
01136
01137 {
01138 Teuchos::TimeMonitor lcltimer( *timerOp_ );
01139 OPT::Apply(*Op_,*X_,*KX_);
01140 count_ApplyOp_ += blockSize_;
01141 }
01142
01143 newstate.R = Teuchos::null;
01144 }
01145
01146
01147 if (hasM_) {
01148 if ( newstate.MX != Teuchos::null ) {
01149 TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.MX) != blockSize_,
01150 std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): vector length of newstate.MX not correct." );
01151 TEST_FOR_EXCEPTION(MVT::GetVecLength(*newstate.MX) != MVT::GetVecLength(*X_),
01152 std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): newstate.MX must have at least block size vectors." );
01153 if (newstate.MX != MX_) {
01154 MVT::SetBlock(*newstate.MX,bsind,*MX_);
01155 }
01156 }
01157 else {
01158
01159 {
01160 Teuchos::TimeMonitor lcltimer( *timerOp_ );
01161 OPT::Apply(*MOp_,*X_,*MX_);
01162 count_ApplyOp_ += blockSize_;
01163 }
01164
01165 newstate.R = Teuchos::null;
01166 }
01167 }
01168 else {
01169
01170 TEST_FOR_EXCEPTION(MX_ != X_, std::logic_error, "Anasazi::BlockDavidson::initialize(): solver invariant not satisfied (MX==X).");
01171 }
01172
01173
01174 if (newstate.R != Teuchos::null) {
01175 TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.R) != blockSize_,
01176 std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): vector length of newstate.R not correct." );
01177 TEST_FOR_EXCEPTION(MVT::GetVecLength(*newstate.R) != MVT::GetVecLength(*X_),
01178 std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): newstate.R must have at least block size vectors." );
01179 if (newstate.R != R_) {
01180 MVT::SetBlock(*newstate.R,bsind,*R_);
01181 }
01182 }
01183 else {
01184 Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
01185
01186
01187 MVT::MvAddMv(ZERO,*KX_,ONE,*KX_,*R_);
01188 Teuchos::SerialDenseMatrix<int,ScalarType> T(blockSize_,blockSize_);
01189 T.putScalar(ZERO);
01190 for (int i=0; i<blockSize_; ++i) T(i,i) = theta_[i];
01191 MVT::MvTimesMatAddMv(-ONE,*MX_,T,ONE,*R_);
01192
01193 }
01194
01195
01196 Rnorms_current_ = false;
01197 R2norms_current_ = false;
01198
01199
01200 initialized_ = true;
01201
01202 if (om_->isVerbosity( Debug ) ) {
01203
01204 CheckList chk;
01205 chk.checkV = true;
01206 chk.checkX = true;
01207 chk.checkKX = true;
01208 chk.checkMX = true;
01209 chk.checkR = true;
01210 chk.checkQ = true;
01211 chk.checkKK = true;
01212 om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
01213 }
01214
01215
01216 if (om_->isVerbosity(Debug)) {
01217 currentStatus( om_->stream(Debug) );
01218 }
01219 else if (om_->isVerbosity(IterationDetails)) {
01220 currentStatus( om_->stream(IterationDetails) );
01221 }
01222 }
01223
01224
01226
01227 template <class ScalarType, class MV, class OP>
01228 void BlockDavidson<ScalarType,MV,OP>::initialize()
01229 {
01230 BlockDavidsonState<ScalarType,MV> empty;
01231 initialize(empty);
01232 }
01233
01234
01236
01237 template <class ScalarType, class MV, class OP>
01238 void BlockDavidson<ScalarType,MV,OP>::iterate ()
01239 {
01240
01241
01242 if (initialized_ == false) {
01243 initialize();
01244 }
01245
01246
01247
01248 const int searchDim = blockSize_*numBlocks_;
01249
01250 Teuchos::BLAS<int,ScalarType> blas;
01251
01252
01253
01254
01255 Teuchos::SerialDenseMatrix<int,ScalarType> S( searchDim, searchDim );
01256
01257
01259
01260
01261 while (tester_->checkStatus(this) != Passed && curDim_ < searchDim) {
01262
01263
01264 if (om_->isVerbosity(Debug)) {
01265 currentStatus( om_->stream(Debug) );
01266 }
01267 else if (om_->isVerbosity(IterationDetails)) {
01268 currentStatus( om_->stream(IterationDetails) );
01269 }
01270
01271 ++iter_;
01272
01273
01274 std::vector<int> curind(blockSize_);
01275 for (int i=0; i<blockSize_; ++i) curind[i] = curDim_ + i;
01276 H_ = MVT::CloneView(*V_,curind);
01277
01278
01279
01280 if (Prec_ != Teuchos::null) {
01281 Teuchos::TimeMonitor lcltimer( *timerPrec_ );
01282 OPT::Apply( *Prec_, *R_, *H_ );
01283 count_ApplyPrec_ += blockSize_;
01284 }
01285 else {
01286 std::vector<int> bsind(blockSize_);
01287 for (int i=0; i<blockSize_; ++i) { bsind[i] = i; }
01288 MVT::SetBlock(*R_,bsind,*H_);
01289 }
01290
01291
01292 if (hasM_) {
01293
01294 MH_ = MX_;
01295 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
01296 OPT::Apply( *MOp_, *H_, *MH_);
01297 count_ApplyM_ += blockSize_;
01298 }
01299 else {
01300 MH_ = H_;
01301 }
01302
01303
01304
01305 std::vector<int> prevind(curDim_);
01306 for (int i=0; i<curDim_; ++i) prevind[i] = i;
01307 Teuchos::RCP<MV> Vprev = MVT::CloneView(*V_,prevind);
01308
01309
01310 {
01311 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
01312
01313 Teuchos::Array<Teuchos::RCP<const MV> > against = auxVecs_;
01314 against.push_back(Vprev);
01315 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
01316 int rank = orthman_->projectAndNormalizeMat(*H_,against,
01317 dummyC,
01318 Teuchos::null,MH_);
01319 TEST_FOR_EXCEPTION(rank != blockSize_,BlockDavidsonOrthoFailure,
01320 "Anasazi::BlockDavidson::iterate(): unable to compute orthonormal basis for H.");
01321 }
01322
01323
01324 {
01325
01326 KH_ = KX_;
01327 Teuchos::TimeMonitor lcltimer( *timerOp_ );
01328 OPT::Apply( *Op_, *H_, *KH_);
01329 count_ApplyOp_ += blockSize_;
01330 }
01331
01332 if (om_->isVerbosity( Debug ) ) {
01333 CheckList chk;
01334 chk.checkH = true;
01335 chk.checkMH = true;
01336 chk.checkKH = true;
01337 om_->print( Debug, accuracyCheck(chk, ": after ortho H") );
01338 }
01339 else if (om_->isVerbosity( OrthoDetails ) ) {
01340 CheckList chk;
01341 chk.checkH = true;
01342 chk.checkMH = true;
01343 chk.checkKH = true;
01344 om_->print( OrthoDetails, accuracyCheck(chk,": after ortho H") );
01345 }
01346
01347
01348
01349 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > nextKK;
01350
01351 nextKK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,blockSize_,0,curDim_) );
01352 MVT::MvTransMv(ONE,*Vprev,*KH_,*nextKK);
01353
01354 nextKK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,blockSize_,blockSize_,curDim_,curDim_) );
01355 MVT::MvTransMv(ONE,*H_,*KH_,*nextKK);
01356
01357
01358 nextKK = Teuchos::null;
01359 for (int i=curDim_; i<curDim_+blockSize_; ++i) {
01360 for (int j=0; j<i; ++j) {
01361 (*KK_)(i,j) = SCT::conjugate((*KK_)(j,i));
01362 }
01363 }
01364
01365
01366 curDim_ += blockSize_;
01367 H_ = KH_ = MH_ = Teuchos::null;
01368 Vprev = Teuchos::null;
01369
01370 if (om_->isVerbosity( Debug ) ) {
01371 CheckList chk;
01372 chk.checkKK = true;
01373 om_->print( Debug, accuracyCheck(chk, ": after expanding KK") );
01374 }
01375
01376
01377 curind.resize(curDim_);
01378 for (int i=0; i<curDim_; ++i) curind[i] = i;
01379 Teuchos::RCP<const MV> curV = MVT::CloneView(*V_,curind);
01380
01381
01382 {
01383 Teuchos::TimeMonitor lcltimer(*timerDS_);
01384 int nevlocal = curDim_;
01385 int info = Utils::directSolver(curDim_,*KK_,Teuchos::null,S,theta_,nevlocal,10);
01386 TEST_FOR_EXCEPTION(info != 0,std::logic_error,"Anasazi::BlockDavidson::iterate(): direct solve returned error code.");
01387
01388 TEST_FOR_EXCEPTION(nevlocal != curDim_,std::logic_error,"Anasazi::BlockDavidson::iterate(): direct solve did not compute all eigenvectors.");
01389 }
01390
01391
01392 {
01393 Teuchos::TimeMonitor lcltimer( *timerSortEval_ );
01394
01395 std::vector<int> order(curDim_);
01396
01397
01398 sm_->sort(theta_, Teuchos::rcp(&order,false), curDim_);
01399
01400
01401 Teuchos::SerialDenseMatrix<int,ScalarType> curS(Teuchos::View,S,curDim_,curDim_);
01402 Utils::permuteVectors(order,curS);
01403 }
01404
01405
01406 Teuchos::SerialDenseMatrix<int,ScalarType> S1( Teuchos::View, S, curDim_, blockSize_ );
01407
01408
01409 {
01410 Teuchos::TimeMonitor lcltimer( *timerLocal_ );
01411 MVT::MvTimesMatAddMv(ONE,*curV,S1,ZERO,*X_);
01412 }
01413
01414
01415 {
01416 Teuchos::TimeMonitor lcltimer( *timerOp_ );
01417 OPT::Apply( *Op_, *X_, *KX_);
01418 count_ApplyOp_ += blockSize_;
01419 }
01420
01421 if (hasM_) {
01422 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
01423 OPT::Apply(*MOp_,*X_,*MX_);
01424 count_ApplyM_ += blockSize_;
01425 }
01426 else {
01427 MX_ = X_;
01428 }
01429
01430
01431
01432 {
01433 Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
01434
01435 MVT::MvAddMv( ONE, *KX_, ZERO, *KX_, *R_ );
01436 Teuchos::SerialDenseMatrix<int,ScalarType> T( blockSize_, blockSize_ );
01437 for (int i = 0; i < blockSize_; ++i) {
01438 T(i,i) = theta_[i];
01439 }
01440 MVT::MvTimesMatAddMv( -ONE, *MX_, T, ONE, *R_ );
01441 }
01442
01443
01444 Rnorms_current_ = false;
01445 R2norms_current_ = false;
01446
01447
01448
01449 if (om_->isVerbosity( Debug ) ) {
01450
01451 CheckList chk;
01452 chk.checkV = true;
01453 chk.checkX = true;
01454 chk.checkKX = true;
01455 chk.checkMX = true;
01456 chk.checkR = true;
01457 om_->print( Debug, accuracyCheck(chk, ": after local update") );
01458 }
01459 else if (om_->isVerbosity( OrthoDetails )) {
01460 CheckList chk;
01461 chk.checkX = true;
01462 chk.checkKX = true;
01463 chk.checkMX = true;
01464 chk.checkR = true;
01465 om_->print( OrthoDetails, accuracyCheck(chk, ": after local update") );
01466 }
01467 }
01468
01469 }
01470
01471
01472
01474
01475 template <class ScalarType, class MV, class OP>
01476 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
01477 BlockDavidson<ScalarType,MV,OP>::getResNorms() {
01478 if (Rnorms_current_ == false) {
01479
01480 orthman_->norm(*R_,Rnorms_);
01481 Rnorms_current_ = true;
01482 }
01483 return Rnorms_;
01484 }
01485
01486
01488
01489 template <class ScalarType, class MV, class OP>
01490 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
01491 BlockDavidson<ScalarType,MV,OP>::getRes2Norms() {
01492 if (R2norms_current_ == false) {
01493
01494 MVT::MvNorm(*R_,R2norms_);
01495 R2norms_current_ = true;
01496 }
01497 return R2norms_;
01498 }
01499
01501
01502 template <class ScalarType, class MV, class OP>
01503 void BlockDavidson<ScalarType,MV,OP>::setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test) {
01504 TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument,
01505 "Anasazi::BlockDavidson::setStatusTest() was passed a null StatusTest.");
01506 tester_ = test;
01507 }
01508
01510
01511 template <class ScalarType, class MV, class OP>
01512 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > BlockDavidson<ScalarType,MV,OP>::getStatusTest() const {
01513 return tester_;
01514 }
01515
01516
01518
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529
01530
01531
01532
01533
01534
01535
01536
01537
01538
01539
01540
01541
01542
01543 template <class ScalarType, class MV, class OP>
01544 std::string BlockDavidson<ScalarType,MV,OP>::accuracyCheck( const CheckList &chk, const std::string &where ) const
01545 {
01546 using std::endl;
01547
01548 std::stringstream os;
01549 os.precision(2);
01550 os.setf(std::ios::scientific, std::ios::floatfield);
01551
01552 os << " Debugging checks: iteration " << iter_ << where << endl;
01553
01554
01555 std::vector<int> lclind(curDim_);
01556 for (int i=0; i<curDim_; ++i) lclind[i] = i;
01557 Teuchos::RCP<MV> lclV;
01558 if (initialized_) {
01559 lclV = MVT::CloneView(*V_,lclind);
01560 }
01561 if (chk.checkV && initialized_) {
01562 MagnitudeType err = orthman_->orthonormError(*lclV);
01563 os << " >> Error in V^H M V == I : " << err << endl;
01564 for (unsigned int i=0; i<auxVecs_.size(); ++i) {
01565 err = orthman_->orthogError(*lclV,*auxVecs_[i]);
01566 os << " >> Error in V^H M Q[" << i << "] == 0 : " << err << endl;
01567 }
01568 Teuchos::SerialDenseMatrix<int,ScalarType> curKK(curDim_,curDim_);
01569 Teuchos::RCP<MV> lclKV = MVT::Clone(*V_,curDim_);
01570 OPT::Apply(*Op_,*lclV,*lclKV);
01571 MVT::MvTransMv(ONE,*lclV,*lclKV,curKK);
01572 Teuchos::SerialDenseMatrix<int,ScalarType> subKK(Teuchos::View,*KK_,curDim_,curDim_);
01573 curKK -= subKK;
01574
01575 for (int j=0; j<curDim_; ++j) {
01576 for (int i=j+1; i<curDim_; ++i) {
01577 curKK(i,j) = curKK(j,i);
01578 }
01579 }
01580 os << " >> Error in V^H K V == KK : " << curKK.normFrobenius() << endl;
01581 }
01582
01583
01584 if (chk.checkX && initialized_) {
01585 MagnitudeType err = orthman_->orthonormError(*X_);
01586 os << " >> Error in X^H M X == I : " << err << endl;
01587 for (unsigned int i=0; i<auxVecs_.size(); ++i) {
01588 err = orthman_->orthogError(*X_,*auxVecs_[i]);
01589 os << " >> Error in X^H M Q[" << i << "] == 0 : " << err << endl;
01590 }
01591 }
01592 if (chk.checkMX && hasM_ && initialized_) {
01593 MagnitudeType err = Utils::errorEquality(*X_, *MX_, MOp_);
01594 os << " >> Error in MX == M*X : " << err << endl;
01595 }
01596 if (chk.checkKX && initialized_) {
01597 MagnitudeType err = Utils::errorEquality(*X_, *KX_, Op_);
01598 os << " >> Error in KX == K*X : " << err << endl;
01599 }
01600
01601
01602 if (chk.checkH && initialized_) {
01603 MagnitudeType err = orthman_->orthonormError(*H_);
01604 os << " >> Error in H^H M H == I : " << err << endl;
01605 err = orthman_->orthogError(*H_,*lclV);
01606 os << " >> Error in H^H M V == 0 : " << err << endl;
01607 err = orthman_->orthogError(*H_,*X_);
01608 os << " >> Error in H^H M X == 0 : " << err << endl;
01609 for (unsigned int i=0; i<auxVecs_.size(); ++i) {
01610 err = orthman_->orthogError(*H_,*auxVecs_[i]);
01611 os << " >> Error in H^H M Q[" << i << "] == 0 : " << err << endl;
01612 }
01613 }
01614 if (chk.checkKH && initialized_) {
01615 MagnitudeType err = Utils::errorEquality(*H_, *KH_, Op_);
01616 os << " >> Error in KH == K*H : " << err << endl;
01617 }
01618 if (chk.checkMH && hasM_ && initialized_) {
01619 MagnitudeType err = Utils::errorEquality(*H_, *MH_, MOp_);
01620 os << " >> Error in MH == M*H : " << err << endl;
01621 }
01622
01623
01624 if (chk.checkR && initialized_) {
01625 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
01626 MVT::MvTransMv(ONE,*X_,*R_,xTx);
01627 MagnitudeType err = xTx.normFrobenius();
01628 os << " >> Error in X^H R == 0 : " << err << endl;
01629 }
01630
01631
01632 if (chk.checkKK && initialized_) {
01633 Teuchos::SerialDenseMatrix<int,ScalarType> SDMerr(curDim_,curDim_), lclKK(Teuchos::View,*KK_,curDim_,curDim_);
01634 for (int j=0; j<curDim_; ++j) {
01635 for (int i=0; i<curDim_; ++i) {
01636 SDMerr(i,j) = lclKK(i,j) - SCT::conjugate(lclKK(j,i));
01637 }
01638 }
01639 os << " >> Error in KK - KK^H == 0 : " << SDMerr.normFrobenius() << endl;
01640 }
01641
01642
01643 if (chk.checkQ) {
01644 for (unsigned int i=0; i<auxVecs_.size(); ++i) {
01645 MagnitudeType err = orthman_->orthonormError(*auxVecs_[i]);
01646 os << " >> Error in Q[" << i << "]^H M Q[" << i << "] == I : " << err << endl;
01647 for (unsigned int j=i+1; j<auxVecs_.size(); ++j) {
01648 err = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
01649 os << " >> Error in Q[" << i << "]^H M Q[" << j << "] == 0 : " << err << endl;
01650 }
01651 }
01652 }
01653
01654 os << endl;
01655
01656 return os.str();
01657 }
01658
01659
01661
01662 template <class ScalarType, class MV, class OP>
01663 void
01664 BlockDavidson<ScalarType,MV,OP>::currentStatus(std::ostream &os)
01665 {
01666 using std::endl;
01667
01668 os.setf(std::ios::scientific, std::ios::floatfield);
01669 os.precision(6);
01670 os <<endl;
01671 os <<"================================================================================" << endl;
01672 os << endl;
01673 os <<" BlockDavidson Solver Status" << endl;
01674 os << endl;
01675 os <<"The solver is "<<(initialized_ ? "initialized." : "not initialized.") << endl;
01676 os <<"The number of iterations performed is " <<iter_<<endl;
01677 os <<"The block size is " << blockSize_<<endl;
01678 os <<"The number of blocks is " << numBlocks_<<endl;
01679 os <<"The current basis size is " << curDim_<<endl;
01680 os <<"The number of auxiliary vectors is "<< numAuxVecs_ << endl;
01681 os <<"The number of operations Op*x is "<<count_ApplyOp_<<endl;
01682 os <<"The number of operations M*x is "<<count_ApplyM_<<endl;
01683 os <<"The number of operations Prec*x is "<<count_ApplyPrec_<<endl;
01684
01685 os.setf(std::ios_base::right, std::ios_base::adjustfield);
01686
01687 if (initialized_) {
01688 os << endl;
01689 os <<"CURRENT EIGENVALUE ESTIMATES "<<endl;
01690 os << std::setw(20) << "Eigenvalue"
01691 << std::setw(20) << "Residual(M)"
01692 << std::setw(20) << "Residual(2)"
01693 << endl;
01694 os <<"--------------------------------------------------------------------------------"<<endl;
01695 for (int i=0; i<blockSize_; ++i) {
01696 os << std::setw(20) << theta_[i];
01697 if (Rnorms_current_) os << std::setw(20) << Rnorms_[i];
01698 else os << std::setw(20) << "not current";
01699 if (R2norms_current_) os << std::setw(20) << R2norms_[i];
01700 else os << std::setw(20) << "not current";
01701 os << endl;
01702 }
01703 }
01704 os <<"================================================================================" << endl;
01705 os << endl;
01706 }
01707
01708 }
01709
01710 #endif
01711
01712