Anasazi  Version of the Day
AnasaziBlockDavidsonSolMgr.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Anasazi: Block Eigensolvers Package
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25 //
26 // ***********************************************************************
27 // @HEADER
28 
29 #ifndef ANASAZI_BLOCKDAVIDSON_SOLMGR_HPP
30 #define ANASAZI_BLOCKDAVIDSON_SOLMGR_HPP
31 
36 #include "AnasaziConfigDefs.hpp"
37 #include "AnasaziTypes.hpp"
38 
39 #include "AnasaziEigenproblem.hpp"
40 #include "AnasaziSolverManager.hpp"
41 #include "AnasaziSolverUtils.hpp"
42 
43 #include "AnasaziBlockDavidson.hpp"
44 #include "AnasaziBasicSort.hpp"
52 #include "Teuchos_BLAS.hpp"
53 #include "Teuchos_LAPACK.hpp"
54 #include "Teuchos_TimeMonitor.hpp"
55 #ifdef TEUCHOS_DEBUG
56 # include <Teuchos_FancyOStream.hpp>
57 #endif
58 #ifdef HAVE_MPI
59 #include <mpi.h>
60 #endif
61 
62 
63 
77 namespace Anasazi {
78 
79 
112 template<class ScalarType, class MV, class OP>
113 class BlockDavidsonSolMgr : public SolverManager<ScalarType,MV,OP> {
114 
115  private:
118  typedef Teuchos::ScalarTraits<ScalarType> SCT;
119  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
120  typedef Teuchos::ScalarTraits<MagnitudeType> MT;
121 
122  public:
123 
125 
126 
149  BlockDavidsonSolMgr( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
150  Teuchos::ParameterList &pl );
151 
153  virtual ~BlockDavidsonSolMgr() {};
155 
157 
158 
161  return *problem_;
162  }
163 
165  int getNumIters() const {
166  return numIters_;
167  }
168 
176  Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
177  return Teuchos::tuple(_timerSolve, _timerRestarting, _timerLocking);
178  }
179 
181 
183 
184 
206  ReturnType solve();
207 
209  void setGlobalStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global);
210 
212  const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getGlobalStatusTest() const;
213 
215  void setLockingStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &locking);
216 
218  const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getLockingStatusTest() const;
219 
221  void setDebugStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &debug);
222 
224  const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getDebugStatusTest() const;
225 
227 
228  private:
229  Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
230 
231  std::string whch_, ortho_;
232 
233  MagnitudeType convtol_, locktol_;
234  int maxRestarts_;
235  bool useLocking_;
236  bool relconvtol_, rellocktol_;
237  int blockSize_, numBlocks_, numIters_;
238  int maxLocked_;
239  int lockQuorum_;
240  bool inSituRestart_;
241  int numRestartBlocks_;
242  enum ResType convNorm_, lockNorm_;
243 
244  Teuchos::RCP<Teuchos::Time> _timerSolve, _timerRestarting, _timerLocking;
245 
246  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > globalTest_;
247  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > lockingTest_;
248  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > debugTest_;
249 
250  Teuchos::RCP<BasicOutputManager<ScalarType> > printer_;
251 };
252 
253 
254 // Constructor
255 template<class ScalarType, class MV, class OP>
257  const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
258  Teuchos::ParameterList &pl ) :
259  problem_(problem),
260  whch_("SR"),
261  ortho_("SVQB"),
262  convtol_(MT::prec()),
263  maxRestarts_(20),
264  useLocking_(false),
265  relconvtol_(true),
266  rellocktol_(true),
267  blockSize_(0),
268  numBlocks_(0),
269  numIters_(0),
270  maxLocked_(0),
271  lockQuorum_(1),
272  inSituRestart_(false),
273  numRestartBlocks_(1)
274 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
275  , _timerSolve(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidsonSolMgr::solve()")),
276  _timerRestarting(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidsonSolMgr restarting")),
277  _timerLocking(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidsonSolMgr locking"))
278 #endif
279 {
280  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
281  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument, "Problem not set.");
282  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isHermitian(), std::invalid_argument, "Problem not symmetric.");
283  TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument, "Problem does not contain initial vectors to clone from.");
284 
285  std::string strtmp;
286 
287  // which values to solve for
288  whch_ = pl.get("Which",whch_);
289  TEUCHOS_TEST_FOR_EXCEPTION(whch_ != "SM" && whch_ != "LM" && whch_ != "SR" && whch_ != "LR",std::invalid_argument, "Invalid sorting string.");
290 
291  // which orthogonalization to use
292  ortho_ = pl.get("Orthogonalization",ortho_);
293  if (ortho_ != "DGKS" && ortho_ != "SVQB") {
294  ortho_ = "SVQB";
295  }
296 
297  // convergence tolerance
298  convtol_ = pl.get("Convergence Tolerance",convtol_);
299  relconvtol_ = pl.get("Relative Convergence Tolerance",relconvtol_);
300  strtmp = pl.get("Convergence Norm",std::string("2"));
301  if (strtmp == "2") {
302  convNorm_ = RES_2NORM;
303  }
304  else if (strtmp == "M") {
305  convNorm_ = RES_ORTH;
306  }
307  else {
308  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
309  "Anasazi::BlockDavidsonSolMgr: Invalid Convergence Norm.");
310  }
311 
312  // locking tolerance
313  useLocking_ = pl.get("Use Locking",useLocking_);
314  rellocktol_ = pl.get("Relative Locking Tolerance",rellocktol_);
315  // default: should be less than convtol_
316  locktol_ = convtol_/10;
317  locktol_ = pl.get("Locking Tolerance",locktol_);
318  strtmp = pl.get("Locking Norm",std::string("2"));
319  if (strtmp == "2") {
320  lockNorm_ = RES_2NORM;
321  }
322  else if (strtmp == "M") {
323  lockNorm_ = RES_ORTH;
324  }
325  else {
326  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
327  "Anasazi::BlockDavidsonSolMgr: Invalid Locking Norm.");
328  }
329 
330  // maximum number of restarts
331  maxRestarts_ = pl.get("Maximum Restarts",maxRestarts_);
332 
333  // block size: default is nev()
334  blockSize_ = pl.get("Block Size",problem_->getNEV());
335  TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
336  "Anasazi::BlockDavidsonSolMgr: \"Block Size\" must be strictly positive.");
337  numBlocks_ = pl.get("Num Blocks",2);
338  TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 1, std::invalid_argument,
339  "Anasazi::BlockDavidsonSolMgr: \"Num Blocks\" must be >= 1.");
340 
341  // max locked: default is nev(), must satisfy maxLocked_ + blockSize_ >= nev
342  if (useLocking_) {
343  maxLocked_ = pl.get("Max Locked",problem_->getNEV());
344  }
345  else {
346  maxLocked_ = 0;
347  }
348  if (maxLocked_ == 0) {
349  useLocking_ = false;
350  }
351  TEUCHOS_TEST_FOR_EXCEPTION(maxLocked_ < 0, std::invalid_argument,
352  "Anasazi::BlockDavidsonSolMgr: \"Max Locked\" must be positive.");
353  TEUCHOS_TEST_FOR_EXCEPTION(maxLocked_ + blockSize_ < problem_->getNEV(),
354  std::invalid_argument,
355  "Anasazi::BlockDavidsonSolMgr: Not enough storage space for requested number of eigenpairs.");
356  TEUCHOS_TEST_FOR_EXCEPTION(static_cast<ptrdiff_t>(numBlocks_)*blockSize_ + maxLocked_ > MVT::GetGlobalLength(*problem_->getInitVec()),
357  std::invalid_argument,
358  "Anasazi::BlockDavidsonSolMgr: Potentially impossible orthogonality requests. Reduce basis size or locking size.");
359 
360  if (useLocking_) {
361  lockQuorum_ = pl.get("Locking Quorum",lockQuorum_);
362  TEUCHOS_TEST_FOR_EXCEPTION(lockQuorum_ <= 0,
363  std::invalid_argument,
364  "Anasazi::BlockDavidsonSolMgr: \"Locking Quorum\" must be strictly positive.");
365  }
366 
367  // restart size
368  numRestartBlocks_ = pl.get("Num Restart Blocks",numRestartBlocks_);
369  TEUCHOS_TEST_FOR_EXCEPTION(numRestartBlocks_ <= 0, std::invalid_argument,
370  "Anasazi::BlockDavidsonSolMgr: \"Num Restart Blocks\" must be strictly positive.");
371  TEUCHOS_TEST_FOR_EXCEPTION(numRestartBlocks_ >= numBlocks_, std::invalid_argument,
372  "Anasazi::BlockDavidsonSolMgr: \"Num Restart Blocks\" must be strictly less than \"Num Blocks\".");
373 
374  // restarting technique: V*Q or applyHouse(V,H,tau)
375  if (pl.isParameter("In Situ Restarting")) {
376  if (Teuchos::isParameterType<bool>(pl,"In Situ Restarting")) {
377  inSituRestart_ = pl.get("In Situ Restarting",inSituRestart_);
378  } else {
379  inSituRestart_ = ( Teuchos::getParameter<int>(pl,"In Situ Restarting") != 0 );
380  }
381  }
382 
383  // output stream
384  std::string fntemplate = "";
385  bool allProcs = false;
386  if (pl.isParameter("Output on all processors")) {
387  if (Teuchos::isParameterType<bool>(pl,"Output on all processors")) {
388  allProcs = pl.get("Output on all processors",allProcs);
389  } else {
390  allProcs = ( Teuchos::getParameter<int>(pl,"Output on all processors") != 0 );
391  }
392  }
393  fntemplate = pl.get("Output filename template",fntemplate);
394  int MyPID;
395 # ifdef HAVE_MPI
396  // Initialize MPI
397  int mpiStarted = 0;
398  MPI_Initialized(&mpiStarted);
399  if (mpiStarted) MPI_Comm_rank(MPI_COMM_WORLD, &MyPID);
400  else MyPID=0;
401 # else
402  MyPID = 0;
403 # endif
404  if (fntemplate != "") {
405  std::ostringstream MyPIDstr;
406  MyPIDstr << MyPID;
407  // replace %d in fntemplate with MyPID
408  int pos, start=0;
409  while ( (pos = fntemplate.find("%d",start)) != -1 ) {
410  fntemplate.replace(pos,2,MyPIDstr.str());
411  start = pos+2;
412  }
413  }
414  Teuchos::RCP<ostream> osp;
415  if (fntemplate != "") {
416  osp = Teuchos::rcp( new std::ofstream(fntemplate.c_str(),std::ios::out | std::ios::app) );
417  if (!*osp) {
418  osp = Teuchos::rcpFromRef(std::cout);
419  std::cout << "Anasazi::BlockDavidsonSolMgr::constructor(): Could not open file for write: " << fntemplate << std::endl;
420  }
421  }
422  else {
423  osp = Teuchos::rcpFromRef(std::cout);
424  }
425  // Output manager
426  int verbosity = Anasazi::Errors;
427  if (pl.isParameter("Verbosity")) {
428  if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
429  verbosity = pl.get("Verbosity", verbosity);
430  } else {
431  verbosity = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
432  }
433  }
434  if (allProcs) {
435  // print on all procs
436  printer_ = Teuchos::rcp( new BasicOutputManager<ScalarType>(verbosity,osp,MyPID) );
437  }
438  else {
439  // print only on proc 0
440  printer_ = Teuchos::rcp( new BasicOutputManager<ScalarType>(verbosity,osp,0) );
441  }
442 }
443 
444 
445 // solve()
446 template<class ScalarType, class MV, class OP>
447 ReturnType
449 
450  typedef SolverUtils<ScalarType,MV,OP> msutils;
451 
452  const int nev = problem_->getNEV();
453 
454 #ifdef TEUCHOS_DEBUG
455  Teuchos::RCP<Teuchos::FancyOStream>
456  out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(printer_->stream(Debug)));
457  out->setShowAllFrontMatter(false).setShowProcRank(true);
458  *out << "Entering Anasazi::BlockDavidsonSolMgr::solve()\n";
459 #endif
460 
462  // Sort manager
463  Teuchos::RCP<BasicSort<MagnitudeType> > sorter = Teuchos::rcp( new BasicSort<MagnitudeType>(whch_) );
464 
466  // Status tests
467  //
468  // convergence
469  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convtest;
470  if (globalTest_ == Teuchos::null) {
471  convtest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(convtol_,nev,convNorm_,relconvtol_) );
472  }
473  else {
474  convtest = globalTest_;
475  }
476  Teuchos::RCP<StatusTestWithOrdering<ScalarType,MV,OP> > ordertest
477  = Teuchos::rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sorter,nev) );
478  // locking
479  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > locktest;
480  if (useLocking_) {
481  if (lockingTest_ == Teuchos::null) {
482  locktest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(locktol_,lockQuorum_,lockNorm_,rellocktol_) );
483  }
484  else {
485  locktest = lockingTest_;
486  }
487  }
488  // for a non-short-circuited OR test, the order doesn't matter
489  Teuchos::Array<Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
490  alltests.push_back(ordertest);
491  if (locktest != Teuchos::null) alltests.push_back(locktest);
492  if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
493 
494  Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > combotest
496  // printing StatusTest
497  Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
498  if ( printer_->isVerbosity(Debug) ) {
499  outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,Passed+Failed+Undefined ) );
500  }
501  else {
502  outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,Passed ) );
503  }
504 
506  // Orthomanager
507  Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho;
508  if (ortho_=="SVQB") {
509  ortho = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
510  } else if (ortho_=="DGKS") {
511  ortho = Teuchos::rcp( new BasicOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
512  } else {
513  TEUCHOS_TEST_FOR_EXCEPTION(ortho_!="SVQB"&&ortho_!="DGKS",std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): Invalid orthogonalization type.");
514  }
515 
517  // Parameter list
518  Teuchos::ParameterList plist;
519  plist.set("Block Size",blockSize_);
520  plist.set("Num Blocks",numBlocks_);
521 
523  // BlockDavidson solver
524  Teuchos::RCP<BlockDavidson<ScalarType,MV,OP> > bd_solver
525  = Teuchos::rcp( new BlockDavidson<ScalarType,MV,OP>(problem_,sorter,printer_,outputtest,ortho,plist) );
526  // set any auxiliary vectors defined in the problem
527  Teuchos::RCP< const MV > probauxvecs = problem_->getAuxVecs();
528  if (probauxvecs != Teuchos::null) {
529  bd_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
530  }
531 
533  // Storage
534  //
535  // lockvecs will contain eigenvectors that have been determined "locked" by the status test
536  int curNumLocked = 0;
537  Teuchos::RCP<MV> lockvecs;
538  // lockvecs is used to hold the locked eigenvectors, as well as for temporary storage when locking.
539  // when locking, we will lock some number of vectors numnew, where numnew <= maxlocked - curlocked
540  // we will produce numnew random vectors, which will go into the space with the new basis.
541  // we will also need numnew storage for the image of these random vectors under A and M;
542  // columns [curlocked+1,curlocked+numnew] will be used for this storage
543  if (maxLocked_ > 0) {
544  lockvecs = MVT::Clone(*problem_->getInitVec(),maxLocked_);
545  }
546  std::vector<MagnitudeType> lockvals;
547  //
548  // Restarting occurs under two scenarios: when the basis is full and after locking.
549  //
550  // For the former, a new basis of size blockSize*numRestartBlocks is generated using the current basis
551  // and the most significant primitive Ritz vectors (projected eigenvectors).
552  // [S,L] = eig(KK)
553  // S = [Sr St] // some for "r"estarting, some are "t"runcated
554  // newV = V*Sr
555  // KK_new = newV'*K*newV = Sr'*V'*K*V*Sr = Sr'*KK*Sr
556  // Therefore, the only multivector operation needed is for the generation of newV.
557  //
558  // * If the multiplication is explicit, it requires a workspace of blockSize*numRestartBlocks vectors.
559  // This space must be specifically allocated for that task, as we don't have any space of that size.
560  // It (workMV) will be allocated at the beginning of solve()
561  // * Optionally, the multiplication can be performed implicitly, via a Householder QR factorization of
562  // Sr. This can be done in situ, using the basis multivector contained in the solver. This requires
563  // that we cast away the const on the multivector returned from getState(). Workspace for this approach
564  // is a single vector. the solver's internal storage must be preserved (X,MX,KX,R), requiring us to
565  // allocate this vector.
566  //
567  // For the latter (restarting after locking), the new basis is the same size as existing basis. If numnew
568  // vectors are locked, they are deflated from the current basis and replaced with randomly generated
569  // vectors.
570  // [S,L] = eig(KK)
571  // S = [Sl Su] // partitioned: "l"ocked and "u"nlocked
572  // newL = V*Sl = X(locked)
573  // defV = V*Su
574  // augV = rand(numnew) // orthogonal to oldL,newL,defV,auxvecs
575  // newV = [defV augV]
576  // Kknew = newV'*K*newV = [Su'*KK*Su defV'*K*augV]
577  // [augV'*K*defV augV'*K*augV]
578  // locked = [oldL newL]
579  // Clearly, this operation is more complicated than the previous.
580  // Here is a list of the significant computations that need to be performed:
581  // - newL will be put into space in lockvecs, but will be copied from getState().X at the end
582  // - defV,augV will be stored in workspace the size of the current basis.
583  // - If inSituRestart==true, we compute defV in situ in bd_solver::V_ and
584  // put augV at the end of bd_solver::V_
585  // - If inSituRestart==false, we must have curDim vectors available for
586  // defV and augV; we will allocate a multivector (workMV) at the beginning of solve()
587  // for this purpose.
588  // - M*augV and K*augV are needed; they will be stored in lockvecs. As a result, newL will
589  // not be put into lockvecs until the end.
590  //
591  // Therefore, we must allocate workMV when ((maxRestarts_ > 0) || (useLocking_ == true)) && inSituRestart == false
592  // It will be allocated to size (numBlocks-1)*blockSize
593  //
594  Teuchos::RCP<MV> workMV;
595  if (inSituRestart_ == false) {
596  // we need storage space to restart, either if we may lock or if may restart after a full basis
597  if (useLocking_==true || maxRestarts_ > 0) {
598  workMV = MVT::Clone(*problem_->getInitVec(),(numBlocks_-1)*blockSize_);
599  }
600  else {
601  // we will never need to restart.
602  workMV = Teuchos::null;
603  }
604  }
605  else { // inSituRestart_ == true
606  // we will restart in situ, if we need to restart
607  // three situation remain:
608  // - never restart => no space needed
609  // - only restart for locking (i.e., never restart full) => no space needed
610  // - restart for full basis => need one vector
611  if (maxRestarts_ > 0) {
612  workMV = MVT::Clone(*problem_->getInitVec(),1);
613  }
614  else {
615  workMV = Teuchos::null;
616  }
617  }
618 
619  // some consts and utils
620  const ScalarType ONE = SCT::one();
621  const ScalarType ZERO = SCT::zero();
622  Teuchos::LAPACK<int,ScalarType> lapack;
623  Teuchos::BLAS<int,ScalarType> blas;
624 
625  // go ahead and initialize the solution to nothing in case we throw an exception
627  sol.numVecs = 0;
628  problem_->setSolution(sol);
629 
630  int numRestarts = 0;
631 
632  // enter solve() iterations
633  {
634 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
635  Teuchos::TimeMonitor slvtimer(*_timerSolve);
636 #endif
637 
638  // tell bd_solver to iterate
639  while (1) {
640  try {
641  bd_solver->iterate();
642 
644  //
645  // check user-specified debug test; if it passed, return an exception
646  //
648  if (debugTest_ != Teuchos::null && debugTest_->getStatus() == Passed) {
649  throw AnasaziError("Anasazi::BlockDavidsonSolMgr::solve(): User-specified debug status test returned Passed.");
650  }
652  //
653  // check convergence next
654  //
656  else if (ordertest->getStatus() == Passed ) {
657  // we have convergence
658  // ordertest->whichVecs() tells us which vectors from lockvecs and solver state are the ones we want
659  // ordertest->howMany() will tell us how many
660  break;
661  }
663  //
664  // check for restarting before locking: if we need to lock, it will happen after the restart
665  //
667  else if ( bd_solver->getCurSubspaceDim() == bd_solver->getMaxSubspaceDim() ) {
668 
669 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
670  Teuchos::TimeMonitor restimer(*_timerRestarting);
671 #endif
672 
673  if ( numRestarts >= maxRestarts_ ) {
674  break; // break from while(1){bd_solver->iterate()}
675  }
676  numRestarts++;
677 
678  printer_->stream(IterationDetails) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << std::endl << std::endl;
679 
680  BlockDavidsonState<ScalarType,MV> state = bd_solver->getState();
681  int curdim = state.curDim;
682  int newdim = numRestartBlocks_*blockSize_;
683 
684 # ifdef TEUCHOS_DEBUG
685  {
686  std::vector<Value<ScalarType> > ritzvalues = bd_solver->getRitzValues();
687  *out << "Ritz values from solver:\n";
688  for (unsigned int i=0; i<ritzvalues.size(); i++) {
689  *out << ritzvalues[i].realpart << " ";
690  }
691  *out << "\n";
692  }
693 # endif
694 
695  //
696  // compute eigenvectors of the projected problem
697  Teuchos::SerialDenseMatrix<int,ScalarType> S(curdim,curdim);
698  std::vector<MagnitudeType> theta(curdim);
699  int rank = curdim;
700 # ifdef TEUCHOS_DEBUG
701  {
702  std::stringstream os;
703  os << "KK before HEEV...\n"
704  << *state.KK << "\n";
705  *out << os.str();
706  }
707 # endif
708  int info = msutils::directSolver(curdim,*state.KK,Teuchos::null,S,theta,rank,10);
709  TEUCHOS_TEST_FOR_EXCEPTION(info != 0 ,std::logic_error,
710  "Anasazi::BlockDavidsonSolMgr::solve(): error calling SolverUtils::directSolver."); // this should never happen
711  TEUCHOS_TEST_FOR_EXCEPTION(rank != curdim,std::logic_error,
712  "Anasazi::BlockDavidsonSolMgr::solve(): direct solve did not compute all eigenvectors."); // this should never happen
713 
714 # ifdef TEUCHOS_DEBUG
715  {
716  std::stringstream os;
717  *out << "Ritz values from heev(KK):\n";
718  for (unsigned int i=0; i<theta.size(); i++) *out << theta[i] << " ";
719  os << "\nRitz vectors from heev(KK):\n"
720  << S << "\n";
721  *out << os.str();
722  }
723 # endif
724 
725  //
726  // sort the eigenvalues (so that we can order the eigenvectors)
727  {
728  std::vector<int> order(curdim);
729  sorter->sort(theta,Teuchos::rcpFromRef(order),curdim);
730  //
731  // apply the same ordering to the primitive ritz vectors
732  msutils::permuteVectors(order,S);
733  }
734 # ifdef TEUCHOS_DEBUG
735  {
736  std::stringstream os;
737  *out << "Ritz values from heev(KK) after sorting:\n";
738  std::copy(theta.begin(), theta.end(), std::ostream_iterator<ScalarType>(*out, " "));
739  os << "\nRitz vectors from heev(KK) after sorting:\n"
740  << S << "\n";
741  *out << os.str();
742  }
743 # endif
744  //
745  // select the significant primitive ritz vectors
746  Teuchos::SerialDenseMatrix<int,ScalarType> Sr(Teuchos::View,S,curdim,newdim);
747 # ifdef TEUCHOS_DEBUG
748  {
749  std::stringstream os;
750  os << "Significant primitive Ritz vectors:\n"
751  << Sr << "\n";
752  *out << os.str();
753  }
754 # endif
755  //
756  // generate newKK = Sr'*KKold*Sr
757  Teuchos::SerialDenseMatrix<int,ScalarType> newKK(newdim,newdim);
758  {
759  Teuchos::SerialDenseMatrix<int,ScalarType> KKtmp(curdim,newdim),
760  KKold(Teuchos::View,*state.KK,curdim,curdim);
761  int teuchosRet;
762  // KKtmp = KKold*Sr
763  teuchosRet = KKtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,KKold,Sr,ZERO);
764  TEUCHOS_TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
765  "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
766  // newKK = Sr'*KKtmp = Sr'*KKold*Sr
767  teuchosRet = newKK.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,Sr,KKtmp,ZERO);
768  TEUCHOS_TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
769  "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
770  // make it Hermitian in memory
771  for (int j=0; j<newdim-1; ++j) {
772  for (int i=j+1; i<newdim; ++i) {
773  newKK(i,j) = SCT::conjugate(newKK(j,i));
774  }
775  }
776  }
777 # ifdef TEUCHOS_DEBUG
778  {
779  std::stringstream os;
780  os << "Sr'*KK*Sr:\n"
781  << newKK << "\n";
782  *out << os.str();
783  }
784 # endif
785 
786  // prepare new state
788  rstate.curDim = newdim;
789  rstate.KK = Teuchos::rcpFromRef(newKK);
790  //
791  // we know that newX = newV*Sr(:,1:bS) = oldV*S(:1:bS) = oldX
792  // the restarting preserves the Ritz vectors and residual
793  // for the Ritz values, we want all of the values associated with newV.
794  // these have already been placed at the beginning of theta
795  rstate.X = state.X;
796  rstate.KX = state.KX;
797  rstate.MX = state.MX;
798  rstate.R = state.R;
799  rstate.T = Teuchos::rcp( new std::vector<MagnitudeType>(theta.begin(),theta.begin()+newdim) );
800 
801  if (inSituRestart_ == true) {
802 # ifdef TEUCHOS_DEBUG
803  *out << "Beginning in-situ restart...\n";
804 # endif
805  //
806  // get non-const pointer to solver's basis so we can work in situ
807  Teuchos::RCP<MV> solverbasis = Teuchos::rcp_const_cast<MV>(state.V);
808  //
809  // perform Householder QR of Sr = Q [D;0], where D is unit diag.
810  // WARNING: this will overwrite Sr; however, we do not need Sr anymore after this
811  std::vector<ScalarType> tau(newdim), work(newdim);
812  int geqrf_info;
813  lapack.GEQRF(curdim,newdim,Sr.values(),Sr.stride(),&tau[0],&work[0],work.size(),&geqrf_info);
814  TEUCHOS_TEST_FOR_EXCEPTION(geqrf_info != 0,std::logic_error,
815  "Anasazi::BlockDavidsonSolMgr::solve(): error calling GEQRF during restarting.");
816  if (printer_->isVerbosity(Debug)) {
817  Teuchos::SerialDenseMatrix<int,ScalarType> R(Teuchos::Copy,Sr,newdim,newdim);
818  for (int j=0; j<newdim; j++) {
819  R(j,j) = SCT::magnitude(R(j,j)) - 1.0;
820  for (int i=j+1; i<newdim; i++) {
821  R(i,j) = ZERO;
822  }
823  }
824  printer_->stream(Debug) << "||Triangular factor of Sr - I||: " << R.normFrobenius() << std::endl;
825  }
826  //
827  // perform implicit oldV*Sr
828  // this actually performs oldV*[Sr Su*M] = [newV truncV], for some unitary M
829  // we are actually interested in only the first newdim vectors of the result
830  {
831  std::vector<int> curind(curdim);
832  for (int i=0; i<curdim; i++) curind[i] = i;
833  Teuchos::RCP<MV> oldV = MVT::CloneViewNonConst(*solverbasis,curind);
834  msutils::applyHouse(newdim,*oldV,Sr,tau,workMV);
835  }
836  //
837  // put the new basis into the state for initialize()
838  // the new basis is contained in the the first newdim columns of solverbasis
839  // initialize() will recognize that pointer bd_solver.V_ == pointer rstate.V, and will neglect the copy.
840  rstate.V = solverbasis;
841  }
842  else { // inSituRestart == false)
843 # ifdef TEUCHOS_DEBUG
844  *out << "Beginning ex-situ restart...\n";
845 # endif
846  // newV = oldV*Sr, explicitly. workspace is in workMV
847  std::vector<int> curind(curdim), newind(newdim);
848  for (int i=0; i<curdim; i++) curind[i] = i;
849  for (int i=0; i<newdim; i++) newind[i] = i;
850  Teuchos::RCP<const MV> oldV = MVT::CloneView(*state.V,curind);
851  Teuchos::RCP<MV> newV = MVT::CloneViewNonConst(*workMV ,newind);
852 
853  MVT::MvTimesMatAddMv(ONE,*oldV,Sr,ZERO,*newV);
854  //
855  // put the new basis into the state for initialize()
856  rstate.V = newV;
857  }
858 
859  //
860  // send the new state to the solver
861  bd_solver->initialize(rstate);
862  } // end of restarting
864  //
865  // check locking if we didn't converge or restart
866  //
868  else if (locktest != Teuchos::null && locktest->getStatus() == Passed) {
869 
870 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
871  Teuchos::TimeMonitor lcktimer(*_timerLocking);
872 #endif
873 
874  //
875  // get current state
876  BlockDavidsonState<ScalarType,MV> state = bd_solver->getState();
877  const int curdim = state.curDim;
878 
879  //
880  // get number,indices of vectors to be locked
881  TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() <= 0,std::logic_error,
882  "Anasazi::BlockDavidsonSolMgr::solve(): status test mistake: howMany() non-positive.");
883  TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() != (int)locktest->whichVecs().size(),std::logic_error,
884  "Anasazi::BlockDavidsonSolMgr::solve(): status test mistake: howMany() not consistent with whichVecs().");
885  // we should have room to lock vectors; otherwise, locking should have been deactivated
886  TEUCHOS_TEST_FOR_EXCEPTION(curNumLocked == maxLocked_,std::logic_error,
887  "Anasazi::BlockDavidsonSolMgr::solve(): status test mistake: locking not deactivated.");
888  //
889  // don't lock more than maxLocked_; we didn't allocate enough space.
890  std::vector<int> tmp_vector_int;
891  if (curNumLocked + locktest->howMany() > maxLocked_) {
892  // just use the first of them
893  tmp_vector_int.insert(tmp_vector_int.begin(),locktest->whichVecs().begin(),locktest->whichVecs().begin()+maxLocked_-curNumLocked);
894  }
895  else {
896  tmp_vector_int = locktest->whichVecs();
897  }
898  const std::vector<int> lockind(tmp_vector_int);
899  const int numNewLocked = lockind.size();
900  //
901  // generate indices of vectors left unlocked
902  // curind = [0,...,curdim-1] = UNION( lockind, unlockind )
903  const int numUnlocked = curdim-numNewLocked;
904  tmp_vector_int.resize(curdim);
905  for (int i=0; i<curdim; i++) tmp_vector_int[i] = i;
906  const std::vector<int> curind(tmp_vector_int); // curind = [0 ... curdim-1]
907  tmp_vector_int.resize(numUnlocked);
908  std::set_difference(curind.begin(),curind.end(),lockind.begin(),lockind.end(),tmp_vector_int.begin());
909  const std::vector<int> unlockind(tmp_vector_int); // unlockind = [0 ... curdim-1] - lockind
910  tmp_vector_int.clear();
911 
912  //
913  // debug printing
914  if (printer_->isVerbosity(Debug)) {
915  printer_->print(Debug,"Locking vectors: ");
916  for (unsigned int i=0; i<lockind.size(); i++) {printer_->stream(Debug) << " " << lockind[i];}
917  printer_->print(Debug,"\n");
918  printer_->print(Debug,"Not locking vectors: ");
919  for (unsigned int i=0; i<unlockind.size(); i++) {printer_->stream(Debug) << " " << unlockind[i];}
920  printer_->print(Debug,"\n");
921  }
922 
923  //
924  // we need primitive ritz vectors/values:
925  // [S,L] = eig(oldKK)
926  //
927  // this will be partitioned as follows:
928  // locked: Sl = S(lockind) // we won't actually need Sl
929  // unlocked: Su = S(unlockind)
930  //
931  Teuchos::SerialDenseMatrix<int,ScalarType> S(curdim,curdim);
932  std::vector<MagnitudeType> theta(curdim);
933  {
934  int rank = curdim;
935  int info = msutils::directSolver(curdim,*state.KK,Teuchos::null,S,theta,rank,10);
936  TEUCHOS_TEST_FOR_EXCEPTION(info != 0 ,std::logic_error,
937  "Anasazi::BlockDavidsonSolMgr::solve(): error calling SolverUtils::directSolver."); // this should never happen
938  TEUCHOS_TEST_FOR_EXCEPTION(rank != curdim,std::logic_error,
939  "Anasazi::BlockDavidsonSolMgr::solve(): direct solve did not compute all eigenvectors."); // this should never happen
940  //
941  // sort the eigenvalues (so that we can order the eigenvectors)
942  std::vector<int> order(curdim);
943  sorter->sort(theta,Teuchos::rcpFromRef(order),curdim);
944  //
945  // apply the same ordering to the primitive ritz vectors
946  msutils::permuteVectors(order,S);
947  }
948  //
949  // select the unlocked ritz vectors
950  // the indexing in unlockind is relative to the ordered primitive ritz vectors
951  // (this is why we ordered theta,S above)
952  Teuchos::SerialDenseMatrix<int,ScalarType> Su(curdim,numUnlocked);
953  for (int i=0; i<numUnlocked; i++) {
954  blas.COPY(curdim, S[unlockind[i]], 1, Su[i], 1);
955  }
956 
957 
958  //
959  // newV has the following form:
960  // newV = [defV augV]
961  // - defV will be of size curdim - numNewLocked, and contain the generated basis: defV = oldV*Su
962  // - augV will be of size numNewLocked, and contain random directions to make up for the lost space
963  //
964  // we will need a pointer to defV below to generate the off-diagonal block of newKK
965  // go ahead and setup pointer to augV
966  //
967  Teuchos::RCP<MV> defV, augV;
968  if (inSituRestart_ == true) {
969  //
970  // get non-const pointer to solver's basis so we can work in situ
971  Teuchos::RCP<MV> solverbasis = Teuchos::rcp_const_cast<MV>(state.V);
972  //
973  // perform Householder QR of Su = Q [D;0], where D is unit diag.
974  // work on a copy of Su, since we need Su below to build newKK
975  Teuchos::SerialDenseMatrix<int,ScalarType> copySu(Su);
976  std::vector<ScalarType> tau(numUnlocked), work(numUnlocked);
977  int info;
978  lapack.GEQRF(curdim,numUnlocked,copySu.values(),copySu.stride(),&tau[0],&work[0],work.size(),&info);
979  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,std::logic_error,
980  "Anasazi::BlockDavidsonSolMgr::solve(): error calling GEQRF during restarting.");
981  if (printer_->isVerbosity(Debug)) {
982  Teuchos::SerialDenseMatrix<int,ScalarType> R(Teuchos::Copy,copySu,numUnlocked,numUnlocked);
983  for (int j=0; j<numUnlocked; j++) {
984  R(j,j) = SCT::magnitude(R(j,j)) - 1.0;
985  for (int i=j+1; i<numUnlocked; i++) {
986  R(i,j) = ZERO;
987  }
988  }
989  printer_->stream(Debug) << "||Triangular factor of Su - I||: " << R.normFrobenius() << std::endl;
990  }
991  //
992  // perform implicit oldV*Su
993  // this actually performs oldV*[Su Sl*M] = [defV lockV], for some unitary M
994  // we are actually interested in only the first numUnlocked vectors of the result
995  {
996  Teuchos::RCP<MV> oldV = MVT::CloneViewNonConst(*solverbasis,curind);
997  msutils::applyHouse(numUnlocked,*oldV,copySu,tau,workMV);
998  }
999  std::vector<int> defind(numUnlocked), augind(numNewLocked);
1000  for (int i=0; i<numUnlocked ; i++) defind[i] = i;
1001  for (int i=0; i<numNewLocked; i++) augind[i] = numUnlocked+i;
1002  defV = MVT::CloneViewNonConst(*solverbasis,defind);
1003  augV = MVT::CloneViewNonConst(*solverbasis,augind);
1004  }
1005  else { // inSituRestart == false)
1006  // defV = oldV*Su, explicitly. workspace is in workMV
1007  std::vector<int> defind(numUnlocked), augind(numNewLocked);
1008  for (int i=0; i<numUnlocked ; i++) defind[i] = i;
1009  for (int i=0; i<numNewLocked; i++) augind[i] = numUnlocked+i;
1010  Teuchos::RCP<const MV> oldV = MVT::CloneView(*state.V,curind);
1011  defV = MVT::CloneViewNonConst(*workMV,defind);
1012  augV = MVT::CloneViewNonConst(*workMV,augind);
1013 
1014  MVT::MvTimesMatAddMv(ONE,*oldV,Su,ZERO,*defV);
1015  }
1016 
1017  //
1018  // lockvecs will be partitioned as follows:
1019  // lockvecs = [curlocked augTmp ...]
1020  // - augTmp will be used for the storage of M*augV and K*augV
1021  // later, the locked vectors (stored in state.X and referenced via const MV view newLocked)
1022  // will be moved into lockvecs on top of augTmp when it is no longer needed as workspace.
1023  // - curlocked will be used in orthogonalization of augV
1024  //
1025  // newL is the new locked vectors; newL = oldV*Sl = RitzVectors(lockind)
1026  // we will not produce them, but instead retrieve them from RitzVectors
1027  //
1028  Teuchos::RCP<const MV> curlocked, newLocked;
1029  Teuchos::RCP<MV> augTmp;
1030  {
1031  // setup curlocked
1032  if (curNumLocked > 0) {
1033  std::vector<int> curlockind(curNumLocked);
1034  for (int i=0; i<curNumLocked; i++) curlockind[i] = i;
1035  curlocked = MVT::CloneView(*lockvecs,curlockind);
1036  }
1037  else {
1038  curlocked = Teuchos::null;
1039  }
1040  // setup augTmp
1041  std::vector<int> augtmpind(numNewLocked);
1042  for (int i=0; i<numNewLocked; i++) augtmpind[i] = curNumLocked+i;
1043  augTmp = MVT::CloneViewNonConst(*lockvecs,augtmpind);
1044  // setup newLocked
1045  newLocked = MVT::CloneView(*bd_solver->getRitzVectors(),lockind);
1046  }
1047 
1048  //
1049  // generate augV and perform orthogonalization
1050  //
1051  MVT::MvRandom(*augV);
1052  //
1053  // orthogonalize it against auxvecs, defV, and all locked vectors (new and current)
1054  // use augTmp as storage for M*augV, if hasM
1055  {
1056  Teuchos::Array<Teuchos::RCP<const MV> > against;
1057  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
1058  if (probauxvecs != Teuchos::null) against.push_back(probauxvecs);
1059  if (curlocked != Teuchos::null) against.push_back(curlocked);
1060  against.push_back(newLocked);
1061  against.push_back(defV);
1062  if (problem_->getM() != Teuchos::null) {
1063  OPT::Apply(*problem_->getM(),*augV,*augTmp);
1064  }
1065  ortho->projectAndNormalizeMat(*augV,against,dummyC,Teuchos::null,augTmp);
1066  }
1067 
1068  //
1069  // form newKK
1070  //
1071  // newKK = newV'*K*newV = [Su'*KK*Su defV'*K*augV]
1072  // [augV'*K*defV augV'*K*augV]
1073  //
1074  // first, generate the principal submatrix, the projection of K onto the unlocked portion of oldV
1075  //
1076  Teuchos::SerialDenseMatrix<int,ScalarType> newKK(curdim,curdim);
1077  {
1078  Teuchos::SerialDenseMatrix<int,ScalarType> KKtmp(curdim,numUnlocked),
1079  KKold(Teuchos::View,*state.KK,curdim,curdim),
1080  KK11(Teuchos::View,newKK,numUnlocked,numUnlocked);
1081  int teuchosRet;
1082  // KKtmp = KKold*Su
1083  teuchosRet = KKtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,KKold,Su,ZERO);
1084  TEUCHOS_TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
1085  "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
1086  // KK11 = Su'*KKtmp = Su'*KKold*Su
1087  teuchosRet = KK11.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,Su,KKtmp,ZERO);
1088  TEUCHOS_TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
1089  "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
1090  }
1091  //
1092  // project the stiffness matrix on augV
1093  {
1094  OPT::Apply(*problem_->getOperator(),*augV,*augTmp);
1095  Teuchos::SerialDenseMatrix<int,ScalarType> KK12(Teuchos::View,newKK,numUnlocked,numNewLocked,0,numUnlocked),
1096  KK22(Teuchos::View,newKK,numNewLocked,numNewLocked,numUnlocked,numUnlocked);
1097  MVT::MvTransMv(ONE,*defV,*augTmp,KK12);
1098  MVT::MvTransMv(ONE,*augV,*augTmp,KK22);
1099  }
1100  //
1101  // done with defV,augV
1102  defV = Teuchos::null;
1103  augV = Teuchos::null;
1104  //
1105  // make it hermitian in memory (fill in KK21)
1106  for (int j=0; j<curdim; ++j) {
1107  for (int i=j+1; i<curdim; ++i) {
1108  newKK(i,j) = SCT::conjugate(newKK(j,i));
1109  }
1110  }
1111  //
1112  // we are done using augTmp as storage
1113  // put newLocked into lockvecs, new values into lockvals
1114  augTmp = Teuchos::null;
1115  {
1116  std::vector<Value<ScalarType> > allvals = bd_solver->getRitzValues();
1117  for (int i=0; i<numNewLocked; i++) {
1118  lockvals.push_back(allvals[lockind[i]].realpart);
1119  }
1120 
1121  std::vector<int> indlock(numNewLocked);
1122  for (int i=0; i<numNewLocked; i++) indlock[i] = curNumLocked+i;
1123  MVT::SetBlock(*newLocked,indlock,*lockvecs);
1124  newLocked = Teuchos::null;
1125 
1126  curNumLocked += numNewLocked;
1127  std::vector<int> curlockind(curNumLocked);
1128  for (int i=0; i<curNumLocked; i++) curlockind[i] = i;
1129  curlocked = MVT::CloneView(*lockvecs,curlockind);
1130  }
1131  // add locked vecs as aux vecs, along with aux vecs from problem
1132  // add lockvals to ordertest
1133  // disable locktest if curNumLocked == maxLocked
1134  {
1135  ordertest->setAuxVals(lockvals);
1136 
1137  Teuchos::Array< Teuchos::RCP<const MV> > aux;
1138  if (probauxvecs != Teuchos::null) aux.push_back(probauxvecs);
1139  aux.push_back(curlocked);
1140  bd_solver->setAuxVecs(aux);
1141 
1142  if (curNumLocked == maxLocked_) {
1143  // disabled locking now
1144  combotest->removeTest(locktest);
1145  }
1146  }
1147 
1148  //
1149  // prepare new state
1151  rstate.curDim = curdim;
1152  if (inSituRestart_) {
1153  // data is already in the solver's memory
1154  rstate.V = state.V;
1155  }
1156  else {
1157  // data is in workspace and will be copied to solver memory
1158  rstate.V = workMV;
1159  }
1160  rstate.KK = Teuchos::rcpFromRef(newKK);
1161  //
1162  // pass new state to the solver
1163  bd_solver->initialize(rstate);
1164  } // end of locking
1166  //
1167  // we returned from iterate(), but none of our status tests Passed.
1168  // something is wrong, and it is probably our fault.
1169  //
1171  else {
1172  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): Invalid return from bd_solver::iterate().");
1173  }
1174  }
1175  catch (const AnasaziError &err) {
1176  printer_->stream(Errors)
1177  << "Anasazi::BlockDavidsonSolMgr::solve() caught unexpected exception from Anasazi::BlockDavidson::iterate() at iteration " << bd_solver->getNumIters() << std::endl
1178  << err.what() << std::endl
1179  << "Anasazi::BlockDavidsonSolMgr::solve() returning Unconverged with no solutions." << std::endl;
1180  return Unconverged;
1181  }
1182  }
1183 
1184  // clear temp space
1185  workMV = Teuchos::null;
1186 
1187  sol.numVecs = ordertest->howMany();
1188  if (sol.numVecs > 0) {
1189  sol.Evecs = MVT::Clone(*problem_->getInitVec(),sol.numVecs);
1190  sol.Espace = sol.Evecs;
1191  sol.Evals.resize(sol.numVecs);
1192  std::vector<MagnitudeType> vals(sol.numVecs);
1193 
1194  // copy them into the solution
1195  std::vector<int> which = ordertest->whichVecs();
1196  // indices between [0,blockSize) refer to vectors/values in the solver
1197  // indices between [-curNumLocked,-1] refer to locked vectors/values [0,curNumLocked)
1198  // everything has already been ordered by the solver; we just have to partition the two references
1199  std::vector<int> inlocked(0), insolver(0);
1200  for (unsigned int i=0; i<which.size(); i++) {
1201  if (which[i] >= 0) {
1202  TEUCHOS_TEST_FOR_EXCEPTION(which[i] >= blockSize_,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): positive indexing mistake from ordertest.");
1203  insolver.push_back(which[i]);
1204  }
1205  else {
1206  // sanity check
1207  TEUCHOS_TEST_FOR_EXCEPTION(which[i] < -curNumLocked,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): negative indexing mistake from ordertest.");
1208  inlocked.push_back(which[i] + curNumLocked);
1209  }
1210  }
1211 
1212  TEUCHOS_TEST_FOR_EXCEPTION(insolver.size() + inlocked.size() != (unsigned int)sol.numVecs,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): indexing mistake.");
1213 
1214  // set the vecs,vals in the solution
1215  if (insolver.size() > 0) {
1216  // set vecs
1217  int lclnum = insolver.size();
1218  std::vector<int> tosol(lclnum);
1219  for (int i=0; i<lclnum; i++) tosol[i] = i;
1220  Teuchos::RCP<const MV> v = MVT::CloneView(*bd_solver->getRitzVectors(),insolver);
1221  MVT::SetBlock(*v,tosol,*sol.Evecs);
1222  // set vals
1223  std::vector<Value<ScalarType> > fromsolver = bd_solver->getRitzValues();
1224  for (unsigned int i=0; i<insolver.size(); i++) {
1225  vals[i] = fromsolver[insolver[i]].realpart;
1226  }
1227  }
1228 
1229  // get the vecs,vals from locked storage
1230  if (inlocked.size() > 0) {
1231  int solnum = insolver.size();
1232  // set vecs
1233  int lclnum = inlocked.size();
1234  std::vector<int> tosol(lclnum);
1235  for (int i=0; i<lclnum; i++) tosol[i] = solnum + i;
1236  Teuchos::RCP<const MV> v = MVT::CloneView(*lockvecs,inlocked);
1237  MVT::SetBlock(*v,tosol,*sol.Evecs);
1238  // set vals
1239  for (unsigned int i=0; i<inlocked.size(); i++) {
1240  vals[i+solnum] = lockvals[inlocked[i]];
1241  }
1242  }
1243 
1244  // sort the eigenvalues and permute the eigenvectors appropriately
1245  {
1246  std::vector<int> order(sol.numVecs);
1247  sorter->sort(vals,Teuchos::rcpFromRef(order),sol.numVecs);
1248  // store the values in the Eigensolution
1249  for (int i=0; i<sol.numVecs; i++) {
1250  sol.Evals[i].realpart = vals[i];
1251  sol.Evals[i].imagpart = MT::zero();
1252  }
1253  // now permute the eigenvectors according to order
1254  msutils::permuteVectors(sol.numVecs,order,*sol.Evecs);
1255  }
1256 
1257  // setup sol.index, remembering that all eigenvalues are real so that index = {0,...,0}
1258  sol.index.resize(sol.numVecs,0);
1259  }
1260  }
1261 
1262  // print final summary
1263  bd_solver->currentStatus(printer_->stream(FinalSummary));
1264 
1265  // print timing information
1266 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1267  if ( printer_->isVerbosity( TimingDetails ) ) {
1268  Teuchos::TimeMonitor::summarize( printer_->stream( TimingDetails ) );
1269  }
1270 #endif
1271 
1272  problem_->setSolution(sol);
1273  printer_->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
1274 
1275  // get the number of iterations taken for this call to solve().
1276  numIters_ = bd_solver->getNumIters();
1277 
1278  if (sol.numVecs < nev) {
1279  return Unconverged; // return from BlockDavidsonSolMgr::solve()
1280  }
1281  return Converged; // return from BlockDavidsonSolMgr::solve()
1282 }
1283 
1284 
1285 template <class ScalarType, class MV, class OP>
1286 void
1288  const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global)
1289 {
1290  globalTest_ = global;
1291 }
1292 
1293 template <class ScalarType, class MV, class OP>
1294 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
1296 {
1297  return globalTest_;
1298 }
1299 
1300 template <class ScalarType, class MV, class OP>
1301 void
1303  const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &debug)
1304 {
1305  debugTest_ = debug;
1306 }
1307 
1308 template <class ScalarType, class MV, class OP>
1309 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
1311 {
1312  return debugTest_;
1313 }
1314 
1315 template <class ScalarType, class MV, class OP>
1316 void
1318  const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &locking)
1319 {
1320  lockingTest_ = locking;
1321 }
1322 
1323 template <class ScalarType, class MV, class OP>
1324 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
1326 {
1327  return lockingTest_;
1328 }
1329 
1330 } // end Anasazi namespace
1331 
1332 #endif /* ANASAZI_BLOCKDAVIDSON_SOLMGR_HPP */
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
Pure virtual base class which describes the basic interface for a solver manager. ...
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
void setLockingStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &locking)
Set the status test defining locking.
ResType
Enumerated type used to specify which residual norm used by residual norm status tests.
A special StatusTest for printing other status tests.
BlockDavidsonSolMgr(const Teuchos::RCP< Eigenproblem< ScalarType, MV, OP > > &problem, Teuchos::ParameterList &pl)
Basic constructor for BlockDavidsonSolMgr.
This class defines the interface required by an eigensolver and status test class to compute solution...
Teuchos::RCP< const MV > R
The current residual vectors.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > KK
The current projected K matrix.
An implementation of the Anasazi::SortManager that performs a collection of common sorting techniques...
Virtual base class which defines basic traits for the operator type.
Teuchos::RCP< MV > Evecs
The computed eigenvectors.
Status test for forming logical combinations of other status tests.
Teuchos::RCP< const MV > KX
The image of the current eigenvectors under K.
The Anasazi::SolverManager is a templated virtual base class that defines the basic interface that an...
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getDebugStatusTest() const
Get the status test for debugging.
static void Apply(const OP &Op, const MV &x, MV &y)
Application method which performs operation y = Op*x. An OperatorError exception is thrown if there i...
Basic implementation of the Anasazi::SortManager class.
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using the SVQB iter...
An exception class parent to all Anasazi exceptions.
This class implements a Block Davidson iteration, a preconditioned iteration for solving linear Hermi...
Implementation of the block Davidson method.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
Anasazi&#39;s templated, static class providing utilities for the solvers.
int numVecs
The number of computed eigenpairs.
int getNumIters() const
Get the iteration count for the most recent call to solve().
Basic output manager for sending information of select verbosity levels to the appropriate output str...
Teuchos::RCP< MV > Espace
An orthonormal basis for the computed eigenspace.
Anasazi&#39;s basic output manager for sending information of select verbosity levels to the appropriate ...
Abstract base class which defines the interface required by an eigensolver and status test class to c...
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
ReturnType
Enumerated type used to pass back information from a solver manager.
Teuchos::RCP< const MV > V
The basis for the Krylov space.
A status test for testing the norm of the eigenvectors residuals.
Traits class which defines basic operations on multivectors.
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &B, Teuchos::SerialDenseMatrix< int, ScalarType > &C)
Compute C := alpha * A^H B.
std::vector< int > index
An index into Evecs to allow compressed storage of eigenvectors for real, non-Hermitian problems...
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getLockingStatusTest() const
Get the status test defining locking.
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
Orthogonalization manager based on the SVQB technique described in "A Block Orthogonalization Procedu...
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
Struct for storing an eigenproblem solution.
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index...
int curDim
The current dimension of the solver.
Structure to contain pointers to BlockDavidson state variables.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getGlobalStatusTest() const
Get the status test defining global convergence.
Special StatusTest for printing status tests.
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
void setDebugStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &debug)
Set the status test for debugging.
Teuchos::RCP< const MV > X
The current eigenvectors.
Teuchos::RCP< const MV > MX
The image of the current eigenvectors under M, or Teuchos::null if M was not specified.
Status test for forming logical combinations of other status tests.
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using (potentially)...
Types and exceptions used within Anasazi solvers and interfaces.
The BlockDavidsonSolMgr provides a powerful solver manager over the BlockDavidson eigensolver...
Common interface of stopping criteria for Anasazi&#39;s solvers.
A status test for testing the norm of the eigenvectors residuals.
Basic implementation of the Anasazi::OrthoManager class.
void setGlobalStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &global)
Set the status test defining global convergence.
ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver&#39;s iterate() routine until ...
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
Teuchos::RCP< const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > T
The current Ritz values. This vector is a copy of the internal data.
Class which provides internal utilities for the Anasazi solvers.