AnasaziLOBPCGSolMgr.hpp

Go to the documentation of this file.
00001 
00002 // @HEADER
00003 // ***********************************************************************
00004 //
00005 //                 Anasazi: Block Eigensolvers Package
00006 //                 Copyright (2004) Sandia Corporation
00007 //
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 //
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA
00025 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00026 //
00027 // ***********************************************************************
00028 // @HEADER
00029 
00030 #ifndef ANASAZI_LOBPCG_SOLMGR_HPP
00031 #define ANASAZI_LOBPCG_SOLMGR_HPP
00032 
00037 #include "AnasaziConfigDefs.hpp"
00038 #include "AnasaziTypes.hpp"
00039 
00040 #include "AnasaziEigenproblem.hpp"
00041 #include "AnasaziSolverManager.hpp"
00042 
00043 #include "AnasaziLOBPCG.hpp"
00044 #include "AnasaziBasicSort.hpp"
00045 #include "AnasaziSVQBOrthoManager.hpp"
00046 #include "AnasaziStatusTestMaxIters.hpp"
00047 #include "AnasaziStatusTestResNorm.hpp"
00048 #include "AnasaziStatusTestOrderedResNorm.hpp"
00049 #include "AnasaziStatusTestCombo.hpp"
00050 #include "AnasaziStatusTestOutput.hpp"
00051 #include "AnasaziBasicOutputManager.hpp"
00052 
00053 
00081 namespace Anasazi {
00082 
00083 template<class ScalarType, class MV, class OP>
00084 class LOBPCGSolMgr : public SolverManager<ScalarType,MV,OP> {
00085 
00086   private:
00087     typedef MultiVecTraits<ScalarType,MV> MVT;
00088     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00089     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00090     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00091     typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00092     
00093   public:
00094 
00096 
00097 
00117   LOBPCGSolMgr( const Teuchos::RefCountPtr<Eigenproblem<ScalarType,MV,OP> > &problem,
00118                              Teuchos::ParameterList &pl );
00119 
00121   virtual ~LOBPCGSolMgr() {};
00123   
00125 
00126 
00127   const Eigenproblem<ScalarType,MV,OP>& getProblem() const {
00128     return *problem_;
00129   }
00130 
00132 
00134 
00135     
00162   ReturnType solve();
00164 
00165   private:
00166   Teuchos::RefCountPtr<Eigenproblem<ScalarType,MV,OP> > problem_;
00167 
00168   string whch_; 
00169 
00170   MagnitudeType convtol_, locktol_;
00171   int maxIters_;
00172   bool useLocking_;
00173   bool relconvtol_, rellocktol_;
00174   int blockSize_;
00175   bool fullOrtho_;
00176   int maxLocked_;
00177   int verbosity_;
00178   int lockQuorum_;
00179   bool recover_;
00180   Teuchos::RefCountPtr<LOBPCGState<ScalarType,MV> > state_;
00181 };
00182 
00183 
00184 // Constructor
00185 template<class ScalarType, class MV, class OP>
00186 LOBPCGSolMgr<ScalarType,MV,OP>::LOBPCGSolMgr( 
00187         const Teuchos::RefCountPtr<Eigenproblem<ScalarType,MV,OP> > &problem,
00188         Teuchos::ParameterList &pl ) : 
00189   problem_(problem),
00190   whch_("SR"),
00191   convtol_(MT::prec()),
00192   maxIters_(100),
00193   useLocking_(false),
00194   relconvtol_(true),
00195   rellocktol_(true),
00196   blockSize_(0),
00197   fullOrtho_(true),
00198   maxLocked_(0),
00199   verbosity_(Anasazi::Errors),
00200   lockQuorum_(1),
00201   recover_(true)
00202 {
00203   TEST_FOR_EXCEPTION(problem_ == Teuchos::null,              std::invalid_argument, "Problem not given to solver manager.");
00204   TEST_FOR_EXCEPTION(!problem_->isProblemSet(),              std::invalid_argument, "Problem not set.");
00205   TEST_FOR_EXCEPTION(!problem_->isHermitian(),               std::invalid_argument, "Problem not symmetric.");
00206   TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument, "Problem does not contain initial vectors to clone from.");
00207 
00208 
00209   // which values to solve for
00210   whch_ = pl.get("Which",whch_);
00211   TEST_FOR_EXCEPTION(whch_ != "SM" && whch_ != "LM" && whch_ != "SR" && whch_ != "LR",std::invalid_argument, "Invalid sorting string.");
00212 
00213   // convergence tolerance
00214   convtol_ = pl.get("Convergence Tolerance",convtol_);
00215   relconvtol_ = pl.get("Relative Convergence Tolerance",relconvtol_);
00216   
00217   // locking tolerance
00218   useLocking_ = pl.get("Use Locking",useLocking_);
00219   rellocktol_ = pl.get("Relative Locking Tolerance",rellocktol_);
00220   // default: should be less than convtol_
00221   locktol_ = convtol_/10;
00222   locktol_ = pl.get("Locking Tolerance",locktol_);
00223 
00224   // maximum number of iterations
00225   maxIters_ = pl.get("Maximum Iterations",maxIters_);
00226 
00227   // block size: default is nev()
00228   blockSize_ = pl.get("Block Size",problem_->getNEV());
00229   TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
00230                      "Anasazi::LOBPCGSolMgr: \"Block Size\" must be strictly positive.");
00231 
00232   // max locked: default is nev(), must satisfy maxLocked_ + blockSize_ >= nev
00233   if (useLocking_) {
00234     maxLocked_ = pl.get("Max Locked",problem_->getNEV());
00235   }
00236   else {
00237     maxLocked_ = 0;
00238   }
00239   if (maxLocked_ == 0) {
00240     useLocking_ = false;
00241   }
00242   TEST_FOR_EXCEPTION(maxLocked_ < 0, std::invalid_argument,
00243                      "Anasazi::LOBPCGSolMgr: \"Max Locked\" must be positive.");
00244   TEST_FOR_EXCEPTION(maxLocked_ + blockSize_ < problem_->getNEV(), 
00245                      std::invalid_argument,
00246                      "Anasazi::LOBPCGSolMgr: Not enough storage space for requested number of eigenpairs.");
00247 
00248   if (useLocking_) {
00249     lockQuorum_ = pl.get("Locking Quorum",lockQuorum_);
00250     TEST_FOR_EXCEPTION(lockQuorum_ <= 0,
00251                        std::invalid_argument,
00252                        "Anasazi::LOBPCGSolMgr: \"Locking Quorum\" must be strictly positive.");
00253   }
00254 
00255   // full orthogonalization: default true
00256   fullOrtho_ = pl.get("Full Ortho",fullOrtho_);
00257 
00258   // verbosity level
00259   if (pl.isParameter("Verbosity")) {
00260     if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
00261       verbosity_ = pl.get("Verbosity", verbosity_);
00262     } else {
00263       verbosity_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
00264     }
00265   }
00266 
00267   // recover from LOBPCGRitzFailure
00268   recover_ = pl.get("Recover",recover_);
00269 
00270   // get (optionally) an initial state
00271   if (pl.isParameter("Init")) {
00272     state_ = Teuchos::getParameter<Teuchos::RefCountPtr<Anasazi::LOBPCGState<ScalarType,MV> > >(pl,"Init");
00273   }
00274 }
00275 
00276 
00277 // solve()
00278 template<class ScalarType, class MV, class OP>
00279 ReturnType 
00280 LOBPCGSolMgr<ScalarType,MV,OP>::solve() {
00281 
00282   const int nev = problem_->getNEV();
00283 
00284 
00285 
00287   // Sort manager
00288   Teuchos::RefCountPtr<BasicSort<ScalarType,MV,OP> > sorter = Teuchos::rcp( new BasicSort<ScalarType,MV,OP>(whch_) );
00289 
00291   // Output manager
00292   Teuchos::RefCountPtr<BasicOutputManager<ScalarType> > printer = Teuchos::rcp( new BasicOutputManager<ScalarType>(verbosity_) );
00293 
00295   // Status tests
00296   //
00297   // maximum number of iterations: optional test
00298   Teuchos::RefCountPtr<StatusTestMaxIters<ScalarType,MV,OP> > maxtest;
00299   if (maxIters_ > 0) {
00300     maxtest = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>(maxIters_) );
00301   }
00302   // convergence
00303   Teuchos::RefCountPtr<StatusTestOrderedResNorm<ScalarType,MV,OP> > convtest 
00304       = Teuchos::rcp( new StatusTestOrderedResNorm<ScalarType,MV,OP>(sorter,convtol_,nev,StatusTestOrderedResNorm<ScalarType,MV,OP>::RES_ORTH,relconvtol_) );
00305   // locking
00306   Teuchos::RefCountPtr<StatusTestResNorm<ScalarType,MV,OP> > locktest;
00307   if (useLocking_) {
00308     locktest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(locktol_,lockQuorum_,StatusTestResNorm<ScalarType,MV,OP>::RES_ORTH,rellocktol_) );
00309   }
00310   Teuchos::Array<Teuchos::RefCountPtr<StatusTest<ScalarType,MV,OP> > > alltests;
00311   // for an OR test, the order doesn't matter
00312   alltests.push_back(convtest);
00313   if (maxtest != Teuchos::null) alltests.push_back(maxtest);
00314   if (locktest != Teuchos::null)   alltests.push_back(locktest);
00315   // combo: convergence || locking || max iters
00316   Teuchos::RefCountPtr<StatusTestCombo<ScalarType,MV,OP> > combotest
00317     = Teuchos::rcp( new StatusTestCombo<ScalarType,MV,OP>( StatusTestCombo<ScalarType,MV,OP>::OR, alltests) );
00318   // printing StatusTest
00319   Teuchos::RefCountPtr<StatusTestOutput<ScalarType,MV,OP> > outputtest;
00320   if ( printer->isVerbosity(Debug) ) {
00321     outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,Passed+Failed+Undefined ) );
00322   }
00323   else {
00324     outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,Passed ) );
00325   }
00326 
00328   // Orthomanager
00329   Teuchos::RefCountPtr<SVQBOrthoManager<ScalarType,MV,OP> > ortho 
00330     = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
00331 
00333   // Parameter list
00334   Teuchos::ParameterList plist;
00335   plist.set("Block Size",blockSize_);
00336   plist.set("Full Ortho",fullOrtho_);
00337 
00338   // utils
00339   ModalSolverUtils<ScalarType,MV,OP> msutils(printer);
00340 
00342   // LOBPCG solver
00343   Teuchos::RefCountPtr<LOBPCG<ScalarType,MV,OP> > lobpcg_solver 
00344     = Teuchos::rcp( new LOBPCG<ScalarType,MV,OP>(problem_,sorter,printer,outputtest,ortho,plist) );
00345   // set any auxiliary vectors defined in the problem
00346   Teuchos::RefCountPtr< const MV > probauxvecs = problem_->getAuxVecs();
00347   if (probauxvecs != Teuchos::null) {
00348     lobpcg_solver->setAuxVecs( Teuchos::tuple< Teuchos::RefCountPtr<const MV> >(probauxvecs) );
00349   }
00350 
00352   // Storage
00353   // 
00354   // lockvecs will contain eigenvectors that have been determined "locked" by the status test
00355   int numlocked = 0;
00356   Teuchos::RefCountPtr<MV> lockvecs;
00357   if (useLocking_) {
00358     lockvecs = MVT::Clone(*problem_->getInitVec(),maxLocked_);
00359   }
00360   std::vector<MagnitudeType> lockvals;
00361   // workMV will be used as work space for LOBPCGRitzFailure recovery and locking
00362   // it will be partitioned in these cases as follows:
00363   // for LOBPCGRitzFailure recovery:
00364   // workMV = [X H P OpX OpH OpP], where OpX OpH OpP will be used for K and M
00365   // total size: 2*3*blocksize
00366   // for locking
00367   // workMV = [X P MX MP], with MX,MP needing storage only if hasM==true
00368   // total size: 2*blocksize or 4*blocksize
00369   Teuchos::RefCountPtr<MV> workMV;
00370   if (fullOrtho_ == false && recover_ == true) {
00371     workMV = MVT::Clone(*problem_->getInitVec(),2*3*blockSize_);
00372   }
00373   else if (useLocking_) {
00374     if (problem_->getM() != Teuchos::null) {
00375       workMV = MVT::Clone(*problem_->getInitVec(),4*blockSize_);
00376     }
00377     else {
00378       workMV = MVT::Clone(*problem_->getInitVec(),2*blockSize_);
00379     }
00380   }
00381 
00382   // initialize the solution to nothing in case we throw an exception
00383   Eigensolution<ScalarType,MV> sol;
00384   sol.numVecs = 0;
00385   problem_->setSolution(sol);
00386 
00387   // initialize the solver if the user specified a state
00388   if (state_ != Teuchos::null) {
00389     lobpcg_solver->initialize(*state_);
00390   }
00391 
00392   // tell the lobpcg_solver to iterate
00393   while (1) {
00394     try {
00395       lobpcg_solver->iterate();
00396 
00397       // check convergence first
00398       if (convtest->getStatus() == Passed || (maxtest != Teuchos::null && maxtest->getStatus() == Passed) ) {
00399         // we have convergence or not
00400         // convtest->whichVecs() tells us which vectors from lockvecs and solver->getRitzVectors() are the ones we want
00401         // convtest->howMany() will tell us how many
00402         break;
00403       }
00404       // check locking if we didn't converge
00405       else if (locktest != Teuchos::null && locktest->getStatus() == Passed) {
00406 
00407         // remove the locked vectors,values from lobpcg_solver: put them in newvecs, newvals
00408         int numnew = locktest->howMany();
00409         TEST_FOR_EXCEPTION(numnew <= 0,std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): status test mistake.");
00410         // get the indices
00411         std::vector<int> indnew = locktest->whichVecs();
00412 
00413         // don't lock more than maxLocked_; we didn't allocate enough space.
00414         if (numlocked + numnew > maxLocked_) {
00415           numnew = maxLocked_ - numlocked;
00416           indnew.resize(numnew);
00417         }
00418 
00419         // the call below to lobpcg_solver->setAuxVecs() will reset the solver to unitialized with hasP() == false
00420         // store the hasP() state for use below
00421         bool hadP = lobpcg_solver->hasP();
00422 
00423         {
00424           // debug printing
00425           printer->print(Debug,"Locking vectors: ");
00426           for (unsigned int i=0; i<indnew.size(); i++) {printer->stream(Debug) << " " << indnew[i];}
00427           printer->print(Debug,"\n");
00428         }
00429         std::vector<MagnitudeType> newvals(numnew);
00430         Teuchos::RefCountPtr<const MV> newvecs;
00431         {
00432           // work in a local scope, to hide the variabes needed for extracting this info
00433           // get the vectors
00434           newvecs = MVT::CloneView(*lobpcg_solver->getRitzVectors(),indnew);
00435           // get the values
00436           std::vector<Value<ScalarType> > allvals = lobpcg_solver->getRitzValues();
00437           for (int i=0; i<numnew; i++) {
00438             newvals[i] = allvals[indnew[i]].realpart;
00439           }
00440         }
00441         // put newvecs into lockvecs
00442         {
00443           std::vector<int> indlock(numnew);
00444           for (int i=0; i<numnew; i++) indlock[i] = numlocked+i;
00445           MVT::SetBlock(*newvecs,indlock,*lockvecs);
00446           newvecs = Teuchos::null;
00447         }
00448         // put newvals into lockvals
00449         lockvals.insert(lockvals.end(),newvals.begin(),newvals.end());
00450         numlocked += numnew;
00451         // add locked vecs as aux vecs, along with aux vecs from problem
00452         {
00453           std::vector<int> indlock(numlocked);
00454           for (int i=0; i<numlocked; i++) indlock[i] = i;
00455           Teuchos::RefCountPtr<const MV> curlocked = MVT::CloneView(*lockvecs,indlock);
00456           if (probauxvecs != Teuchos::null) {
00457             lobpcg_solver->setAuxVecs( Teuchos::tuple< Teuchos::RefCountPtr<const MV> >(probauxvecs,curlocked) );
00458           }
00459           else {
00460             lobpcg_solver->setAuxVecs( Teuchos::tuple< Teuchos::RefCountPtr<const MV> >(curlocked) );
00461           }
00462         }
00463         // add locked vals to convtest
00464         convtest->setAuxVals(lockvals);
00465         // fill out the empty state in the solver
00466         {
00467           LOBPCGState<ScalarType,MV> state = lobpcg_solver->getState();
00468           Teuchos::RefCountPtr<MV> newstateX, newstateMX, newstateP, newstateMP;
00469           //
00470           // workMV will be partitioned as follows: workMV = [X P MX MP], 
00471           //
00472           // make a copy of the current X,MX state
00473           std::vector<int> bsind(blockSize_); 
00474           for (int i=0; i<blockSize_; i++) bsind[i] = i;
00475           newstateX = MVT::CloneView(*workMV,bsind);
00476           MVT::SetBlock(*state.X,bsind,*newstateX);
00477 
00478           if (state.MX != Teuchos::null) {
00479             std::vector<int> block3(blockSize_);
00480             for (int i=0; i<blockSize_; i++) block3[i] = 2*blockSize_+i;
00481             newstateMX = MVT::CloneView(*workMV,block3);
00482             MVT::SetBlock(*state.MX,bsind,*newstateMX);
00483           }
00484           //
00485           // get select part, set to random, apply M
00486           {
00487             Teuchos::RefCountPtr<MV> newX = MVT::CloneView(*newstateX,indnew);
00488             MVT::MvRandom(*newX);
00489 
00490             if (newstateMX != Teuchos::null) {
00491               Teuchos::RefCountPtr<MV> newMX = MVT::CloneView(*newstateMX,indnew);
00492               OPT::Apply(*problem_->getM(),*newX,*newMX);
00493             }
00494           }
00495 
00496           Teuchos::Array<Teuchos::RefCountPtr<const MV> > curauxvecs = lobpcg_solver->getAuxVecs();
00497           Teuchos::Array<Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummy;
00498           // ortho X against the aux vectors
00499           ortho->projectAndNormalize(*newstateX,newstateMX,dummy,Teuchos::null,curauxvecs);
00500 
00501           if (hadP) {
00502             //
00503             // get P and optionally MP, orthogonalize against X and auxiliary vectors
00504             std::vector<int> block2(blockSize_);
00505             for (int i=0; i<blockSize_; i++) block2[i] = blockSize_+i;
00506             newstateP = MVT::CloneView(*workMV,block2);
00507             MVT::SetBlock(*state.P,bsind,*newstateP);
00508 
00509             if (state.MP != Teuchos::null) {
00510               std::vector<int> block4(blockSize_);
00511               for (int i=0; i<blockSize_; i++) block4[i] = 3*blockSize_+i;
00512               newstateMP = MVT::CloneView(*workMV,block4);
00513               MVT::SetBlock(*state.MP,bsind,*newstateMP);
00514             }
00515 
00516             if (fullOrtho_) {
00517               // ortho P against the new aux vectors and new X
00518               curauxvecs.push_back(newstateX);
00519               ortho->projectAndNormalize(*newstateP,newstateMP,dummy,Teuchos::null,curauxvecs);
00520             }
00521             else {
00522               // ortho P against the new aux vectors
00523               ortho->projectAndNormalize(*newstateP,newstateMP,dummy,Teuchos::null,curauxvecs);
00524             }
00525           }
00526           // set the new state
00527           LOBPCGState<ScalarType,MV> newstate;
00528           newstate.X  = newstateX;
00529           newstate.MX = newstateMX;
00530           newstate.P  = newstateP;
00531           newstate.MP = newstateMP;
00532           lobpcg_solver->initialize(newstate);
00533         }
00534 
00535         if (numlocked == maxLocked_) {
00536           // disabled locking now
00537           locktest->setQuorum(blockSize_+1);
00538         }
00539       }
00540       else {
00541         TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): Invalid return from lobpcg_solver::iterate().");
00542       }
00543     }
00544     catch (LOBPCGRitzFailure re) {
00545       if (fullOrtho_==true || recover_==false) {
00546         // if we are already using full orthogonalization, there isn't much we can do here. 
00547         // the most recent information in the status tests is still valid, and can be used to extract/return the 
00548         // eigenpairs that have converged.
00549         printer->stream(Warnings) << "Error! Caught LOBPCGRitzFailure at iteration " << lobpcg_solver->getNumIters() << endl
00550                                 << "Will not try to recover." << endl;
00551         break; // while(1)
00552       }
00553       printer->stream(Warnings) << "Error! Caught LOBPCGRitzFailure at iteration " << lobpcg_solver->getNumIters() << endl
00554                               << "Full orthogonalization is off; will try to recover." << endl;
00555       // get the current "basis" from the solver, orthonormalize it, do a rayleigh-ritz, and restart with the ritz vectors
00556       // if there aren't enough, break and quit with what we have
00557       //
00558       // workMV = [X H P OpX OpH OpP], where OpX OpH OpP will be used for K and M
00559       LOBPCGState<ScalarType,MV> curstate = lobpcg_solver->getState();
00560       Teuchos::RefCountPtr<MV> restart, Krestart, Mrestart;
00561       int localsize = lobpcg_solver->hasP() ? 3*blockSize_ : 2*blockSize_;
00562       bool hasM = problem_->getM() != Teuchos::null;
00563       {
00564         std::vector<int> recind(localsize);
00565         for (int i=0; i<localsize; i++) recind[i] = i;
00566         restart = MVT::CloneView(*workMV,recind);
00567       }
00568       {
00569         std::vector<int> recind(localsize);
00570         for (int i=0; i<localsize; i++) recind[i] = localsize+i;
00571         Krestart = MVT::CloneView(*workMV,recind);
00572       }
00573       if (hasM) {
00574         Mrestart = Krestart;
00575       }
00576       else {
00577         Mrestart = restart;
00578       }
00579       //
00580       // set restart = [X H P] and Mrestart = M*[X H P]
00581       //
00582       // put X into [0 , blockSize)
00583       {
00584         std::vector<int> blk1(blockSize_);
00585         for (int i=0; i < blockSize_; i++) blk1[i] = i;
00586         MVT::SetBlock(*curstate.X,blk1,*restart);
00587 
00588         // put MX into [0 , blockSize)
00589         if (hasM) {
00590           MVT::SetBlock(*curstate.MX,blk1,*Mrestart);
00591         }
00592       }
00593       //
00594       // put H into [blockSize_ , 2*blockSize)
00595       {
00596         std::vector<int> blk2(blockSize_);
00597         for (int i=0; i < blockSize_; i++) blk2[i] = blockSize_+i;
00598         MVT::SetBlock(*curstate.H,blk2,*restart);
00599 
00600         // put MH into [blockSize_ , 2*blockSize)
00601         if (hasM) {
00602           MVT::SetBlock(*curstate.MH,blk2,*Mrestart);
00603         }
00604       }
00605       // optionally, put P into [2*blockSize,3*blockSize)
00606       if (localsize == 3*blockSize_) {
00607         std::vector<int> blk3(blockSize_);
00608         for (int i=0; i < blockSize_; i++) blk3[i] = 2*blockSize_+i;
00609         MVT::SetBlock(*curstate.P,blk3,*restart);
00610 
00611         // put MP into [2*blockSize,3*blockSize)
00612         if (hasM) {
00613           MVT::SetBlock(*curstate.MP,blk3,*Mrestart);
00614         }
00615       }
00616       // project against auxvecs and locked vecs, and orthonormalize the basis
00617       Teuchos::Array<Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummy;
00618       Teuchos::Array<Teuchos::RefCountPtr<const MV> > Q;
00619       {
00620         if (numlocked > 0) {
00621           std::vector<int> indlock(numlocked);
00622           for (int i=0; i<numlocked; i++) indlock[i] = i;
00623           Teuchos::RefCountPtr<const MV> curlocked = MVT::CloneView(*lockvecs,indlock);
00624           Q.push_back(curlocked);
00625         }
00626         if (probauxvecs != Teuchos::null) {
00627           Q.push_back(probauxvecs);
00628         }
00629       }
00630       int rank = ortho->projectAndNormalize(*restart,Mrestart,dummy,Teuchos::null,Q);
00631       if (rank < blockSize_) {
00632         // quit
00633         printer->stream(Errors) << "Error! Recovered basis only rank " << rank << ". Block size is " << blockSize_ << ".\n"
00634                                 << "Recovery failed." << endl;
00635         break;
00636       }
00637       // reduce multivec size if necessary
00638       if (rank < localsize) {
00639         localsize = rank;
00640         std::vector<int> redind(localsize);
00641         for (int i=0; i<localsize; i++) redind[i] = i;
00642         // grab the first part of restart,Krestart
00643         restart = MVT::CloneView(*restart,redind);
00644         Krestart = MVT::CloneView(*Krestart,redind);
00645         if (hasM) {
00646           Mrestart = Krestart;
00647         }
00648         else {
00649           Mrestart = restart;
00650         }
00651       }
00652       Teuchos::SerialDenseMatrix<int,ScalarType> KK(localsize,localsize), MM(localsize,localsize), S(localsize,localsize);
00653       std::vector<MagnitudeType> theta(localsize);
00654       // project the matrices
00655       //
00656       // MM = restart^H M restart
00657       MVT::MvTransMv(1.0,*restart,*Mrestart,MM);
00658       // 
00659       // compute Krestart = K*restart
00660       OPT::Apply(*problem_->getOperator(),*restart,*Krestart);
00661       //
00662       // KK = restart^H K restart
00663       MVT::MvTransMv(1.0,*restart,*Krestart,KK);
00664       rank = localsize;
00665       msutils.directSolver(localsize,KK,&MM,&S,&theta,&rank,1);
00666       if (rank < blockSize_) {
00667         printer->stream(Errors) << "Error! Recovered basis of rank " << rank << " produced only " << rank << "ritz vectors.\n"
00668                                 << "Block size is " << blockSize_ << ".\n"
00669                                 << "Recovery failed." << endl;
00670         break;
00671       }
00672       theta.resize(rank);
00673       //
00674       // sort the ritz values using the sort manager
00675       {
00676         Teuchos::BLAS<int,ScalarType> blas;
00677         std::vector<int> order(rank);
00678         // sort
00679         sorter->sort( lobpcg_solver.get(), rank, theta, &order );   // don't catch exception
00680         // Sort the primitive ritz vectors
00681         Teuchos::SerialDenseMatrix<int,ScalarType> curS(Teuchos::View,S,rank,rank);
00682         msutils.permuteVectors(order,curS);
00683       }
00684       //
00685       Teuchos::SerialDenseMatrix<int,ScalarType> S1(Teuchos::View,S,localsize,blockSize_);
00686       //
00687       // compute the ritz vectors: store them in Krestart
00688       LOBPCGState<ScalarType,MV> newstate;
00689       Teuchos::RefCountPtr<MV> newX; 
00690       {
00691         std::vector<int> bsind(blockSize_);
00692         for (int i=0; i<blockSize_; i++) bsind[i] = i;
00693         newX = MVT::CloneView(*Krestart,bsind);
00694       }
00695       MVT::MvTimesMatAddMv(1.0,*restart,S1,0.0,*newX);
00696       // send X and theta into the solver
00697       newstate.X = newX;
00698       theta.resize(blockSize_);
00699       newstate.T = Teuchos::rcp( &theta, false );
00700       // initialize
00701       lobpcg_solver->initialize(newstate);
00702     }
00703     // don't catch any other exceptions
00704   }
00705 
00706   sol.numVecs = convtest->howMany();
00707   if (sol.numVecs > 0) {
00708     sol.Evecs = MVT::Clone(*problem_->getInitVec(),sol.numVecs);
00709     sol.Espace = sol.Evecs;
00710     sol.Evals.resize(sol.numVecs);
00711     std::vector<MagnitudeType> vals(sol.numVecs);
00712 
00713     // copy them into the solution
00714     std::vector<int> which = convtest->whichVecs();
00715     // indices between [0,blockSize) refer to vectors/values in the solver
00716     // indices between [blockSize,blocksize+numlocked) refer to locked vectors/values
00717     // everything has already been ordered by the solver; we just have to partition the two references
00718     std::vector<int> inlocked(0), insolver(0);
00719     for (unsigned int i=0; i<which.size(); i++) {
00720       if (which[i] < blockSize_) {
00721         insolver.push_back(which[i]);
00722       }
00723       else {
00724         // sanity check
00725         TEST_FOR_EXCEPTION(which[i] >= numlocked+blockSize_,std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): indexing mistake.");
00726         inlocked.push_back(which[i] - blockSize_);
00727       }
00728     }
00729 
00730     TEST_FOR_EXCEPTION(insolver.size() + inlocked.size() != (unsigned int)sol.numVecs,std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): indexing mistake.");
00731 
00732     // set the vecs,vals in the solution
00733     if (insolver.size() > 0) {
00734       // set vecs
00735       int lclnum = insolver.size();
00736       std::vector<int> tosol(lclnum);
00737       for (int i=0; i<lclnum; i++) tosol[i] = i;
00738       Teuchos::RefCountPtr<const MV> v = MVT::CloneView(*lobpcg_solver->getRitzVectors(),insolver);
00739       MVT::SetBlock(*v,tosol,*sol.Evecs);
00740       // set vals
00741       std::vector<Value<ScalarType> > fromsolver = lobpcg_solver->getRitzValues();
00742       for (unsigned int i=0; i<insolver.size(); i++) {
00743         vals[i] = fromsolver[insolver[i]].realpart;
00744       }
00745     }
00746 
00747     // get the vecs,vals from locked storage
00748     if (inlocked.size() > 0) {
00749       int solnum = insolver.size();
00750       // set vecs
00751       int lclnum = inlocked.size();
00752       std::vector<int> tosol(lclnum);
00753       for (int i=0; i<lclnum; i++) tosol[i] = solnum + i;
00754       Teuchos::RefCountPtr<const MV> v = MVT::CloneView(*lockvecs,inlocked);
00755       MVT::SetBlock(*v,tosol,*sol.Evecs);
00756       // set vals
00757       for (unsigned int i=0; i<inlocked.size(); i++) {
00758         vals[i+solnum] = lockvals[inlocked[i]];
00759       }
00760     }
00761 
00762     // sort the eigenvalues and permute the eigenvectors appropriately
00763     {
00764       std::vector<int> order(sol.numVecs);
00765       sorter->sort( lobpcg_solver.get(), sol.numVecs, vals, &order );
00766       // store the values in the Eigensolution
00767       for (int i=0; i<sol.numVecs; i++) {
00768         sol.Evals[i].realpart = vals[i];
00769         sol.Evals[i].imagpart = MT::zero();
00770       }
00771       // now permute the eigenvectors according to order
00772       msutils.permuteVectors(sol.numVecs,order,*sol.Evecs);
00773     }
00774 
00775     // setup sol.index, remembering that all eigenvalues are real so that index = {0,...,0}
00776     sol.index.resize(sol.numVecs,0);
00777   }
00778 
00779   // print final summary
00780   lobpcg_solver->currentStatus(printer->stream(FinalSummary));
00781 
00782   // print timing information
00783   Teuchos::TimeMonitor::summarize(printer->stream(TimingDetails));
00784 
00785   problem_->setSolution(sol);
00786   printer->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << endl;
00787 
00788   if (sol.numVecs < nev) return Unconverged; // return from LOBPCGSolMgr::solve() 
00789   return Converged; // return from LOBPCGSolMgr::solve() 
00790 }
00791 
00792 
00793 } // end Anasazi namespace
00794 
00795 #endif /* ANASAZI_LOBPCG_SOLMGR_HPP */

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