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