AnasaziBlockDavidsonSolMgr.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_BLOCKDAVIDSON_SOLMGR_HPP
00031 #define ANASAZI_BLOCKDAVIDSON_SOLMGR_HPP
00032 
00037 #include "AnasaziConfigDefs.hpp"
00038 #include "AnasaziTypes.hpp"
00039 
00040 #include "AnasaziEigenproblem.hpp"
00041 #include "AnasaziSolverManager.hpp"
00042 #include "AnasaziModalSolverUtils.hpp"
00043 
00044 #include "AnasaziBlockDavidson.hpp"
00045 #include "AnasaziBasicSort.hpp"
00046 #include "AnasaziSVQBOrthoManager.hpp"
00047 #include "AnasaziStatusTestMaxIters.hpp"
00048 #include "AnasaziStatusTestResNorm.hpp"
00049 #include "AnasaziStatusTestOrderedResNorm.hpp"
00050 #include "AnasaziStatusTestCombo.hpp"
00051 #include "AnasaziStatusTestOutput.hpp"
00052 #include "AnasaziBasicOutputManager.hpp"
00053 #include "Teuchos_BLAS.hpp"
00054 
00055 
00076 namespace Anasazi {
00077 
00078 template<class ScalarType, class MV, class OP>
00079 class BlockDavidsonSolMgr : public SolverManager<ScalarType,MV,OP> {
00080 
00081   private:
00082     typedef MultiVecTraits<ScalarType,MV> MVT;
00083     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00084     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00085     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00086     typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00087     
00088   public:
00089 
00091 
00092 
00110   BlockDavidsonSolMgr( const Teuchos::RefCountPtr<Eigenproblem<ScalarType,MV,OP> > &problem,
00111                              Teuchos::ParameterList &pl );
00112 
00114   virtual ~BlockDavidsonSolMgr() {};
00116   
00118 
00119 
00120   const Eigenproblem<ScalarType,MV,OP>& getProblem() const {
00121     return *problem_;
00122   }
00123 
00125 
00127 
00128     
00150   ReturnType solve();
00152 
00153   private:
00154   Teuchos::RefCountPtr<Eigenproblem<ScalarType,MV,OP> > problem_;
00155 
00156   string whch_; 
00157 
00158   MagnitudeType convtol_, locktol_;
00159   int maxRestarts_;
00160   bool useLocking_;
00161   bool relconvtol_, rellocktol_;
00162   int blockSize_, numBlocks_;
00163   int maxLocked_;
00164   int verbosity_;
00165   int lockQuorum_;
00166 };
00167 
00168 
00169 // Constructor
00170 template<class ScalarType, class MV, class OP>
00171 BlockDavidsonSolMgr<ScalarType,MV,OP>::BlockDavidsonSolMgr( 
00172         const Teuchos::RefCountPtr<Eigenproblem<ScalarType,MV,OP> > &problem,
00173         Teuchos::ParameterList &pl ) : 
00174   problem_(problem),
00175   whch_("SR"),
00176   convtol_(0),
00177   locktol_(0),
00178   maxRestarts_(20),
00179   useLocking_(false),
00180   relconvtol_(true),
00181   rellocktol_(true),
00182   blockSize_(0),
00183   numBlocks_(0),
00184   maxLocked_(0),
00185   verbosity_(Anasazi::Errors),
00186   lockQuorum_(1)
00187 {
00188   TEST_FOR_EXCEPTION(problem_ == Teuchos::null,               std::invalid_argument, "Problem not given to solver manager.");
00189   TEST_FOR_EXCEPTION(!problem_->isProblemSet(),               std::invalid_argument, "Problem not set.");
00190   TEST_FOR_EXCEPTION(!problem_->isHermitian(),                std::invalid_argument, "Problem not symmetric.");
00191   TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null, std::invalid_argument, "Problem does not contain initial vectors to clone from.");
00192 
00193   // which values to solve for
00194   whch_ = pl.get("Which",whch_);
00195   TEST_FOR_EXCEPTION(whch_ != "SM" && whch_ != "LM" && whch_ != "SR" && whch_ != "LR",std::invalid_argument, "Invalid sorting string.");
00196 
00197   // convergence tolerance
00198   convtol_ = pl.get("Convergence Tolerance",MT::prec());
00199   relconvtol_ = pl.get("Relative Convergence Tolerance",relconvtol_);
00200   
00201   // locking tolerance
00202   useLocking_ = pl.get("Use Locking",useLocking_);
00203   rellocktol_ = pl.get("Relative Locking Tolerance",rellocktol_);
00204   locktol_ = pl.get("Locking Tolerance",convtol_/10.0);
00205 
00206   // maximum number of restarts
00207   maxRestarts_ = pl.get("Maximum Restarts",maxRestarts_);
00208 
00209   // block size: default is nev()
00210   blockSize_ = pl.get("Block Size",problem_->getNEV());
00211   TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
00212                      "Anasazi::BlockDavidsonSolMgr: \"Block Size\" must be strictly positive.");
00213   numBlocks_ = pl.get("Num Blocks",2);
00214   TEST_FOR_EXCEPTION(numBlocks_ <= 1, std::invalid_argument,
00215                      "Anasazi::BlockDavidsonSolMgr: \"Num Blocks\" must be strictly positive.");
00216 
00217   // max locked: default is nev(), must satisfy maxLocked_ + blockSize_ >= nev
00218   if (useLocking_) {
00219     maxLocked_ = pl.get("Max Locked",problem_->getNEV());
00220   }
00221   else {
00222     maxLocked_ = 0;
00223   }
00224   if (maxLocked_ == 0) {
00225     useLocking_ = false;
00226   }
00227   TEST_FOR_EXCEPTION(maxLocked_ < 0, std::invalid_argument,
00228                      "Anasazi::BlockDavidsonSolMgr: \"Max Locked\" must be positive.");
00229   TEST_FOR_EXCEPTION(maxLocked_ + blockSize_ < problem_->getNEV(), 
00230                      std::invalid_argument,
00231                      "Anasazi::BlockDavidsonSolMgr: Not enough storage space for requested number of eigenpairs.");
00232   TEST_FOR_EXCEPTION(numBlocks_*blockSize_ + maxLocked_ > MVT::GetVecLength(*problem_->getInitVec()),
00233                      std::invalid_argument,
00234                      "Anasazi::BlockDavidsonSolMgr: Potentially impossible orthogonality requests. Reduce basis size or locking size.");
00235 
00236   if (useLocking_) {
00237     lockQuorum_ = pl.get("Locking Quorum",lockQuorum_);
00238     TEST_FOR_EXCEPTION(lockQuorum_ <= 0,
00239                        std::invalid_argument,
00240                        "Anasazi::BlockDavidsonSolMgr: \"Locking Quorum\" must be strictly positive.");
00241   }
00242 
00243   // verbosity level
00244   if (pl.isParameter("Verbosity")) {
00245     if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
00246       verbosity_ = pl.get("Verbosity", verbosity_);
00247     } else {
00248       verbosity_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
00249     }
00250   }
00251 }
00252 
00253 
00254 // solve()
00255 template<class ScalarType, class MV, class OP>
00256 ReturnType 
00257 BlockDavidsonSolMgr<ScalarType,MV,OP>::solve() {
00258 
00259   const int nev = problem_->getNEV();
00260 
00262   // Sort manager
00263   Teuchos::RefCountPtr<BasicSort<ScalarType,MV,OP> > sorter = Teuchos::rcp( new BasicSort<ScalarType,MV,OP>(whch_) );
00264 
00266   // Output manager
00267   Teuchos::RefCountPtr<BasicOutputManager<ScalarType> > printer = Teuchos::rcp( new BasicOutputManager<ScalarType>(verbosity_) );
00268 
00270   // Status tests
00271   //
00272   // convergence
00273   Teuchos::RefCountPtr<StatusTestOrderedResNorm<ScalarType,MV,OP> > convtest 
00274       = Teuchos::rcp( new StatusTestOrderedResNorm<ScalarType,MV,OP>(sorter,convtol_,nev,StatusTestOrderedResNorm<ScalarType,MV,OP>::RES_ORTH,relconvtol_) );
00275   // locking
00276   Teuchos::RefCountPtr<StatusTestResNorm<ScalarType,MV,OP> > locktest;
00277   if (useLocking_) {
00278     locktest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(locktol_,lockQuorum_,StatusTestResNorm<ScalarType,MV,OP>::RES_ORTH,rellocktol_) );
00279   }
00280   // combo class
00281   Teuchos::Array<Teuchos::RefCountPtr<StatusTest<ScalarType,MV,OP> > > alltests;
00282   // for an OR test, the order doesn't matter
00283   alltests.push_back(convtest);
00284   if (locktest != Teuchos::null)   alltests.push_back(locktest);
00285   // combo: convergence || locking 
00286   Teuchos::RefCountPtr<StatusTestCombo<ScalarType,MV,OP> > combotest
00287     = Teuchos::rcp( new StatusTestCombo<ScalarType,MV,OP>( StatusTestCombo<ScalarType,MV,OP>::OR, alltests) );
00288   // printing StatusTest
00289   Teuchos::RefCountPtr<StatusTestOutput<ScalarType,MV,OP> > outputtest
00290     = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,Passed ) );
00291 
00293   // Orthomanager
00294   Teuchos::RefCountPtr<SVQBOrthoManager<ScalarType,MV,OP> > ortho 
00295     = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
00296 
00297   // utils
00298   ModalSolverUtils<ScalarType,MV,OP> msutils(printer);
00299 
00301   // Parameter list
00302   Teuchos::ParameterList plist;
00303   plist.set("Block Size",blockSize_);
00304   plist.set("Num Blocks",numBlocks_);
00305 
00307   // BlockDavidson solver
00308   Teuchos::RefCountPtr<BlockDavidson<ScalarType,MV,OP> > bd_solver 
00309     = Teuchos::rcp( new BlockDavidson<ScalarType,MV,OP>(problem_,sorter,printer,outputtest,ortho,plist) );
00310   // set any auxiliary vectors defined in the problem
00311   Teuchos::RefCountPtr< const MV > probauxvecs = problem_->getAuxVecs();
00312   if (probauxvecs != Teuchos::null) {
00313     bd_solver->setAuxVecs( Teuchos::tuple< Teuchos::RefCountPtr<const MV> >(probauxvecs) );
00314   }
00315 
00317   // Storage
00318   // for locked vectors
00319   int numlocked = 0;
00320   Teuchos::RefCountPtr<MV> lockvecs;
00321   // lockvecs is used to hold the locked eigenvectors, as well as for temporary storage when locking
00322   // when locking, we will lock some number of vectors numnew, where numnew <= maxlocked - curlocked
00323   // we will allocated numnew random vectors, which will go into workMV (see below)
00324   // we will also need numnew storage for the image of these random vectors under A and M; 
00325   // columns [curlocked+1,curlocked+numnew] will be used for this storage
00326   if (maxLocked_ > 0) {
00327     lockvecs = MVT::Clone(*problem_->getInitVec(),maxLocked_);
00328   }
00329   std::vector<MagnitudeType> lockvals;
00330   // workMV will be used when restarting because of 
00331   // a) full basis, in which case we need 2*blockSize, for X and KX
00332   // b) locking, in which case we will need as many vectors as in the current basis, 
00333   //    a number which will be always <= (numblocks-1)*blocksize
00334   //    [note: this is because we never lock with curdim == numblocks*blocksize]
00335   Teuchos::RefCountPtr<MV> workMV;
00336   if (useLocking_) {
00337     workMV = MVT::Clone(*problem_->getInitVec(),(numBlocks_-1)*blockSize_);
00338   }
00339   else {
00340     workMV = MVT::Clone(*problem_->getInitVec(),2*blockSize_);
00341   }
00342 
00343   // go ahead and initialize the solution to nothing in case we throw an exception
00344   Eigensolution<ScalarType,MV> sol;
00345   sol.numVecs = 0;
00346   problem_->setSolution(sol);
00347 
00348   int numRestarts = 0;
00349 
00350   // tell bd_solver to iterate
00351   while (1) {
00352     try {
00353       bd_solver->iterate();
00354 
00355       // check convergence first
00356       if (convtest->getStatus() == Passed ) {
00357         // we have convergence
00358         // convtest->whichVecs() tells us which vectors from lockvecs and solver state are the ones we want
00359         // convtest->howMany() will tell us how many
00360         break;
00361       }
00362       // check for restarting before locking: if we need to lock, it will happen after the restart
00363       else if ( bd_solver->getCurSubspaceDim() == bd_solver->getMaxSubspaceDim() ) {
00364 
00365         if ( numRestarts >= maxRestarts_ ) {
00366           break; // break from while(1){bd_solver->iterate()}
00367         }
00368         numRestarts++;
00369 
00370         printer->stream(IterationDetails) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << endl << endl;
00371 
00372         // the solver has filled its basis. 
00373         // the current eigenvectors will be used to restart the basis.
00374         std::vector<int> b1ind(blockSize_), b2ind(blockSize_);
00375         for (int i=0;i<blockSize_;i++) {
00376           b1ind[i] = i;
00377           b2ind[i] = blockSize_+i;
00378         }
00379 
00380         // these will be pointers into workMV
00381         Teuchos::RefCountPtr<MV> newV, newKV, newMV;
00382         { 
00383           newV = MVT::CloneView(*workMV,b1ind);
00384           MVT::SetBlock(*bd_solver->getRitzVectors(),b1ind,*newV);
00385         }
00386 
00387         if (problem_->getM() != Teuchos::null) {
00388           newMV = MVT::CloneView(*workMV,b2ind);
00389           OPT::Apply(*problem_->getM(),*newV,*newMV);
00390         }
00391         else {
00392           newMV = Teuchos::null;
00393         }
00394 
00395         // send this basis to the orthomanager to ensure orthonormality
00396         // we don't want anything in our Krylov basis that hasn't been sent through the ortho manager
00397         Teuchos::Array<Teuchos::RefCountPtr<const MV> > curauxvecs = bd_solver->getAuxVecs();
00398         Teuchos::Array<Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummy;
00399         ortho->projectAndNormalize(*newV,newMV,dummy,Teuchos::null,curauxvecs);
00400         // don't need newMV anymore, and the storage it points to will be needed for newKV
00401         newMV = Teuchos::null;
00402 
00403         // compute K*newV
00404         newKV = MVT::CloneView(*workMV,b2ind);
00405         OPT::Apply(*problem_->getOperator(),*newV,*newKV);
00406 
00407         // compute projected stiffness matrix
00408         Teuchos::RefCountPtr< Teuchos::SerialDenseMatrix<int,ScalarType> > 
00409             newKK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(blockSize_,blockSize_) );
00410         MVT::MvTransMv(SCT::one(),*newV,*newKV,*newKK);
00411 
00412         // initialize() will do the rest
00413         BlockDavidsonState<ScalarType,MV> newstate;
00414         newstate.V = newV;
00415         newstate.KK = newKK;
00416         bd_solver->initialize(newstate);
00417 
00418         // don't need either of these anymore
00419         newKV = Teuchos::null;
00420         newV = Teuchos::null;
00421       }
00422       // check locking if we didn't converge or restart
00423       else if (locktest != Teuchos::null && locktest->getStatus() == Passed) {
00424 
00425         // get number,indices of vectors to be locked
00426         int numnew = locktest->howMany();
00427         TEST_FOR_EXCEPTION(numnew <= 0,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): status test mistake.");
00428         std::vector<int> newind = locktest->whichVecs();
00429 
00430         // don't lock more than maxLocked_; we didn't allocate enough space.
00431         if (numlocked + numnew > maxLocked_) {
00432           numnew = maxLocked_ - numlocked;
00433           // just use the first of them
00434           newind.resize(numnew);
00435         }
00436 
00437         {
00438           // debug printing
00439           printer->print(Debug,"Locking vectors: ");
00440           for (unsigned int i=0; i<newind.size(); i++) {printer->stream(Debug) << " " << newind[i];}
00441           printer->print(Debug,"\n");
00442         }
00443 
00444         BlockDavidsonState<ScalarType,MV> state = bd_solver->getState();
00445         //
00446         // get current size of basis, build index, and get a view of the current basis and projected stiffness matrix
00447         int curdim = state.curDim;
00448 
00449         // workMV will be partitioned as follows:
00450         // workMV = [genV augV ...]
00451         // genV will be of size curdim - numnew, and contain the generated basis
00452         // augV will be of size numnew, and contain random directions to make up for
00453         //      the lost space
00454         //
00455         // lockvecs will be partitioned as follows:
00456         // lockvecs = [curlocked augTmp ...]
00457         // augTmp will be used for the storage of M*augV and K*augV
00458         // later, the locked vectors (storeg inside the eigensolver and referenced via 
00459         // an MV view) will be moved into lockvecs on top of augTmp when the space is no
00460         // longer needed.
00461         Teuchos::RefCountPtr<const MV> newvecs, curlocked;
00462         Teuchos::RefCountPtr<MV> genV, augV, augTmp;
00463         Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > newKK;
00464         newKK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(curdim,curdim) );
00465         // newvecs,newvals will contain the to-be-locked eigenpairs
00466         std::vector<MagnitudeType> newvals(numnew);
00467         {
00468           // setup genV
00469           std::vector<int> genind(curdim-numnew);
00470           for (int i=0; i<curdim-numnew; i++) genind[i] = i;
00471           genV = MVT::CloneView(*workMV,genind);
00472           // setup augV
00473           std::vector<int> augind(numnew);
00474           for (int i=0; i<numnew; i++) augind[i] = curdim-numnew+i;
00475           augV = MVT::CloneView(*workMV,augind);
00476           // setup curlocked
00477           if (numlocked > 0) {
00478             std::vector<int> lockind(numlocked);
00479             for (int i=0; i<numlocked; i++) lockind[i] = i;
00480             curlocked = MVT::CloneView(*lockvecs,lockind);
00481           }
00482           else {
00483             curlocked = Teuchos::null;
00484           }
00485           // setup augTmp
00486           std::vector<int> augtmpind(numnew); 
00487           for (int i=0; i<numnew; i++) augtmpind[i] = numlocked+i;
00488           augTmp = MVT::CloneView(*lockvecs,augtmpind);
00489           // setup newvecs
00490           newvecs = MVT::CloneView(*bd_solver->getRitzVectors(),newind);
00491           // setup newvals
00492           std::vector<Value<ScalarType> > allvals = bd_solver->getRitzValues();
00493           for (int i=0; i<numnew; i++) {
00494             newvals[i] = allvals[newind[i]].realpart;
00495           }
00496         }
00497 
00498         //
00499         // restart the solver
00500         //
00501         // we have to get the projected stiffness matrix from the solver, compute the projected eigenvectors and sort them
00502         // then all of the ones that don't correspond to the ones we wish to lock get orthogonalized using QR factorization
00503         // then we generate new, slightly smaller basis (in genV) and augment this basis with random information (in augV)
00504         // to preserve rank (which must be a muliple of blockSize)
00505         Teuchos::BLAS<int,ScalarType> blas;
00506         {
00507           const ScalarType ONE = SCT::one();
00508           const ScalarType ZERO = SCT::zero();
00509 
00510           std::vector<int> curind(curdim);
00511           for (int i=0; i<curdim; i++) curind[i] = i;
00512           Teuchos::RefCountPtr<const MV> curV = MVT::CloneView(*state.V,curind);
00513           Teuchos::RefCountPtr<const Teuchos::SerialDenseMatrix<int,ScalarType> > curKK;
00514           curKK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*state.KK,curdim,curdim) );
00515           //
00516           // compute eigenvectors of the projected stiffness matrix
00517           Teuchos::SerialDenseMatrix<int,ScalarType> S(curdim,curdim);
00518           std::vector<MagnitudeType> theta(curdim);
00519           int rank = curdim;
00520           msutils.directSolver(curdim,*curKK,0,&S,&theta,&rank,10);
00521           TEST_FOR_EXCEPTION(rank != curdim,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): direct solve did not compute all eigenvectors."); // this should never happen
00522           // 
00523           // sort the eigenvalues (so that we can order the eigenvectors)
00524           {
00525             std::vector<int> order(curdim);
00526             // make a ScalarType copy of theta
00527             sorter->sort(bd_solver.get(),curdim,theta,&order);
00528             //
00529             // apply the same ordering to the primitive ritz vectors
00530             msutils.permuteVectors(order,S);
00531           }
00532           // select the non-locked eigenvectors
00533           std::vector<int> unlockind(curdim-numnew);
00534           set_difference(curind.begin(),curind.end(),newind.begin(),newind.end(),unlockind.begin());
00535           Teuchos::SerialDenseMatrix<int,ScalarType> Sunlocked(curdim,curdim-numnew);
00536           for (int i=0; i<curdim-numnew; i++) {
00537             blas.COPY(curdim, S[unlockind[i]], 1, Sunlocked[i], 1);
00538           }
00539           // 
00540           // compute the qr factorization
00541           {
00542             Teuchos::LAPACK<int,ScalarType> lapack;
00543             int info, lwork = (curdim*numnew)*lapack.ILAENV( 1, "geqrf", "", curdim, curdim-numnew );
00544             std::vector<ScalarType> tau(curdim-numnew), work(lwork);
00545             lapack.GEQRF(curdim,curdim-numnew,Sunlocked.values(),Sunlocked.stride(),&tau[0],&work[0],lwork,&info);
00546             TEST_FOR_EXCEPTION(info != 0,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): Error calling GEQRF in locking code.");
00547             lapack.UNGQR(curdim,curdim-numnew,curdim-numnew,Sunlocked.values(),Sunlocked.stride(),&tau[0],&work[0],lwork,&info);
00548             TEST_FOR_EXCEPTION(info != 0,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): Error calling UNGQR in locking code.");
00549 
00550             if (printer->isVerbosity(Debug)) {
00551               Teuchos::SerialDenseMatrix<int,ScalarType> StS(curdim-numnew,curdim-numnew);
00552               int info = StS.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,Sunlocked,Sunlocked,ZERO);
00553               TEST_FOR_EXCEPTION(info != 0, std::logic_error, "Anasazi::BlockDavidsonSolMgr::solve(): Input error to SerialDenseMatrix::multiply.");
00554               for (int i=0; i<curdim-numnew; i++) {
00555                 StS(i,i) -= ONE;
00556               }
00557               printer->stream(Debug) << "Locking: Error in Snew^T Snew == I : " << StS.normFrobenius() << endl;
00558             }
00559           }
00560 
00561           // compute the first part of the new basis: genV = curV*Sunlocked
00562           MVT::MvTimesMatAddMv(ONE,*curV,Sunlocked,ZERO,*genV);
00563           //
00564           // compute leading block of the new projected stiffness matrix
00565           {
00566             Teuchos::SerialDenseMatrix<int,ScalarType> tmpKK(curdim,curdim-numnew),
00567                                                        newKK11(Teuchos::View,*newKK,curdim-numnew,curdim-numnew),
00568                                                        curKKsym(*curKK);
00569             for (int j=0; j<curdim; j++) {
00570               for (int i=j+1; i<curdim; i++) {
00571                 curKKsym(i,j) = curKKsym(j,i);
00572               }
00573             }
00574             int info = tmpKK.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,curKKsym,Sunlocked,ZERO);
00575             TEST_FOR_EXCEPTION(info != 0, std::logic_error, "Anasazi::BlockDavidsonSolMgr::solve(): Input error to SerialDenseMatrix::multiply.");
00576             info = newKK11.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,Sunlocked,tmpKK,ZERO);
00577             TEST_FOR_EXCEPTION(info != 0, std::logic_error, "Anasazi::BlockDavidsonSolMgr::solve(): Input error to SerialDenseMatrix::multiply.");
00578           }
00579           //
00580           // generate random data to fill the rest of the basis back to curdim
00581           MVT::MvRandom(*augV);
00582           // 
00583           // orthogonalize it against auxvecs, the current basis, and all locked vectors
00584           {
00585             Teuchos::Array<Teuchos::RefCountPtr<const MV> > against;
00586             Teuchos::Array<Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummy;
00587             if (probauxvecs != Teuchos::null) against.push_back(probauxvecs);
00588             if (curlocked != Teuchos::null)   against.push_back(curlocked);
00589             against.push_back(newvecs);
00590             against.push_back(genV);
00591             ortho->projectAndNormalize(*augV,Teuchos::null,dummy,Teuchos::null,against);
00592           }
00593           //
00594           // project the stiffness matrix on the new part
00595           {
00596             OPT::Apply(*problem_->getOperator(),*augV,*augTmp);
00597             Teuchos::SerialDenseMatrix<int,ScalarType> KK12(Teuchos::View,*newKK,curdim-numnew,numnew,0,curdim-numnew),
00598                                                        KK22(Teuchos::View,*newKK,numnew,numnew,curdim-numnew,curdim-numnew);
00599             MVT::MvTransMv(ONE,*genV,*augTmp,KK12);
00600             MVT::MvTransMv(ONE,*augV,*augTmp,KK22);
00601           }
00602         }
00603 
00604         // done with pointers
00605         augV = genV = augTmp = Teuchos::null;
00606         curlocked = Teuchos::null;
00607 
00608         // put newvecs into lockvecs, newvals into lockvals
00609         {
00610           lockvals.insert(lockvals.end(),newvals.begin(),newvals.end());
00611 
00612           std::vector<int> indlock(numnew);
00613           for (int i=0; i<numnew; i++) indlock[i] = numlocked+i;
00614           MVT::SetBlock(*newvecs,indlock,*lockvecs);
00615           newvecs = Teuchos::null;
00616 
00617           numlocked += numnew;
00618           std::vector<int> curind(numlocked);
00619           for (int i=0; i<numlocked; i++) curind[i] = i;
00620           curlocked = MVT::CloneView(*lockvecs,curind);
00621         }
00622         // add locked vecs as aux vecs, along with aux vecs from problem
00623         // add lockvals to convtest
00624         {
00625           convtest->setAuxVals(lockvals);
00626 
00627           Teuchos::Array< Teuchos::RefCountPtr<const MV> > aux;
00628           if (probauxvecs != Teuchos::null) aux.push_back(probauxvecs);
00629           aux.push_back(curlocked);
00630           bd_solver->setAuxVecs(aux);
00631         }
00632 
00633         // get pointer to new basis, KK
00634         {
00635           std::vector<int> curind(curdim);
00636           for (int i=0; i<curdim; i++) curind[i] = i;
00637           state.V = MVT::CloneView(*workMV,curind);
00638           state.KK = newKK;
00639           state.curDim = curdim;
00640           // clear the rest
00641           state.X = Teuchos::null;
00642           state.KX = Teuchos::null;
00643           state.MX = Teuchos::null;
00644           state.R = Teuchos::null;
00645           state.H = Teuchos::null;
00646           state.T = Teuchos::null;
00647           //
00648           // pass new state to the solver
00649           bd_solver->initialize(state);
00650         }
00651 
00652         if (numlocked == maxLocked_) {
00653           // disabled locking now by setting quorum to unreachable number
00654           locktest->setQuorum(blockSize_+1);
00655         }
00656       }
00657       else {
00658         TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): Invalid return from bd_solver::iterate().");
00659       }
00660     }
00661     catch (std::exception e) {
00662       printer->stream(Errors) << "Error! Caught exception in BlockDavidson::iterate() at iteration " << bd_solver->getNumIters() << endl 
00663                               << e.what() << endl;
00664       throw;
00665     }
00666   }
00667 
00668   // clear temp space
00669   workMV = Teuchos::null;
00670 
00671   sol.numVecs = convtest->howMany();
00672   if (sol.numVecs > 0) {
00673     sol.Evecs = MVT::Clone(*problem_->getInitVec(),sol.numVecs);
00674     sol.Espace = sol.Evecs;
00675     sol.Evals.resize(sol.numVecs);
00676     std::vector<MagnitudeType> vals(sol.numVecs);
00677 
00678     // copy them into the solution
00679     std::vector<int> which = convtest->whichVecs();
00680     // indices between [0,blockSize) refer to vectors/values in the solver
00681     // indices between [blockSize,blocksize+numlocked) refer to locked vectors/values
00682     // everything has already been ordered by the solver; we just have to partition the two references
00683     std::vector<int> inlocked(0), insolver(0);
00684     for (unsigned int i=0; i<which.size(); i++) {
00685       if (which[i] < blockSize_) {
00686         insolver.push_back(which[i]);
00687       }
00688       else {
00689         // sanity check
00690         TEST_FOR_EXCEPTION(which[i] >= numlocked+blockSize_,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): indexing mistake.");
00691         inlocked.push_back(which[i] - blockSize_);
00692       }
00693     }
00694 
00695     TEST_FOR_EXCEPTION(insolver.size() + inlocked.size() != (unsigned int)sol.numVecs,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): indexing mistake.");
00696 
00697     // set the vecs,vals in the solution
00698     if (insolver.size() > 0) {
00699       // set vecs
00700       int lclnum = insolver.size();
00701       std::vector<int> tosol(lclnum);
00702       for (int i=0; i<lclnum; i++) tosol[i] = i;
00703       Teuchos::RefCountPtr<const MV> v = MVT::CloneView(*bd_solver->getRitzVectors(),insolver);
00704       MVT::SetBlock(*v,tosol,*sol.Evecs);
00705       // set vals
00706       std::vector<Value<ScalarType> > fromsolver = bd_solver->getRitzValues();
00707       for (unsigned int i=0; i<insolver.size(); i++) {
00708         vals[i] = fromsolver[insolver[i]].realpart;
00709       }
00710     }
00711 
00712     // get the vecs,vals from locked storage
00713     if (inlocked.size() > 0) {
00714       int solnum = insolver.size();
00715       // set vecs
00716       int lclnum = inlocked.size();
00717       std::vector<int> tosol(lclnum);
00718       for (int i=0; i<lclnum; i++) tosol[i] = solnum + i;
00719       Teuchos::RefCountPtr<const MV> v = MVT::CloneView(*lockvecs,inlocked);
00720       MVT::SetBlock(*v,tosol,*sol.Evecs);
00721       // set vals
00722       for (unsigned int i=0; i<inlocked.size(); i++) {
00723         vals[i+solnum] = lockvals[inlocked[i]];
00724       }
00725     }
00726 
00727     // sort the eigenvalues and permute the eigenvectors appropriately
00728     {
00729       std::vector<int> order(sol.numVecs);
00730       sorter->sort(bd_solver.get(), sol.numVecs, vals, &order );
00731       // store the values in the Eigensolution
00732       for (int i=0; i<sol.numVecs; i++) {
00733         sol.Evals[i].realpart = vals[i];
00734         sol.Evals[i].imagpart = MT::zero();
00735       }
00736       // now permute the eigenvectors according to order
00737       msutils.permuteVectors(sol.numVecs,order,*sol.Evecs);
00738     }
00739 
00740     // setup sol.index, remembering that all eigenvalues are real so that index = {0,...,0}
00741     sol.index.resize(sol.numVecs,0);
00742   }
00743 
00744   // print final summary
00745   bd_solver->currentStatus(printer->stream(FinalSummary));
00746 
00747   // print timing information
00748   Teuchos::TimeMonitor::summarize(printer->stream(TimingDetails));
00749 
00750   problem_->setSolution(sol);
00751   printer->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << endl;
00752 
00753   if (sol.numVecs < nev) {
00754     return Unconverged; // return from BlockDavidsonSolMgr::solve() 
00755   }
00756   return Converged; // return from BlockDavidsonSolMgr::solve() 
00757 }
00758 
00759 
00760 } // end Anasazi namespace
00761 
00762 #endif /* ANASAZI_BLOCKDAVIDSON_SOLMGR_HPP */

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