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