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 "AnasaziModalSolverUtils.hpp"
00042 
00043 #include "AnasaziBlockKrylovSchur.hpp"
00044 #include "AnasaziBasicSort.hpp"
00045 #include "AnasaziSVQBOrthoManager.hpp"
00046 #include "AnasaziBasicOrthoManager.hpp"
00047 #include "AnasaziStatusTestMaxIters.hpp"
00048 #include "AnasaziStatusTestResNorm.hpp"
00049 #include "AnasaziStatusTestOrderedResNorm.hpp"
00050 #include "AnasaziStatusTestCombo.hpp"
00051 #include "AnasaziStatusTestOutput.hpp"
00052 #include "AnasaziBasicOutputManager.hpp"
00053 #include "Teuchos_BLAS.hpp"
00054 #include "Teuchos_TimeMonitor.hpp"
00055 
00069 namespace Anasazi {
00070 
00071 template<class ScalarType, class MV, class OP>
00072 class BlockKrylovSchurSolMgr : public SolverManager<ScalarType,MV,OP> {
00073 
00074   private:
00075     typedef MultiVecTraits<ScalarType,MV> MVT;
00076     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00077     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00078     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00079     typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00080     
00081   public:
00082 
00084 
00085 
00101   BlockKrylovSchurSolMgr( const Teuchos::RefCountPtr<Eigenproblem<ScalarType,MV,OP> > &problem,
00102                              Teuchos::ParameterList &pl );
00103 
00105   virtual ~BlockKrylovSchurSolMgr() {};
00107   
00109 
00110 
00111   const Eigenproblem<ScalarType,MV,OP>& getProblem() const {
00112     return *_problem;
00113   }
00114 
00117   std::vector<Value<ScalarType> > getRitzValues() const {
00118     std::vector<Value<ScalarType> > ret( _ritzValues );
00119     return ret;
00120   }
00121 
00128    Teuchos::Array<Teuchos::RefCountPtr<Teuchos::Time> > getTimers() const {
00129      return tuple(_timerSolve, _timerRestarting);
00130    }
00131 
00133 
00135 
00136     
00155   ReturnType solve();
00157 
00158   private:
00159   Teuchos::RefCountPtr<Eigenproblem<ScalarType,MV,OP> > _problem;
00160   Teuchos::RefCountPtr<SortManager<ScalarType,MV,OP> > _sort;
00161 
00162   string _whch, _ortho; 
00163   MagnitudeType _ortho_kappa;
00164 
00165   MagnitudeType _convtol;
00166   int _maxRestarts;
00167   bool _relconvtol,_conjSplit;
00168   int _blockSize, _numBlocks, _stepSize, _nevBlocks, _xtra_nevBlocks;
00169   int _verbosity;
00170 
00171   std::vector<Value<ScalarType> > _ritzValues;
00172 
00173   Teuchos::RefCountPtr<Teuchos::Time> _timerSolve, _timerRestarting;
00174 };
00175 
00176 
00177 // Constructor
00178 template<class ScalarType, class MV, class OP>
00179 BlockKrylovSchurSolMgr<ScalarType,MV,OP>::BlockKrylovSchurSolMgr( 
00180         const Teuchos::RefCountPtr<Eigenproblem<ScalarType,MV,OP> > &problem,
00181         Teuchos::ParameterList &pl ) : 
00182   _problem(problem),
00183   _whch("LM"),
00184   _ortho("SVQB"),
00185   _ortho_kappa(-1.0),
00186   _convtol(0),
00187   _maxRestarts(20),
00188   _relconvtol(true),
00189   _conjSplit(false),
00190   _blockSize(0),
00191   _numBlocks(0),
00192   _stepSize(0),
00193   _nevBlocks(0),
00194   _xtra_nevBlocks(0),
00195   _verbosity(Anasazi::Errors),
00196   _timerSolve(Teuchos::TimeMonitor::getNewTimer("BKSSolMgr::solve()")),
00197   _timerRestarting(Teuchos::TimeMonitor::getNewTimer("BKSSolMgr restarting"))
00198 {
00199   TEST_FOR_EXCEPTION(_problem == Teuchos::null,               std::invalid_argument, "Problem not given to solver manager.");
00200   TEST_FOR_EXCEPTION(!_problem->isProblemSet(),               std::invalid_argument, "Problem not set.");
00201   TEST_FOR_EXCEPTION(_problem->getInitVec() == Teuchos::null, std::invalid_argument, "Problem does not contain initial vectors to clone from.");
00202 
00203   const int nev = _problem->getNEV();
00204 
00205   // convergence tolerance
00206   _convtol = pl.get("Convergence Tolerance",MT::prec());
00207   _relconvtol = pl.get("Relative Convergence Tolerance",_relconvtol);
00208   
00209   // maximum number of restarts
00210   _maxRestarts = pl.get("Maximum Restarts",_maxRestarts);
00211 
00212   // block size: default is 1
00213   _blockSize = pl.get("Block Size",1);
00214   TEST_FOR_EXCEPTION(_blockSize <= 0, std::invalid_argument,
00215                      "Anasazi::BlockKrylovSchurSolMgr: \"Block Size\" must be strictly positive.");
00216 
00217   // set the number of blocks we need to save to compute the nev eigenvalues of interest.
00218   _xtra_nevBlocks = pl.get("Extra NEV Blocks",0);
00219   if (nev%_blockSize) {
00220     _nevBlocks = nev/_blockSize + _xtra_nevBlocks + 1;
00221   } else {
00222     _nevBlocks = nev/_blockSize + _xtra_nevBlocks;
00223   }
00224 
00225   _numBlocks = pl.get("Num Blocks",3*_nevBlocks);
00226   TEST_FOR_EXCEPTION(_numBlocks <= _nevBlocks, std::invalid_argument,
00227                      "Anasazi::BlockKrylovSchurSolMgr: \"Num Blocks\" must be strictly positive and large enough to compute the requested eigenvalues.");
00228 
00229   TEST_FOR_EXCEPTION(_numBlocks*_blockSize > MVT::GetVecLength(*_problem->getInitVec()),
00230                      std::invalid_argument,
00231                      "Anasazi::BlockKrylovSchurSolMgr: Potentially impossible orthogonality requests. Reduce basis size.");
00232   
00233   // step size: the default is _maxRestarts*_numBlocks, so that Ritz values are only computed every restart.
00234   if (_maxRestarts) {
00235     _stepSize = pl.get("Step Size", (_maxRestarts+1)*(_numBlocks+1));
00236   } else {
00237     _stepSize = pl.get("Step Size", _numBlocks+1);
00238   }
00239   TEST_FOR_EXCEPTION(_stepSize < 1, std::invalid_argument,
00240                      "Anasazi::BlockKrylovSchurSolMgr: \"Step Size\" must be strictly positive.");
00241 
00242   // get the sort manager
00243   if (pl.isParameter("Sort Manager")) {
00244     _sort = Teuchos::getParameter<Teuchos::RefCountPtr<Anasazi::SortManager<ScalarType,MV,OP> > >(pl,"Sort Manager");
00245   } else {
00246     // which values to solve for
00247     _whch = pl.get("Which",_whch);
00248     TEST_FOR_EXCEPTION(_whch != "SM" && _whch != "LM" && _whch != "SR" && _whch != "LR" && _whch != "SI" && _whch != "LI",
00249                        std::invalid_argument, "Invalid sorting string.");
00250     _sort = Teuchos::rcp( new BasicSort<ScalarType,MV,OP>(_whch) );
00251   }
00252 
00253   // which orthogonalization to use
00254   _ortho = pl.get("Orthogonalization",_ortho);
00255   if (_ortho != "DGKS" && _ortho != "SVQB") {
00256     _ortho = "SVQB";
00257   }
00258 
00259   // which orthogonalization constant to use
00260   _ortho_kappa = pl.get("Orthogonalization Constant",_ortho_kappa);
00261 
00262   // verbosity level
00263   if (pl.isParameter("Verbosity")) {
00264     if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
00265       _verbosity = pl.get("Verbosity", _verbosity);
00266     } else {
00267       _verbosity = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
00268     }
00269   }
00270 }
00271 
00272 
00273 // solve()
00274 template<class ScalarType, class MV, class OP>
00275 ReturnType 
00276 BlockKrylovSchurSolMgr<ScalarType,MV,OP>::solve() {
00277 
00278   const int nev = _problem->getNEV();
00279   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00280   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00281 
00282   Teuchos::BLAS<int,ScalarType> blas;
00283 
00285   // Output manager
00286   Teuchos::RefCountPtr<BasicOutputManager<ScalarType> > printer = Teuchos::rcp( new BasicOutputManager<ScalarType>(_verbosity) );
00287 
00289   // Status tests
00290   //
00291   // convergence
00292   Teuchos::RefCountPtr<StatusTestOrderedResNorm<ScalarType,MV,OP> > convtest 
00293     = Teuchos::rcp( new StatusTestOrderedResNorm<ScalarType,MV,OP>(_sort,_convtol,nev,StatusTestOrderedResNorm<ScalarType,MV,OP>::RITZRES_2NORM,_relconvtol) );
00294 
00295   // printing StatusTest
00296   Teuchos::RefCountPtr<StatusTestOutput<ScalarType,MV,OP> > outputtest
00297     = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,convtest,1,Passed ) );
00298 
00300   // Orthomanager
00301   //
00302   Teuchos::RefCountPtr<OrthoManager<ScalarType,MV> > ortho; 
00303   if (_ortho=="SVQB") {
00304     ortho = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>(_problem->getM()) );
00305   } else if (_ortho=="DGKS") {
00306     if (_ortho_kappa <= 0) {
00307       ortho = Teuchos::rcp( new BasicOrthoManager<ScalarType,MV,OP>(_problem->getM()) );
00308     }
00309     else {
00310       ortho = Teuchos::rcp( new BasicOrthoManager<ScalarType,MV,OP>(_problem->getM(),_ortho_kappa) );
00311     }
00312   } else {
00313     TEST_FOR_EXCEPTION(_ortho!="SVQB"&&_ortho!="DGKS",std::logic_error,"Anasazi::BlockKrylovSchurSolMgr::solve(): Invalid orthogonalization type.");
00314   }
00315   
00316   // utils
00317   ModalSolverUtils<ScalarType,MV,OP> msutils(printer);
00318 
00320   // Parameter list
00321   Teuchos::ParameterList plist;
00322   plist.set("Block Size",_blockSize);
00323   plist.set("Num Blocks",_numBlocks);
00324   plist.set("Step Size",_stepSize);
00325   plist.set("Print Number of Ritz Values",_nevBlocks*_blockSize);
00326 
00328   // BlockKrylovSchur solver
00329   Teuchos::RefCountPtr<BlockKrylovSchur<ScalarType,MV,OP> > bks_solver 
00330     = Teuchos::rcp( new BlockKrylovSchur<ScalarType,MV,OP>(_problem,_sort,printer,outputtest,ortho,plist) );
00331   // set any auxiliary vectors defined in the problem
00332   Teuchos::RefCountPtr< const MV > probauxvecs = _problem->getAuxVecs();
00333   if (probauxvecs != Teuchos::null) {
00334     bks_solver->setAuxVecs( Teuchos::tuple< Teuchos::RefCountPtr<const MV> >(probauxvecs) );
00335   }
00336 
00337   // go ahead and initialize the solution to nothing in case we throw an exception
00338   Eigensolution<ScalarType,MV> sol;
00339   sol.numVecs = 0;
00340   _problem->setSolution(sol);
00341 
00342   int numRestarts = 0;
00343   int cur_nevBlocks = 0;
00344 
00345   // enter solve() iterations
00346   {
00347     Teuchos::TimeMonitor slvtimer(*_timerSolve);
00348   
00349   
00350     // tell bks_solver to iterate
00351     while (1) {
00352       try {
00353         bks_solver->iterate();
00354     
00355         // check convergence first
00356         if (convtest->getStatus() == Passed ) {
00357           // we have convergence
00358           // convtest->whichVecs() tells us which vectors from solver state are the ones we want
00359           // convtest->howMany() will tell us how many
00360           break;
00361         }
00362         // check for restarting:  this is for the Hermitian case, or non-Hermitian conjugate split situation.
00363         // --> for the Hermitian case the current subspace dimension needs to match the maximum subspace dimension
00364         // --> for the non-Hermitican case:
00365         //     --> if a conjugate pair was detected in the previous restart then the current subspace dimension needs to match the
00366         //         maximum subspace dimension (the BKS solver keeps one extra vector if the problem is non-Hermitian).
00367         //     --> if a conjugate pair was not detected in the previous restart then the current subspace dimension will be one less
00368         //         than the maximum subspace dimension.
00369         else if ( (bks_solver->getCurSubspaceDim() == bks_solver->getMaxSubspaceDim()) ||
00370                   (!_problem->isHermitian() && !_conjSplit && (bks_solver->getCurSubspaceDim()+1 == bks_solver->getMaxSubspaceDim())) ) {
00371   
00372           Teuchos::TimeMonitor restimer(*_timerRestarting);
00373   
00374           if ( numRestarts >= _maxRestarts ) {
00375             break; // break from while(1){bks_solver->iterate()}
00376           }
00377           numRestarts++;
00378   
00379           printer->stream(Debug) << " Performing restart number " << numRestarts << " of " << _maxRestarts << endl << endl;
00380   
00381           // Update the Schur form of the projected eigenproblem, then sort it.
00382           if (!bks_solver->isSchurCurrent())
00383             bks_solver->computeSchurForm( true );
00384           
00385           // Get the most current Ritz values before we continue.
00386           _ritzValues = bks_solver->getRitzValues();
00387   
00388           // Get the state.
00389           BlockKrylovSchurState<ScalarType,MV> oldState = bks_solver->getState();
00390           
00391           // Get the current dimension of the factorization
00392           int curDim = oldState.curDim;
00393   
00394           // Determine if the storage for the nev eigenvalues of interest splits a complex conjugate pair.
00395           std::vector<int> ritzIndex = bks_solver->getRitzIndex();
00396           if (ritzIndex[_nevBlocks*_blockSize-1]==1) {
00397             _conjSplit = true;
00398             cur_nevBlocks = _nevBlocks*_blockSize+1;
00399           } else {
00400             _conjSplit = false;
00401             cur_nevBlocks = _nevBlocks*_blockSize;
00402           }
00403           
00404           // Update the Krylov-Schur decomposition
00405   
00406           // Get a view of the Schur vectors of interest.
00407           Teuchos::SerialDenseMatrix<int,ScalarType> Qnev(Teuchos::View, *(oldState.Q), curDim, cur_nevBlocks);
00408   
00409           // Get a view of the current Krylov basis.
00410           std::vector<int> curind( curDim );
00411           for (int i=0; i<curDim; i++) { curind[i] = i; }
00412           Teuchos::RefCountPtr<const MV> basistemp = MVT::CloneView( *(oldState.V), curind );
00413   
00414           // Create a vector for the new Krylov basis.
00415           // Need cur_nevBlocks for the updated factorization and another block for the current factorization residual block (F).
00416           Teuchos::RefCountPtr<MV> newV = MVT::Clone( *(oldState.V), cur_nevBlocks + _blockSize );
00417           curind.resize(cur_nevBlocks);
00418           Teuchos::RefCountPtr<MV> tmp_newV = MVT::CloneView( *newV, curind );
00419           
00420           // Compute the new Krylov basis.
00421           // --> Vnew = V*Qnev
00422           MVT::MvTimesMatAddMv( one, *basistemp, Qnev, zero, *tmp_newV );
00423   
00424           // Move the current factorization residual block (F) to the last block of newV.
00425           curind.resize(_blockSize);
00426           for (int i=0; i<_blockSize; i++) { curind[i] = curDim + i; }
00427           Teuchos::RefCountPtr<const MV> oldF = MVT::CloneView( *(oldState.V), curind );
00428           for (int i=0; i<_blockSize; i++) { curind[i] = cur_nevBlocks + i; }
00429           MVT::SetBlock( *oldF, curind, *newV );
00430           
00431           // Update the Krylov-Schur quasi-triangular matrix.
00432           
00433           // Create storage for the new Schur matrix of the Krylov-Schur factorization
00434           // Copy over the current quasi-triangular factorization of oldState.H which is stored in oldState.S.
00435           Teuchos::SerialDenseMatrix<int,ScalarType> oldS(Teuchos::View, *(oldState.S), cur_nevBlocks+_blockSize, cur_nevBlocks);
00436           Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > newH = 
00437             Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( oldS ) );
00438   
00439           // Get a view of the B block of the current factorization
00440           Teuchos::SerialDenseMatrix<int,ScalarType> oldB(Teuchos::View, *(oldState.H), _blockSize, _blockSize, curDim, curDim-_blockSize);
00441   
00442           // Get a view of the a block row of the Schur vectors.
00443           Teuchos::SerialDenseMatrix<int,ScalarType> subQ(Teuchos::View, *(oldState.Q), _blockSize, cur_nevBlocks, curDim-_blockSize);
00444           
00445           // Get a view of the new B block of the updated Krylov-Schur factorization
00446           Teuchos::SerialDenseMatrix<int,ScalarType> newB(Teuchos::View, *newH,  _blockSize, cur_nevBlocks, cur_nevBlocks);
00447   
00448           // Compute the new B block.
00449           blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, _blockSize, cur_nevBlocks, _blockSize, one, 
00450                      oldB.values(), oldB.stride(), subQ.values(), subQ.stride(), zero, newB.values(), newB.stride() );
00451    
00452           // Set the new state and initialize the solver.
00453           BlockKrylovSchurState<ScalarType,MV> newstate;
00454           newstate.V = newV;
00455           newstate.H = newH;
00456           newstate.curDim = cur_nevBlocks;
00457           bks_solver->initialize(newstate);
00458   
00459         }
00460         else {
00461           TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::BlockKrylovSchurSolMgr::solve(): Invalid return from bks_solver::iterate().");
00462         }
00463       }
00464       catch (std::exception e) {
00465         printer->stream(Errors) << "Error! Caught exception in BlockKrylovSchur::iterate() at iteration " << bks_solver->getNumIters() << endl 
00466                                 << e.what() << endl;
00467         throw;
00468       }
00469     }
00470   
00471     // Get the most current Ritz values before we return
00472     _ritzValues = bks_solver->getRitzValues();
00473     
00474     sol.numVecs = convtest->howMany();
00475     if (sol.numVecs > 0) {
00476       sol.index = bks_solver->getRitzIndex();
00477       sol.Evals = bks_solver->getRitzValues();
00478       // Check to see if conjugate pair is on the boundary.
00479       if (sol.index[sol.numVecs-1]==1) {
00480         sol.numVecs++;
00481         sol.Evals.resize(sol.numVecs);
00482         sol.index.resize(sol.numVecs);
00483         bks_solver->setNumRitzVectors(sol.numVecs);
00484       } else {
00485         sol.Evals.resize(sol.numVecs);
00486         sol.index.resize(sol.numVecs);
00487         bks_solver->setNumRitzVectors(sol.numVecs);
00488       }
00489       bks_solver->computeRitzVectors();
00490       sol.Evecs = MVT::CloneCopy( *(bks_solver->getRitzVectors()) );
00491       sol.Espace = sol.Evecs;
00492     } 
00493   }
00494 
00495   // print final summary
00496   bks_solver->currentStatus(printer->stream(FinalSummary));
00497 
00498   // print timing information
00499   Teuchos::TimeMonitor::summarize(printer->stream(TimingDetails));
00500 
00501   _problem->setSolution(sol);
00502   printer->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << endl;
00503 
00504   if (sol.numVecs < nev) {
00505     return Unconverged; // return from BlockKrylovSchurSolMgr::solve() 
00506   }
00507   return Converged; // return from BlockKrylovSchurSolMgr::solve() 
00508 }
00509 
00510 
00511 } // end Anasazi namespace
00512 
00513 #endif /* ANASAZI_BLOCK_KRYLOV_SCHUR_SOLMGR_HPP */

Generated on Thu Sep 18 12:31:36 2008 for Anasazi by doxygen 1.3.9.1