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 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 Teuchos::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("Anasazi: BlockKrylovSchurSolMgr::solve()")),
00262   _timerRestarting(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchurSolMgr restarting")),
00263   _printNum(-1)
00264 {
00265   TEUCHOS_TEST_FOR_EXCEPTION(_problem == Teuchos::null,               std::invalid_argument, "Problem not given to solver manager.");
00266   TEUCHOS_TEST_FOR_EXCEPTION(!_problem->isProblemSet(),               std::invalid_argument, "Problem not set.");
00267   TEUCHOS_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   TEUCHOS_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   TEUCHOS_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   TEUCHOS_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   TEUCHOS_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     TEUCHOS_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     TEUCHOS_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             TEUCHOS_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             {
00592               Teuchos::RCP<MV> oldV = MVT::CloneViewNonConst(*solverbasis,curind);
00593               msutils::applyHouse(cur_nevBlocks,*oldV,copyQnev,tau,workMV);
00594             }
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             {
00600               Teuchos::RCP<MV> newV = MVT::CloneViewNonConst( *solverbasis, curind );
00601               MVT::MvScale(*newV,d);
00602             }
00603             // get pointer to new location for F
00604             curind.resize(_blockSize);
00605             for (int i=0; i<_blockSize; i++) { curind[i] = cur_nevBlocks + i; }
00606             newF = MVT::CloneViewNonConst( *solverbasis, curind );
00607           }
00608           else {
00609             // get pointer to first part of work space
00610             curind.resize(cur_nevBlocks);
00611             for (int i=0; i<cur_nevBlocks; i++) { curind[i] = i; }
00612             Teuchos::RCP<MV> tmp_newV = MVT::CloneViewNonConst(*workMV, curind );
00613             // perform V*Qnev
00614             MVT::MvTimesMatAddMv( one, *basistemp, Qnev, zero, *tmp_newV );
00615             tmp_newV = Teuchos::null;
00616             // get pointer to new location for F
00617             curind.resize(_blockSize);
00618             for (int i=0; i<_blockSize; i++) { curind[i] = cur_nevBlocks + i; }
00619             newF = MVT::CloneViewNonConst( *workMV, curind );
00620           }
00621 
00622           // Move the current factorization residual block (F) to the last block of newV.
00623           curind.resize(_blockSize);
00624           for (int i=0; i<_blockSize; i++) { curind[i] = curDim + i; }
00625           Teuchos::RCP<const MV> oldF = MVT::CloneView( *(oldState.V), curind );
00626           for (int i=0; i<_blockSize; i++) { curind[i] = i; }
00627           MVT::SetBlock( *oldF, curind, *newF );
00628           newF = Teuchos::null;
00629 
00630           // Update the Krylov-Schur quasi-triangular matrix.
00631           //
00632           // Create storage for the new Schur matrix of the Krylov-Schur factorization
00633           // Copy over the current quasi-triangular factorization of oldState.H which is stored in oldState.S.
00634           Teuchos::SerialDenseMatrix<int,ScalarType> oldS(Teuchos::View, *(oldState.S), cur_nevBlocks+_blockSize, cur_nevBlocks);
00635           Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > newH = 
00636             Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( oldS ) );
00637           //
00638           // Get a view of the B block of the current factorization
00639           Teuchos::SerialDenseMatrix<int,ScalarType> oldB(Teuchos::View, *(oldState.H), _blockSize, _blockSize, curDim, curDim-_blockSize);
00640           //
00641           // Get a view of the a block row of the Schur vectors.
00642           Teuchos::SerialDenseMatrix<int,ScalarType> subQ(Teuchos::View, *(oldState.Q), _blockSize, cur_nevBlocks, curDim-_blockSize);
00643           //
00644           // Get a view of the new B block of the updated Krylov-Schur factorization
00645           Teuchos::SerialDenseMatrix<int,ScalarType> newB(Teuchos::View, *newH,  _blockSize, cur_nevBlocks, cur_nevBlocks);
00646           //
00647           // Compute the new B block.
00648           blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, _blockSize, cur_nevBlocks, _blockSize, one, 
00649                      oldB.values(), oldB.stride(), subQ.values(), subQ.stride(), zero, newB.values(), newB.stride() );
00650 
00651 
00652           //
00653           // Set the new state and initialize the solver.
00654           BlockKrylovSchurState<ScalarType,MV> newstate;
00655           if (_inSituRestart) {
00656             newstate.V = oldState.V;
00657           } else {
00658             newstate.V = workMV;
00659           }
00660           newstate.H = newH;
00661           newstate.curDim = cur_nevBlocks;
00662           bks_solver->initialize(newstate);
00663   
00664         } // end of restarting
00666         //
00667         // we returned from iterate(), but none of our status tests Passed.
00668         // something is wrong, and it is probably our fault.
00669         //
00671         else {
00672           TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::BlockKrylovSchurSolMgr::solve(): Invalid return from bks_solver::iterate().");
00673         }
00674       }
00675       catch (const AnasaziError &err) {
00676         printer->stream(Errors) 
00677           << "Anasazi::BlockKrylovSchurSolMgr::solve() caught unexpected exception from Anasazi::BlockKrylovSchur::iterate() at iteration " << bks_solver->getNumIters() << std::endl
00678           << err.what() << std::endl
00679           << "Anasazi::BlockKrylovSchurSolMgr::solve() returning Unconverged with no solutions." << std::endl;
00680         return Unconverged;
00681       }
00682     }
00683 
00684     //
00685     // free temporary space
00686     workMV = Teuchos::null;
00687 
00688     // Get the most current Ritz values before we return
00689     _ritzValues = bks_solver->getRitzValues();
00690 
00691     sol.numVecs = ordertest->howMany();
00692     std::vector<int> whichVecs = ordertest->whichVecs();
00693 
00694     // Place any converged eigenpairs in the solution container.
00695     if (sol.numVecs > 0) {
00696 
00697       // Next determine if there is a conjugate pair on the boundary and resize.
00698       std::vector<int> tmpIndex = bks_solver->getRitzIndex();
00699       printer->stream(Debug) << "Number of converged eigenpairs (before) = " << sol.numVecs << std::endl;
00700       for (int i=0; i<sol.numVecs; ++i) {
00701         printer->stream(Debug) << "whichVecs[" << i << "] = " << whichVecs[i] << ", tmpIndex[" << whichVecs[i] << "] = " << tmpIndex[whichVecs[i]] << std::endl;
00702       }
00703       if (tmpIndex[whichVecs[sol.numVecs-1]]==1) {
00704         printer->stream(Debug) << "There is a conjugate pair on the boundary, resizing sol.numVecs" << std::endl;
00705   sol.numVecs++;
00706       }
00707 
00708       bool keepMore = false;
00709       int numEvecs = sol.numVecs;
00710       printer->stream(Debug) << "Number of converged eigenpairs (after) = " << sol.numVecs << std::endl;
00711       if (whichVecs[sol.numVecs-1] > (sol.numVecs-1)) {
00712   keepMore = true;
00713   numEvecs = whichVecs[sol.numVecs-1]+1;  // Add 1 to fix zero-based indexing
00714         printer->stream(Debug) << "keepMore = true; numEvecs = " << numEvecs << std::endl;
00715       }
00716 
00717       // Next set the number of Ritz vectors that the iteration must compute and compute them.
00718       bks_solver->setNumRitzVectors(numEvecs);
00719       bks_solver->computeRitzVectors();
00720 
00721       // If the leading Ritz pairs are the converged ones, get the information 
00722       // from the iteration to the solution container. Otherwise copy the necessary
00723       // information using 'whichVecs'.
00724       if (!keepMore) {
00725   sol.index = bks_solver->getRitzIndex();
00726   sol.Evals = bks_solver->getRitzValues();
00727   sol.Evecs = MVT::CloneCopy( *(bks_solver->getRitzVectors()) );
00728       }
00729 
00730       // Resize based on the number of solutions being returned and set the number of Ritz
00731       // vectors for the iteration to compute.
00732       sol.Evals.resize(sol.numVecs);
00733       sol.index.resize(sol.numVecs);
00734  
00735       // If the converged Ritz pairs are not the leading ones, copy over the information directly.      
00736       if (keepMore) {
00737   std::vector<Anasazi::Value<ScalarType> > tmpEvals = bks_solver->getRitzValues();
00738   for (int vec_i=0; vec_i<sol.numVecs; ++vec_i) {
00739     sol.index[vec_i] = tmpIndex[whichVecs[vec_i]];
00740     sol.Evals[vec_i] = tmpEvals[whichVecs[vec_i]];
00741   }
00742   sol.Evecs = MVT::CloneCopy( *(bks_solver->getRitzVectors()), whichVecs );
00743       }
00744 
00745       // Set the solution space to be the Ritz vectors at this time.
00746       sol.Espace = sol.Evecs;
00747     } 
00748   }
00749 
00750   // print final summary
00751   bks_solver->currentStatus(printer->stream(FinalSummary));
00752 
00753   // print timing information
00754   Teuchos::TimeMonitor::summarize(printer->stream(TimingDetails));
00755 
00756   _problem->setSolution(sol);
00757   printer->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
00758 
00759   // get the number of iterations performed during this solve.
00760   _numIters = bks_solver->getNumIters();
00761 
00762   if (sol.numVecs < nev) {
00763     return Unconverged; // return from BlockKrylovSchurSolMgr::solve() 
00764   }
00765   return Converged; // return from BlockKrylovSchurSolMgr::solve() 
00766 }
00767 
00768 
00769 template <class ScalarType, class MV, class OP>
00770 void 
00771 BlockKrylovSchurSolMgr<ScalarType,MV,OP>::setGlobalStatusTest(
00772     const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global) 
00773 {
00774   globalTest_ = global;
00775 }
00776 
00777 template <class ScalarType, class MV, class OP>
00778 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & 
00779 BlockKrylovSchurSolMgr<ScalarType,MV,OP>::getGlobalStatusTest() const 
00780 {
00781   return globalTest_;
00782 }
00783 
00784 template <class ScalarType, class MV, class OP>
00785 void 
00786 BlockKrylovSchurSolMgr<ScalarType,MV,OP>::setDebugStatusTest(
00787     const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &debug)
00788 {
00789   debugTest_ = debug;
00790 }
00791 
00792 template <class ScalarType, class MV, class OP>
00793 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & 
00794 BlockKrylovSchurSolMgr<ScalarType,MV,OP>::getDebugStatusTest() const
00795 {
00796   return debugTest_;
00797 }
00798 
00799 } // end Anasazi namespace
00800 
00801 #endif /* ANASAZI_BLOCK_KRYLOV_SCHUR_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends