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

Generated on Tue Jul 13 09:22:46 2010 for Anasazi by  doxygen 1.4.7