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