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