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