|
Anasazi Version of the Day
|
|
Anasazi Version of the Day
|
00001 // @HEADER 00002 // *********************************************************************** 00003 // 00004 // Anasazi: Block Eigensolvers Package 00005 // Copyright (2004) Sandia Corporation 00006 // 00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive 00008 // license for use of this work by or on behalf of the U.S. Government. 00009 // 00010 // This library is free software; you can redistribute it and/or modify 00011 // it under the terms of the GNU Lesser General Public License as 00012 // published by the Free Software Foundation; either version 2.1 of the 00013 // License, or (at your option) any later version. 00014 // 00015 // This library is distributed in the hope that it will be useful, but 00016 // WITHOUT ANY WARRANTY; without even the implied warranty of 00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00018 // Lesser General Public License for more details. 00019 // 00020 // You should have received a copy of the GNU Lesser General Public 00021 // License along with this library; if not, write to the Free Software 00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 00023 // USA 00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 00025 // 00026 // *********************************************************************** 00027 // @HEADER 00028 00033 #ifndef ANASAZI_RTRBASE_HPP 00034 #define ANASAZI_RTRBASE_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 "AnasaziGenOrthoManager.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 00101 namespace Anasazi { 00102 00104 00105 00110 template <class ScalarType, class MV> 00111 struct RTRState { 00113 Teuchos::RCP<const MV> X; 00115 Teuchos::RCP<const MV> AX; 00117 Teuchos::RCP<const MV> BX; 00119 Teuchos::RCP<const MV> R; 00121 Teuchos::RCP<const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > T; 00125 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType rho; 00126 RTRState() : X(Teuchos::null),AX(Teuchos::null),BX(Teuchos::null), 00127 R(Teuchos::null),T(Teuchos::null),rho(0) {}; 00128 }; 00129 00131 00133 00134 00148 class RTRRitzFailure : public AnasaziError {public: 00149 RTRRitzFailure(const std::string& what_arg) : AnasaziError(what_arg) 00150 {}}; 00151 00160 class RTRInitFailure : public AnasaziError {public: 00161 RTRInitFailure(const std::string& what_arg) : AnasaziError(what_arg) 00162 {}}; 00163 00180 class RTROrthoFailure : public AnasaziError {public: 00181 RTROrthoFailure(const std::string& what_arg) : AnasaziError(what_arg) 00182 {}}; 00183 00184 00186 00187 00188 template <class ScalarType, class MV, class OP> 00189 class RTRBase : public Eigensolver<ScalarType,MV,OP> { 00190 public: 00191 00193 00194 00200 RTRBase(const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem, 00201 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter, 00202 const Teuchos::RCP<OutputManager<ScalarType> > &printer, 00203 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester, 00204 const Teuchos::RCP<GenOrthoManager<ScalarType,MV,OP> > &ortho, 00205 Teuchos::ParameterList ¶ms, 00206 const std::string &solverLabel, bool skinnySolver 00207 ); 00208 00210 virtual ~RTRBase() {}; 00211 00213 00215 00216 00238 virtual void iterate() = 0; 00239 00264 void initialize(RTRState<ScalarType,MV> newstate); 00265 00269 void initialize(); 00270 00283 bool isInitialized() const; 00284 00292 RTRState<ScalarType,MV> getState() const; 00293 00295 00297 00298 00300 int getNumIters() const; 00301 00303 void resetNumIters(); 00304 00312 Teuchos::RCP<const MV> getRitzVectors(); 00313 00319 std::vector<Value<ScalarType> > getRitzValues(); 00320 00328 std::vector<int> getRitzIndex(); 00329 00335 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getResNorms(); 00336 00337 00343 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRes2Norms(); 00344 00345 00350 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRitzRes2Norms(); 00351 00352 00361 int getCurSubspaceDim() const; 00362 00365 int getMaxSubspaceDim() const; 00366 00368 00370 00371 00373 void setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test); 00374 00376 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > getStatusTest() const; 00377 00379 const Eigenproblem<ScalarType,MV,OP>& getProblem() const; 00380 00381 00390 void setBlockSize(int blockSize); 00391 00392 00394 int getBlockSize() const; 00395 00396 00417 void setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs); 00418 00420 Teuchos::Array<Teuchos::RCP<const MV> > getAuxVecs() const; 00421 00423 00425 00426 00428 virtual void currentStatus(std::ostream &os); 00429 00431 00432 protected: 00433 // 00434 // Convenience typedefs 00435 // 00436 typedef SolverUtils<ScalarType,MV,OP> Utils; 00437 typedef MultiVecTraits<ScalarType,MV> MVT; 00438 typedef OperatorTraits<ScalarType,MV,OP> OPT; 00439 typedef Teuchos::ScalarTraits<ScalarType> SCT; 00440 typedef typename SCT::magnitudeType MagnitudeType; 00441 typedef Teuchos::ScalarTraits<MagnitudeType> MAT; 00442 const MagnitudeType ONE; 00443 const MagnitudeType ZERO; 00444 const MagnitudeType NANVAL; 00445 typedef typename std::vector<MagnitudeType>::iterator vecMTiter; 00446 typedef typename std::vector<ScalarType>::iterator vecSTiter; 00447 // 00448 // Internal structs 00449 // 00450 struct CheckList { 00451 bool checkX, checkAX, checkBX; 00452 bool checkEta, checkAEta, checkBEta; 00453 bool checkR, checkQ, checkBR; 00454 bool checkZ, checkPBX; 00455 CheckList() : checkX(false),checkAX(false),checkBX(false), 00456 checkEta(false),checkAEta(false),checkBEta(false), 00457 checkR(false),checkQ(false),checkBR(false), 00458 checkZ(false), checkPBX(false) {}; 00459 }; 00460 // 00461 // Internal methods 00462 // 00463 std::string accuracyCheck(const CheckList &chk, const std::string &where) const; 00464 // solves the model minimization 00465 virtual void solveTRSubproblem() = 0; 00466 // Riemannian metric 00467 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType ginner(const MV &xi) const; 00468 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType ginner(const MV &xi, const MV &zeta) const; 00469 void ginnersep(const MV &xi, std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &d) const; 00470 void ginnersep(const MV &xi, const MV &zeta, std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &d) const; 00471 // 00472 // Classes input through constructor that define the eigenproblem to be solved. 00473 // 00474 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_; 00475 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sm_; 00476 const Teuchos::RCP<OutputManager<ScalarType> > om_; 00477 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > tester_; 00478 const Teuchos::RCP<GenOrthoManager<ScalarType,MV,OP> > orthman_; 00479 // 00480 // Information obtained from the eigenproblem 00481 // 00482 Teuchos::RCP<const OP> AOp_; 00483 Teuchos::RCP<const OP> BOp_; 00484 Teuchos::RCP<const OP> Prec_; 00485 bool hasBOp_, hasPrec_, olsenPrec_; 00486 // 00487 // Internal timers 00488 // 00489 Teuchos::RCP<Teuchos::Time> timerAOp_, timerBOp_, timerPrec_, 00490 timerSort_, 00491 timerLocalProj_, timerDS_, 00492 timerLocalUpdate_, timerCompRes_, 00493 timerOrtho_, timerInit_; 00494 // 00495 // Counters 00496 // 00497 // Number of operator applications 00498 int counterAOp_, counterBOp_, counterPrec_; 00499 00500 // 00501 // Algorithmic parameters. 00502 // 00503 // blockSize_ is the solver block size 00504 int blockSize_; 00505 // 00506 // Current solver state 00507 // 00508 // initialized_ specifies that the basis vectors have been initialized and the iterate() routine 00509 // is capable of running; _initialize is controlled by the initialize() member method 00510 // For the implications of the state of initialized_, please see documentation for initialize() 00511 bool initialized_; 00512 // 00513 // nevLocal_ reflects how much of the current basis is valid (0 <= nevLocal_ <= blockSize_) 00514 // this tells us how many of the values in theta_ are valid Ritz values 00515 int nevLocal_; 00516 // 00517 // are we implementing a skinny solver? (SkinnyIRTR) 00518 bool isSkinny_; 00519 // 00520 // are we computing leftmost or rightmost eigenvalue? 00521 bool leftMost_; 00522 // 00523 // State Multivecs 00524 // 00525 // if we are implementing a skinny solver (SkinnyIRTR), 00526 // then some of these will never be allocated 00527 // 00528 // In order to handle auxiliary vectors, we need to handle the projector 00529 // P_{[BQ BX],[BQ BX]} 00530 // Using an orthomanager with B-inner product, this requires calling with multivectors 00531 // [BQ,BX] and [Q,X]. 00532 // These multivectors must be combined because <[BQ,BX],[Q,X]>_B != I 00533 // Hence, we will create two multivectors V and BV, which store 00534 // V = [Q,X] 00535 // BV = [BQ,BX] 00536 // 00537 // In the context of preconditioning, we may need to apply the projector 00538 // P_{prec*[BQ,BX],[BQ,BX]} 00539 // Because [BQ,BX] do not change during the supproblem solver, we will apply 00540 // the preconditioner to [BQ,BX] only once. This result is stored in PBV. 00541 // 00542 // X,BX are views into V,BV 00543 // We don't need views for Q,BQ 00544 // Inside the subproblem solver, X,BX are constant, so we can allow these 00545 // views to exist alongside the full view of V,BV without violating 00546 // view semantics. 00547 // 00548 // Skinny solver allocates 6/7/8 multivectors: 00549 // V_, BV_ (if hasB) 00550 // PBV_ (if hasPrec and olsenPrec) 00551 // R_, Z_ (regardless of hasPrec) 00552 // eta_, delta_, Hdelta_ 00553 // 00554 // Hefty solver allocates 10/11/12/13 multivectors: 00555 // V_, AX_, BV_ (if hasB) 00556 // PBV_ (if hasPrec and olsenPrec) 00557 // R_, Z_ (if hasPrec) 00558 // eta_, Aeta_, Beta_ 00559 // delta_, Adelta_, Bdelta_, Hdelta_ 00560 // 00561 Teuchos::RCP<MV> V_, BV_, PBV_, // V = [Q,X]; B*V; Prec*B*V 00562 AX_, R_, // A*X_; residual, gradient, and residual of model minimization 00563 eta_, Aeta_, Beta_, // update vector from model minimization 00564 delta_, Adelta_, Bdelta_, Hdelta_, // search direction in model minimization 00565 Z_; // preconditioned residual 00566 Teuchos::RCP<const MV> X_, BX_; 00567 // 00568 // auxiliary vectors 00569 Teuchos::Array<Teuchos::RCP<const MV> > auxVecs_; 00570 int numAuxVecs_; 00571 // 00572 // Number of iterations that have been performed. 00573 int iter_; 00574 // 00575 // Current eigenvalues, residual norms 00576 std::vector<MagnitudeType> theta_, Rnorms_, R2norms_, ritz2norms_; 00577 // 00578 // are the residual norms current with the residual? 00579 bool Rnorms_current_, R2norms_current_; 00580 // 00581 // parameters solver and inner solver 00582 MagnitudeType conv_kappa_, conv_theta_; 00583 MagnitudeType rho_; 00584 // 00585 // current objective function value 00586 MagnitudeType fx_; 00587 }; 00588 00589 00590 00591 00593 // Constructor 00594 template <class ScalarType, class MV, class OP> 00595 RTRBase<ScalarType,MV,OP>::RTRBase( 00596 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem, 00597 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter, 00598 const Teuchos::RCP<OutputManager<ScalarType> > &printer, 00599 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester, 00600 const Teuchos::RCP<GenOrthoManager<ScalarType,MV,OP> > &ortho, 00601 Teuchos::ParameterList ¶ms, 00602 const std::string &solverLabel, bool skinnySolver 00603 ) : 00604 ONE(Teuchos::ScalarTraits<MagnitudeType>::one()), 00605 ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()), 00606 NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan()), 00607 // problem, tools 00608 problem_(problem), 00609 sm_(sorter), 00610 om_(printer), 00611 tester_(tester), 00612 orthman_(ortho), 00613 // timers, counters 00614 timerAOp_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Operation A*x")), 00615 timerBOp_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Operation B*x")), 00616 timerPrec_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Operation Prec*x")), 00617 timerSort_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Sorting eigenvalues")), 00618 timerLocalProj_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Local projection")), 00619 timerDS_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Direct solve")), 00620 timerLocalUpdate_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Local update")), 00621 timerCompRes_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Computing residuals")), 00622 timerOrtho_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Orthogonalization")), 00623 timerInit_(Teuchos::TimeMonitor::getNewTimer(solverLabel+"::Initialization")), 00624 counterAOp_(0), 00625 counterBOp_(0), 00626 counterPrec_(0), 00627 // internal data 00628 blockSize_(0), 00629 initialized_(false), 00630 nevLocal_(0), 00631 isSkinny_(skinnySolver), 00632 leftMost_(true), 00633 auxVecs_( Teuchos::Array<Teuchos::RCP<const MV> >(0) ), 00634 numAuxVecs_(0), 00635 iter_(0), 00636 Rnorms_current_(false), 00637 R2norms_current_(false), 00638 conv_kappa_(.1), 00639 conv_theta_(1), 00640 rho_( MAT::nan() ), 00641 fx_( MAT::nan() ) 00642 { 00643 TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument, 00644 "Anasazi::RTRBase::constructor: user passed null problem pointer."); 00645 TEST_FOR_EXCEPTION(sm_ == Teuchos::null,std::invalid_argument, 00646 "Anasazi::RTRBase::constructor: user passed null sort manager pointer."); 00647 TEST_FOR_EXCEPTION(om_ == Teuchos::null,std::invalid_argument, 00648 "Anasazi::RTRBase::constructor: user passed null output manager pointer."); 00649 TEST_FOR_EXCEPTION(tester_ == Teuchos::null,std::invalid_argument, 00650 "Anasazi::RTRBase::constructor: user passed null status test pointer."); 00651 TEST_FOR_EXCEPTION(orthman_ == Teuchos::null,std::invalid_argument, 00652 "Anasazi::RTRBase::constructor: user passed null orthogonalization manager pointer."); 00653 TEST_FOR_EXCEPTION(problem_->isProblemSet() == false, std::invalid_argument, 00654 "Anasazi::RTRBase::constructor: problem is not set."); 00655 TEST_FOR_EXCEPTION(problem_->isHermitian() == false, std::invalid_argument, 00656 "Anasazi::RTRBase::constructor: problem is not Hermitian."); 00657 00658 // get the problem operators 00659 AOp_ = problem_->getOperator(); 00660 TEST_FOR_EXCEPTION(AOp_ == Teuchos::null, std::invalid_argument, 00661 "Anasazi::RTRBase::constructor: problem provides no A matrix."); 00662 BOp_ = problem_->getM(); 00663 Prec_ = problem_->getPrec(); 00664 hasBOp_ = (BOp_ != Teuchos::null); 00665 hasPrec_ = (Prec_ != Teuchos::null); 00666 olsenPrec_ = params.get<bool>("Olsen Prec", true); 00667 00668 TEST_FOR_EXCEPTION(orthman_->getOp() != BOp_,std::invalid_argument, 00669 "Anasazi::RTRBase::constructor: orthogonalization manager must use mass matrix for inner product."); 00670 00671 // set the block size and allocate data 00672 int bs = params.get("Block Size", problem_->getNEV()); 00673 setBlockSize(bs); 00674 00675 // leftmost or rightmost? 00676 leftMost_ = params.get("Leftmost",leftMost_); 00677 00678 conv_kappa_ = params.get("Kappa Convergence",conv_kappa_); 00679 TEST_FOR_EXCEPTION(conv_kappa_ <= 0 || conv_kappa_ >= 1,std::invalid_argument, 00680 "Anasazi::RTRBase::constructor: kappa must be in (0,1)."); 00681 conv_theta_ = params.get("Theta Convergence",conv_theta_); 00682 TEST_FOR_EXCEPTION(conv_theta_ <= 0,std::invalid_argument, 00683 "Anasazi::RTRBase::constructor: theta must be strictly postitive."); 00684 } 00685 00686 00688 // Set the block size and make necessary adjustments. 00689 template <class ScalarType, class MV, class OP> 00690 void RTRBase<ScalarType,MV,OP>::setBlockSize (int blockSize) 00691 { 00692 // time spent here counts towards timerInit_ 00693 Teuchos::TimeMonitor lcltimer( *timerInit_ ); 00694 00695 // This routine only allocates space; it doesn't not perform any computation 00696 // if solver is initialized and size is to be decreased, take the first blockSize vectors of all to preserve state 00697 // otherwise, shrink/grow/allocate space and set solver to unitialized 00698 00699 Teuchos::RCP<const MV> tmp; 00700 // grab some Multivector to Clone 00701 // in practice, getInitVec() should always provide this, but it is possible to use a 00702 // Eigenproblem with nothing in getInitVec() by manually initializing with initialize(); 00703 // in case of that strange scenario, we will try to Clone from R_ 00704 // we like R_ for this, because it has minimal size (blockSize_), as opposed to V_ (blockSize_+numAuxVecs_) 00705 if (blockSize_ > 0) { 00706 tmp = R_; 00707 } 00708 else { 00709 tmp = problem_->getInitVec(); 00710 TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::logic_error, 00711 "Anasazi::RTRBase::setBlockSize(): Eigenproblem did not specify initial vectors to clone from"); 00712 } 00713 00714 TEST_FOR_EXCEPTION(blockSize <= 0 || blockSize > MVT::GetVecLength(*tmp), std::invalid_argument, 00715 "Anasazi::RTRBase::setBlockSize was passed a non-positive block size"); 00716 00717 // last chance to quit before causing side-effects 00718 if (blockSize == blockSize_) { 00719 // do nothing 00720 return; 00721 } 00722 00723 // clear views 00724 X_ = Teuchos::null; 00725 BX_ = Teuchos::null; 00726 00727 // regardless of whether we preserve any data, any content in V, BV and PBV corresponding to the 00728 // auxiliary vectors must be retained 00729 // go ahead and do these first 00730 // 00731 // two cases here: 00732 // * we are growing (possibly, from empty) 00733 // any aux data must be copied over, nothing from X need copying 00734 // * we are shrinking 00735 // any aux data must be copied over, go ahead and copy X material (initialized or not) 00736 // 00737 if (blockSize > blockSize_) 00738 { 00739 // GROWING 00740 // get a pointer for Q-related material, and an index for extracting/setting it 00741 Teuchos::RCP<const MV> Q; 00742 std::vector<int> indQ(numAuxVecs_); 00743 for (int i=0; i<numAuxVecs_; i++) indQ[i] = i; 00744 // if numAuxVecs_ > 0, then necessarily blockSize_ > 0 (we have already been allocated once) 00745 TEST_FOR_EXCEPTION(numAuxVecs_ > 0 && blockSize_ == 0, std::logic_error, 00746 "Anasazi::RTRBase::setSize(): logic error. Please contact Anasazi team."); 00747 // V 00748 if (numAuxVecs_ > 0) Q = MVT::CloneView(*V_,indQ); 00749 V_ = MVT::Clone(*tmp,numAuxVecs_ + blockSize); 00750 if (numAuxVecs_ > 0) MVT::SetBlock(*Q,indQ,*V_); 00751 // BV 00752 if (hasBOp_) { 00753 if (numAuxVecs_ > 0) Q = MVT::CloneView(*BV_,indQ); 00754 BV_ = MVT::Clone(*tmp,numAuxVecs_ + blockSize); 00755 if (numAuxVecs_ > 0) MVT::SetBlock(*Q,indQ,*BV_); 00756 } 00757 else { 00758 BV_ = V_; 00759 } 00760 // PBV 00761 if (hasPrec_ && olsenPrec_) { 00762 if (numAuxVecs_ > 0) Q = MVT::CloneView(*PBV_,indQ); 00763 PBV_ = MVT::Clone(*tmp,numAuxVecs_ + blockSize); 00764 if (numAuxVecs_ > 0) MVT::SetBlock(*Q,indQ,*PBV_); 00765 } 00766 } 00767 else 00768 { 00769 // SHRINKING 00770 std::vector<int> indV(numAuxVecs_+blockSize); 00771 for (int i=0; i<numAuxVecs_+blockSize; i++) indV[i] = i; 00772 // V 00773 V_ = MVT::CloneCopy(*V_,indV); 00774 // BV 00775 if (hasBOp_) { 00776 BV_ = MVT::CloneCopy(*BV_,indV); 00777 } 00778 else { 00779 BV_ = V_; 00780 } 00781 // PBV 00782 if (hasPrec_ && olsenPrec_) { 00783 PBV_ = MVT::CloneCopy(*PBV_,indV); 00784 } 00785 } 00786 00787 if (blockSize < blockSize_) { 00788 // shrink vectors 00789 blockSize_ = blockSize; 00790 00791 theta_.resize(blockSize_); 00792 ritz2norms_.resize(blockSize_); 00793 Rnorms_.resize(blockSize_); 00794 R2norms_.resize(blockSize_); 00795 00796 if (initialized_) { 00797 // shrink multivectors with copy 00798 std::vector<int> ind(blockSize_); 00799 for (int i=0; i<blockSize_; i++) ind[i] = i; 00800 00801 // Z can be deleted, no important info there 00802 Z_ = Teuchos::null; 00803 00804 // we will not use tmp for cloning, clear it and free some space 00805 tmp = Teuchos::null; 00806 00807 R_ = MVT::CloneCopy(*R_ ,ind); 00808 eta_ = MVT::CloneCopy(*eta_ ,ind); 00809 delta_ = MVT::CloneCopy(*delta_ ,ind); 00810 Hdelta_ = MVT::CloneCopy(*Hdelta_,ind); 00811 if (!isSkinny_) { 00812 AX_ = MVT::CloneCopy(*AX_ ,ind); 00813 Aeta_ = MVT::CloneCopy(*Aeta_ ,ind); 00814 Adelta_ = MVT::CloneCopy(*Adelta_,ind); 00815 } 00816 else { 00817 AX_ = Teuchos::null; 00818 Aeta_ = Teuchos::null; 00819 Adelta_ = Teuchos::null; 00820 } 00821 00822 if (hasBOp_) { 00823 if (!isSkinny_) { 00824 Beta_ = MVT::CloneCopy(*Beta_,ind); 00825 Bdelta_ = MVT::CloneCopy(*Bdelta_,ind); 00826 } 00827 else { 00828 Beta_ = Teuchos::null; 00829 Bdelta_ = Teuchos::null; 00830 } 00831 } 00832 else { 00833 Beta_ = eta_; 00834 Bdelta_ = delta_; 00835 } 00836 00837 // we need Z if we have a preconditioner 00838 // we also use Z for temp storage in the skinny solvers 00839 if (hasPrec_ || isSkinny_) { 00840 Z_ = MVT::Clone(*V_,blockSize_); 00841 } 00842 else { 00843 Z_ = R_; 00844 } 00845 00846 } 00847 else { 00848 // release previous multivectors 00849 // shrink multivectors without copying 00850 AX_ = Teuchos::null; 00851 R_ = Teuchos::null; 00852 eta_ = Teuchos::null; 00853 Aeta_ = Teuchos::null; 00854 delta_ = Teuchos::null; 00855 Adelta_ = Teuchos::null; 00856 Hdelta_ = Teuchos::null; 00857 Beta_ = Teuchos::null; 00858 Bdelta_ = Teuchos::null; 00859 Z_ = Teuchos::null; 00860 00861 R_ = MVT::Clone(*tmp,blockSize_); 00862 eta_ = MVT::Clone(*tmp,blockSize_); 00863 delta_ = MVT::Clone(*tmp,blockSize_); 00864 Hdelta_ = MVT::Clone(*tmp,blockSize_); 00865 if (!isSkinny_) { 00866 AX_ = MVT::Clone(*tmp,blockSize_); 00867 Aeta_ = MVT::Clone(*tmp,blockSize_); 00868 Adelta_ = MVT::Clone(*tmp,blockSize_); 00869 } 00870 00871 if (hasBOp_) { 00872 if (!isSkinny_) { 00873 Beta_ = MVT::Clone(*tmp,blockSize_); 00874 Bdelta_ = MVT::Clone(*tmp,blockSize_); 00875 } 00876 } 00877 else { 00878 Beta_ = eta_; 00879 Bdelta_ = delta_; 00880 } 00881 00882 // we need Z if we have a preconditioner 00883 // we also use Z for temp storage in the skinny solvers 00884 if (hasPrec_ || isSkinny_) { 00885 Z_ = MVT::Clone(*tmp,blockSize_); 00886 } 00887 else { 00888 Z_ = R_; 00889 } 00890 } // if initialized_ 00891 } // if blockSize is shrinking 00892 else { // if blockSize > blockSize_ 00893 // this is also the scenario for our initial call to setBlockSize(), in the constructor 00894 initialized_ = false; 00895 00896 // grow/allocate vectors 00897 theta_.resize(blockSize,NANVAL); 00898 ritz2norms_.resize(blockSize,NANVAL); 00899 Rnorms_.resize(blockSize,NANVAL); 00900 R2norms_.resize(blockSize,NANVAL); 00901 00902 // deallocate old multivectors 00903 AX_ = Teuchos::null; 00904 R_ = Teuchos::null; 00905 eta_ = Teuchos::null; 00906 Aeta_ = Teuchos::null; 00907 delta_ = Teuchos::null; 00908 Adelta_ = Teuchos::null; 00909 Hdelta_ = Teuchos::null; 00910 Beta_ = Teuchos::null; 00911 Bdelta_ = Teuchos::null; 00912 Z_ = Teuchos::null; 00913 00914 // clone multivectors off of tmp 00915 R_ = MVT::Clone(*tmp,blockSize); 00916 eta_ = MVT::Clone(*tmp,blockSize); 00917 delta_ = MVT::Clone(*tmp,blockSize); 00918 Hdelta_ = MVT::Clone(*tmp,blockSize); 00919 if (!isSkinny_) { 00920 AX_ = MVT::Clone(*tmp,blockSize); 00921 Aeta_ = MVT::Clone(*tmp,blockSize); 00922 Adelta_ = MVT::Clone(*tmp,blockSize); 00923 } 00924 00925 if (hasBOp_) { 00926 if (!isSkinny_) { 00927 Beta_ = MVT::Clone(*tmp,blockSize); 00928 Bdelta_ = MVT::Clone(*tmp,blockSize); 00929 } 00930 } 00931 else { 00932 Beta_ = eta_; 00933 Bdelta_ = delta_; 00934 } 00935 if (hasPrec_ || isSkinny_) { 00936 Z_ = MVT::Clone(*tmp,blockSize); 00937 } 00938 else { 00939 Z_ = R_; 00940 } 00941 blockSize_ = blockSize; 00942 } 00943 00944 // get view of X from V, BX from BV 00945 // these are located after the first numAuxVecs columns 00946 { 00947 std::vector<int> indX(blockSize_); 00948 for (int i=0; i<blockSize_; i++) indX[i] = numAuxVecs_+i; 00949 X_ = MVT::CloneView(*V_,indX); 00950 if (hasBOp_) { 00951 BX_ = MVT::CloneView(*BV_,indX); 00952 } 00953 else { 00954 BX_ = X_; 00955 } 00956 } 00957 } 00958 00959 00961 // Set a new StatusTest for the solver. 00962 template <class ScalarType, class MV, class OP> 00963 void RTRBase<ScalarType,MV,OP>::setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test) { 00964 TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument, 00965 "Anasazi::RTRBase::setStatusTest() was passed a null StatusTest."); 00966 tester_ = test; 00967 } 00968 00969 00971 // Get the current StatusTest used by the solver. 00972 template <class ScalarType, class MV, class OP> 00973 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > RTRBase<ScalarType,MV,OP>::getStatusTest() const { 00974 return tester_; 00975 } 00976 00977 00979 // Set the auxiliary vectors 00980 template <class ScalarType, class MV, class OP> 00981 void RTRBase<ScalarType,MV,OP>::setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs) { 00982 typedef typename Teuchos::Array<Teuchos::RCP<const MV> >::const_iterator tarcpmv; 00983 00984 // set new auxiliary vectors 00985 auxVecs_.resize(0); 00986 auxVecs_.reserve(auxvecs.size()); 00987 00988 numAuxVecs_ = 0; 00989 for (tarcpmv v=auxvecs.begin(); v != auxvecs.end(); ++v) { 00990 numAuxVecs_ += MVT::GetNumberVecs(**v); 00991 } 00992 00993 // If the solver has been initialized, X is not necessarily orthogonal to new auxiliary vectors 00994 if (numAuxVecs_ > 0 && initialized_) { 00995 initialized_ = false; 00996 } 00997 00998 // clear X,BX views 00999 X_ = Teuchos::null; 01000 BX_ = Teuchos::null; 01001 // deallocate, we'll clone off R_ below 01002 V_ = Teuchos::null; 01003 BV_ = Teuchos::null; 01004 PBV_ = Teuchos::null; 01005 01006 // put auxvecs into V, update BV and PBV if necessary 01007 if (numAuxVecs_ > 0) { 01008 V_ = MVT::Clone(*R_,numAuxVecs_ + blockSize_); 01009 int numsofar = 0; 01010 for (tarcpmv v=auxvecs.begin(); v != auxvecs.end(); ++v) { 01011 std::vector<int> ind(MVT::GetNumberVecs(**v)); 01012 for (size_t j=0; j<ind.size(); j++) ind[j] = numsofar++; 01013 MVT::SetBlock(**v,ind,*V_); 01014 auxVecs_.push_back(MVT::CloneView(*Teuchos::rcp_static_cast<const MV>(V_),ind)); 01015 } 01016 TEST_FOR_EXCEPTION(numsofar != numAuxVecs_, std::logic_error, 01017 "Anasazi::RTRBase::setAuxVecs(): Logic error. Please contact Anasazi team."); 01018 // compute B*V, Prec*B*V 01019 if (hasBOp_) { 01020 BV_ = MVT::Clone(*R_,numAuxVecs_ + blockSize_); 01021 OPT::Apply(*BOp_,*V_,*BV_); 01022 } 01023 else { 01024 BV_ = V_; 01025 } 01026 if (hasPrec_ && olsenPrec_) { 01027 PBV_ = MVT::Clone(*R_,numAuxVecs_ + blockSize_); 01028 OPT::Apply(*Prec_,*BV_,*V_); 01029 } 01030 } 01031 // 01032 01033 if (om_->isVerbosity( Debug ) ) { 01034 // Check almost everything here 01035 CheckList chk; 01036 chk.checkQ = true; 01037 om_->print( Debug, accuracyCheck(chk, "in setAuxVecs()") ); 01038 } 01039 } 01040 01041 01043 /* Initialize the state of the solver 01044 * 01045 * POST-CONDITIONS: 01046 * 01047 * initialized_ == true 01048 * X is orthonormal, orthogonal to auxVecs_ 01049 * AX = A*X if not skinnny 01050 * BX = B*X if hasBOp_ 01051 * theta_ contains Ritz values of X 01052 * R = AX - BX*diag(theta_) 01053 */ 01054 template <class ScalarType, class MV, class OP> 01055 void RTRBase<ScalarType,MV,OP>::initialize(RTRState<ScalarType,MV> newstate) 01056 { 01057 // NOTE: memory has been allocated by setBlockSize(). Use SetBlock below; do not Clone 01058 // NOTE: Overall time spent in this routine is counted to timerInit_; portions will also be counted towards other primitives 01059 01060 // clear const views to X,BX in V,BV 01061 // work with temporary non-const views 01062 X_ = Teuchos::null; 01063 BX_ = Teuchos::null; 01064 Teuchos::RCP<MV> X, BX; 01065 { 01066 std::vector<int> ind(blockSize_); 01067 for (int i=0; i<blockSize_; ++i) ind[i] = numAuxVecs_+i; 01068 X = MVT::CloneViewNonConst(*V_,ind); 01069 if (hasBOp_) { 01070 BX = MVT::CloneViewNonConst(*BV_,ind); 01071 } 01072 else { 01073 BX = X; 01074 } 01075 } 01076 01077 Teuchos::TimeMonitor inittimer( *timerInit_ ); 01078 01079 std::vector<int> bsind(blockSize_); 01080 for (int i=0; i<blockSize_; i++) bsind[i] = i; 01081 01082 // in RTR, X (the subspace iterate) is primary 01083 // the order of dependence follows like so. 01084 // --init-> X 01085 // --op apply-> AX,BX 01086 // --ritz analysis-> theta 01087 // 01088 // if the user specifies all data for a level, we will accept it. 01089 // otherwise, we will generate the whole level, and all subsequent levels. 01090 // 01091 // the data members are ordered based on dependence, and the levels are 01092 // partitioned according to the amount of work required to produce the 01093 // items in a level. 01094 // 01095 // inconsitent multivectors widths and lengths will not be tolerated, and 01096 // will be treated with exceptions. 01097 01098 // set up X, AX, BX: get them from "state" if user specified them 01099 if (newstate.X != Teuchos::null) { 01100 TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.X) != MVT::GetVecLength(*X), 01101 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): vector length of newstate.X not correct." ); 01102 // newstate.X must have blockSize_ vectors; any more will be ignored 01103 TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.X) < blockSize_, 01104 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.X must have at least block size vectors."); 01105 01106 // put data in X 01107 MVT::SetBlock(*newstate.X,bsind,*X); 01108 01109 // put data in AX 01110 // if we are implementing a skinny solver, then we don't have memory allocated for AX 01111 // in this case, point AX at Z (skinny solvers allocate Z) and use it for temporary storage 01112 // we will clear the pointer later 01113 if (isSkinny_) { 01114 AX_ = Z_; 01115 } 01116 if (newstate.AX != Teuchos::null) { 01117 TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.AX) != MVT::GetVecLength(*X), 01118 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): vector length of newstate.AX not correct." ); 01119 // newstate.AX must have blockSize_ vectors; any more will be ignored 01120 TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.AX) < blockSize_, 01121 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.AX must have at least block size vectors."); 01122 MVT::SetBlock(*newstate.AX,bsind,*AX_); 01123 } 01124 else { 01125 { 01126 Teuchos::TimeMonitor lcltimer( *timerAOp_ ); 01127 OPT::Apply(*AOp_,*X,*AX_); 01128 counterAOp_ += blockSize_; 01129 } 01130 // we generated AX; we will generate R as well 01131 newstate.R = Teuchos::null; 01132 } 01133 01134 // put data in BX 01135 // skinny solvers always allocate BX if hasB, so this is unconditionally appropriate 01136 if (hasBOp_) { 01137 if (newstate.BX != Teuchos::null) { 01138 TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.BX) != MVT::GetVecLength(*X), 01139 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): vector length of newstate.BX not correct." ); 01140 // newstate.BX must have blockSize_ vectors; any more will be ignored 01141 TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.BX) < blockSize_, 01142 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.BX must have at least block size vectors."); 01143 MVT::SetBlock(*newstate.BX,bsind,*BX); 01144 } 01145 else { 01146 { 01147 Teuchos::TimeMonitor lcltimer( *timerBOp_ ); 01148 OPT::Apply(*BOp_,*X,*BX); 01149 counterBOp_ += blockSize_; 01150 } 01151 // we generated BX; we will generate R as well 01152 newstate.R = Teuchos::null; 01153 } 01154 } 01155 else { 01156 // the assignment BX_==X_ would be redundant; take advantage of this opportunity to debug a little 01157 TEST_FOR_EXCEPTION(BX != X, std::logic_error, "Anasazi::RTRBase::initialize(): solver invariant not satisfied (BX==X)."); 01158 } 01159 01160 } 01161 else { 01162 // user did not specify X 01163 01164 // clear state so we won't use any data from it below 01165 newstate.R = Teuchos::null; 01166 newstate.T = Teuchos::null; 01167 01168 // generate something and projectAndNormalize 01169 Teuchos::RCP<const MV> ivec = problem_->getInitVec(); 01170 TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::logic_error, 01171 "Anasazi::RTRBase::initialize(): Eigenproblem did not specify initial vectors to clone from."); 01172 01173 int initSize = MVT::GetNumberVecs(*ivec); 01174 if (initSize > blockSize_) { 01175 // we need only the first blockSize_ vectors from ivec; get a view of them 01176 initSize = blockSize_; 01177 std::vector<int> ind(blockSize_); 01178 for (int i=0; i<blockSize_; i++) ind[i] = i; 01179 ivec = MVT::CloneView(*ivec,ind); 01180 } 01181 01182 // assign ivec to first part of X 01183 if (initSize > 0) { 01184 std::vector<int> ind(initSize); 01185 for (int i=0; i<initSize; i++) ind[i] = i; 01186 MVT::SetBlock(*ivec,ind,*X); 01187 } 01188 // fill the rest of X with random 01189 if (blockSize_ > initSize) { 01190 std::vector<int> ind(blockSize_ - initSize); 01191 for (int i=0; i<blockSize_ - initSize; i++) ind[i] = initSize + i; 01192 Teuchos::RCP<MV> rX = MVT::CloneViewNonConst(*X,ind); 01193 MVT::MvRandom(*rX); 01194 rX = Teuchos::null; 01195 } 01196 01197 // put data in BX 01198 if (hasBOp_) { 01199 Teuchos::TimeMonitor lcltimer( *timerBOp_ ); 01200 OPT::Apply(*BOp_,*X,*BX); 01201 counterBOp_ += blockSize_; 01202 } 01203 else { 01204 // the assignment BX==X would be redundant; take advantage of this opportunity to debug a little 01205 TEST_FOR_EXCEPTION(BX != X, std::logic_error, "Anasazi::RTRBase::initialize(): solver invariant not satisfied (BX==X)."); 01206 } 01207 01208 // remove auxVecs from X and normalize it 01209 if (numAuxVecs_ > 0) { 01210 Teuchos::TimeMonitor lcltimer( *timerOrtho_ ); 01211 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC; 01212 int rank = orthman_->projectAndNormalizeMat(*X,auxVecs_,dummyC,Teuchos::null,BX); 01213 TEST_FOR_EXCEPTION(rank != blockSize_, RTRInitFailure, 01214 "Anasazi::RTRBase::initialize(): Couldn't generate initial basis of full rank."); 01215 } 01216 else { 01217 Teuchos::TimeMonitor lcltimer( *timerOrtho_ ); 01218 int rank = orthman_->normalizeMat(*X,Teuchos::null,BX); 01219 TEST_FOR_EXCEPTION(rank != blockSize_, RTRInitFailure, 01220 "Anasazi::RTRBase::initialize(): Couldn't generate initial basis of full rank."); 01221 } 01222 01223 // put data in AX 01224 if (isSkinny_) { 01225 AX_ = Z_; 01226 } 01227 { 01228 Teuchos::TimeMonitor lcltimer( *timerAOp_ ); 01229 OPT::Apply(*AOp_,*X,*AX_); 01230 counterAOp_ += blockSize_; 01231 } 01232 01233 } // end if (newstate.X != Teuchos::null) 01234 01235 01236 // set up Ritz values 01237 if (newstate.T != Teuchos::null) { 01238 TEST_FOR_EXCEPTION( (signed int)(newstate.T->size()) < blockSize_, 01239 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.T must contain at least block size Ritz values."); 01240 for (int i=0; i<blockSize_; i++) { 01241 theta_[i] = (*newstate.T)[i]; 01242 } 01243 } 01244 else { 01245 // get ritz vecs/vals 01246 Teuchos::SerialDenseMatrix<int,ScalarType> AA(blockSize_,blockSize_), 01247 BB(blockSize_,blockSize_), 01248 S(blockSize_,blockSize_); 01249 // project A 01250 { 01251 Teuchos::TimeMonitor lcltimer( *timerLocalProj_ ); 01252 MVT::MvTransMv(ONE,*X,*AX_,AA); 01253 if (hasBOp_) { 01254 MVT::MvTransMv(ONE,*X,*BX,BB); 01255 } 01256 } 01257 nevLocal_ = blockSize_; 01258 01259 // solve the projected problem 01260 // project B 01261 int ret; 01262 if (hasBOp_) { 01263 Teuchos::TimeMonitor lcltimer( *timerDS_ ); 01264 ret = Utils::directSolver(blockSize_,AA,Teuchos::rcpFromRef(BB),S,theta_,nevLocal_,1); 01265 } 01266 else { 01267 Teuchos::TimeMonitor lcltimer( *timerDS_ ); 01268 ret = Utils::directSolver(blockSize_,AA,Teuchos::null,S,theta_,nevLocal_,10); 01269 } 01270 TEST_FOR_EXCEPTION(ret != 0,RTRInitFailure, 01271 "Anasazi::RTRBase::initialize(): failure solving projected eigenproblem after retraction. LAPACK returns " << ret); 01272 TEST_FOR_EXCEPTION(nevLocal_ != blockSize_,RTRInitFailure,"Anasazi::RTRBase::initialize(): retracted iterate failed in Ritz analysis."); 01273 01274 // We only have blockSize_ ritz pairs, ergo we do not need to select. 01275 // However, we still require them to be ordered correctly 01276 { 01277 Teuchos::TimeMonitor lcltimer( *timerSort_ ); 01278 01279 std::vector<int> order(blockSize_); 01280 // 01281 // sort the first blockSize_ values in theta_ 01282 sm_->sort(theta_, Teuchos::rcpFromRef(order), blockSize_); // don't catch exception 01283 // 01284 // apply the same ordering to the primitive ritz vectors 01285 Utils::permuteVectors(order,S); 01286 } 01287 01288 // compute ritz residual norms 01289 { 01290 Teuchos::BLAS<int,ScalarType> blas; 01291 Teuchos::SerialDenseMatrix<int,ScalarType> RR(blockSize_,blockSize_); 01292 // RR = AA*S - BB*S*diag(theta) 01293 int info; 01294 if (hasBOp_) { 01295 info = RR.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,BB,S,ZERO); 01296 TEST_FOR_EXCEPTION(info != 0, std::logic_error, "Anasazi::RTRBase::initialize(): Logic error calling SerialDenseMatrix::multiply."); 01297 } 01298 else { 01299 RR.assign(S); 01300 } 01301 for (int i=0; i<blockSize_; i++) { 01302 blas.SCAL(blockSize_,theta_[i],RR[i],1); 01303 } 01304 info = RR.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,AA,S,-ONE); 01305 TEST_FOR_EXCEPTION(info != 0, std::logic_error, "Anasazi::RTRBase::initialize(): Logic error calling SerialDenseMatrix::multiply."); 01306 for (int i=0; i<blockSize_; i++) { 01307 ritz2norms_[i] = blas.NRM2(blockSize_,RR[i],1); 01308 } 01309 } 01310 01311 // update the solution 01312 // use R as local work space 01313 // Z may already be in use as work space (holding AX if isSkinny) 01314 { 01315 Teuchos::TimeMonitor lcltimer( *timerLocalUpdate_ ); 01316 // X <- X*S 01317 MVT::MvAddMv( ONE, *X, ZERO, *X, *R_ ); 01318 MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *X ); 01319 // AX <- AX*S 01320 MVT::MvAddMv( ONE, *AX_, ZERO, *AX_, *R_ ); 01321 MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *AX_ ); 01322 if (hasBOp_) { 01323 // BX <- BX*S 01324 MVT::MvAddMv( ONE, *BX, ZERO, *BX, *R_ ); 01325 MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *BX ); 01326 } 01327 } 01328 } 01329 01330 // done modifying X,BX; clear temp views and setup const views 01331 X = Teuchos::null; 01332 BX = Teuchos::null; 01333 { 01334 std::vector<int> ind(blockSize_); 01335 for (int i=0; i<blockSize_; ++i) ind[i] = numAuxVecs_+i; 01336 this->X_ = MVT::CloneView(static_cast<const MV&>(*V_),ind); 01337 if (hasBOp_) { 01338 this->BX_ = MVT::CloneView(static_cast<const MV&>(*BV_),ind); 01339 } 01340 else { 01341 this->BX_ = this->X_; 01342 } 01343 } 01344 01345 01346 // get objective function value 01347 fx_ = std::accumulate(theta_.begin(),theta_.end(),ZERO); 01348 01349 // set up R 01350 if (newstate.R != Teuchos::null) { 01351 TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) < blockSize_, 01352 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.R must have blockSize number of vectors." ); 01353 TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.R) != MVT::GetVecLength(*R_), 01354 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): vector length of newstate.R not correct." ); 01355 MVT::SetBlock(*newstate.R,bsind,*R_); 01356 } 01357 else { 01358 Teuchos::TimeMonitor lcltimer( *timerCompRes_ ); 01359 // form R <- AX - BX*T 01360 MVT::MvAddMv(ZERO,*AX_,ONE,*AX_,*R_); 01361 Teuchos::SerialDenseMatrix<int,ScalarType> T(blockSize_,blockSize_); 01362 T.putScalar(ZERO); 01363 for (int i=0; i<blockSize_; i++) T(i,i) = theta_[i]; 01364 MVT::MvTimesMatAddMv(-ONE,*BX_,T,ONE,*R_); 01365 } 01366 01367 // R has been updated; mark the norms as out-of-date 01368 Rnorms_current_ = false; 01369 R2norms_current_ = false; 01370 01371 // if isSkinny, then AX currently points to Z for temp storage 01372 // set AX back to null 01373 if (isSkinny_) { 01374 AX_ = Teuchos::null; 01375 } 01376 01377 // finally, we are initialized 01378 initialized_ = true; 01379 01380 if (om_->isVerbosity( Debug ) ) { 01381 // Check almost everything here 01382 CheckList chk; 01383 chk.checkX = true; 01384 chk.checkAX = true; 01385 chk.checkBX = true; 01386 chk.checkR = true; 01387 chk.checkQ = true; 01388 om_->print( Debug, accuracyCheck(chk, "after initialize()") ); 01389 } 01390 } 01391 01392 template <class ScalarType, class MV, class OP> 01393 void RTRBase<ScalarType,MV,OP>::initialize() 01394 { 01395 RTRState<ScalarType,MV> empty; 01396 initialize(empty); 01397 } 01398 01399 01400 01401 01403 // compute/return residual B-norms 01404 template <class ScalarType, class MV, class OP> 01405 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> 01406 RTRBase<ScalarType,MV,OP>::getResNorms() { 01407 if (Rnorms_current_ == false) { 01408 // Update the residual norms 01409 orthman_->norm(*R_,Rnorms_); 01410 Rnorms_current_ = true; 01411 } 01412 return Rnorms_; 01413 } 01414 01415 01417 // compute/return residual 2-norms 01418 template <class ScalarType, class MV, class OP> 01419 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> 01420 RTRBase<ScalarType,MV,OP>::getRes2Norms() { 01421 if (R2norms_current_ == false) { 01422 // Update the residual 2-norms 01423 MVT::MvNorm(*R_,R2norms_); 01424 R2norms_current_ = true; 01425 } 01426 return R2norms_; 01427 } 01428 01429 01430 01431 01433 // Check accuracy, orthogonality, and other debugging stuff 01434 // 01435 // bools specify which tests we want to run (instead of running more than we actually care about) 01436 // 01437 // we don't bother checking the following because they are computed explicitly: 01438 // AH == A*H 01439 // 01440 // 01441 // checkX : X orthonormal 01442 // orthogonal to auxvecs 01443 // checkAX : check AX == A*X 01444 // checkBX : check BX == B*X 01445 // checkEta : check that Eta'*B*X == 0 01446 // orthogonal to auxvecs 01447 // checkAEta : check that AEta = A*Eta 01448 // checkBEta : check that BEta = B*Eta 01449 // checkR : check R orthogonal to X 01450 // checkBR : check R B-orthogonal to X, valid in and immediately after solveTRSubproblem 01451 // checkQ : check that auxiliary vectors are actually orthonormal 01452 // 01453 // TODO: 01454 // add checkTheta 01455 // 01456 template <class ScalarType, class MV, class OP> 01457 std::string RTRBase<ScalarType,MV,OP>::accuracyCheck( const CheckList &chk, const std::string &where ) const 01458 { 01459 using std::setprecision; 01460 using std::scientific; 01461 using std::setw; 01462 using std::endl; 01463 std::stringstream os; 01464 MagnitudeType tmp; 01465 01466 os << " Debugging checks: " << where << endl; 01467 01468 // X and friends 01469 if (chk.checkX && initialized_) { 01470 tmp = orthman_->orthonormError(*X_); 01471 os << " >> Error in X^H B X == I : " << scientific << setprecision(10) << tmp << endl; 01472 for (Array_size_type i=0; i<auxVecs_.size(); i++) { 01473 tmp = orthman_->orthogError(*X_,*auxVecs_[i]); 01474 os << " >> Error in X^H B Q[" << i << "] == 0 : " << scientific << setprecision(10) << tmp << endl; 01475 } 01476 } 01477 if (chk.checkAX && !isSkinny_ && initialized_) { 01478 tmp = Utils::errorEquality(*X_, *AX_, AOp_); 01479 os << " >> Error in AX == A*X : " << scientific << setprecision(10) << tmp << endl; 01480 } 01481 if (chk.checkBX && hasBOp_ && initialized_) { 01482 Teuchos::RCP<MV> BX = MVT::Clone(*this->X_,this->blockSize_); 01483 OPT::Apply(*BOp_,*this->X_,*BX); 01484 tmp = Utils::errorEquality(*BX, *BX_); 01485 os << " >> Error in BX == B*X : " << scientific << setprecision(10) << tmp << endl; 01486 } 01487 01488 // Eta and friends 01489 if (chk.checkEta && initialized_) { 01490 tmp = orthman_->orthogError(*X_,*eta_); 01491 os << " >> Error in X^H B Eta == 0 : " << scientific << setprecision(10) << tmp << endl; 01492 for (Array_size_type i=0; i<auxVecs_.size(); i++) { 01493 tmp = orthman_->orthogError(*eta_,*auxVecs_[i]); 01494 os << " >> Error in Eta^H B Q[" << i << "]==0 : " << scientific << setprecision(10) << tmp << endl; 01495 } 01496 } 01497 if (chk.checkAEta && !isSkinny_ && initialized_) { 01498 tmp = Utils::errorEquality(*eta_, *Aeta_, AOp_); 01499 os << " >> Error in AEta == A*Eta : " << scientific << setprecision(10) << tmp << endl; 01500 } 01501 if (chk.checkBEta && !isSkinny_ && hasBOp_ && initialized_) { 01502 tmp = Utils::errorEquality(*eta_, *Beta_, BOp_); 01503 os << " >> Error in BEta == B*Eta : " << scientific << setprecision(10) << tmp << endl; 01504 } 01505 01506 // R: this is not B-orthogonality, but standard euclidean orthogonality 01507 if (chk.checkR && initialized_) { 01508 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_); 01509 MVT::MvTransMv(ONE,*X_,*R_,xTx); 01510 tmp = xTx.normFrobenius(); 01511 os << " >> Error in X^H R == 0 : " << scientific << setprecision(10) << tmp << endl; 01512 } 01513 01514 // BR: this is B-orthogonality: this is only valid inside and immediately after solveTRSubproblem 01515 // only check if B != I, otherwise, it is equivalent to the above test 01516 if (chk.checkBR && hasBOp_ && initialized_) { 01517 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_); 01518 MVT::MvTransMv(ONE,*BX_,*R_,xTx); 01519 tmp = xTx.normFrobenius(); 01520 os << " >> Error in X^H B R == 0 : " << scientific << setprecision(10) << tmp << endl; 01521 } 01522 01523 // Z: Z is preconditioned R, should be on tangent plane 01524 if (chk.checkZ && initialized_) { 01525 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_); 01526 MVT::MvTransMv(ONE,*BX_,*Z_,xTx); 01527 tmp = xTx.normFrobenius(); 01528 os << " >> Error in X^H B Z == 0 : " << scientific << setprecision(10) << tmp << endl; 01529 } 01530 01531 // Q 01532 if (chk.checkQ) { 01533 for (Array_size_type i=0; i<auxVecs_.size(); i++) { 01534 tmp = orthman_->orthonormError(*auxVecs_[i]); 01535 os << " >> Error in Q[" << i << "]^H B Q[" << i << "]==I: " << scientific << setprecision(10) << tmp << endl; 01536 for (Array_size_type j=i+1; j<auxVecs_.size(); j++) { 01537 tmp = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]); 01538 os << " >> Error in Q[" << i << "]^H B Q[" << j << "]==0: " << scientific << setprecision(10) << tmp << endl; 01539 } 01540 } 01541 } 01542 os << endl; 01543 return os.str(); 01544 } 01545 01546 01548 // Print the current status of the solver 01549 template <class ScalarType, class MV, class OP> 01550 void 01551 RTRBase<ScalarType,MV,OP>::currentStatus(std::ostream &os) 01552 { 01553 using std::setprecision; 01554 using std::scientific; 01555 using std::setw; 01556 using std::endl; 01557 01558 os <<endl; 01559 os <<"================================================================================" << endl; 01560 os << endl; 01561 os <<" RTRBase Solver Status" << endl; 01562 os << endl; 01563 os <<"The solver is "<<(initialized_ ? "initialized." : "not initialized.") << endl; 01564 os <<"The number of iterations performed is " << iter_ << endl; 01565 os <<"The current block size is " << blockSize_ << endl; 01566 os <<"The number of auxiliary vectors is " << numAuxVecs_ << endl; 01567 os <<"The number of operations A*x is " << counterAOp_ << endl; 01568 os <<"The number of operations B*x is " << counterBOp_ << endl; 01569 os <<"The number of operations Prec*x is " << counterPrec_ << endl; 01570 os <<"The most recent rho was " << scientific << setprecision(10) << rho_ << endl; 01571 os <<"The current value of f(x) is " << scientific << setprecision(10) << fx_ << endl; 01572 01573 if (initialized_) { 01574 os << endl; 01575 os <<"CURRENT EIGENVALUE ESTIMATES "<<endl; 01576 os << setw(20) << "Eigenvalue" 01577 << setw(20) << "Residual(B)" 01578 << setw(20) << "Residual(2)" 01579 << endl; 01580 os <<"--------------------------------------------------------------------------------"<<endl; 01581 for (int i=0; i<blockSize_; i++) { 01582 os << scientific << setprecision(10) << setw(20) << theta_[i]; 01583 if (Rnorms_current_) os << scientific << setprecision(10) << setw(20) << Rnorms_[i]; 01584 else os << scientific << setprecision(10) << setw(20) << "not current"; 01585 if (R2norms_current_) os << scientific << setprecision(10) << setw(20) << R2norms_[i]; 01586 else os << scientific << setprecision(10) << setw(20) << "not current"; 01587 os << endl; 01588 } 01589 } 01590 os <<"================================================================================" << endl; 01591 os << endl; 01592 } 01593 01594 01596 // Inner product 1 01597 template <class ScalarType, class MV, class OP> 01598 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 01599 RTRBase<ScalarType,MV,OP>::ginner(const MV &xi) const 01600 { 01601 std::vector<MagnitudeType> d(MVT::GetNumberVecs(xi)); 01602 MVT::MvNorm(xi,d); 01603 MagnitudeType ret = 0; 01604 for (vecMTiter i=d.begin(); i != d.end(); ++i) { 01605 ret += (*i)*(*i); 01606 } 01607 return ret; 01608 } 01609 01610 01612 // Inner product 2 01613 template <class ScalarType, class MV, class OP> 01614 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 01615 RTRBase<ScalarType,MV,OP>::ginner(const MV &xi, const MV &zeta) const 01616 { 01617 std::vector<ScalarType> d(MVT::GetNumberVecs(xi)); 01618 MVT::MvDot(xi,zeta,d); 01619 return SCT::real(std::accumulate(d.begin(),d.end(),SCT::zero())); 01620 } 01621 01622 01624 // Inner product 1 without trace accumulation 01625 template <class ScalarType, class MV, class OP> 01626 void RTRBase<ScalarType,MV,OP>::ginnersep( 01627 const MV &xi, 01628 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &d) const 01629 { 01630 MVT::MvNorm(xi,d); 01631 for (vecMTiter i=d.begin(); i != d.end(); ++i) { 01632 (*i) = (*i)*(*i); 01633 } 01634 } 01635 01636 01638 // Inner product 2 without trace accumulation 01639 template <class ScalarType, class MV, class OP> 01640 void RTRBase<ScalarType,MV,OP>::ginnersep( 01641 const MV &xi, const MV &zeta, 01642 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &d) const 01643 { 01644 std::vector<ScalarType> dC(MVT::GetNumberVecs(xi)); 01645 MVT::MvDot(xi,zeta,dC); 01646 vecMTiter iR=d.begin(); 01647 vecSTiter iS=dC.begin(); 01648 for (; iR != d.end(); iR++, iS++) { 01649 (*iR) = SCT::real(*iS); 01650 } 01651 } 01652 01653 template <class ScalarType, class MV, class OP> 01654 Teuchos::Array<Teuchos::RCP<const MV> > RTRBase<ScalarType,MV,OP>::getAuxVecs() const { 01655 return auxVecs_; 01656 } 01657 01658 template <class ScalarType, class MV, class OP> 01659 int RTRBase<ScalarType,MV,OP>::getBlockSize() const { 01660 return(blockSize_); 01661 } 01662 01663 template <class ScalarType, class MV, class OP> 01664 const Eigenproblem<ScalarType,MV,OP>& RTRBase<ScalarType,MV,OP>::getProblem() const { 01665 return(*problem_); 01666 } 01667 01668 template <class ScalarType, class MV, class OP> 01669 int RTRBase<ScalarType,MV,OP>::getMaxSubspaceDim() const { 01670 return blockSize_; 01671 } 01672 01673 template <class ScalarType, class MV, class OP> 01674 int RTRBase<ScalarType,MV,OP>::getCurSubspaceDim() const 01675 { 01676 if (!initialized_) return 0; 01677 return nevLocal_; 01678 } 01679 01680 template <class ScalarType, class MV, class OP> 01681 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> 01682 RTRBase<ScalarType,MV,OP>::getRitzRes2Norms() 01683 { 01684 std::vector<MagnitudeType> ret = ritz2norms_; 01685 ret.resize(nevLocal_); 01686 return ret; 01687 } 01688 01689 template <class ScalarType, class MV, class OP> 01690 std::vector<Value<ScalarType> > 01691 RTRBase<ScalarType,MV,OP>::getRitzValues() 01692 { 01693 std::vector<Value<ScalarType> > ret(nevLocal_); 01694 for (int i=0; i<nevLocal_; i++) { 01695 ret[i].realpart = theta_[i]; 01696 ret[i].imagpart = ZERO; 01697 } 01698 return ret; 01699 } 01700 01701 template <class ScalarType, class MV, class OP> 01702 Teuchos::RCP<const MV> 01703 RTRBase<ScalarType,MV,OP>::getRitzVectors() 01704 { 01705 return X_; 01706 } 01707 01708 template <class ScalarType, class MV, class OP> 01709 void RTRBase<ScalarType,MV,OP>::resetNumIters() 01710 { 01711 iter_=0; 01712 } 01713 01714 template <class ScalarType, class MV, class OP> 01715 int RTRBase<ScalarType,MV,OP>::getNumIters() const 01716 { 01717 return iter_; 01718 } 01719 01720 template <class ScalarType, class MV, class OP> 01721 RTRState<ScalarType,MV> RTRBase<ScalarType,MV,OP>::getState() const 01722 { 01723 RTRState<ScalarType,MV> state; 01724 state.X = X_; 01725 state.AX = AX_; 01726 if (hasBOp_) { 01727 state.BX = BX_; 01728 } 01729 else { 01730 state.BX = Teuchos::null; 01731 } 01732 state.rho = rho_; 01733 state.R = R_; 01734 state.T = Teuchos::rcp(new std::vector<MagnitudeType>(theta_)); 01735 return state; 01736 } 01737 01738 template <class ScalarType, class MV, class OP> 01739 bool RTRBase<ScalarType,MV,OP>::isInitialized() const 01740 { 01741 return initialized_; 01742 } 01743 01744 template <class ScalarType, class MV, class OP> 01745 std::vector<int> RTRBase<ScalarType,MV,OP>::getRitzIndex() 01746 { 01747 std::vector<int> ret(nevLocal_,0); 01748 return ret; 01749 } 01750 01751 01752 } // end Anasazi namespace 01753 01754 #endif // ANASAZI_RTRBASE_HPP
1.7.4
