AnasaziBlockDavidsonSolMgr.hpp

Go to the documentation of this file.
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 
00029 #ifndef ANASAZI_BLOCKDAVIDSON_SOLMGR_HPP
00030 #define ANASAZI_BLOCKDAVIDSON_SOLMGR_HPP
00031 
00036 #include "AnasaziConfigDefs.hpp"
00037 #include "AnasaziTypes.hpp"
00038 
00039 #include "AnasaziEigenproblem.hpp"
00040 #include "AnasaziSolverManager.hpp"
00041 #include "AnasaziSolverUtils.hpp"
00042 
00043 #include "AnasaziBlockDavidson.hpp"
00044 #include "AnasaziBasicSort.hpp"
00045 #include "AnasaziSVQBOrthoManager.hpp"
00046 #include "AnasaziBasicOrthoManager.hpp"
00047 #include "AnasaziStatusTestResNorm.hpp"
00048 #include "AnasaziStatusTestWithOrdering.hpp"
00049 #include "AnasaziStatusTestCombo.hpp"
00050 #include "AnasaziStatusTestOutput.hpp"
00051 #include "AnasaziBasicOutputManager.hpp"
00052 #include "Teuchos_BLAS.hpp"
00053 #include "Teuchos_LAPACK.hpp"
00054 #include "Teuchos_TimeMonitor.hpp"
00055 #ifdef TEUCHOS_DEBUG
00056 #  include <Teuchos_FancyOStream.hpp>
00057 #endif
00058 #ifdef HAVE_MPI
00059 #include <mpi.h>
00060 #endif
00061 
00062 
00063 
00077 namespace Anasazi {
00078 
00079 
00112 template<class ScalarType, class MV, class OP>
00113 class BlockDavidsonSolMgr : public SolverManager<ScalarType,MV,OP> {
00114 
00115   private:
00116     typedef MultiVecTraits<ScalarType,MV> MVT;
00117     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00118     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00119     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00120     typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00121     
00122   public:
00123 
00125 
00126 
00149   BlockDavidsonSolMgr( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00150                              Teuchos::ParameterList &pl );
00151 
00153   virtual ~BlockDavidsonSolMgr() {};
00155   
00157 
00158 
00160   const Eigenproblem<ScalarType,MV,OP>& getProblem() const {
00161     return *problem_;
00162   }
00163 
00165   int getNumIters() const { 
00166     return numIters_; 
00167   }
00168 
00176    Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00177      return tuple(_timerSolve, _timerRestarting, _timerLocking);
00178    }
00179 
00181 
00183 
00184     
00206   ReturnType solve();
00207 
00209   void setGlobalStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global);
00210 
00212   const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getGlobalStatusTest() const;
00213 
00215   void setLockingStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &locking);
00216 
00218   const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getLockingStatusTest() const;
00219 
00221   void setDebugStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &debug);
00222 
00224   const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getDebugStatusTest() const;
00225 
00227 
00228   private:
00229   Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
00230 
00231   std::string whch_, ortho_;
00232 
00233   MagnitudeType convtol_, locktol_;
00234   int maxRestarts_;
00235   bool useLocking_;
00236   bool relconvtol_, rellocktol_;
00237   int blockSize_, numBlocks_, numIters_;
00238   int maxLocked_;
00239   int lockQuorum_;
00240   bool inSituRestart_;
00241   int numRestartBlocks_;
00242   enum StatusTestResNorm<ScalarType,MV,OP>::ResType convNorm_, lockNorm_;
00243 
00244   Teuchos::RCP<Teuchos::Time> _timerSolve, _timerRestarting, _timerLocking;
00245 
00246   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > globalTest_;
00247   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > lockingTest_; 
00248   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > debugTest_;
00249   
00250   Teuchos::RCP<BasicOutputManager<ScalarType> > printer_;
00251 };
00252 
00253 
00254 // Constructor
00255 template<class ScalarType, class MV, class OP>
00256 BlockDavidsonSolMgr<ScalarType,MV,OP>::BlockDavidsonSolMgr( 
00257         const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00258         Teuchos::ParameterList &pl ) : 
00259   problem_(problem),
00260   whch_("SR"),
00261   ortho_("SVQB"),
00262   convtol_(MT::prec()),
00263   maxRestarts_(20),
00264   useLocking_(false),
00265   relconvtol_(true),
00266   rellocktol_(true),
00267   blockSize_(0),
00268   numBlocks_(0),
00269   numIters_(0),
00270   maxLocked_(0),
00271   lockQuorum_(1),
00272   inSituRestart_(false),
00273   numRestartBlocks_(1),
00274   _timerSolve(Teuchos::TimeMonitor::getNewTimer("BlockDavidsonSolMgr::solve()")),
00275   _timerRestarting(Teuchos::TimeMonitor::getNewTimer("BlockDavidsonSolMgr restarting")),
00276   _timerLocking(Teuchos::TimeMonitor::getNewTimer("BlockDavidsonSolMgr locking"))
00277 {
00278   TEST_FOR_EXCEPTION(problem_ == Teuchos::null,              std::invalid_argument, "Problem not given to solver manager.");
00279   TEST_FOR_EXCEPTION(!problem_->isProblemSet(),              std::invalid_argument, "Problem not set.");
00280   TEST_FOR_EXCEPTION(!problem_->isHermitian(),               std::invalid_argument, "Problem not symmetric.");
00281   TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument, "Problem does not contain initial vectors to clone from.");
00282 
00283   std::string strtmp;
00284 
00285   // which values to solve for
00286   whch_ = pl.get("Which",whch_);
00287   TEST_FOR_EXCEPTION(whch_ != "SM" && whch_ != "LM" && whch_ != "SR" && whch_ != "LR",std::invalid_argument, "Invalid sorting string.");
00288 
00289   // which orthogonalization to use
00290   ortho_ = pl.get("Orthogonalization",ortho_);
00291   if (ortho_ != "DGKS" && ortho_ != "SVQB") {
00292     ortho_ = "SVQB";
00293   }
00294 
00295   // convergence tolerance
00296   convtol_ = pl.get("Convergence Tolerance",convtol_);
00297   relconvtol_ = pl.get("Relative Convergence Tolerance",relconvtol_);
00298   strtmp = pl.get("Convergence Norm",std::string("2"));
00299   if (strtmp == "2") {
00300     convNorm_ = StatusTestResNorm<ScalarType,MV,OP>::RES_2NORM;
00301   }
00302   else if (strtmp == "M") {
00303     convNorm_ = StatusTestResNorm<ScalarType,MV,OP>::RES_ORTH;
00304   }
00305   else {
00306     TEST_FOR_EXCEPTION(true, std::invalid_argument, 
00307         "Anasazi::BlockDavidsonSolMgr: Invalid Convergence Norm.");
00308   }
00309   
00310   // locking tolerance
00311   useLocking_ = pl.get("Use Locking",useLocking_);
00312   rellocktol_ = pl.get("Relative Locking Tolerance",rellocktol_);
00313   // default: should be less than convtol_
00314   locktol_ = convtol_/10;
00315   locktol_ = pl.get("Locking Tolerance",locktol_);
00316   strtmp = pl.get("Locking Norm",std::string("2"));
00317   if (strtmp == "2") {
00318     lockNorm_ = StatusTestResNorm<ScalarType,MV,OP>::RES_2NORM;
00319   }
00320   else if (strtmp == "M") {
00321     lockNorm_ = StatusTestResNorm<ScalarType,MV,OP>::RES_ORTH;
00322   }
00323   else {
00324     TEST_FOR_EXCEPTION(true, std::invalid_argument, 
00325         "Anasazi::BlockDavidsonSolMgr: Invalid Locking Norm.");
00326   }
00327 
00328   // maximum number of restarts
00329   maxRestarts_ = pl.get("Maximum Restarts",maxRestarts_);
00330 
00331   // block size: default is nev()
00332   blockSize_ = pl.get("Block Size",problem_->getNEV());
00333   TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
00334                      "Anasazi::BlockDavidsonSolMgr: \"Block Size\" must be strictly positive.");
00335   numBlocks_ = pl.get("Num Blocks",2);
00336   TEST_FOR_EXCEPTION(numBlocks_ <= 1, std::invalid_argument,
00337                      "Anasazi::BlockDavidsonSolMgr: \"Num Blocks\" must be >= 1.");
00338 
00339   // max locked: default is nev(), must satisfy maxLocked_ + blockSize_ >= nev
00340   if (useLocking_) {
00341     maxLocked_ = pl.get("Max Locked",problem_->getNEV());
00342   }
00343   else {
00344     maxLocked_ = 0;
00345   }
00346   if (maxLocked_ == 0) {
00347     useLocking_ = false;
00348   }
00349   TEST_FOR_EXCEPTION(maxLocked_ < 0, std::invalid_argument,
00350                      "Anasazi::BlockDavidsonSolMgr: \"Max Locked\" must be positive.");
00351   TEST_FOR_EXCEPTION(maxLocked_ + blockSize_ < problem_->getNEV(), 
00352                      std::invalid_argument,
00353                      "Anasazi::BlockDavidsonSolMgr: Not enough storage space for requested number of eigenpairs.");
00354   TEST_FOR_EXCEPTION(numBlocks_*blockSize_ + maxLocked_ > MVT::GetVecLength(*problem_->getInitVec()),
00355                      std::invalid_argument,
00356                      "Anasazi::BlockDavidsonSolMgr: Potentially impossible orthogonality requests. Reduce basis size or locking size.");
00357 
00358   if (useLocking_) {
00359     lockQuorum_ = pl.get("Locking Quorum",lockQuorum_);
00360     TEST_FOR_EXCEPTION(lockQuorum_ <= 0,
00361                        std::invalid_argument,
00362                        "Anasazi::BlockDavidsonSolMgr: \"Locking Quorum\" must be strictly positive.");
00363   }
00364 
00365   // restart size
00366   numRestartBlocks_ = pl.get("Num Restart Blocks",numRestartBlocks_);
00367   TEST_FOR_EXCEPTION(numRestartBlocks_ <= 0, std::invalid_argument,
00368                      "Anasazi::BlockDavidsonSolMgr: \"Num Restart Blocks\" must be strictly positive.");
00369   TEST_FOR_EXCEPTION(numRestartBlocks_ >= numBlocks_, std::invalid_argument,
00370                      "Anasazi::BlockDavidsonSolMgr: \"Num Restart Blocks\" must be strictly less than \"Num Blocks\".");
00371 
00372   // restarting technique: V*Q or applyHouse(V,H,tau)
00373   if (pl.isParameter("In Situ Restarting")) {
00374     if (Teuchos::isParameterType<bool>(pl,"In Situ Restarting")) {
00375       inSituRestart_ = pl.get("In Situ Restarting",inSituRestart_);
00376     } else {
00377       inSituRestart_ = ( Teuchos::getParameter<int>(pl,"In Situ Restarting") != 0 );
00378     }
00379   }
00380 
00381   // output stream
00382   std::string fntemplate = "";
00383   bool allProcs = false;
00384   if (pl.isParameter("Output on all processors")) {
00385     if (Teuchos::isParameterType<bool>(pl,"Output on all processors")) {
00386       allProcs = pl.get("Output on all processors",allProcs);
00387     } else {
00388       allProcs = ( Teuchos::getParameter<int>(pl,"Output on all processors") != 0 );
00389     }
00390   }
00391   fntemplate = pl.get("Output filename template",fntemplate);
00392   int MyPID;
00393 # ifdef HAVE_MPI
00394     // Initialize MPI
00395     int mpiStarted = 0;
00396     MPI_Initialized(&mpiStarted);
00397     if (mpiStarted) MPI_Comm_rank(MPI_COMM_WORLD, &MyPID);
00398     else MyPID=0;
00399 # else 
00400     MyPID = 0;
00401 # endif
00402   if (fntemplate != "") {
00403     std::ostringstream MyPIDstr;
00404     MyPIDstr << MyPID;
00405     // replace %d in fntemplate with MyPID
00406     int pos, start=0;
00407     while ( (pos = fntemplate.find("%d",start)) != -1 ) {
00408       fntemplate.replace(pos,2,MyPIDstr.str());
00409       start = pos+2;
00410     }
00411   }
00412   Teuchos::RCP<ostream> osp;
00413   if (fntemplate != "") {
00414     osp = Teuchos::rcp( new std::ofstream(fntemplate.c_str(),std::ios::out | std::ios::app) );
00415     if (!*osp) {
00416       osp = Teuchos::rcpFromRef(std::cout);
00417       std::cout << "Anasazi::BlockDavidsonSolMgr::constructor(): Could not open file for write: " << fntemplate << std::endl;
00418     }
00419   }
00420   else {
00421     osp = Teuchos::rcpFromRef(std::cout);
00422   }
00423   // Output manager
00424   int verbosity = Anasazi::Errors;
00425   if (pl.isParameter("Verbosity")) {
00426     if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
00427       verbosity = pl.get("Verbosity", verbosity);
00428     } else {
00429       verbosity = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
00430     }
00431   }
00432   if (allProcs) {
00433     // print on all procs
00434     printer_ = Teuchos::rcp( new BasicOutputManager<ScalarType>(verbosity,osp,MyPID) );
00435   }
00436   else {
00437     // print only on proc 0
00438     printer_ = Teuchos::rcp( new BasicOutputManager<ScalarType>(verbosity,osp,0) );
00439   }
00440 }
00441 
00442 
00443 // solve()
00444 template<class ScalarType, class MV, class OP>
00445 ReturnType 
00446 BlockDavidsonSolMgr<ScalarType,MV,OP>::solve() {
00447 
00448   typedef SolverUtils<ScalarType,MV,OP> msutils;
00449 
00450   const int nev = problem_->getNEV();
00451 
00452 #ifdef TEUCHOS_DEBUG
00453     Teuchos::RCP<Teuchos::FancyOStream>
00454       out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(printer_->stream(Debug)));
00455     out->setShowAllFrontMatter(false).setShowProcRank(true);
00456     *out << "Entering Anasazi::BlockDavidsonSolMgr::solve()\n";
00457 #endif
00458 
00460   // Sort manager
00461   Teuchos::RCP<BasicSort<MagnitudeType> > sorter = Teuchos::rcp( new BasicSort<MagnitudeType>(whch_) );
00462 
00464   // Status tests
00465   //
00466   // convergence
00467   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convtest;
00468   if (globalTest_ == Teuchos::null) {
00469     convtest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(convtol_,nev,convNorm_,relconvtol_) );
00470   }
00471   else {
00472     convtest = globalTest_;
00473   }
00474   Teuchos::RCP<StatusTestWithOrdering<ScalarType,MV,OP> > ordertest 
00475     = Teuchos::rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sorter,nev) );
00476   // locking
00477   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > locktest;
00478   if (useLocking_) {
00479     if (lockingTest_ == Teuchos::null) {
00480       locktest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(locktol_,lockQuorum_,lockNorm_,rellocktol_) );
00481     }
00482     else {
00483       locktest = lockingTest_;
00484     }
00485   }
00486   // for a non-short-circuited OR test, the order doesn't matter
00487   Teuchos::Array<Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
00488   alltests.push_back(ordertest);
00489   if (locktest != Teuchos::null) alltests.push_back(locktest);
00490   if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
00491 
00492   Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > combotest
00493     = Teuchos::rcp( new StatusTestCombo<ScalarType,MV,OP>( StatusTestCombo<ScalarType,MV,OP>::OR, alltests) );
00494   // printing StatusTest
00495   Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
00496   if ( printer_->isVerbosity(Debug) ) {
00497     outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,Passed+Failed+Undefined ) );
00498   }
00499   else {
00500     outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,Passed ) );
00501   }
00502 
00504   // Orthomanager
00505   Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho; 
00506   if (ortho_=="SVQB") {
00507     ortho = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
00508   } else if (ortho_=="DGKS") {
00509     ortho = Teuchos::rcp( new BasicOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
00510   } else {
00511     TEST_FOR_EXCEPTION(ortho_!="SVQB"&&ortho_!="DGKS",std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): Invalid orthogonalization type.");
00512   }
00513 
00515   // Parameter list
00516   Teuchos::ParameterList plist;
00517   plist.set("Block Size",blockSize_);
00518   plist.set("Num Blocks",numBlocks_);
00519 
00521   // BlockDavidson solver
00522   Teuchos::RCP<BlockDavidson<ScalarType,MV,OP> > bd_solver 
00523     = Teuchos::rcp( new BlockDavidson<ScalarType,MV,OP>(problem_,sorter,printer_,outputtest,ortho,plist) );
00524   // set any auxiliary vectors defined in the problem
00525   Teuchos::RCP< const MV > probauxvecs = problem_->getAuxVecs();
00526   if (probauxvecs != Teuchos::null) {
00527     bd_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
00528   }
00529 
00531   // Storage
00532   // 
00533   // lockvecs will contain eigenvectors that have been determined "locked" by the status test
00534   int curNumLocked = 0;
00535   Teuchos::RCP<MV> lockvecs;
00536   // lockvecs is used to hold the locked eigenvectors, as well as for temporary storage when locking.
00537   // when locking, we will lock some number of vectors numnew, where numnew <= maxlocked - curlocked
00538   // we will produce numnew random vectors, which will go into the space with the new basis.
00539   // we will also need numnew storage for the image of these random vectors under A and M; 
00540   // columns [curlocked+1,curlocked+numnew] will be used for this storage
00541   if (maxLocked_ > 0) {
00542     lockvecs = MVT::Clone(*problem_->getInitVec(),maxLocked_);
00543   }
00544   std::vector<MagnitudeType> lockvals;
00545   //
00546   // Restarting occurs under two scenarios: when the basis is full and after locking.
00547   //
00548   // For the former, a new basis of size blockSize*numRestartBlocks is generated using the current basis
00549   // and the most significant primitive Ritz vectors (projected eigenvectors).
00550   //     [S,L] = eig(KK)
00551   //     S = [Sr St]   // some for "r"estarting, some are "t"runcated
00552   //     newV = V*Sr
00553   //     KK_new = newV'*K*newV = Sr'*V'*K*V*Sr = Sr'*KK*Sr
00554   //  Therefore, the only multivector operation needed is for the generation of newV.
00555   //
00556   //  * If the multiplication is explicit, it requires a workspace of blockSize*numRestartBlocks vectors. 
00557   //    This space must be specifically allocated for that task, as we don't have any space of that size.
00558   //    It (workMV) will be allocated at the beginning of solve()
00559   //  * Optionally, the multiplication can be performed implicitly, via a Householder QR factorization of 
00560   //    Sr. This can be done in situ, using the basis multivector contained in the solver. This requires
00561   //    that we cast away the const on the multivector returned from getState(). Workspace for this approach
00562   //    is a single vector. the solver's internal storage must be preserved (X,MX,KX,R), requiring us to 
00563   //    allocate this vector.
00564   //
00565   // For the latter (restarting after locking), the new basis is the same size as existing basis. If numnew
00566   // vectors are locked, they are deflated from the current basis and replaced with randomly generated 
00567   // vectors.
00568   //     [S,L] = eig(KK)
00569   //     S = [Sl Su]  // partitioned: "l"ocked and "u"nlocked
00570   //     newL = V*Sl = X(locked)
00571   //     defV = V*Su
00572   //     augV = rand(numnew)  // orthogonal to oldL,newL,defV,auxvecs
00573   //     newV = [defV augV]
00574   //     Kknew = newV'*K*newV = [Su'*KK*Su    defV'*K*augV]
00575   //                            [augV'*K*defV augV'*K*augV]
00576   //     locked = [oldL newL]
00577   // Clearly, this operation is more complicated than the previous.
00578   // Here is a list of the significant computations that need to be performed:
00579   // - newL will be put into space in lockvecs, but will be copied from getState().X at the end
00580   // - defV,augV will be stored in workspace the size of the current basis.
00581   //   - If inSituRestart==true, we compute defV in situ in bd_solver::V_ and
00582   //     put augV at the end of bd_solver::V_
00583   //   - If inSituRestart==false, we must have curDim vectors available for 
00584   //     defV and augV; we will allocate a multivector (workMV) at the beginning of solve()
00585   //     for this purpose.
00586   // - M*augV and K*augV are needed; they will be stored in lockvecs. As a result, newL will
00587   //   not be put into lockvecs until the end.
00588   //
00589   // Therefore, we must allocate workMV when ((maxRestarts_ > 0) || (useLocking_ == true)) && inSituRestart == false
00590   // It will be allocated to size (numBlocks-1)*blockSize
00591   //
00592   Teuchos::RCP<MV> workMV;
00593   if (inSituRestart_ == false) {
00594     // we need storage space to restart, either if we may lock or if may restart after a full basis
00595     if (useLocking_==true || maxRestarts_ > 0) {
00596       workMV = MVT::Clone(*problem_->getInitVec(),(numBlocks_-1)*blockSize_);
00597     }
00598     else {
00599       // we will never need to restart.
00600       workMV = Teuchos::null;
00601     }
00602   }
00603   else { // inSituRestart_ == true 
00604     // we will restart in situ, if we need to restart
00605     // three situation remain: 
00606     // - never restart                                       => no space needed
00607     // - only restart for locking (i.e., never restart full) => no space needed
00608     // - restart for full basis                              => need one vector
00609     if (maxRestarts_ > 0) {
00610       workMV = MVT::Clone(*problem_->getInitVec(),1);
00611     }
00612     else {
00613       workMV = Teuchos::null;
00614     }
00615   }
00616 
00617   // some consts and utils
00618   const ScalarType ONE = SCT::one();
00619   const ScalarType ZERO = SCT::zero();
00620   Teuchos::LAPACK<int,ScalarType> lapack;
00621   Teuchos::BLAS<int,ScalarType> blas;
00622 
00623   // go ahead and initialize the solution to nothing in case we throw an exception
00624   Eigensolution<ScalarType,MV> sol;
00625   sol.numVecs = 0;
00626   problem_->setSolution(sol);
00627 
00628   int numRestarts = 0;
00629 
00630   // enter solve() iterations
00631   {
00632     Teuchos::TimeMonitor slvtimer(*_timerSolve);
00633 
00634     // tell bd_solver to iterate
00635     while (1) {
00636       try {
00637         bd_solver->iterate();
00638 
00640         //
00641         // check user-specified debug test; if it passed, return an exception
00642         //
00644         if (debugTest_ != Teuchos::null && debugTest_->getStatus() == Passed) {
00645           throw AnasaziError("Anasazi::BlockDavidsonSolMgr::solve(): User-specified debug status test returned Passed.");
00646         }
00648         //
00649         // check convergence next
00650         //
00652         else if (ordertest->getStatus() == Passed ) {
00653           // we have convergence
00654           // ordertest->whichVecs() tells us which vectors from lockvecs and solver state are the ones we want
00655           // ordertest->howMany() will tell us how many
00656           break;
00657         }
00659         //
00660         // check for restarting before locking: if we need to lock, it will happen after the restart
00661         //
00663         else if ( bd_solver->getCurSubspaceDim() == bd_solver->getMaxSubspaceDim() ) {
00664 
00665           Teuchos::TimeMonitor restimer(*_timerRestarting);
00666 
00667           if ( numRestarts >= maxRestarts_ ) {
00668             break; // break from while(1){bd_solver->iterate()}
00669           }
00670           numRestarts++;
00671 
00672           printer_->stream(IterationDetails) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << std::endl << std::endl;
00673 
00674           BlockDavidsonState<ScalarType,MV> state = bd_solver->getState();
00675           int curdim = state.curDim;
00676           int newdim = numRestartBlocks_*blockSize_;
00677 
00678 #       ifdef TEUCHOS_DEBUG
00679           {
00680             std::vector<Value<ScalarType> > ritzvalues = bd_solver->getRitzValues();
00681             *out << "Ritz values from solver:\n";
00682             for (unsigned int i=0; i<ritzvalues.size(); i++) {
00683               *out << ritzvalues[i].realpart << " ";
00684             }
00685             *out << "\n";
00686           }
00687 #       endif
00688 
00689           //
00690           // compute eigenvectors of the projected problem
00691           Teuchos::SerialDenseMatrix<int,ScalarType> S(curdim,curdim);
00692           std::vector<MagnitudeType> theta(curdim);
00693           int rank = curdim;
00694 #       ifdef TEUCHOS_DEBUG
00695           {
00696             std::stringstream os;
00697             os << "KK before HEEV...\n"
00698               << *state.KK << "\n";
00699             *out << os.str();
00700           }
00701 #       endif
00702           int info = msutils::directSolver(curdim,*state.KK,Teuchos::null,S,theta,rank,10);
00703           TEST_FOR_EXCEPTION(info != 0     ,std::logic_error,
00704               "Anasazi::BlockDavidsonSolMgr::solve(): error calling SolverUtils::directSolver.");       // this should never happen
00705           TEST_FOR_EXCEPTION(rank != curdim,std::logic_error,
00706               "Anasazi::BlockDavidsonSolMgr::solve(): direct solve did not compute all eigenvectors."); // this should never happen
00707 
00708 #       ifdef TEUCHOS_DEBUG
00709           {
00710             std::stringstream os;
00711             *out << "Ritz values from heev(KK):\n";
00712             for (unsigned int i=0; i<theta.size(); i++) *out << theta[i] << " ";
00713             os << "\nRitz vectors from heev(KK):\n" 
00714                                                  << S << "\n";
00715             *out << os.str();
00716           }
00717 #       endif
00718 
00719           // 
00720           // sort the eigenvalues (so that we can order the eigenvectors)
00721           {
00722             std::vector<int> order(curdim);
00723             sorter->sort(theta,Teuchos::rcpFromRef(order),curdim);
00724             //
00725             // apply the same ordering to the primitive ritz vectors
00726             msutils::permuteVectors(order,S);
00727           }
00728 #       ifdef TEUCHOS_DEBUG
00729           {
00730             std::stringstream os;
00731             *out << "Ritz values from heev(KK) after sorting:\n";
00732             std::copy(theta.begin(), theta.end(), std::ostream_iterator<ScalarType>(*out, " "));
00733             os << "\nRitz vectors from heev(KK) after sorting:\n" 
00734               << S << "\n";
00735             *out << os.str();
00736           }
00737 #       endif
00738           //
00739           // select the significant primitive ritz vectors
00740           Teuchos::SerialDenseMatrix<int,ScalarType> Sr(Teuchos::View,S,curdim,newdim);
00741 #       ifdef TEUCHOS_DEBUG
00742           {
00743             std::stringstream os;
00744             os << "Significant primitive Ritz vectors:\n"
00745               << Sr << "\n";
00746             *out << os.str();
00747           }
00748 #       endif
00749           //
00750           // generate newKK = Sr'*KKold*Sr
00751           Teuchos::SerialDenseMatrix<int,ScalarType> newKK(newdim,newdim);
00752           {
00753             Teuchos::SerialDenseMatrix<int,ScalarType> KKtmp(curdim,newdim), 
00754               KKold(Teuchos::View,*state.KK,curdim,curdim);
00755             int teuchosRet;
00756             // KKtmp = KKold*Sr
00757             teuchosRet = KKtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,KKold,Sr,ZERO);
00758             TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
00759                 "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
00760             // newKK = Sr'*KKtmp = Sr'*KKold*Sr
00761             teuchosRet = newKK.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,Sr,KKtmp,ZERO);
00762             TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
00763                 "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
00764             // make it Hermitian in memory
00765             for (int j=0; j<newdim-1; ++j) {
00766               for (int i=j+1; i<newdim; ++i) {
00767                 newKK(i,j) = SCT::conjugate(newKK(j,i));
00768               }
00769             }
00770           }
00771 #       ifdef TEUCHOS_DEBUG
00772           {
00773             std::stringstream os;
00774             os << "Sr'*KK*Sr:\n"
00775               << newKK << "\n";
00776             *out << os.str();
00777           }
00778 #       endif
00779 
00780           // prepare new state
00781           BlockDavidsonState<ScalarType,MV> rstate;
00782           rstate.curDim = newdim;
00783           rstate.KK = Teuchos::rcpFromRef(newKK);
00784           // 
00785           // we know that newX = newV*Sr(:,1:bS) = oldV*S(:1:bS) = oldX
00786           // the restarting preserves the Ritz vectors and residual
00787           // for the Ritz values, we want all of the values associated with newV. 
00788           // these have already been placed at the beginning of theta
00789           rstate.X  = state.X;
00790           rstate.KX = state.KX;
00791           rstate.MX = state.MX;
00792           rstate.R  = state.R;
00793           rstate.T  = Teuchos::rcp( new std::vector<MagnitudeType>(theta.begin(),theta.begin()+newdim) );
00794 
00795           if (inSituRestart_ == true) {
00796 #         ifdef TEUCHOS_DEBUG
00797             *out << "Beginning in-situ restart...\n";
00798 #         endif
00799             //
00800             // get non-const pointer to solver's basis so we can work in situ
00801             Teuchos::RCP<MV> solverbasis = Teuchos::rcp_const_cast<MV>(state.V);
00802             // 
00803             // perform Householder QR of Sr = Q [D;0], where D is unit diag.
00804             // WARNING: this will overwrite Sr; however, we do not need Sr anymore after this
00805             std::vector<ScalarType> tau(newdim), work(newdim);
00806             int geqrf_info;
00807             lapack.GEQRF(curdim,newdim,Sr.values(),Sr.stride(),&tau[0],&work[0],work.size(),&geqrf_info);
00808             TEST_FOR_EXCEPTION(geqrf_info != 0,std::logic_error,
00809                 "Anasazi::BlockDavidsonSolMgr::solve(): error calling GEQRF during restarting.");
00810             if (printer_->isVerbosity(Debug)) {
00811               Teuchos::SerialDenseMatrix<int,ScalarType> R(Teuchos::Copy,Sr,newdim,newdim);
00812               for (int j=0; j<newdim; j++) {
00813                 R(j,j) = SCT::magnitude(R(j,j)) - 1.0;
00814                 for (int i=j+1; i<newdim; i++) {
00815                   R(i,j) = ZERO;
00816                 }
00817               }
00818               printer_->stream(Debug) << "||Triangular factor of Sr - I||: " << R.normFrobenius() << std::endl;
00819             }
00820             // 
00821             // perform implicit oldV*Sr
00822             // this actually performs oldV*[Sr Su*M] = [newV truncV], for some unitary M
00823             // we are actually interested in only the first newdim vectors of the result
00824             {
00825               std::vector<int> curind(curdim);
00826               for (int i=0; i<curdim; i++) curind[i] = i;
00827               Teuchos::RCP<MV> oldV = MVT::CloneViewNonConst(*solverbasis,curind);
00828               msutils::applyHouse(newdim,*oldV,Sr,tau,workMV);
00829             }
00830             // 
00831             // put the new basis into the state for initialize()
00832             // the new basis is contained in the the first newdim columns of solverbasis
00833             // initialize() will recognize that pointer bd_solver.V_ == pointer rstate.V, and will neglect the copy.
00834             rstate.V = solverbasis;
00835           }
00836           else { // inSituRestart == false)
00837 #         ifdef TEUCHOS_DEBUG
00838             *out << "Beginning ex-situ restart...\n";
00839 #         endif
00840             // newV = oldV*Sr, explicitly. workspace is in workMV
00841             std::vector<int> curind(curdim), newind(newdim);
00842             for (int i=0; i<curdim; i++) curind[i] = i;
00843             for (int i=0; i<newdim; i++) newind[i] = i;
00844             Teuchos::RCP<const MV> oldV = MVT::CloneView(*state.V,curind);
00845             Teuchos::RCP<MV>       newV = MVT::CloneViewNonConst(*workMV ,newind);
00846 
00847             MVT::MvTimesMatAddMv(ONE,*oldV,Sr,ZERO,*newV);
00848             // 
00849             // put the new basis into the state for initialize()
00850             rstate.V = newV;
00851           }
00852 
00853           //
00854           // send the new state to the solver
00855           bd_solver->initialize(rstate);
00856         } // end of restarting
00858         //
00859         // check locking if we didn't converge or restart
00860         //
00862         else if (locktest != Teuchos::null && locktest->getStatus() == Passed) {
00863 
00864           Teuchos::TimeMonitor lcktimer(*_timerLocking);
00865 
00866           // 
00867           // get current state
00868           BlockDavidsonState<ScalarType,MV> state = bd_solver->getState();
00869           const int curdim = state.curDim;
00870 
00871           //
00872           // get number,indices of vectors to be locked
00873           TEST_FOR_EXCEPTION(locktest->howMany() <= 0,std::logic_error,
00874               "Anasazi::BlockDavidsonSolMgr::solve(): status test mistake: howMany() non-positive.");
00875           TEST_FOR_EXCEPTION(locktest->howMany() != (int)locktest->whichVecs().size(),std::logic_error,
00876               "Anasazi::BlockDavidsonSolMgr::solve(): status test mistake: howMany() not consistent with whichVecs().");
00877           // we should have room to lock vectors; otherwise, locking should have been deactivated
00878           TEST_FOR_EXCEPTION(curNumLocked == maxLocked_,std::logic_error,
00879               "Anasazi::BlockDavidsonSolMgr::solve(): status test mistake: locking not deactivated.");
00880           //
00881           // don't lock more than maxLocked_; we didn't allocate enough space.
00882           std::vector<int> tmp_vector_int;
00883           if (curNumLocked + locktest->howMany() > maxLocked_) {
00884             // just use the first of them
00885             tmp_vector_int.insert(tmp_vector_int.begin(),locktest->whichVecs().begin(),locktest->whichVecs().begin()+maxLocked_-curNumLocked);
00886           }
00887           else {
00888             tmp_vector_int = locktest->whichVecs();
00889           }
00890           const std::vector<int> lockind(tmp_vector_int);
00891           const int numNewLocked = lockind.size();
00892           //
00893           // generate indices of vectors left unlocked
00894           // curind = [0,...,curdim-1] = UNION( lockind, unlockind )
00895           const int numUnlocked = curdim-numNewLocked;
00896           tmp_vector_int.resize(curdim);
00897           for (int i=0; i<curdim; i++) tmp_vector_int[i] = i;
00898           const std::vector<int> curind(tmp_vector_int);       // curind = [0 ... curdim-1]
00899           tmp_vector_int.resize(numUnlocked); 
00900           std::set_difference(curind.begin(),curind.end(),lockind.begin(),lockind.end(),tmp_vector_int.begin());
00901           const std::vector<int> unlockind(tmp_vector_int);    // unlockind = [0 ... curdim-1] - lockind
00902           tmp_vector_int.clear();
00903 
00904           //
00905           // debug printing
00906           if (printer_->isVerbosity(Debug)) {
00907             printer_->print(Debug,"Locking vectors: ");
00908             for (unsigned int i=0; i<lockind.size(); i++) {printer_->stream(Debug) << " " << lockind[i];}
00909             printer_->print(Debug,"\n");
00910             printer_->print(Debug,"Not locking vectors: ");
00911             for (unsigned int i=0; i<unlockind.size(); i++) {printer_->stream(Debug) << " " << unlockind[i];}
00912             printer_->print(Debug,"\n");
00913           }
00914 
00915           //
00916           // we need primitive ritz vectors/values:
00917           // [S,L] = eig(oldKK)
00918           //
00919           // this will be partitioned as follows:
00920           //   locked: Sl = S(lockind)      // we won't actually need Sl
00921           // unlocked: Su = S(unlockind)
00922           //
00923           Teuchos::SerialDenseMatrix<int,ScalarType> S(curdim,curdim);
00924           std::vector<MagnitudeType> theta(curdim);
00925           {
00926             int rank = curdim;
00927             int info = msutils::directSolver(curdim,*state.KK,Teuchos::null,S,theta,rank,10);
00928             TEST_FOR_EXCEPTION(info != 0     ,std::logic_error,
00929                 "Anasazi::BlockDavidsonSolMgr::solve(): error calling SolverUtils::directSolver.");       // this should never happen
00930             TEST_FOR_EXCEPTION(rank != curdim,std::logic_error,
00931                 "Anasazi::BlockDavidsonSolMgr::solve(): direct solve did not compute all eigenvectors."); // this should never happen
00932             // 
00933             // sort the eigenvalues (so that we can order the eigenvectors)
00934             std::vector<int> order(curdim);
00935             sorter->sort(theta,Teuchos::rcpFromRef(order),curdim);
00936             //
00937             // apply the same ordering to the primitive ritz vectors
00938             msutils::permuteVectors(order,S);
00939           }
00940           //
00941           // select the unlocked ritz vectors
00942           // the indexing in unlockind is relative to the ordered primitive ritz vectors
00943           // (this is why we ordered theta,S above)
00944           Teuchos::SerialDenseMatrix<int,ScalarType> Su(curdim,numUnlocked);
00945           for (int i=0; i<numUnlocked; i++) {
00946             blas.COPY(curdim, S[unlockind[i]], 1, Su[i], 1);
00947           }
00948 
00949 
00950           //
00951           // newV has the following form:
00952           // newV = [defV augV]
00953           // - defV will be of size curdim - numNewLocked, and contain the generated basis: defV = oldV*Su
00954           // - augV will be of size numNewLocked, and contain random directions to make up for the lost space
00955           //
00956           // we will need a pointer to defV below to generate the off-diagonal block of newKK
00957           // go ahead and setup pointer to augV
00958           //
00959           Teuchos::RCP<MV> defV, augV;
00960           if (inSituRestart_ == true) {
00961             //
00962             // get non-const pointer to solver's basis so we can work in situ
00963             Teuchos::RCP<MV> solverbasis = Teuchos::rcp_const_cast<MV>(state.V);
00964             // 
00965             // perform Householder QR of Su = Q [D;0], where D is unit diag.
00966             // work on a copy of Su, since we need Su below to build newKK
00967             Teuchos::SerialDenseMatrix<int,ScalarType> copySu(Su);
00968             std::vector<ScalarType> tau(numUnlocked), work(numUnlocked);
00969             int info;
00970             lapack.GEQRF(curdim,numUnlocked,copySu.values(),copySu.stride(),&tau[0],&work[0],work.size(),&info);
00971             TEST_FOR_EXCEPTION(info != 0,std::logic_error,
00972                 "Anasazi::BlockDavidsonSolMgr::solve(): error calling GEQRF during restarting.");
00973             if (printer_->isVerbosity(Debug)) {
00974               Teuchos::SerialDenseMatrix<int,ScalarType> R(Teuchos::Copy,copySu,numUnlocked,numUnlocked);
00975               for (int j=0; j<numUnlocked; j++) {
00976                 R(j,j) = SCT::magnitude(R(j,j)) - 1.0;
00977                 for (int i=j+1; i<numUnlocked; i++) {
00978                   R(i,j) = ZERO;
00979                 }
00980               }
00981               printer_->stream(Debug) << "||Triangular factor of Su - I||: " << R.normFrobenius() << std::endl;
00982             }
00983             // 
00984             // perform implicit oldV*Su
00985             // this actually performs oldV*[Su Sl*M] = [defV lockV], for some unitary M
00986             // we are actually interested in only the first numUnlocked vectors of the result
00987             {
00988               Teuchos::RCP<MV> oldV = MVT::CloneViewNonConst(*solverbasis,curind);
00989               msutils::applyHouse(numUnlocked,*oldV,copySu,tau,workMV);
00990             }
00991             std::vector<int> defind(numUnlocked), augind(numNewLocked);
00992             for (int i=0; i<numUnlocked ; i++) defind[i] = i;
00993             for (int i=0; i<numNewLocked; i++) augind[i] = numUnlocked+i;
00994             defV = MVT::CloneViewNonConst(*solverbasis,defind);
00995             augV = MVT::CloneViewNonConst(*solverbasis,augind);
00996           }
00997           else { // inSituRestart == false)
00998             // defV = oldV*Su, explicitly. workspace is in workMV
00999             std::vector<int> defind(numUnlocked), augind(numNewLocked);
01000             for (int i=0; i<numUnlocked ; i++) defind[i] = i;
01001             for (int i=0; i<numNewLocked; i++) augind[i] = numUnlocked+i;
01002             Teuchos::RCP<const MV> oldV = MVT::CloneView(*state.V,curind);
01003             defV = MVT::CloneViewNonConst(*workMV,defind);
01004             augV = MVT::CloneViewNonConst(*workMV,augind);
01005 
01006             MVT::MvTimesMatAddMv(ONE,*oldV,Su,ZERO,*defV);
01007           }
01008 
01009           //
01010           // lockvecs will be partitioned as follows:
01011           // lockvecs = [curlocked augTmp ...]
01012           // - augTmp will be used for the storage of M*augV and K*augV
01013           //   later, the locked vectors (stored in state.X and referenced via const MV view newLocked) 
01014           //   will be moved into lockvecs on top of augTmp when it is no longer needed as workspace.
01015           // - curlocked will be used in orthogonalization of augV
01016           //
01017           // newL is the new locked vectors; newL = oldV*Sl = RitzVectors(lockind)
01018           // we will not produce them, but instead retrieve them from RitzVectors
01019           //
01020           Teuchos::RCP<const MV> curlocked, newLocked;
01021           Teuchos::RCP<MV> augTmp;
01022           {
01023             // setup curlocked
01024             if (curNumLocked > 0) {
01025               std::vector<int> curlockind(curNumLocked);
01026               for (int i=0; i<curNumLocked; i++) curlockind[i] = i;
01027               curlocked = MVT::CloneView(*lockvecs,curlockind);
01028             }
01029             else {
01030               curlocked = Teuchos::null;
01031             }
01032             // setup augTmp
01033             std::vector<int> augtmpind(numNewLocked); 
01034             for (int i=0; i<numNewLocked; i++) augtmpind[i] = curNumLocked+i;
01035             augTmp = MVT::CloneViewNonConst(*lockvecs,augtmpind);
01036             // setup newLocked
01037             newLocked = MVT::CloneView(*bd_solver->getRitzVectors(),lockind);
01038           }
01039 
01040           // 
01041           // generate augV and perform orthogonalization
01042           //
01043           MVT::MvRandom(*augV);
01044           // 
01045           // orthogonalize it against auxvecs, defV, and all locked vectors (new and current)
01046           // use augTmp as storage for M*augV, if hasM
01047           {
01048             Teuchos::Array<Teuchos::RCP<const MV> > against;
01049             Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
01050             if (probauxvecs != Teuchos::null) against.push_back(probauxvecs);
01051             if (curlocked != Teuchos::null)   against.push_back(curlocked);
01052             against.push_back(newLocked);
01053             against.push_back(defV);
01054             if (problem_->getM() != Teuchos::null) {
01055               OPT::Apply(*problem_->getM(),*augV,*augTmp);
01056             }
01057             ortho->projectAndNormalizeMat(*augV,against,dummyC,Teuchos::null,augTmp);
01058           }
01059 
01060           //
01061           // form newKK
01062           //
01063           // newKK = newV'*K*newV = [Su'*KK*Su    defV'*K*augV]
01064           //                        [augV'*K*defV augV'*K*augV]
01065           //
01066           // first, generate the principal submatrix, the projection of K onto the unlocked portion of oldV
01067           //
01068           Teuchos::SerialDenseMatrix<int,ScalarType> newKK(curdim,curdim);
01069           {
01070             Teuchos::SerialDenseMatrix<int,ScalarType> KKtmp(curdim,numUnlocked), 
01071               KKold(Teuchos::View,*state.KK,curdim,curdim),
01072               KK11(Teuchos::View,newKK,numUnlocked,numUnlocked);
01073             int teuchosRet;
01074             // KKtmp = KKold*Su
01075             teuchosRet = KKtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,KKold,Su,ZERO);
01076             TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
01077                 "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
01078             // KK11 = Su'*KKtmp = Su'*KKold*Su
01079             teuchosRet = KK11.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,Su,KKtmp,ZERO);
01080             TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
01081                 "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
01082           }
01083           //
01084           // project the stiffness matrix on augV
01085           {
01086             OPT::Apply(*problem_->getOperator(),*augV,*augTmp);
01087             Teuchos::SerialDenseMatrix<int,ScalarType> KK12(Teuchos::View,newKK,numUnlocked,numNewLocked,0,numUnlocked),
01088               KK22(Teuchos::View,newKK,numNewLocked,numNewLocked,numUnlocked,numUnlocked);
01089             MVT::MvTransMv(ONE,*defV,*augTmp,KK12);
01090             MVT::MvTransMv(ONE,*augV,*augTmp,KK22);
01091           }
01092           // 
01093           // done with defV,augV
01094           defV = Teuchos::null;
01095           augV = Teuchos::null;
01096           //
01097           // make it hermitian in memory (fill in KK21)
01098           for (int j=0; j<curdim; ++j) {
01099             for (int i=j+1; i<curdim; ++i) {
01100               newKK(i,j) = SCT::conjugate(newKK(j,i));
01101             }
01102           }
01103           //
01104           // we are done using augTmp as storage
01105           // put newLocked into lockvecs, new values into lockvals
01106           augTmp = Teuchos::null;
01107           {
01108             std::vector<Value<ScalarType> > allvals = bd_solver->getRitzValues();
01109             for (int i=0; i<numNewLocked; i++) {
01110               lockvals.push_back(allvals[lockind[i]].realpart);
01111             }
01112 
01113             std::vector<int> indlock(numNewLocked);
01114             for (int i=0; i<numNewLocked; i++) indlock[i] = curNumLocked+i;
01115             MVT::SetBlock(*newLocked,indlock,*lockvecs);
01116             newLocked = Teuchos::null;
01117 
01118             curNumLocked += numNewLocked;
01119             std::vector<int> curlockind(curNumLocked);
01120             for (int i=0; i<curNumLocked; i++) curlockind[i] = i;
01121             curlocked = MVT::CloneView(*lockvecs,curlockind);
01122           }
01123           // add locked vecs as aux vecs, along with aux vecs from problem
01124           // add lockvals to ordertest
01125           // disable locktest if curNumLocked == maxLocked
01126           {
01127             ordertest->setAuxVals(lockvals);
01128 
01129             Teuchos::Array< Teuchos::RCP<const MV> > aux;
01130             if (probauxvecs != Teuchos::null) aux.push_back(probauxvecs);
01131             aux.push_back(curlocked);
01132             bd_solver->setAuxVecs(aux);
01133 
01134             if (curNumLocked == maxLocked_) {
01135               // disabled locking now
01136               combotest->removeTest(locktest);
01137             }
01138           }
01139 
01140           //
01141           // prepare new state
01142           BlockDavidsonState<ScalarType,MV> rstate;
01143           rstate.curDim = curdim;
01144           if (inSituRestart_) {
01145             // data is already in the solver's memory
01146             rstate.V = state.V;
01147           }
01148           else {
01149             // data is in workspace and will be copied to solver memory
01150             rstate.V = workMV;
01151           }
01152           rstate.KK = Teuchos::rcpFromRef(newKK);
01153           //
01154           // pass new state to the solver
01155           bd_solver->initialize(rstate);
01156         } // end of locking
01158         //
01159         // we returned from iterate(), but none of our status tests Passed.
01160         // something is wrong, and it is probably our fault.
01161         //
01163         else {
01164           TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): Invalid return from bd_solver::iterate().");
01165         }
01166       }
01167       catch (const AnasaziError &err) {
01168         printer_->stream(Errors) 
01169           << "Anasazi::BlockDavidsonSolMgr::solve() caught unexpected exception from Anasazi::BlockDavidson::iterate() at iteration " << bd_solver->getNumIters() << std::endl
01170           << err.what() << std::endl
01171           << "Anasazi::BlockDavidsonSolMgr::solve() returning Unconverged with no solutions." << std::endl;
01172         return Unconverged;
01173       }
01174     }
01175 
01176     // clear temp space
01177     workMV = Teuchos::null;
01178 
01179     sol.numVecs = ordertest->howMany();
01180     if (sol.numVecs > 0) {
01181       sol.Evecs = MVT::Clone(*problem_->getInitVec(),sol.numVecs);
01182       sol.Espace = sol.Evecs;
01183       sol.Evals.resize(sol.numVecs);
01184       std::vector<MagnitudeType> vals(sol.numVecs);
01185 
01186       // copy them into the solution
01187       std::vector<int> which = ordertest->whichVecs();
01188       // indices between [0,blockSize) refer to vectors/values in the solver
01189       // indices between [-curNumLocked,-1] refer to locked vectors/values [0,curNumLocked)
01190       // everything has already been ordered by the solver; we just have to partition the two references
01191       std::vector<int> inlocked(0), insolver(0);
01192       for (unsigned int i=0; i<which.size(); i++) {
01193         if (which[i] >= 0) {
01194           TEST_FOR_EXCEPTION(which[i] >= blockSize_,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): positive indexing mistake from ordertest.");
01195           insolver.push_back(which[i]);
01196         }
01197         else {
01198           // sanity check
01199           TEST_FOR_EXCEPTION(which[i] < -curNumLocked,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): negative indexing mistake from ordertest.");
01200           inlocked.push_back(which[i] + curNumLocked);
01201         }
01202       }
01203 
01204       TEST_FOR_EXCEPTION(insolver.size() + inlocked.size() != (unsigned int)sol.numVecs,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): indexing mistake.");
01205 
01206       // set the vecs,vals in the solution
01207       if (insolver.size() > 0) {
01208         // set vecs
01209         int lclnum = insolver.size();
01210         std::vector<int> tosol(lclnum);
01211         for (int i=0; i<lclnum; i++) tosol[i] = i;
01212         Teuchos::RCP<const MV> v = MVT::CloneView(*bd_solver->getRitzVectors(),insolver);
01213         MVT::SetBlock(*v,tosol,*sol.Evecs);
01214         // set vals
01215         std::vector<Value<ScalarType> > fromsolver = bd_solver->getRitzValues();
01216         for (unsigned int i=0; i<insolver.size(); i++) {
01217           vals[i] = fromsolver[insolver[i]].realpart;
01218         }
01219       }
01220 
01221       // get the vecs,vals from locked storage
01222       if (inlocked.size() > 0) {
01223         int solnum = insolver.size();
01224         // set vecs
01225         int lclnum = inlocked.size();
01226         std::vector<int> tosol(lclnum);
01227         for (int i=0; i<lclnum; i++) tosol[i] = solnum + i;
01228         Teuchos::RCP<const MV> v = MVT::CloneView(*lockvecs,inlocked);
01229         MVT::SetBlock(*v,tosol,*sol.Evecs);
01230         // set vals
01231         for (unsigned int i=0; i<inlocked.size(); i++) {
01232           vals[i+solnum] = lockvals[inlocked[i]];
01233         }
01234       }
01235 
01236       // sort the eigenvalues and permute the eigenvectors appropriately
01237       {
01238         std::vector<int> order(sol.numVecs);
01239         sorter->sort(vals,Teuchos::rcpFromRef(order),sol.numVecs);
01240         // store the values in the Eigensolution
01241         for (int i=0; i<sol.numVecs; i++) {
01242           sol.Evals[i].realpart = vals[i];
01243           sol.Evals[i].imagpart = MT::zero();
01244         }
01245         // now permute the eigenvectors according to order
01246         msutils::permuteVectors(sol.numVecs,order,*sol.Evecs);
01247       }
01248 
01249       // setup sol.index, remembering that all eigenvalues are real so that index = {0,...,0}
01250       sol.index.resize(sol.numVecs,0);
01251     }
01252   }
01253 
01254   // print final summary
01255   bd_solver->currentStatus(printer_->stream(FinalSummary));
01256 
01257   // print timing information
01258   Teuchos::TimeMonitor::summarize(printer_->stream(TimingDetails));
01259 
01260   problem_->setSolution(sol);
01261   printer_->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
01262 
01263   // get the number of iterations taken for this call to solve().
01264   numIters_ = bd_solver->getNumIters();
01265 
01266   if (sol.numVecs < nev) {
01267     return Unconverged; // return from BlockDavidsonSolMgr::solve() 
01268   }
01269   return Converged; // return from BlockDavidsonSolMgr::solve() 
01270 }
01271 
01272 
01273 template <class ScalarType, class MV, class OP>
01274 void 
01275 BlockDavidsonSolMgr<ScalarType,MV,OP>::setGlobalStatusTest(
01276     const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global) 
01277 {
01278   globalTest_ = global;
01279 }
01280 
01281 template <class ScalarType, class MV, class OP>
01282 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & 
01283 BlockDavidsonSolMgr<ScalarType,MV,OP>::getGlobalStatusTest() const 
01284 {
01285   return globalTest_;
01286 }
01287 
01288 template <class ScalarType, class MV, class OP>
01289 void 
01290 BlockDavidsonSolMgr<ScalarType,MV,OP>::setDebugStatusTest(
01291     const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &debug)
01292 {
01293   debugTest_ = debug;
01294 }
01295 
01296 template <class ScalarType, class MV, class OP>
01297 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & 
01298 BlockDavidsonSolMgr<ScalarType,MV,OP>::getDebugStatusTest() const
01299 {
01300   return debugTest_;
01301 }
01302 
01303 template <class ScalarType, class MV, class OP>
01304 void 
01305 BlockDavidsonSolMgr<ScalarType,MV,OP>::setLockingStatusTest(
01306     const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &locking) 
01307 {
01308   lockingTest_ = locking;
01309 }
01310 
01311 template <class ScalarType, class MV, class OP>
01312 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & 
01313 BlockDavidsonSolMgr<ScalarType,MV,OP>::getLockingStatusTest() const 
01314 {
01315   return lockingTest_;
01316 }
01317 
01318 } // end Anasazi namespace
01319 
01320 #endif /* ANASAZI_BLOCKDAVIDSON_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends
Generated on Wed Apr 13 09:56:58 2011 for Anasazi by  doxygen 1.6.3