Anasazi Version of the Day
AnasaziBlockKrylovSchurSolMgr.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_BLOCK_KRYLOV_SCHUR_SOLMGR_HPP
00030 #define ANASAZI_BLOCK_KRYLOV_SCHUR_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 "AnasaziBlockKrylovSchur.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 
00056 
00061 
00062 
00063 
00064 
00065 
00066 
00067 
00068 
00069 
00070 
00071 
00072 
00073 
00074 
00075 
00076 
00077 
00078 
00079 
00080 
00081 
00082 
00083 
00084 
00085 
00086 
00103 namespace Anasazi {
00104 
00105 
00132 template<class ScalarType, class MV, class OP>
00133 class BlockKrylovSchurSolMgr : public SolverManager<ScalarType,MV,OP> {
00134 
00135   private:
00136     typedef MultiVecTraits<ScalarType,MV> MVT;
00137     typedef MultiVecTraitsExt<ScalarType,MV> MVText;
00138     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00139     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00140     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00141     typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00142 
00143   public:
00144 
00146 
00147 
00165   BlockKrylovSchurSolMgr( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00166                              Teuchos::ParameterList &pl );
00167 
00169   virtual ~BlockKrylovSchurSolMgr() {};
00171 
00173 
00174 
00176   const Eigenproblem<ScalarType,MV,OP>& getProblem() const {
00177     return *_problem;
00178   }
00179 
00181   int getNumIters() const {
00182     return _numIters;
00183   }
00184 
00187   std::vector<Value<ScalarType> > getRitzValues() const {
00188     std::vector<Value<ScalarType> > ret( _ritzValues );
00189     return ret;
00190   }
00191 
00198    Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00199      return Teuchos::tuple(_timerSolve, _timerRestarting);
00200    }
00201 
00203 
00205 
00206 
00225   ReturnType solve();
00226 
00228   void setGlobalStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global);
00229 
00231   const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getGlobalStatusTest() const;
00232 
00234   void setDebugStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &debug);
00235 
00237   const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getDebugStatusTest() const;
00238 
00240 
00241   private:
00242   Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > _problem;
00243   Teuchos::RCP<SortManager<MagnitudeType> > _sort;
00244 
00245   std::string _whch, _ortho;
00246   MagnitudeType _ortho_kappa;
00247 
00248   MagnitudeType _convtol;
00249   int _maxRestarts;
00250   bool _relconvtol,_conjSplit;
00251   int _blockSize, _numBlocks, _stepSize, _nevBlocks, _xtra_nevBlocks;
00252   int _numIters;
00253   int _verbosity;
00254   bool _inSituRestart, _dynXtraNev;
00255 
00256   std::vector<Value<ScalarType> > _ritzValues;
00257 
00258   int _printNum;
00259   Teuchos::RCP<Teuchos::Time> _timerSolve, _timerRestarting;
00260 
00261   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > globalTest_;
00262   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > debugTest_;
00263 
00264 };
00265 
00266 
00267 // Constructor
00268 template<class ScalarType, class MV, class OP>
00269 BlockKrylovSchurSolMgr<ScalarType,MV,OP>::BlockKrylovSchurSolMgr(
00270         const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00271         Teuchos::ParameterList &pl ) :
00272   _problem(problem),
00273   _whch("LM"),
00274   _ortho("SVQB"),
00275   _ortho_kappa(-1.0),
00276   _convtol(0),
00277   _maxRestarts(20),
00278   _relconvtol(true),
00279   _conjSplit(false),
00280   _blockSize(0),
00281   _numBlocks(0),
00282   _stepSize(0),
00283   _nevBlocks(0),
00284   _xtra_nevBlocks(0),
00285   _numIters(0),
00286   _verbosity(Anasazi::Errors),
00287   _inSituRestart(false),
00288   _dynXtraNev(false),
00289   _printNum(-1)
00290 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
00291   ,_timerSolve(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchurSolMgr::solve()")),
00292   _timerRestarting(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchurSolMgr restarting"))
00293 #endif
00294 {
00295   TEUCHOS_TEST_FOR_EXCEPTION(_problem == Teuchos::null,               std::invalid_argument, "Problem not given to solver manager.");
00296   TEUCHOS_TEST_FOR_EXCEPTION(!_problem->isProblemSet(),               std::invalid_argument, "Problem not set.");
00297   TEUCHOS_TEST_FOR_EXCEPTION(_problem->getInitVec() == Teuchos::null, std::invalid_argument, "Problem does not contain initial vectors to clone from.");
00298 
00299   const int nev = _problem->getNEV();
00300 
00301   // convergence tolerance
00302   _convtol = pl.get("Convergence Tolerance",MT::prec());
00303   _relconvtol = pl.get("Relative Convergence Tolerance",_relconvtol);
00304 
00305   // maximum number of restarts
00306   _maxRestarts = pl.get("Maximum Restarts",_maxRestarts);
00307 
00308   // block size: default is 1
00309   _blockSize = pl.get("Block Size",1);
00310   TEUCHOS_TEST_FOR_EXCEPTION(_blockSize <= 0, std::invalid_argument,
00311                      "Anasazi::BlockKrylovSchurSolMgr: \"Block Size\" must be strictly positive.");
00312 
00313   // set the number of blocks we need to save to compute the nev eigenvalues of interest.
00314   _xtra_nevBlocks = pl.get("Extra NEV Blocks",0);
00315   if (nev%_blockSize) {
00316     _nevBlocks = nev/_blockSize + 1;
00317   } else {
00318     _nevBlocks = nev/_blockSize;
00319   }
00320 
00321   // determine if we should use the dynamic scheme for selecting the current number of retained eigenvalues.
00322   // NOTE:  This employs a technique similar to ARPACK in that it increases the number of retained eigenvalues
00323   //        by one for every converged eigenpair to accelerate convergence.
00324   if (pl.isParameter("Dynamic Extra NEV")) {
00325     if (Teuchos::isParameterType<bool>(pl,"Dynamic Extra NEV")) {
00326       _dynXtraNev = pl.get("Dynamic Extra NEV",_dynXtraNev);
00327     } else {
00328       _dynXtraNev = ( Teuchos::getParameter<int>(pl,"Dynamic Extra NEV") != 0 );
00329     }
00330   }
00331 
00332   _numBlocks = pl.get("Num Blocks",3*_nevBlocks);
00333   TEUCHOS_TEST_FOR_EXCEPTION(_numBlocks <= _nevBlocks, std::invalid_argument,
00334                      "Anasazi::BlockKrylovSchurSolMgr: \"Num Blocks\" must be strictly positive and large enough to compute the requested eigenvalues.");
00335 
00336   TEUCHOS_TEST_FOR_EXCEPTION(static_cast<ptrdiff_t>(_numBlocks)*_blockSize > MVText::GetGlobalLength(*_problem->getInitVec()),
00337                      std::invalid_argument,
00338                      "Anasazi::BlockKrylovSchurSolMgr: Potentially impossible orthogonality requests. Reduce basis size.");
00339 
00340   // step size: the default is _maxRestarts*_numBlocks, so that Ritz values are only computed every restart.
00341   if (_maxRestarts) {
00342     _stepSize = pl.get("Step Size", (_maxRestarts+1)*(_numBlocks+1));
00343   } else {
00344     _stepSize = pl.get("Step Size", _numBlocks+1);
00345   }
00346   TEUCHOS_TEST_FOR_EXCEPTION(_stepSize < 1, std::invalid_argument,
00347                      "Anasazi::BlockKrylovSchurSolMgr: \"Step Size\" must be strictly positive.");
00348 
00349   // get the sort manager
00350   if (pl.isParameter("Sort Manager")) {
00351     _sort = Teuchos::getParameter<Teuchos::RCP<Anasazi::SortManager<MagnitudeType> > >(pl,"Sort Manager");
00352   } else {
00353     // which values to solve for
00354     _whch = pl.get("Which",_whch);
00355     TEUCHOS_TEST_FOR_EXCEPTION(_whch != "SM" && _whch != "LM" && _whch != "SR" && _whch != "LR" && _whch != "SI" && _whch != "LI",
00356                        std::invalid_argument, "Invalid sorting string.");
00357     _sort = Teuchos::rcp( new BasicSort<MagnitudeType>(_whch) );
00358   }
00359 
00360   // which orthogonalization to use
00361   _ortho = pl.get("Orthogonalization",_ortho);
00362   if (_ortho != "DGKS" && _ortho != "SVQB") {
00363     _ortho = "SVQB";
00364   }
00365 
00366   // which orthogonalization constant to use
00367   _ortho_kappa = pl.get("Orthogonalization Constant",_ortho_kappa);
00368 
00369   // verbosity level
00370   if (pl.isParameter("Verbosity")) {
00371     if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
00372       _verbosity = pl.get("Verbosity", _verbosity);
00373     } else {
00374       _verbosity = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
00375     }
00376   }
00377 
00378   // restarting technique: V*Q or applyHouse(V,H,tau)
00379   if (pl.isParameter("In Situ Restarting")) {
00380     if (Teuchos::isParameterType<bool>(pl,"In Situ Restarting")) {
00381       _inSituRestart = pl.get("In Situ Restarting",_inSituRestart);
00382     } else {
00383       _inSituRestart = ( Teuchos::getParameter<int>(pl,"In Situ Restarting") != 0 );
00384     }
00385   }
00386 
00387   _printNum = pl.get<int>("Print Number of Ritz Values",-1);
00388 }
00389 
00390 
00391 // solve()
00392 template<class ScalarType, class MV, class OP>
00393 ReturnType
00394 BlockKrylovSchurSolMgr<ScalarType,MV,OP>::solve() {
00395 
00396   const int nev = _problem->getNEV();
00397   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00398   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00399 
00400   Teuchos::BLAS<int,ScalarType> blas;
00401   Teuchos::LAPACK<int,ScalarType> lapack;
00402   typedef SolverUtils<ScalarType,MV,OP> msutils;
00403 
00405   // Output manager
00406   Teuchos::RCP<BasicOutputManager<ScalarType> > printer = Teuchos::rcp( new BasicOutputManager<ScalarType>(_verbosity) );
00407 
00409   // Status tests
00410   //
00411   // convergence
00412   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convtest;
00413   if (globalTest_ == Teuchos::null) {
00414     convtest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(_convtol,nev,StatusTestResNorm<ScalarType,MV,OP>::RITZRES_2NORM,_relconvtol) );
00415   }
00416   else {
00417     convtest = globalTest_;
00418   }
00419   Teuchos::RCP<StatusTestWithOrdering<ScalarType,MV,OP> > ordertest
00420     = Teuchos::rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,_sort,nev) );
00421   // for a non-short-circuited OR test, the order doesn't matter
00422   Teuchos::Array<Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
00423   alltests.push_back(ordertest);
00424 
00425   if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
00426 
00427   Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > combotest
00428     = Teuchos::rcp( new StatusTestCombo<ScalarType,MV,OP>( StatusTestCombo<ScalarType,MV,OP>::OR, alltests) );
00429   // printing StatusTest
00430   Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
00431   if ( printer->isVerbosity(Debug) ) {
00432     outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,Passed+Failed+Undefined ) );
00433   }
00434   else {
00435     outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,Passed ) );
00436   }
00437 
00439   // Orthomanager
00440   Teuchos::RCP<OrthoManager<ScalarType,MV> > ortho;
00441   if (_ortho=="SVQB") {
00442     ortho = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>(_problem->getM()) );
00443   } else if (_ortho=="DGKS") {
00444     if (_ortho_kappa <= 0) {
00445       ortho = Teuchos::rcp( new BasicOrthoManager<ScalarType,MV,OP>(_problem->getM()) );
00446     }
00447     else {
00448       ortho = Teuchos::rcp( new BasicOrthoManager<ScalarType,MV,OP>(_problem->getM(),_ortho_kappa) );
00449     }
00450   } else {
00451     TEUCHOS_TEST_FOR_EXCEPTION(_ortho!="SVQB"&&_ortho!="DGKS",std::logic_error,"Anasazi::BlockKrylovSchurSolMgr::solve(): Invalid orthogonalization type.");
00452   }
00453 
00455   // Parameter list
00456   Teuchos::ParameterList plist;
00457   plist.set("Block Size",_blockSize);
00458   plist.set("Num Blocks",_numBlocks);
00459   plist.set("Step Size",_stepSize);
00460   if (_printNum == -1) {
00461     plist.set("Print Number of Ritz Values",_nevBlocks*_blockSize);
00462   }
00463   else {
00464     plist.set("Print Number of Ritz Values",_printNum);
00465   }
00466 
00468   // BlockKrylovSchur solver
00469   Teuchos::RCP<BlockKrylovSchur<ScalarType,MV,OP> > bks_solver
00470     = Teuchos::rcp( new BlockKrylovSchur<ScalarType,MV,OP>(_problem,_sort,printer,outputtest,ortho,plist) );
00471   // set any auxiliary vectors defined in the problem
00472   Teuchos::RCP< const MV > probauxvecs = _problem->getAuxVecs();
00473   if (probauxvecs != Teuchos::null) {
00474     bks_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
00475   }
00476 
00477   // Create workspace for the Krylov basis generated during a restart
00478   // Need at most (_nevBlocks*_blockSize+1) for the updated factorization and another block for the current factorization residual block (F).
00479   //  ---> (_nevBlocks*_blockSize+1) + _blockSize
00480   // If Hermitian, this becomes _nevBlocks*_blockSize + _blockSize
00481   // we only need this if there is the possibility of restarting, ex situ
00482 
00483   // Maximum allowable extra vectors that BKS can keep to accelerate convergence
00484   int maxXtraBlocks = 0;
00485   if ( _dynXtraNev ) maxXtraBlocks = ( ( bks_solver->getMaxSubspaceDim() - nev ) / _blockSize ) / 2;
00486 
00487   Teuchos::RCP<MV> workMV;
00488   if (_maxRestarts > 0) {
00489     if (_inSituRestart==true) {
00490       // still need one work vector for applyHouse()
00491       workMV = MVT::Clone( *_problem->getInitVec(), 1 );
00492     }
00493     else { // inSituRestart == false
00494       if (_problem->isHermitian()) {
00495         workMV = MVT::Clone( *_problem->getInitVec(), (_nevBlocks+maxXtraBlocks)*_blockSize + _blockSize );
00496       } else {
00497         // Increase workspace by 1 because of potential complex conjugate pairs on the Ritz vector boundary
00498         workMV = MVT::Clone( *_problem->getInitVec(), (_nevBlocks+maxXtraBlocks)*_blockSize + 1 + _blockSize );
00499       }
00500     }
00501   } else {
00502     workMV = Teuchos::null;
00503   }
00504 
00505   // go ahead and initialize the solution to nothing in case we throw an exception
00506   Eigensolution<ScalarType,MV> sol;
00507   sol.numVecs = 0;
00508   _problem->setSolution(sol);
00509 
00510   int numRestarts = 0;
00511   int cur_nevBlocks = 0;
00512 
00513   // enter solve() iterations
00514   {
00515 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
00516     Teuchos::TimeMonitor slvtimer(*_timerSolve);
00517 #endif
00518 
00519     // tell bks_solver to iterate
00520     while (1) {
00521       try {
00522         bks_solver->iterate();
00523 
00525         //
00526         // check convergence first
00527         //
00529         if ( ordertest->getStatus() == Passed ) {
00530           // we have convergence
00531           // ordertest->whichVecs() tells us which vectors from solver state are the ones we want
00532           // ordertest->howMany() will tell us how many
00533           break;
00534         }
00536         //
00537         // check for restarting, i.e. the subspace is full
00538         //
00540         // this is for the Hermitian case, or non-Hermitian conjugate split situation.
00541         // --> for the Hermitian case the current subspace dimension needs to match the maximum subspace dimension
00542         // --> for the non-Hermitian case:
00543         //     --> if a conjugate pair was detected in the previous restart then the current subspace dimension needs to match the
00544         //         maximum subspace dimension (the BKS solver keeps one extra vector if the problem is non-Hermitian).
00545         //     --> if a conjugate pair was not detected in the previous restart then the current subspace dimension will be one less
00546         //         than the maximum subspace dimension.
00547         else if ( (bks_solver->getCurSubspaceDim() == bks_solver->getMaxSubspaceDim()) ||
00548                   (!_problem->isHermitian() && !_conjSplit && (bks_solver->getCurSubspaceDim()+1 == bks_solver->getMaxSubspaceDim())) ) {
00549 
00550           // Update the Schur form of the projected eigenproblem, then sort it.
00551           if (!bks_solver->isSchurCurrent()) {
00552             bks_solver->computeSchurForm( true );
00553 
00554             // Check for convergence, just in case we wait for every restart to check
00555             outputtest->checkStatus( &*bks_solver );
00556           }
00557 
00558           // Don't bother to restart if we've converged or reached the maximum number of restarts
00559           if ( numRestarts >= _maxRestarts || ordertest->getStatus() == Passed) {
00560             break; // break from while(1){bks_solver->iterate()}
00561           }
00562 
00563           // Start restarting timer and increment counter
00564 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
00565           Teuchos::TimeMonitor restimer(*_timerRestarting);
00566 #endif
00567           numRestarts++;
00568 
00569           int numConv = ordertest->howMany();
00570           cur_nevBlocks = _nevBlocks*_blockSize;
00571 
00572           // Add in extra blocks for restarting if either static or dynamic boundaries are being used.
00573           int moreNevBlocks = std::min( maxXtraBlocks, std::max( numConv/_blockSize, _xtra_nevBlocks) );
00574           if ( _dynXtraNev )
00575             cur_nevBlocks += moreNevBlocks * _blockSize;
00576           else if ( _xtra_nevBlocks )
00577             cur_nevBlocks += _xtra_nevBlocks * _blockSize;
00578 /*
00579             int cur_numConv = numConv;
00580             while ( (cur_nevBlocks < (_nevBlocks + maxXtraVecs)) && cur_numConv > 0 ) {
00581               cur_nevBlocks++;
00582               cur_numConv--;
00583 */
00584 
00585           printer->stream(Debug) << " Performing restart number " << numRestarts << " of " << _maxRestarts << std::endl << std::endl;
00586           printer->stream(Debug) << "   - Current NEV blocks is " << cur_nevBlocks << ", the minimum is " << _nevBlocks*_blockSize << std::endl;
00587 
00588           // Get the most current Ritz values before we continue.
00589           _ritzValues = bks_solver->getRitzValues();
00590 
00591           // Get the state.
00592           BlockKrylovSchurState<ScalarType,MV> oldState = bks_solver->getState();
00593 
00594           // Get the current dimension of the factorization
00595           int curDim = oldState.curDim;
00596 
00597           // Determine if the storage for the nev eigenvalues of interest splits a complex conjugate pair.
00598           std::vector<int> ritzIndex = bks_solver->getRitzIndex();
00599           if (ritzIndex[cur_nevBlocks-1]==1) {
00600             _conjSplit = true;
00601             cur_nevBlocks++;
00602           } else {
00603             _conjSplit = false;
00604           }
00605 
00606           // Print out a warning to the user if complex eigenvalues were found on the boundary of the restart subspace
00607           // and the eigenproblem is Hermitian.  This solver is not prepared to handle this situation.
00608           if (_problem->isHermitian() && _conjSplit)
00609           {
00610             printer->stream(Warnings)
00611               << " Eigenproblem is Hermitian, complex eigenvalues have been detected, and eigenvalues of interest split a conjugate pair!!!"
00612               << std::endl
00613               << " Block Krylov-Schur eigensolver cannot guarantee correct behavior in this situation, please turn Hermitian flag off!!!"
00614               << std::endl;
00615           }
00616 
00617           // Update the Krylov-Schur decomposition
00618 
00619           // Get a view of the Schur vectors of interest.
00620           Teuchos::SerialDenseMatrix<int,ScalarType> Qnev(Teuchos::View, *(oldState.Q), curDim, cur_nevBlocks);
00621 
00622           // Get a view of the current Krylov basis.
00623           std::vector<int> curind( curDim );
00624           for (int i=0; i<curDim; i++) { curind[i] = i; }
00625           Teuchos::RCP<const MV> basistemp = MVT::CloneView( *(oldState.V), curind );
00626 
00627           // Compute the new Krylov basis: Vnew = V*Qnev
00628           //
00629           // this will occur ex situ in workspace allocated for this purpose (tmpMV)
00630           // or in situ in the solver's memory space.
00631           //
00632           // we will also set a pointer for the location that the current factorization residual block (F),
00633           // currently located after the current basis in oldstate.V, will be moved to
00634           //
00635           Teuchos::RCP<MV> newF;
00636           if (_inSituRestart) {
00637             //
00638             // get non-const pointer to solver's basis so we can work in situ
00639             Teuchos::RCP<MV> solverbasis = Teuchos::rcp_const_cast<MV>(oldState.V);
00640             Teuchos::SerialDenseMatrix<int,ScalarType> copyQnev(Qnev);
00641             //
00642             // perform Householder QR of copyQnev = Q [D;0], where D is unit diag. We will want D below.
00643             std::vector<ScalarType> tau(cur_nevBlocks), work(cur_nevBlocks);
00644             int info;
00645             lapack.GEQRF(curDim,cur_nevBlocks,copyQnev.values(),copyQnev.stride(),&tau[0],&work[0],work.size(),&info);
00646             TEUCHOS_TEST_FOR_EXCEPTION(info != 0,std::logic_error,
00647                                "Anasazi::BlockKrylovSchurSolMgr::solve(): error calling GEQRF during restarting.");
00648             // we need to get the diagonal of D
00649             std::vector<ScalarType> d(cur_nevBlocks);
00650             for (int j=0; j<copyQnev.numCols(); j++) {
00651               d[j] = copyQnev(j,j);
00652             }
00653             if (printer->isVerbosity(Debug)) {
00654               Teuchos::SerialDenseMatrix<int,ScalarType> R(Teuchos::Copy,copyQnev,cur_nevBlocks,cur_nevBlocks);
00655               for (int j=0; j<R.numCols(); j++) {
00656                 R(j,j) = SCT::magnitude(R(j,j)) - 1.0;
00657                 for (int i=j+1; i<R.numRows(); i++) {
00658                   R(i,j) = zero;
00659                 }
00660               }
00661               printer->stream(Debug) << "||Triangular factor of Su - I||: " << R.normFrobenius() << std::endl;
00662             }
00663             //
00664             // perform implicit V*Qnev
00665             // this actually performs V*[Qnev Qtrunc*M] = [newV truncV], for some unitary M
00666             // we are interested in only the first cur_nevBlocks vectors of the result
00667             curind.resize(curDim);
00668             for (int i=0; i<curDim; i++) curind[i] = i;
00669             {
00670               Teuchos::RCP<MV> oldV = MVT::CloneViewNonConst(*solverbasis,curind);
00671               msutils::applyHouse(cur_nevBlocks,*oldV,copyQnev,tau,workMV);
00672             }
00673             // multiply newV*D
00674             // get pointer to new basis
00675             curind.resize(cur_nevBlocks);
00676             for (int i=0; i<cur_nevBlocks; i++) { curind[i] = i; }
00677             {
00678               Teuchos::RCP<MV> newV = MVT::CloneViewNonConst( *solverbasis, curind );
00679               MVT::MvScale(*newV,d);
00680             }
00681             // get pointer to new location for F
00682             curind.resize(_blockSize);
00683             for (int i=0; i<_blockSize; i++) { curind[i] = cur_nevBlocks + i; }
00684             newF = MVT::CloneViewNonConst( *solverbasis, curind );
00685           }
00686           else {
00687             // get pointer to first part of work space
00688             curind.resize(cur_nevBlocks);
00689             for (int i=0; i<cur_nevBlocks; i++) { curind[i] = i; }
00690             Teuchos::RCP<MV> tmp_newV = MVT::CloneViewNonConst(*workMV, curind );
00691             // perform V*Qnev
00692             MVT::MvTimesMatAddMv( one, *basistemp, Qnev, zero, *tmp_newV );
00693             tmp_newV = Teuchos::null;
00694             // get pointer to new location for F
00695             curind.resize(_blockSize);
00696             for (int i=0; i<_blockSize; i++) { curind[i] = cur_nevBlocks + i; }
00697             newF = MVT::CloneViewNonConst( *workMV, curind );
00698           }
00699 
00700           // Move the current factorization residual block (F) to the last block of newV.
00701           curind.resize(_blockSize);
00702           for (int i=0; i<_blockSize; i++) { curind[i] = curDim + i; }
00703           Teuchos::RCP<const MV> oldF = MVT::CloneView( *(oldState.V), curind );
00704           for (int i=0; i<_blockSize; i++) { curind[i] = i; }
00705           MVT::SetBlock( *oldF, curind, *newF );
00706           newF = Teuchos::null;
00707 
00708           // Update the Krylov-Schur quasi-triangular matrix.
00709           //
00710           // Create storage for the new Schur matrix of the Krylov-Schur factorization
00711           // Copy over the current quasi-triangular factorization of oldState.H which is stored in oldState.S.
00712           Teuchos::SerialDenseMatrix<int,ScalarType> oldS(Teuchos::View, *(oldState.S), cur_nevBlocks+_blockSize, cur_nevBlocks);
00713           Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > newH =
00714             Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( oldS ) );
00715           //
00716           // Get a view of the B block of the current factorization
00717           Teuchos::SerialDenseMatrix<int,ScalarType> oldB(Teuchos::View, *(oldState.H), _blockSize, _blockSize, curDim, curDim-_blockSize);
00718           //
00719           // Get a view of the a block row of the Schur vectors.
00720           Teuchos::SerialDenseMatrix<int,ScalarType> subQ(Teuchos::View, *(oldState.Q), _blockSize, cur_nevBlocks, curDim-_blockSize);
00721           //
00722           // Get a view of the new B block of the updated Krylov-Schur factorization
00723           Teuchos::SerialDenseMatrix<int,ScalarType> newB(Teuchos::View, *newH,  _blockSize, cur_nevBlocks, cur_nevBlocks);
00724           //
00725           // Compute the new B block.
00726           blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, _blockSize, cur_nevBlocks, _blockSize, one,
00727                      oldB.values(), oldB.stride(), subQ.values(), subQ.stride(), zero, newB.values(), newB.stride() );
00728 
00729 
00730           //
00731           // Set the new state and initialize the solver.
00732           BlockKrylovSchurState<ScalarType,MV> newstate;
00733           if (_inSituRestart) {
00734             newstate.V = oldState.V;
00735           } else {
00736             newstate.V = workMV;
00737           }
00738           newstate.H = newH;
00739           newstate.curDim = cur_nevBlocks;
00740           bks_solver->initialize(newstate);
00741 
00742         } // end of restarting
00744         //
00745         // we returned from iterate(), but none of our status tests Passed.
00746         // something is wrong, and it is probably our fault.
00747         //
00749         else {
00750           TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::BlockKrylovSchurSolMgr::solve(): Invalid return from bks_solver::iterate().");
00751         }
00752       }
00753       catch (const AnasaziError &err) {
00754         printer->stream(Errors)
00755           << "Anasazi::BlockKrylovSchurSolMgr::solve() caught unexpected exception from Anasazi::BlockKrylovSchur::iterate() at iteration " << bks_solver->getNumIters() << std::endl
00756           << err.what() << std::endl
00757           << "Anasazi::BlockKrylovSchurSolMgr::solve() returning Unconverged with no solutions." << std::endl;
00758         return Unconverged;
00759       }
00760     }
00761 
00762     //
00763     // free temporary space
00764     workMV = Teuchos::null;
00765 
00766     // Get the most current Ritz values before we return
00767     _ritzValues = bks_solver->getRitzValues();
00768 
00769     sol.numVecs = ordertest->howMany();
00770     printer->stream(Debug) << "ordertest->howMany() : " << sol.numVecs << std::endl;
00771     std::vector<int> whichVecs = ordertest->whichVecs();
00772 
00773     // Place any converged eigenpairs in the solution container.
00774     if (sol.numVecs > 0) {
00775 
00776       // Next determine if there is a conjugate pair on the boundary and resize.
00777       std::vector<int> tmpIndex = bks_solver->getRitzIndex();
00778       for (int i=0; i<(int)_ritzValues.size(); ++i) {
00779         printer->stream(Debug) << _ritzValues[i].realpart << " + i " << _ritzValues[i].imagpart << ", Index = " << tmpIndex[i] << std::endl;
00780       }
00781       printer->stream(Debug) << "Number of converged eigenpairs (before) = " << sol.numVecs << std::endl;
00782       for (int i=0; i<sol.numVecs; ++i) {
00783         printer->stream(Debug) << "whichVecs[" << i << "] = " << whichVecs[i] << ", tmpIndex[" << whichVecs[i] << "] = " << tmpIndex[whichVecs[i]] << std::endl;
00784       }
00785       if (tmpIndex[whichVecs[sol.numVecs-1]]==1) {
00786         printer->stream(Debug) << "There is a conjugate pair on the boundary, resizing sol.numVecs" << std::endl;
00787         whichVecs.push_back(whichVecs[sol.numVecs-1]+1);
00788         sol.numVecs++;
00789         for (int i=0; i<sol.numVecs; ++i) {
00790           printer->stream(Debug) << "whichVecs[" << i << "] = " << whichVecs[i] << ", tmpIndex[" << whichVecs[i] << "] = " << tmpIndex[whichVecs[i]] << std::endl;
00791         }
00792       }
00793 
00794       bool keepMore = false;
00795       int numEvecs = sol.numVecs;
00796       printer->stream(Debug) << "Number of converged eigenpairs (after) = " << sol.numVecs << std::endl;
00797       printer->stream(Debug) << "whichVecs[sol.numVecs-1] > sol.numVecs-1 : " << whichVecs[sol.numVecs-1] << " > " << sol.numVecs-1 << std::endl;
00798       if (whichVecs[sol.numVecs-1] > (sol.numVecs-1)) {
00799         keepMore = true;
00800         numEvecs = whichVecs[sol.numVecs-1]+1;  // Add 1 to fix zero-based indexing
00801         printer->stream(Debug) << "keepMore = true; numEvecs = " << numEvecs << std::endl;
00802       }
00803 
00804       // Next set the number of Ritz vectors that the iteration must compute and compute them.
00805       bks_solver->setNumRitzVectors(numEvecs);
00806       bks_solver->computeRitzVectors();
00807 
00808       // If the leading Ritz pairs are the converged ones, get the information
00809       // from the iteration to the solution container. Otherwise copy the necessary
00810       // information using 'whichVecs'.
00811       if (!keepMore) {
00812         sol.index = bks_solver->getRitzIndex();
00813         sol.Evals = bks_solver->getRitzValues();
00814         sol.Evecs = MVT::CloneCopy( *(bks_solver->getRitzVectors()) );
00815       }
00816 
00817       // Resize based on the number of solutions being returned and set the number of Ritz
00818       // vectors for the iteration to compute.
00819       sol.Evals.resize(sol.numVecs);
00820       sol.index.resize(sol.numVecs);
00821 
00822       // If the converged Ritz pairs are not the leading ones, copy over the information directly.
00823       if (keepMore) {
00824         std::vector<Anasazi::Value<ScalarType> > tmpEvals = bks_solver->getRitzValues();
00825         for (int vec_i=0; vec_i<sol.numVecs; ++vec_i) {
00826           sol.index[vec_i] = tmpIndex[whichVecs[vec_i]];
00827           sol.Evals[vec_i] = tmpEvals[whichVecs[vec_i]];
00828         }
00829         sol.Evecs = MVT::CloneCopy( *(bks_solver->getRitzVectors()), whichVecs );
00830       }
00831 
00832       // Set the solution space to be the Ritz vectors at this time.
00833       sol.Espace = sol.Evecs;
00834     }
00835   }
00836 
00837   // print final summary
00838   bks_solver->currentStatus(printer->stream(FinalSummary));
00839 
00840   // print timing information
00841 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
00842   if ( printer->isVerbosity( TimingDetails ) ) {
00843     Teuchos::TimeMonitor::summarize( printer->stream( TimingDetails ) );
00844   }
00845 #endif
00846 
00847   _problem->setSolution(sol);
00848   printer->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
00849 
00850   // get the number of iterations performed during this solve.
00851   _numIters = bks_solver->getNumIters();
00852 
00853   if (sol.numVecs < nev) {
00854     return Unconverged; // return from BlockKrylovSchurSolMgr::solve()
00855   }
00856   return Converged; // return from BlockKrylovSchurSolMgr::solve()
00857 }
00858 
00859 
00860 template <class ScalarType, class MV, class OP>
00861 void
00862 BlockKrylovSchurSolMgr<ScalarType,MV,OP>::setGlobalStatusTest(
00863     const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global)
00864 {
00865   globalTest_ = global;
00866 }
00867 
00868 template <class ScalarType, class MV, class OP>
00869 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
00870 BlockKrylovSchurSolMgr<ScalarType,MV,OP>::getGlobalStatusTest() const
00871 {
00872   return globalTest_;
00873 }
00874 
00875 template <class ScalarType, class MV, class OP>
00876 void
00877 BlockKrylovSchurSolMgr<ScalarType,MV,OP>::setDebugStatusTest(
00878     const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &debug)
00879 {
00880   debugTest_ = debug;
00881 }
00882 
00883 template <class ScalarType, class MV, class OP>
00884 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
00885 BlockKrylovSchurSolMgr<ScalarType,MV,OP>::getDebugStatusTest() const
00886 {
00887   return debugTest_;
00888 }
00889 
00890 } // end Anasazi namespace
00891 
00892 #endif /* ANASAZI_BLOCK_KRYLOV_SCHUR_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends