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_KRYLOV_SCHUR_HPP
00034 #define ANASAZI_BLOCK_KRYLOV_SCHUR_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 #include "AnasaziHelperTraits.hpp"
00043
00044 #include "AnasaziOrthoManager.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
00066 namespace Anasazi {
00067
00069
00070
00075 template <class ScalarType, class MulVec>
00076 struct BlockKrylovSchurState {
00081 int curDim;
00083 Teuchos::RCP<const MulVec> V;
00089 Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > H;
00091 Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > S;
00093 Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > Q;
00094 BlockKrylovSchurState() : curDim(0), V(Teuchos::null),
00095 H(Teuchos::null), S(Teuchos::null),
00096 Q(Teuchos::null) {}
00097 };
00098
00100
00102
00103
00116 class BlockKrylovSchurInitFailure : public AnasaziError {public:
00117 BlockKrylovSchurInitFailure(const std::string& what_arg) : AnasaziError(what_arg)
00118 {}};
00119
00126 class BlockKrylovSchurOrthoFailure : public AnasaziError {public:
00127 BlockKrylovSchurOrthoFailure(const std::string& what_arg) : AnasaziError(what_arg)
00128 {}};
00129
00131
00132
00133 template <class ScalarType, class MV, class OP>
00134 class BlockKrylovSchur : public Eigensolver<ScalarType,MV,OP> {
00135 public:
00137
00138
00150 BlockKrylovSchur( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00151 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
00152 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00153 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00154 const Teuchos::RCP<OrthoManager<ScalarType,MV> > &ortho,
00155 Teuchos::ParameterList ¶ms
00156 );
00157
00159 virtual ~BlockKrylovSchur() {};
00161
00162
00164
00165
00187 void iterate();
00188
00210 void initialize(BlockKrylovSchurState<ScalarType,MV> state);
00211
00215 void initialize();
00216
00225 bool isInitialized() const { return initialized_; }
00226
00234 BlockKrylovSchurState<ScalarType,MV> getState() const {
00235 BlockKrylovSchurState<ScalarType,MV> state;
00236 state.curDim = curDim_;
00237 state.V = V_;
00238 state.H = H_;
00239 state.Q = Q_;
00240 state.S = schurH_;
00241 return state;
00242 }
00243
00245
00246
00248
00249
00251 int getNumIters() const { return(iter_); }
00252
00254 void resetNumIters() { iter_=0; }
00255
00263 Teuchos::RCP<const MV> getRitzVectors() { return ritzVectors_; }
00264
00272 std::vector<Value<ScalarType> > getRitzValues() {
00273 std::vector<Value<ScalarType> > ret = ritzValues_;
00274 ret.resize(ritzIndex_.size());
00275 return ret;
00276 }
00277
00283 std::vector<int> getRitzIndex() { return ritzIndex_; }
00284
00289 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getResNorms() {
00290 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> ret(0);
00291 return ret;
00292 }
00293
00298 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRes2Norms() {
00299 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> ret(0);
00300 return ret;
00301 }
00302
00307 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRitzRes2Norms() {
00308 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> ret = ritzResiduals_;
00309 ret.resize(ritzIndex_.size());
00310 return ret;
00311 }
00312
00314
00316
00317
00319 void setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test);
00320
00322 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > getStatusTest() const;
00323
00325 const Eigenproblem<ScalarType,MV,OP>& getProblem() const { return(*problem_); };
00326
00333 void setSize(int blockSize, int numBlocks);
00334
00336 void setBlockSize(int blockSize);
00337
00339 void setStepSize(int stepSize);
00340
00342 void setNumRitzVectors(int numRitzVecs);
00343
00345 int getStepSize() const { return(stepSize_); }
00346
00348 int getBlockSize() const { return(blockSize_); }
00349
00351 int getNumRitzVectors() const { return(numRitzVecs_); }
00352
00358 int getCurSubspaceDim() const {
00359 if (!initialized_) return 0;
00360 return curDim_;
00361 }
00362
00364 int getMaxSubspaceDim() const { return (problem_->isHermitian()?blockSize_*numBlocks_:blockSize_*numBlocks_+1); }
00365
00366
00379 void setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs);
00380
00382 Teuchos::Array<Teuchos::RCP<const MV> > getAuxVecs() const {return auxVecs_;}
00383
00385
00387
00388
00390 void currentStatus(std::ostream &os);
00391
00393
00395
00396
00398 bool isRitzVecsCurrent() const { return ritzVecsCurrent_; }
00399
00401 bool isRitzValsCurrent() const { return ritzValsCurrent_; }
00402
00404 bool isSchurCurrent() const { return schurCurrent_; }
00405
00407
00409
00410
00412 void computeRitzVectors();
00413
00415 void computeRitzValues();
00416
00418 void computeSchurForm( const bool sort = true );
00419
00421
00422 private:
00423
00424
00425
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 typedef typename std::vector<ScalarType>::iterator STiter;
00431 typedef typename std::vector<MagnitudeType>::iterator MTiter;
00432 const MagnitudeType MT_ONE;
00433 const MagnitudeType MT_ZERO;
00434 const MagnitudeType NANVAL;
00435 const ScalarType ST_ONE;
00436 const ScalarType ST_ZERO;
00437
00438
00439
00440 struct CheckList {
00441 bool checkV;
00442 bool checkArn;
00443 bool checkAux;
00444 CheckList() : checkV(false), checkArn(false), checkAux(false) {};
00445 };
00446
00447
00448
00449 std::string accuracyCheck(const CheckList &chk, const std::string &where) const;
00450 void sortSchurForm( Teuchos::SerialDenseMatrix<int,ScalarType>& H,
00451 Teuchos::SerialDenseMatrix<int,ScalarType>& Q,
00452 std::vector<int>& order );
00453
00454
00455
00456 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
00457 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sm_;
00458 const Teuchos::RCP<OutputManager<ScalarType> > om_;
00459 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > tester_;
00460 const Teuchos::RCP<OrthoManager<ScalarType,MV> > orthman_;
00461
00462
00463
00464 Teuchos::RCP<const OP> Op_;
00465
00466
00467
00468 Teuchos::RCP<Teuchos::Time> timerOp_, timerSortRitzVal_,
00469 timerCompSF_, timerSortSF_,
00470 timerCompRitzVec_, timerOrtho_;
00471
00472
00473
00474 int count_ApplyOp_;
00475
00476
00477
00478
00479
00480
00481
00482 int blockSize_;
00483
00484 int numBlocks_;
00485
00486
00487 int stepSize_;
00488
00489
00490
00491
00492
00493
00494
00495 bool initialized_;
00496
00497
00498
00499
00500
00501 int curDim_;
00502
00503
00504 Teuchos::RCP<MV> ritzVectors_, V_;
00505 int numRitzVecs_;
00506
00507
00508
00509
00510 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H_;
00511
00512
00513
00514
00515 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > schurH_;
00516 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Q_;
00517
00518
00519 Teuchos::Array<Teuchos::RCP<const MV> > auxVecs_;
00520 int numAuxVecs_;
00521
00522
00523 int iter_;
00524
00525
00526 bool ritzVecsCurrent_, ritzValsCurrent_, schurCurrent_;
00527
00528
00529 std::vector<Value<ScalarType> > ritzValues_;
00530 std::vector<MagnitudeType> ritzResiduals_;
00531
00532
00533 std::vector<int> ritzIndex_;
00534 std::vector<int> ritzOrder_;
00535
00536
00537 int numRitzPrint_;
00538 };
00539
00540
00542
00543 template <class ScalarType, class MV, class OP>
00544 BlockKrylovSchur<ScalarType,MV,OP>::BlockKrylovSchur(
00545 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00546 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
00547 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00548 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00549 const Teuchos::RCP<OrthoManager<ScalarType,MV> > &ortho,
00550 Teuchos::ParameterList ¶ms
00551 ) :
00552 MT_ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
00553 MT_ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
00554 NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan()),
00555 ST_ONE(Teuchos::ScalarTraits<ScalarType>::one()),
00556 ST_ZERO(Teuchos::ScalarTraits<ScalarType>::zero()),
00557
00558 problem_(problem),
00559 sm_(sorter),
00560 om_(printer),
00561 tester_(tester),
00562 orthman_(ortho),
00563
00564 timerOp_(Teuchos::TimeMonitor::getNewTimer("BlockKrylovSchur::Operation Op*x")),
00565 timerSortRitzVal_(Teuchos::TimeMonitor::getNewTimer("BlockKrylovSchur::Sorting Ritz values")),
00566 timerCompSF_(Teuchos::TimeMonitor::getNewTimer("BlockKrylovSchur::Computing Schur form")),
00567 timerSortSF_(Teuchos::TimeMonitor::getNewTimer("BlockKrylovSchur::Sorting Schur form")),
00568 timerCompRitzVec_(Teuchos::TimeMonitor::getNewTimer("BlockKrylovSchur::Computing Ritz vectors")),
00569 timerOrtho_(Teuchos::TimeMonitor::getNewTimer("BlockKrylovSchur::Orthogonalization")),
00570 count_ApplyOp_(0),
00571
00572 blockSize_(0),
00573 numBlocks_(0),
00574 stepSize_(0),
00575 initialized_(false),
00576 curDim_(0),
00577 numRitzVecs_(0),
00578 auxVecs_( Teuchos::Array<Teuchos::RCP<const MV> >(0) ),
00579 numAuxVecs_(0),
00580 iter_(0),
00581 ritzVecsCurrent_(false),
00582 ritzValsCurrent_(false),
00583 schurCurrent_(false),
00584 numRitzPrint_(0)
00585 {
00586 TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument,
00587 "Anasazi::BlockKrylovSchur::constructor: user specified null problem pointer.");
00588 TEST_FOR_EXCEPTION(sm_ == Teuchos::null,std::invalid_argument,
00589 "Anasazi::BlockKrylovSchur::constructor: user passed null sort manager pointer.");
00590 TEST_FOR_EXCEPTION(om_ == Teuchos::null,std::invalid_argument,
00591 "Anasazi::BlockKrylovSchur::constructor: user passed null output manager pointer.");
00592 TEST_FOR_EXCEPTION(tester_ == Teuchos::null,std::invalid_argument,
00593 "Anasazi::BlockKrylovSchur::constructor: user passed null status test pointer.");
00594 TEST_FOR_EXCEPTION(orthman_ == Teuchos::null,std::invalid_argument,
00595 "Anasazi::BlockKrylovSchur::constructor: user passed null orthogonalization manager pointer.");
00596 TEST_FOR_EXCEPTION(problem_->isProblemSet() == false, std::invalid_argument,
00597 "Anasazi::BlockKrylovSchur::constructor: user specified problem is not set.");
00598 TEST_FOR_EXCEPTION(sorter == Teuchos::null,std::invalid_argument,
00599 "Anasazi::BlockKrylovSchur::constructor: user specified null sort manager pointer.");
00600 TEST_FOR_EXCEPTION(printer == Teuchos::null,std::invalid_argument,
00601 "Anasazi::BlockKrylovSchur::constructor: user specified null output manager pointer.");
00602 TEST_FOR_EXCEPTION(tester == Teuchos::null,std::invalid_argument,
00603 "Anasazi::BlockKrylovSchur::constructor: user specified null status test pointer.");
00604 TEST_FOR_EXCEPTION(ortho == Teuchos::null,std::invalid_argument,
00605 "Anasazi::BlockKrylovSchur::constructor: user specified null ortho manager pointer.");
00606
00607
00608 Op_ = problem_->getOperator();
00609
00610
00611 TEST_FOR_EXCEPTION(!params.isParameter("Step Size"), std::invalid_argument,
00612 "Anasazi::BlockKrylovSchur::constructor: mandatory parameter 'Step Size' is not specified.");
00613 int ss = params.get("Step Size",numBlocks_);
00614 setStepSize(ss);
00615
00616
00617 int bs = params.get("Block Size", 1);
00618 int nb = params.get("Num Blocks", 3*problem_->getNEV());
00619 setSize(bs,nb);
00620
00621
00622
00623 int numRitzVecs = params.get("Number of Ritz Vectors", 0);
00624 setNumRitzVectors( numRitzVecs );
00625
00626
00627 numRitzPrint_ = params.get("Print Number of Ritz Values", bs);
00628 }
00629
00630
00632
00633
00634 template <class ScalarType, class MV, class OP>
00635 void BlockKrylovSchur<ScalarType,MV,OP>::setBlockSize (int blockSize)
00636 {
00637 setSize(blockSize,numBlocks_);
00638 }
00639
00640
00642
00643 template <class ScalarType, class MV, class OP>
00644 void BlockKrylovSchur<ScalarType,MV,OP>::setStepSize (int stepSize)
00645 {
00646 TEST_FOR_EXCEPTION(stepSize <= 0, std::invalid_argument, "Anasazi::BlockKrylovSchur::setStepSize(): new step size must be positive and non-zero.");
00647 stepSize_ = stepSize;
00648 }
00649
00651
00652 template <class ScalarType, class MV, class OP>
00653 void BlockKrylovSchur<ScalarType,MV,OP>::setNumRitzVectors (int numRitzVecs)
00654 {
00655
00656
00657
00658 TEST_FOR_EXCEPTION(numRitzVecs < 0, std::invalid_argument, "Anasazi::BlockKrylovSchur::setNumRitzVectors(): number of Ritz vectors to compute must be positive.");
00659
00660
00661 if (numRitzVecs != numRitzVecs_) {
00662 if (numRitzVecs) {
00663 ritzVectors_ = Teuchos::null;
00664 ritzVectors_ = MVT::Clone(*V_, numRitzVecs);
00665 } else {
00666 ritzVectors_ = Teuchos::null;
00667 }
00668 numRitzVecs_ = numRitzVecs;
00669 ritzVecsCurrent_ = false;
00670 }
00671 }
00672
00674
00675 template <class ScalarType, class MV, class OP>
00676 void BlockKrylovSchur<ScalarType,MV,OP>::setSize (int blockSize, int numBlocks)
00677 {
00678
00679
00680
00681 TEST_FOR_EXCEPTION(numBlocks <= 0 || blockSize <= 0, std::invalid_argument, "Anasazi::BlockKrylovSchur::setSize was passed a non-positive argument.");
00682 TEST_FOR_EXCEPTION(numBlocks < 3, std::invalid_argument, "Anasazi::BlockKrylovSchur::setSize(): numBlocks must be at least three.");
00683 if (blockSize == blockSize_ && numBlocks == numBlocks_) {
00684
00685 return;
00686 }
00687
00688 blockSize_ = blockSize;
00689 numBlocks_ = numBlocks;
00690
00691 Teuchos::RCP<const MV> tmp;
00692
00693
00694
00695
00696
00697 if (problem_->getInitVec() != Teuchos::null) {
00698 tmp = problem_->getInitVec();
00699 }
00700 else {
00701 tmp = V_;
00702 TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
00703 "Anasazi::BlockKrylovSchur::setSize(): eigenproblem did not specify initial vectors to clone from.");
00704 }
00705
00706
00708
00709
00710 int newsd;
00711 if (problem_->isHermitian()) {
00712 newsd = blockSize_*numBlocks_;
00713 } else {
00714 newsd = blockSize_*numBlocks_+1;
00715 }
00716
00717 TEST_FOR_EXCEPTION(newsd > MVT::GetVecLength(*tmp),std::invalid_argument,
00718 "Anasazi::BlockKrylovSchur::setSize(): maximum basis size is larger than problem dimension.");
00719
00720 ritzValues_.resize(newsd);
00721 ritzResiduals_.resize(newsd,MT_ONE);
00722 ritzOrder_.resize(newsd);
00723 V_ = Teuchos::null;
00724 V_ = MVT::Clone(*tmp,newsd+blockSize_);
00725 H_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(newsd+blockSize_,newsd) );
00726 Q_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(newsd,newsd) );
00727
00728 initialized_ = false;
00729 curDim_ = 0;
00730 }
00731
00732
00734
00735 template <class ScalarType, class MV, class OP>
00736 void BlockKrylovSchur<ScalarType,MV,OP>::setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs) {
00737 typedef typename Teuchos::Array<Teuchos::RCP<const MV> >::iterator tarcpmv;
00738
00739
00740 auxVecs_ = auxvecs;
00741
00742 if (om_->isVerbosity( Debug ) ) {
00743
00744 CheckList chk;
00745 chk.checkAux = true;
00746 om_->print( Debug, accuracyCheck(chk, ": in setAuxVecs()") );
00747 }
00748
00749 numAuxVecs_ = 0;
00750 for (tarcpmv i=auxVecs_.begin(); i != auxVecs_.end(); i++) {
00751 numAuxVecs_ += MVT::GetNumberVecs(**i);
00752 }
00753
00754
00755 if (numAuxVecs_ > 0 && initialized_) {
00756 initialized_ = false;
00757 }
00758 }
00759
00761
00762
00763
00764
00765
00766
00767
00768
00769 template <class ScalarType, class MV, class OP>
00770 void BlockKrylovSchur<ScalarType,MV,OP>::initialize(BlockKrylovSchurState<ScalarType,MV> newstate)
00771 {
00772
00773
00774 std::vector<int> bsind(blockSize_);
00775 for (int i=0; i<blockSize_; i++) bsind[i] = i;
00776
00777
00778
00779
00780
00781
00782
00783 std::string errstr("Anasazi::BlockKrylovSchur::initialize(): specified multivectors must have a consistent length and width.");
00784
00785
00786
00787 if (newstate.V != Teuchos::null && newstate.H != Teuchos::null) {
00788
00789
00790
00791 TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.V) != MVT::GetVecLength(*V_),
00792 std::invalid_argument, errstr );
00793 if (newstate.V != V_) {
00794 TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < blockSize_,
00795 std::invalid_argument, errstr );
00796 TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) > getMaxSubspaceDim(),
00797 std::invalid_argument, errstr );
00798 }
00799 TEST_FOR_EXCEPTION( newstate.curDim > getMaxSubspaceDim(),
00800 std::invalid_argument, errstr );
00801
00802 curDim_ = newstate.curDim;
00803 int lclDim = MVT::GetNumberVecs(*newstate.V);
00804
00805
00806 TEST_FOR_EXCEPTION(newstate.H->numRows() < curDim_ || newstate.H->numCols() < curDim_, std::invalid_argument, errstr);
00807
00808 if (curDim_ == 0 && lclDim > blockSize_) {
00809 om_->stream(Warnings) << "Anasazi::BlockKrylovSchur::initialize(): the solver was initialized with a kernel of " << lclDim << std::endl
00810 << "The block size however is only " << blockSize_ << std::endl
00811 << "The last " << lclDim - blockSize_ << " vectors of the kernel will be overwritten on the first call to iterate()." << std::endl;
00812 }
00813
00814
00815
00816 if (newstate.V != V_) {
00817 std::vector<int> nevind(lclDim);
00818 for (int i=0; i<lclDim; i++) nevind[i] = i;
00819 MVT::SetBlock(*newstate.V,nevind,*V_);
00820 }
00821
00822
00823 if (newstate.H != H_) {
00824 H_->putScalar( ST_ZERO );
00825 Teuchos::SerialDenseMatrix<int,ScalarType> newH(Teuchos::View,*newstate.H,curDim_+blockSize_,curDim_);
00826 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclH;
00827 lclH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*H_,curDim_+blockSize_,curDim_) );
00828 lclH->assign(newH);
00829
00830
00831 lclH = Teuchos::null;
00832 }
00833
00834 }
00835 else {
00836
00837
00838 Teuchos::RCP<const MV> ivec = problem_->getInitVec();
00839 TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::invalid_argument,
00840 "Anasazi::BlockKrylovSchur::initialize(): eigenproblem did not specify initial vectors to clone from.");
00841
00842 int lclDim = MVT::GetNumberVecs(*ivec);
00843 bool userand = false;
00844 if (lclDim < blockSize_) {
00845
00846
00847 userand = true;
00848 }
00849
00850 if (userand) {
00851
00852 std::vector<int> dimind2(lclDim);
00853 for (int i=0; i<lclDim; i++) { dimind2[i] = i; }
00854
00855
00856 Teuchos::RCP<MV> newV1 = MVT::CloneView(*V_,dimind2);
00857
00858
00859 MVT::SetBlock(*ivec,dimind2,*newV1);
00860
00861
00862 dimind2.resize(blockSize_-lclDim);
00863 for (int i=0; i<blockSize_-lclDim; i++) { dimind2[i] = lclDim + i; }
00864
00865
00866 Teuchos::RCP<MV> newV2 = MVT::CloneView(*V_,dimind2);
00867 MVT::MvRandom(*newV2);
00868 }
00869 else {
00870
00871 Teuchos::RCP<MV> newV1 = MVT::CloneView(*V_,bsind);
00872
00873
00874 Teuchos::RCP<const MV> ivecV = MVT::CloneView(*ivec,bsind);
00875
00876
00877 MVT::SetBlock(*ivecV,bsind,*newV1);
00878 }
00879
00880
00881 Teuchos::RCP<MV> newV = MVT::CloneView(*V_,bsind);
00882
00883
00884 if (auxVecs_.size() > 0) {
00885 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
00886
00887 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummy;
00888 int rank = orthman_->projectAndNormalize(*newV,auxVecs_);
00889 TEST_FOR_EXCEPTION( rank != blockSize_,BlockKrylovSchurInitFailure,
00890 "Anasazi::BlockKrylovSchur::initialize(): couldn't generate initial basis of full rank." );
00891 }
00892 else {
00893 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
00894
00895 int rank = orthman_->normalize(*newV);
00896 TEST_FOR_EXCEPTION( rank != blockSize_,BlockKrylovSchurInitFailure,
00897 "Anasazi::BlockKrylovSchur::initialize(): couldn't generate initial basis of full rank." );
00898 }
00899
00900
00901 curDim_ = 0;
00902
00903
00904 newV = Teuchos::null;
00905 }
00906
00907
00908 ritzVecsCurrent_ = false;
00909 ritzValsCurrent_ = false;
00910 schurCurrent_ = false;
00911
00912
00913 initialized_ = true;
00914
00915 if (om_->isVerbosity( Debug ) ) {
00916
00917 CheckList chk;
00918 chk.checkV = true;
00919 chk.checkArn = true;
00920 chk.checkAux = true;
00921 om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
00922 }
00923
00924
00925 if (om_->isVerbosity(Debug)) {
00926 currentStatus( om_->stream(Debug) );
00927 }
00928 else if (om_->isVerbosity(IterationDetails)) {
00929 currentStatus( om_->stream(IterationDetails) );
00930 }
00931 }
00932
00933
00935
00936 template <class ScalarType, class MV, class OP>
00937 void BlockKrylovSchur<ScalarType,MV,OP>::initialize()
00938 {
00939 BlockKrylovSchurState<ScalarType,MV> empty;
00940 initialize(empty);
00941 }
00942
00943
00945
00946 template <class ScalarType, class MV, class OP>
00947 void BlockKrylovSchur<ScalarType,MV,OP>::iterate()
00948 {
00949
00950
00951
00952 if (initialized_ == false) {
00953 initialize();
00954 }
00955
00956
00957
00958 int searchDim = blockSize_*numBlocks_;
00959 if (problem_->isHermitian() == false) {
00960 searchDim++;
00961 }
00962
00964
00965
00966
00967
00968 while (tester_->checkStatus(this) != Passed && curDim_+blockSize_ <= searchDim) {
00969
00970 iter_++;
00971
00972
00973 int lclDim = curDim_ + blockSize_;
00974
00975
00976 std::vector<int> curind(blockSize_);
00977 for (int i=0; i<blockSize_; i++) { curind[i] = lclDim + i; }
00978 Teuchos::RCP<MV> Vnext = MVT::CloneView(*V_,curind);
00979
00980
00981
00982 for (int i=0; i<blockSize_; i++) { curind[i] = curDim_ + i; }
00983 Teuchos::RCP<MV> Vprev = MVT::CloneView(*V_,curind);
00984
00985
00986 {
00987 Teuchos::TimeMonitor lcltimer( *timerOp_ );
00988 OPT::Apply(*Op_,*Vprev,*Vnext);
00989 count_ApplyOp_ += blockSize_;
00990 }
00991 Vprev = Teuchos::null;
00992
00993
00994 {
00995 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
00996
00997
00998 std::vector<int> prevind(lclDim);
00999 for (int i=0; i<lclDim; i++) { prevind[i] = i; }
01000 Vprev = MVT::CloneView(*V_,prevind);
01001 Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
01002
01003
01004 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
01005 subH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
01006 ( Teuchos::View,*H_,lclDim,blockSize_,0,curDim_ ) );
01007 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubH;
01008 AsubH.append( subH );
01009
01010
01011 if (auxVecs_.size() > 0) {
01012 for (unsigned int i=0; i<auxVecs_.size(); i++) {
01013 AVprev.append( auxVecs_[i] );
01014 AsubH.append( Teuchos::null );
01015 }
01016 }
01017
01018
01019 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
01020 subR = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
01021 ( Teuchos::View,*H_,blockSize_,blockSize_,lclDim,curDim_ ) );
01022 int rank = orthman_->projectAndNormalize(*Vnext,AVprev,AsubH,subR);
01023 TEST_FOR_EXCEPTION(rank != blockSize_,BlockKrylovSchurOrthoFailure,
01024 "Anasazi::BlockKrylovSchur::iterate(): couldn't generate basis of full rank.");
01025 }
01026
01027
01028
01029
01030 Vnext = Teuchos::null;
01031 curDim_ += blockSize_;
01032
01033 ritzVecsCurrent_ = false;
01034 ritzValsCurrent_ = false;
01035 schurCurrent_ = false;
01036
01037
01038 if (!(iter_%stepSize_)) {
01039 computeRitzValues();
01040 }
01041
01042
01043 if (om_->isVerbosity( Debug ) ) {
01044
01045 CheckList chk;
01046 chk.checkV = true;
01047 chk.checkArn = true;
01048 om_->print( Debug, accuracyCheck(chk, ": after local update") );
01049 }
01050 else if (om_->isVerbosity( OrthoDetails ) ) {
01051 CheckList chk;
01052 chk.checkV = true;
01053 om_->print( OrthoDetails, accuracyCheck(chk, ": after local update") );
01054 }
01055
01056
01057 if (om_->isVerbosity(Debug)) {
01058 currentStatus( om_->stream(Debug) );
01059 }
01060 else if (om_->isVerbosity(IterationDetails)) {
01061 currentStatus( om_->stream(IterationDetails) );
01062 }
01063
01064 }
01065
01066 }
01067
01068
01070
01071
01072
01073
01074
01075
01076
01077
01078
01079
01080
01081
01082 template <class ScalarType, class MV, class OP>
01083 std::string BlockKrylovSchur<ScalarType,MV,OP>::accuracyCheck( const CheckList &chk, const std::string &where ) const
01084 {
01085 std::stringstream os;
01086 os.precision(2);
01087 os.setf(std::ios::scientific, std::ios::floatfield);
01088 MagnitudeType tmp;
01089
01090 os << " Debugging checks: iteration " << iter_ << where << std::endl;
01091
01092
01093 std::vector<int> lclind(curDim_);
01094 for (int i=0; i<curDim_; i++) lclind[i] = i;
01095 std::vector<int> bsind(blockSize_);
01096 for (int i=0; i<blockSize_; i++) { bsind[i] = curDim_ + i; }
01097
01098 Teuchos::RCP<MV> lclV,lclF,lclAV;
01099 if (curDim_)
01100 lclV = MVT::CloneView(*V_,lclind);
01101 lclF = MVT::CloneView(*V_,bsind);
01102
01103 if (chk.checkV) {
01104 if (curDim_) {
01105 tmp = orthman_->orthonormError(*lclV);
01106 os << " >> Error in V^H M V == I : " << tmp << std::endl;
01107 }
01108 tmp = orthman_->orthonormError(*lclF);
01109 os << " >> Error in F^H M F == I : " << tmp << std::endl;
01110 if (curDim_) {
01111 tmp = orthman_->orthogError(*lclV,*lclF);
01112 os << " >> Error in V^H M F == 0 : " << tmp << std::endl;
01113 }
01114 for (unsigned int i=0; i<auxVecs_.size(); i++) {
01115 if (curDim_) {
01116 tmp = orthman_->orthogError(*lclV,*auxVecs_[i]);
01117 os << " >> Error in V^H M Aux[" << i << "] == 0 : " << tmp << std::endl;
01118 }
01119 tmp = orthman_->orthogError(*lclF,*auxVecs_[i]);
01120 os << " >> Error in F^H M Aux[" << i << "] == 0 : " << tmp << std::endl;
01121 }
01122 }
01123
01124 if (chk.checkArn) {
01125
01126 if (curDim_) {
01127
01128 lclAV = MVT::Clone(*V_,curDim_);
01129 {
01130 Teuchos::TimeMonitor lcltimer( *timerOp_ );
01131 OPT::Apply(*Op_,*lclV,*lclAV);
01132 }
01133
01134
01135 Teuchos::SerialDenseMatrix<int,ScalarType> subH(Teuchos::View,*H_,curDim_,curDim_);
01136 MVT::MvTimesMatAddMv( -ST_ONE, *lclV, subH, ST_ONE, *lclAV );
01137
01138
01139 Teuchos::SerialDenseMatrix<int,ScalarType> curB(Teuchos::View,*H_,
01140 blockSize_,curDim_, curDim_ );
01141 MVT::MvTimesMatAddMv( -ST_ONE, *lclF, curB, ST_ONE, *lclAV );
01142
01143
01144 std::vector<MagnitudeType> arnNorms( curDim_ );
01145 orthman_->norm( *lclAV, arnNorms );
01146
01147 for (int i=0; i<curDim_; i++) {
01148 os << " >> Error in Krylov-Schur factorization (R = AV-VS-FB^H), ||R[" << i << "]|| : " << arnNorms[i] << std::endl;
01149 }
01150 }
01151 }
01152
01153 if (chk.checkAux) {
01154 for (unsigned int i=0; i<auxVecs_.size(); i++) {
01155 tmp = orthman_->orthonormError(*auxVecs_[i]);
01156 os << " >> Error in Aux[" << i << "]^H M Aux[" << i << "] == I : " << tmp << std::endl;
01157 for (unsigned int j=i+1; j<auxVecs_.size(); j++) {
01158 tmp = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
01159 os << " >> Error in Aux[" << i << "]^H M Aux[" << j << "] == 0 : " << tmp << std::endl;
01160 }
01161 }
01162 }
01163
01164 os << std::endl;
01165
01166 return os.str();
01167 }
01168
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179
01180 template <class ScalarType, class MV, class OP>
01181 void BlockKrylovSchur<ScalarType,MV,OP>::computeRitzValues()
01182 {
01183
01184 if (initialized_) {
01185
01186
01187
01188 if (!ritzValsCurrent_) {
01189
01190
01191 computeSchurForm( false );
01192 }
01193 }
01194 }
01195
01197
01198
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208 template <class ScalarType, class MV, class OP>
01209 void BlockKrylovSchur<ScalarType,MV,OP>::computeRitzVectors()
01210 {
01211 Teuchos::TimeMonitor LocalTimer(*timerCompRitzVec_);
01212
01213 TEST_FOR_EXCEPTION(numRitzVecs_==0, std::invalid_argument,
01214 "Anasazi::BlockKrylovSchur::computeRitzVectors(): no Ritz vectors were required from this solver.");
01215
01216 TEST_FOR_EXCEPTION(curDim_ < numRitzVecs_, std::invalid_argument,
01217 "Anasazi::BlockKrylovSchur::computeRitzVectors(): the current subspace is not large enough to compute the number of requested Ritz vectors.");
01218
01219
01220
01221 if (curDim_ && initialized_) {
01222
01223
01224 if (!ritzVecsCurrent_) {
01225
01226
01227 if (!schurCurrent_) {
01228
01229
01230 computeSchurForm( true );
01231 }
01232
01233
01234
01235 TEST_FOR_EXCEPTION(ritzIndex_[numRitzVecs_-1]==1, std::logic_error,
01236 "Anasazi::BlockKrylovSchur::computeRitzVectors(): the number of required Ritz vectors splits a complex conjugate pair.");
01237
01238 Teuchos::LAPACK<int,ScalarType> lapack;
01239 Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249 std::vector<int> curind( curDim_ );
01250 for (int i=0; i<curDim_; i++) { curind[i] = i; }
01251 Teuchos::RCP<MV> Vtemp = MVT::CloneView( *V_, curind );
01252 if (problem_->isHermitian()) {
01253
01254 Teuchos::SerialDenseMatrix<int,ScalarType> subQ( Teuchos::View, *Q_, curDim_, numRitzVecs_ );
01255
01256
01257 MVT::MvTimesMatAddMv( ST_ONE, *Vtemp, subQ, ST_ZERO, *ritzVectors_ );
01258
01259 } else {
01260
01261
01262 Teuchos::SerialDenseMatrix<int,ScalarType> subQ( Teuchos::View, *Q_, curDim_, curDim_ );
01263
01264
01265 Teuchos::RCP<MV> tmpritzVectors_ = MVT::Clone( *V_, curDim_ );
01266
01267
01268 MVT::MvTimesMatAddMv( ST_ONE, *Vtemp, subQ, ST_ZERO, *tmpritzVectors_ );
01269
01270
01271
01272
01273
01274
01275
01276 int lwork = 3*curDim_;
01277 std::vector<ScalarType> work( lwork );
01278 std::vector<MagnitudeType> rwork( curDim_ );
01279 char side = 'R';
01280 int mm, info = 0;
01281 const int ldvl = 1;
01282 ScalarType vl[ ldvl ];
01283 Teuchos::SerialDenseMatrix<int,ScalarType> copyQ( Teuchos::Copy, *Q_, curDim_, curDim_ );
01284 lapack.TREVC( side, curDim_, schurH_->values(), schurH_->stride(), vl, ldvl,
01285 copyQ.values(), copyQ.stride(), curDim_, &mm, &work[0], &rwork[0], &info );
01286 TEST_FOR_EXCEPTION(info != 0, std::logic_error,
01287 "Anasazi::BlockKrylovSchur::computeRitzVectors(): TREVC(n==" << curDim_ << ") returned info " << info << " != 0.");
01288
01289
01290 Teuchos::SerialDenseMatrix<int,ScalarType> subCopyQ( Teuchos::View, copyQ, curDim_, numRitzVecs_ );
01291
01292
01293 curind.resize( numRitzVecs_ );
01294 Teuchos::RCP<MV> view_ritzVectors = MVT::CloneView( *ritzVectors_, curind );
01295 MVT::MvTimesMatAddMv( ST_ONE, *tmpritzVectors_, subCopyQ, ST_ZERO, *view_ritzVectors );
01296
01297
01298 std::vector<MagnitudeType> ritzNrm( numRitzVecs_ );
01299 MVT::MvNorm( *view_ritzVectors, ritzNrm );
01300
01301
01302 tmpritzVectors_ = Teuchos::null;
01303 view_ritzVectors = Teuchos::null;
01304
01305
01306 ScalarType ritzScale = ST_ONE;
01307 for (int i=0; i<numRitzVecs_; i++) {
01308
01309
01310 if (ritzIndex_[i] == 1 ) {
01311 ritzScale = ST_ONE/lapack_mag.LAPY2(ritzNrm[i],ritzNrm[i+1]);
01312 std::vector<int> newind(2);
01313 newind[0] = i; newind[1] = i+1;
01314 tmpritzVectors_ = MVT::CloneCopy( *ritzVectors_, newind );
01315 view_ritzVectors = MVT::CloneView( *ritzVectors_, newind );
01316 MVT::MvAddMv( ritzScale, *tmpritzVectors_, ST_ZERO, *tmpritzVectors_, *view_ritzVectors );
01317
01318
01319 i++;
01320 } else {
01321
01322
01323 std::vector<int> newind(1);
01324 newind[0] = i;
01325 tmpritzVectors_ = MVT::CloneCopy( *ritzVectors_, newind );
01326 view_ritzVectors = MVT::CloneView( *ritzVectors_, newind );
01327 MVT::MvAddMv( ST_ONE/ritzNrm[i], *tmpritzVectors_, ST_ZERO, *tmpritzVectors_, *view_ritzVectors );
01328 }
01329 }
01330
01331 }
01332
01333
01334 ritzVecsCurrent_ = true;
01335
01336 }
01337 }
01338 }
01339
01340
01342
01343 template <class ScalarType, class MV, class OP>
01344 void BlockKrylovSchur<ScalarType,MV,OP>::setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test) {
01345 TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument,
01346 "Anasazi::BlockKrylovSchur::setStatusTest() was passed a null StatusTest.");
01347 tester_ = test;
01348 }
01349
01351
01352 template <class ScalarType, class MV, class OP>
01353 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > BlockKrylovSchur<ScalarType,MV,OP>::getStatusTest() const {
01354 return tester_;
01355 }
01356
01357
01359
01360
01361
01362
01363
01364
01365
01366
01367
01368
01369
01370 template <class ScalarType, class MV, class OP>
01371 void BlockKrylovSchur<ScalarType,MV,OP>::computeSchurForm( const bool sort )
01372 {
01373
01374 Teuchos::TimeMonitor LocalTimer(*timerCompSF_);
01375
01376
01377 if (curDim_) {
01378
01379
01380 if (!schurCurrent_) {
01381
01382
01383
01384 if (!ritzValsCurrent_) {
01385 Teuchos::LAPACK<int,ScalarType> lapack;
01386 Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
01387 Teuchos::BLAS<int,ScalarType> blas;
01388 Teuchos::BLAS<int,MagnitudeType> blas_mag;
01389
01390
01391 Teuchos::SerialDenseMatrix<int,ScalarType> subQ( Teuchos::View, *Q_, curDim_, curDim_ );
01392
01393
01394 schurH_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *H_, curDim_, curDim_ ) );
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408 int lwork = 3*curDim_;
01409 std::vector<ScalarType> work( lwork );
01410 std::vector<MagnitudeType> rwork( curDim_ );
01411 std::vector<MagnitudeType> tmp_rRitzValues( curDim_ );
01412 std::vector<MagnitudeType> tmp_iRitzValues( curDim_ );
01413 std::vector<int> bwork( curDim_ );
01414 int info = 0, sdim = 0;
01415 char jobvs = 'V';
01416 lapack.GEES( jobvs,curDim_, schurH_->values(), schurH_->stride(), &sdim, &tmp_rRitzValues[0],
01417 &tmp_iRitzValues[0], subQ.values(), subQ.stride(), &work[0], lwork,
01418 &rwork[0], &bwork[0], &info );
01419
01420 TEST_FOR_EXCEPTION(info != 0, std::logic_error,
01421 "Anasazi::BlockKrylovSchur::computeSchurForm(): GEES(n==" << curDim_ << ") returned info " << info << " != 0.");
01422
01423
01424
01425
01426
01427
01428
01429
01430
01431
01432
01433
01434
01435
01436
01437 Teuchos::SerialDenseMatrix<int,ScalarType> curB(Teuchos::View, *H_,
01438 blockSize_, curDim_, curDim_ );
01439
01440
01441 Teuchos::SerialDenseMatrix<int,ScalarType> subB( blockSize_, curDim_ );
01442 blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, blockSize_, curDim_, curDim_, ST_ONE,
01443 curB.values(), curB.stride(), subQ.values(), subQ.stride(),
01444 ST_ZERO, subB.values(), subB.stride() );
01445
01446
01447
01448 ScalarType* b_ptr = subB.values();
01449 if (problem_->isHermitian()) {
01450
01451
01452
01453 for (int i=0; i<curDim_ ; i++) {
01454 ritzResiduals_[i] = blas.NRM2(blockSize_, b_ptr + i*blockSize_, 1);
01455 }
01456 } else {
01457
01458
01459
01460 char side = 'R';
01461 int mm;
01462 const int ldvl = 1;
01463 ScalarType vl[ ldvl ];
01464 Teuchos::SerialDenseMatrix<int,ScalarType> S( curDim_, curDim_ );
01465 lapack.TREVC( side, curDim_, schurH_->values(), schurH_->stride(), vl, ldvl,
01466 S.values(), S.stride(), curDim_, &mm, &work[0], &rwork[0], &info );
01467
01468 TEST_FOR_EXCEPTION(info != 0, std::logic_error,
01469 "Anasazi::BlockKrylovSchur::computeSchurForm(): TREVC(n==" << curDim_ << ") returned info " << info << " != 0.");
01470
01471
01472
01473 HelperTraits<ScalarType>::scaleRitzVectors( tmp_iRitzValues, &S );
01474
01475
01476
01477 Teuchos::SerialDenseMatrix<int,ScalarType> ritzRes( blockSize_, curDim_ );
01478 blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, blockSize_, curDim_, curDim_, ST_ONE,
01479 subB.values(), subB.stride(), S.values(), S.stride(),
01480 ST_ZERO, ritzRes.values(), ritzRes.stride() );
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
01507
01508 HelperTraits<ScalarType>::computeRitzResiduals( tmp_iRitzValues, ritzRes, &ritzResiduals_ );
01509 }
01510
01511
01512
01513 {
01514 Teuchos::TimeMonitor LocalTimer2(*timerSortRitzVal_);
01515 int i=0;
01516 if (problem_->isHermitian()) {
01517
01518
01519 sm_->sort(tmp_rRitzValues, Teuchos::rcpFromRef(ritzOrder_), curDim_);
01520 ritzIndex_.clear();
01521 while ( i < curDim_ ) {
01522
01523 ritzValues_[i].set(tmp_rRitzValues[i], MT_ZERO);
01524 ritzIndex_.push_back(0);
01525 i++;
01526 }
01527 }
01528 else {
01529
01530
01531 sm_->sort(tmp_rRitzValues, tmp_iRitzValues, Teuchos::rcpFromRef(ritzOrder_) , curDim_);
01532 HelperTraits<ScalarType>::sortRitzValues( tmp_rRitzValues, tmp_iRitzValues, &ritzValues_, &ritzOrder_, &ritzIndex_ );
01533 }
01534
01535
01536 std::vector<MagnitudeType> ritz2( curDim_ );
01537 for (i=0; i<curDim_; i++) { ritz2[i] = ritzResiduals_[ ritzOrder_[i] ]; }
01538 blas_mag.COPY( curDim_, &ritz2[0], 1, &ritzResiduals_[0], 1 );
01539
01540
01541 ritzValsCurrent_ = true;
01542 }
01543
01544 }
01545
01546
01547
01548
01549
01550
01551
01552
01553
01554
01555
01556 if (sort) {
01557 sortSchurForm( *schurH_, *Q_, ritzOrder_ );
01558
01559
01560 schurCurrent_ = true;
01561 }
01562 }
01563
01564 }
01565
01566 }
01567
01568
01570
01571 template <class ScalarType, class MV, class OP>
01572 void BlockKrylovSchur<ScalarType,MV,OP>::sortSchurForm( Teuchos::SerialDenseMatrix<int,ScalarType>& H,
01573 Teuchos::SerialDenseMatrix<int,ScalarType>& Q,
01574 std::vector<int>& order )
01575 {
01576
01577 Teuchos::TimeMonitor LocalTimer(*timerSortSF_);
01578
01579
01580
01581
01582
01583
01584
01585
01586
01587 int i = 0, nevtemp = 0;
01588 char compq = 'V';
01589 std::vector<int> offset2( curDim_ );
01590 std::vector<int> order2( curDim_ );
01591
01592
01593 Teuchos::LAPACK<int,ScalarType> lapack;
01594 int lwork = 3*curDim_;
01595 std::vector<ScalarType> work( lwork );
01596
01597 while (i < curDim_) {
01598 if ( ritzIndex_[i] != 0 ) {
01599 offset2[nevtemp] = 0;
01600 for (int j=i; j<curDim_; j++) {
01601 if (order[j] > order[i]) { offset2[nevtemp]++; }
01602 }
01603 order2[nevtemp] = order[i];
01604 i = i+2;
01605 } else {
01606 offset2[nevtemp] = 0;
01607 for (int j=i; j<curDim_; j++) {
01608 if (order[j] > order[i]) { offset2[nevtemp]++; }
01609 }
01610 order2[nevtemp] = order[i];
01611 i++;
01612 }
01613 nevtemp++;
01614 }
01615 ScalarType *ptr_h = H.values();
01616 ScalarType *ptr_q = Q.values();
01617 int ldh = H.stride(), ldq = Q.stride();
01618 int info = 0;
01619 for (i=nevtemp-1; i>=0; i--) {
01620 lapack.TREXC( compq, curDim_, ptr_h, ldh, ptr_q, ldq, order2[i]+1+offset2[i],
01621 1, &work[0], &info );
01622 TEST_FOR_EXCEPTION(info != 0, std::logic_error,
01623 "Anasazi::BlockKrylovSchur::computeSchurForm(): TREXC(n==" << curDim_ << ") returned info " << info << " != 0.");
01624 }
01625 }
01626
01628
01629 template <class ScalarType, class MV, class OP>
01630 void BlockKrylovSchur<ScalarType,MV,OP>::currentStatus(std::ostream &os)
01631 {
01632 using std::endl;
01633
01634 os.setf(std::ios::scientific, std::ios::floatfield);
01635 os.precision(6);
01636 os <<"================================================================================" << endl;
01637 os << endl;
01638 os <<" BlockKrylovSchur Solver Status" << endl;
01639 os << endl;
01640 os <<"The solver is "<<(initialized_ ? "initialized." : "not initialized.") << endl;
01641 os <<"The number of iterations performed is " <<iter_<<endl;
01642 os <<"The block size is " << blockSize_<<endl;
01643 os <<"The number of blocks is " << numBlocks_<<endl;
01644 os <<"The current basis size is " << curDim_<<endl;
01645 os <<"The number of auxiliary vectors is " << numAuxVecs_ << endl;
01646 os <<"The number of operations Op*x is "<<count_ApplyOp_<<endl;
01647
01648 os.setf(std::ios_base::right, std::ios_base::adjustfield);
01649
01650 os << endl;
01651 if (initialized_) {
01652 os <<"CURRENT RITZ VALUES "<<endl;
01653 if (ritzIndex_.size() != 0) {
01654 int numPrint = (curDim_ < numRitzPrint_? curDim_: numRitzPrint_);
01655 if (problem_->isHermitian()) {
01656 os << std::setw(20) << "Ritz Value"
01657 << std::setw(20) << "Ritz Residual"
01658 << endl;
01659 os <<"--------------------------------------------------------------------------------"<<endl;
01660 for (int i=0; i<numPrint; i++) {
01661 os << std::setw(20) << ritzValues_[i].realpart
01662 << std::setw(20) << ritzResiduals_[i]
01663 << endl;
01664 }
01665 } else {
01666 os << std::setw(24) << "Ritz Value"
01667 << std::setw(30) << "Ritz Residual"
01668 << endl;
01669 os <<"--------------------------------------------------------------------------------"<<endl;
01670 for (int i=0; i<numPrint; i++) {
01671
01672 os << std::setw(15) << ritzValues_[i].realpart;
01673 if (ritzValues_[i].imagpart < MT_ZERO) {
01674 os << " - i" << std::setw(15) << Teuchos::ScalarTraits<MagnitudeType>::magnitude(ritzValues_[i].imagpart);
01675 } else {
01676 os << " + i" << std::setw(15) << ritzValues_[i].imagpart;
01677 }
01678 os << std::setw(20) << ritzResiduals_[i] << endl;
01679 }
01680 }
01681 } else {
01682 os << std::setw(20) << "[ NONE COMPUTED ]" << endl;
01683 }
01684 }
01685 os << endl;
01686 os <<"================================================================================" << endl;
01687 os << endl;
01688 }
01689
01690 }
01691
01692 #endif
01693
01694