Anasazi  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #ifndef ANASAZI_BLOCKDAVIDSON_SOLMGR_HPP
43 #define ANASAZI_BLOCKDAVIDSON_SOLMGR_HPP
44 
49 #include "AnasaziConfigDefs.hpp"
50 #include "AnasaziTypes.hpp"
51 
52 #include "AnasaziEigenproblem.hpp"
53 #include "AnasaziSolverManager.hpp"
54 #include "AnasaziSolverUtils.hpp"
55 
56 #include "AnasaziBlockDavidson.hpp"
57 #include "AnasaziBasicSort.hpp"
64 #include "AnasaziOutputManager.hpp"
66 #include "Teuchos_BLAS.hpp"
67 #include "Teuchos_LAPACK.hpp"
68 #include "Teuchos_TimeMonitor.hpp"
69 #include "Teuchos_FancyOStream.hpp"
70 
71 
85 namespace Anasazi {
86 
87 
120 template<class ScalarType, class MV, class OP>
121 class BlockDavidsonSolMgr : public SolverManager<ScalarType,MV,OP> {
122 
123  private:
127  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
129 
130  public:
131 
133 
134 
162 
164  virtual ~BlockDavidsonSolMgr() {};
166 
168 
169 
172  return *problem_;
173  }
174 
176  int getNumIters() const {
177  return numIters_;
178  }
179 
188  return Teuchos::tuple(_timerSolve, _timerRestarting, _timerLocking);
189  }
190 
192 
194 
195 
217  ReturnType solve();
218 
221 
224 
227 
230 
233 
236 
238 
239  private:
241 
242  std::string whch_, ortho_;
243 
244  MagnitudeType convtol_, locktol_;
245  int maxRestarts_;
246  bool useLocking_;
247  bool relconvtol_, rellocktol_;
248  int blockSize_, numBlocks_, numIters_;
249  int maxLocked_;
250  int lockQuorum_;
251  bool inSituRestart_;
252  int numRestartBlocks_;
253  enum ResType convNorm_, lockNorm_;
254 
255  Teuchos::RCP<Teuchos::Time> _timerSolve, _timerRestarting, _timerLocking;
256 
260 
262 };
263 
264 
265 // Constructor
266 template<class ScalarType, class MV, class OP>
269  Teuchos::ParameterList &pl ) :
270  problem_(problem),
271  whch_("SR"),
272  ortho_("SVQB"),
273  convtol_(MT::prec()),
274  maxRestarts_(20),
275  useLocking_(false),
276  relconvtol_(true),
277  rellocktol_(true),
278  blockSize_(0),
279  numBlocks_(0),
280  numIters_(0),
281  maxLocked_(0),
282  lockQuorum_(1),
283  inSituRestart_(false),
284  numRestartBlocks_(1)
285 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
286  , _timerSolve(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidsonSolMgr::solve()")),
287  _timerRestarting(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidsonSolMgr restarting")),
288  _timerLocking(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidsonSolMgr locking"))
289 #endif
290 {
291  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
292  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument, "Problem not set.");
293  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isHermitian(), std::invalid_argument, "Problem not symmetric.");
294  TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument, "Problem does not contain initial vectors to clone from.");
295 
296  std::string strtmp;
297 
298  // which values to solve for
299  whch_ = pl.get("Which",whch_);
300  TEUCHOS_TEST_FOR_EXCEPTION(whch_ != "SM" && whch_ != "LM" && whch_ != "SR" && whch_ != "LR",std::invalid_argument, "Invalid sorting string.");
301 
302  // which orthogonalization to use
303  ortho_ = pl.get("Orthogonalization",ortho_);
304  if (ortho_ != "DGKS" && ortho_ != "SVQB") {
305  ortho_ = "SVQB";
306  }
307 
308  // convergence tolerance
309  convtol_ = pl.get("Convergence Tolerance",convtol_);
310  relconvtol_ = pl.get("Relative Convergence Tolerance",relconvtol_);
311  strtmp = pl.get("Convergence Norm",std::string("2"));
312  if (strtmp == "2") {
313  convNorm_ = RES_2NORM;
314  }
315  else if (strtmp == "M") {
316  convNorm_ = RES_ORTH;
317  }
318  else {
319  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
320  "Anasazi::BlockDavidsonSolMgr: Invalid Convergence Norm.");
321  }
322 
323  // locking tolerance
324  useLocking_ = pl.get("Use Locking",useLocking_);
325  rellocktol_ = pl.get("Relative Locking Tolerance",rellocktol_);
326  // default: should be less than convtol_
327  locktol_ = convtol_/10;
328  locktol_ = pl.get("Locking Tolerance",locktol_);
329  strtmp = pl.get("Locking Norm",std::string("2"));
330  if (strtmp == "2") {
331  lockNorm_ = RES_2NORM;
332  }
333  else if (strtmp == "M") {
334  lockNorm_ = RES_ORTH;
335  }
336  else {
337  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
338  "Anasazi::BlockDavidsonSolMgr: Invalid Locking Norm.");
339  }
340 
341  // maximum number of restarts
342  maxRestarts_ = pl.get("Maximum Restarts",maxRestarts_);
343 
344  // block size: default is nev()
345  blockSize_ = pl.get("Block Size",problem_->getNEV());
346  TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
347  "Anasazi::BlockDavidsonSolMgr: \"Block Size\" must be strictly positive.");
348  numBlocks_ = pl.get("Num Blocks",2);
349  TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 1, std::invalid_argument,
350  "Anasazi::BlockDavidsonSolMgr: \"Num Blocks\" must be >= 1.");
351 
352  // max locked: default is nev(), must satisfy maxLocked_ + blockSize_ >= nev
353  if (useLocking_) {
354  maxLocked_ = pl.get("Max Locked",problem_->getNEV());
355  }
356  else {
357  maxLocked_ = 0;
358  }
359  if (maxLocked_ == 0) {
360  useLocking_ = false;
361  }
362  TEUCHOS_TEST_FOR_EXCEPTION(maxLocked_ < 0, std::invalid_argument,
363  "Anasazi::BlockDavidsonSolMgr: \"Max Locked\" must be positive.");
364  TEUCHOS_TEST_FOR_EXCEPTION(maxLocked_ + blockSize_ < problem_->getNEV(),
365  std::invalid_argument,
366  "Anasazi::BlockDavidsonSolMgr: Not enough storage space for requested number of eigenpairs.");
367  TEUCHOS_TEST_FOR_EXCEPTION(static_cast<ptrdiff_t>(numBlocks_)*blockSize_ + maxLocked_ > MVT::GetGlobalLength(*problem_->getInitVec()),
368  std::invalid_argument,
369  "Anasazi::BlockDavidsonSolMgr: Potentially impossible orthogonality requests. Reduce basis size or locking size.");
370 
371  if (useLocking_) {
372  lockQuorum_ = pl.get("Locking Quorum",lockQuorum_);
373  TEUCHOS_TEST_FOR_EXCEPTION(lockQuorum_ <= 0,
374  std::invalid_argument,
375  "Anasazi::BlockDavidsonSolMgr: \"Locking Quorum\" must be strictly positive.");
376  }
377 
378  // restart size
379  numRestartBlocks_ = pl.get("Num Restart Blocks",numRestartBlocks_);
380  TEUCHOS_TEST_FOR_EXCEPTION(numRestartBlocks_ <= 0, std::invalid_argument,
381  "Anasazi::BlockDavidsonSolMgr: \"Num Restart Blocks\" must be strictly positive.");
382  TEUCHOS_TEST_FOR_EXCEPTION(numRestartBlocks_ >= numBlocks_, std::invalid_argument,
383  "Anasazi::BlockDavidsonSolMgr: \"Num Restart Blocks\" must be strictly less than \"Num Blocks\".");
384 
385  // restarting technique: V*Q or applyHouse(V,H,tau)
386  if (pl.isParameter("In Situ Restarting")) {
387  if (Teuchos::isParameterType<bool>(pl,"In Situ Restarting")) {
388  inSituRestart_ = pl.get("In Situ Restarting",inSituRestart_);
389  } else {
390  inSituRestart_ = ( Teuchos::getParameter<int>(pl,"In Situ Restarting") != 0 );
391  }
392  }
393 
394  // Output manager
395  // Create a formatted output stream to print to.
396  // See if user requests output processor.
397  int osProc = pl.get("Output Processor", 0);
398 
399  // If not passed in by user, it will be chosen based upon operator type.
401 
402  // output stream
403  if (pl.isParameter("Output Stream")) {
404  osp = Teuchos::getParameter<Teuchos::RCP<Teuchos::FancyOStream> >(pl,"Output Stream");
405  }
406  else {
407  osp = OutputStreamTraits<OP>::getOutputStream (*problem_->getOperator(), osProc);
408  }
409 
410  // verbosity
411  int verbosity = Anasazi::Errors;
412  if (pl.isParameter("Verbosity")) {
413  if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
414  verbosity = pl.get("Verbosity", verbosity);
415  } else {
416  verbosity = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
417  }
418  }
419  printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity,osp) );
420 
421 }
422 
423 
424 // solve()
425 template<class ScalarType, class MV, class OP>
426 ReturnType
428 
429  typedef SolverUtils<ScalarType,MV,OP> msutils;
430 
431  const int nev = problem_->getNEV();
432 
433 #ifdef TEUCHOS_DEBUG
435  out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(printer_->stream(Debug)));
436  out->setShowAllFrontMatter(false).setShowProcRank(true);
437  *out << "Entering Anasazi::BlockDavidsonSolMgr::solve()\n";
438 #endif
439 
441  // Sort manager
443 
445  // Status tests
446  //
447  // convergence
449  if (globalTest_ == Teuchos::null) {
450  convtest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(convtol_,nev,convNorm_,relconvtol_) );
451  }
452  else {
453  convtest = globalTest_;
454  }
456  = Teuchos::rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sorter,nev) );
457  // locking
459  if (useLocking_) {
460  if (lockingTest_ == Teuchos::null) {
461  locktest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(locktol_,lockQuorum_,lockNorm_,rellocktol_) );
462  }
463  else {
464  locktest = lockingTest_;
465  }
466  }
467  // for a non-short-circuited OR test, the order doesn't matter
469  alltests.push_back(ordertest);
470  if (locktest != Teuchos::null) alltests.push_back(locktest);
471  if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
472 
475  // printing StatusTest
477  if ( printer_->isVerbosity(Debug) ) {
478  outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,Passed+Failed+Undefined ) );
479  }
480  else {
481  outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,Passed ) );
482  }
483 
485  // Orthomanager
487  if (ortho_=="SVQB") {
488  ortho = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
489  } else if (ortho_=="DGKS") {
490  ortho = Teuchos::rcp( new BasicOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
491  } else {
492  TEUCHOS_TEST_FOR_EXCEPTION(ortho_!="SVQB"&&ortho_!="DGKS",std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): Invalid orthogonalization type.");
493  }
494 
496  // Parameter list
498  plist.set("Block Size",blockSize_);
499  plist.set("Num Blocks",numBlocks_);
500 
502  // BlockDavidson solver
504  = Teuchos::rcp( new BlockDavidson<ScalarType,MV,OP>(problem_,sorter,printer_,outputtest,ortho,plist) );
505  // set any auxiliary vectors defined in the problem
506  Teuchos::RCP< const MV > probauxvecs = problem_->getAuxVecs();
507  if (probauxvecs != Teuchos::null) {
508  bd_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
509  }
510 
512  // Storage
513  //
514  // lockvecs will contain eigenvectors that have been determined "locked" by the status test
515  int curNumLocked = 0;
516  Teuchos::RCP<MV> lockvecs;
517  // lockvecs is used to hold the locked eigenvectors, as well as for temporary storage when locking.
518  // when locking, we will lock some number of vectors numnew, where numnew <= maxlocked - curlocked
519  // we will produce numnew random vectors, which will go into the space with the new basis.
520  // we will also need numnew storage for the image of these random vectors under A and M;
521  // columns [curlocked+1,curlocked+numnew] will be used for this storage
522  if (maxLocked_ > 0) {
523  lockvecs = MVT::Clone(*problem_->getInitVec(),maxLocked_);
524  }
525  std::vector<MagnitudeType> lockvals;
526  //
527  // Restarting occurs under two scenarios: when the basis is full and after locking.
528  //
529  // For the former, a new basis of size blockSize*numRestartBlocks is generated using the current basis
530  // and the most significant primitive Ritz vectors (projected eigenvectors).
531  // [S,L] = eig(KK)
532  // S = [Sr St] // some for "r"estarting, some are "t"runcated
533  // newV = V*Sr
534  // KK_new = newV'*K*newV = Sr'*V'*K*V*Sr = Sr'*KK*Sr
535  // Therefore, the only multivector operation needed is for the generation of newV.
536  //
537  // * If the multiplication is explicit, it requires a workspace of blockSize*numRestartBlocks vectors.
538  // This space must be specifically allocated for that task, as we don't have any space of that size.
539  // It (workMV) will be allocated at the beginning of solve()
540  // * Optionally, the multiplication can be performed implicitly, via a Householder QR factorization of
541  // Sr. This can be done in situ, using the basis multivector contained in the solver. This requires
542  // that we cast away the const on the multivector returned from getState(). Workspace for this approach
543  // is a single vector. the solver's internal storage must be preserved (X,MX,KX,R), requiring us to
544  // allocate this vector.
545  //
546  // For the latter (restarting after locking), the new basis is the same size as existing basis. If numnew
547  // vectors are locked, they are deflated from the current basis and replaced with randomly generated
548  // vectors.
549  // [S,L] = eig(KK)
550  // S = [Sl Su] // partitioned: "l"ocked and "u"nlocked
551  // newL = V*Sl = X(locked)
552  // defV = V*Su
553  // augV = rand(numnew) // orthogonal to oldL,newL,defV,auxvecs
554  // newV = [defV augV]
555  // Kknew = newV'*K*newV = [Su'*KK*Su defV'*K*augV]
556  // [augV'*K*defV augV'*K*augV]
557  // locked = [oldL newL]
558  // Clearly, this operation is more complicated than the previous.
559  // Here is a list of the significant computations that need to be performed:
560  // - newL will be put into space in lockvecs, but will be copied from getState().X at the end
561  // - defV,augV will be stored in workspace the size of the current basis.
562  // - If inSituRestart==true, we compute defV in situ in bd_solver::V_ and
563  // put augV at the end of bd_solver::V_
564  // - If inSituRestart==false, we must have curDim vectors available for
565  // defV and augV; we will allocate a multivector (workMV) at the beginning of solve()
566  // for this purpose.
567  // - M*augV and K*augV are needed; they will be stored in lockvecs. As a result, newL will
568  // not be put into lockvecs until the end.
569  //
570  // Therefore, we must allocate workMV when ((maxRestarts_ > 0) || (useLocking_ == true)) && inSituRestart == false
571  // It will be allocated to size (numBlocks-1)*blockSize
572  //
573  Teuchos::RCP<MV> workMV;
574  if (inSituRestart_ == false) {
575  // we need storage space to restart, either if we may lock or if may restart after a full basis
576  if (useLocking_==true || maxRestarts_ > 0) {
577  workMV = MVT::Clone(*problem_->getInitVec(),(numBlocks_-1)*blockSize_);
578  }
579  else {
580  // we will never need to restart.
581  workMV = Teuchos::null;
582  }
583  }
584  else { // inSituRestart_ == true
585  // we will restart in situ, if we need to restart
586  // three situation remain:
587  // - never restart => no space needed
588  // - only restart for locking (i.e., never restart full) => no space needed
589  // - restart for full basis => need one vector
590  if (maxRestarts_ > 0) {
591  workMV = MVT::Clone(*problem_->getInitVec(),1);
592  }
593  else {
594  workMV = Teuchos::null;
595  }
596  }
597 
598  // some consts and utils
599  const ScalarType ONE = SCT::one();
600  const ScalarType ZERO = SCT::zero();
603 
604  // go ahead and initialize the solution to nothing in case we throw an exception
606  sol.numVecs = 0;
607  problem_->setSolution(sol);
608 
609  int numRestarts = 0;
610 
611  // enter solve() iterations
612  {
613 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
614  Teuchos::TimeMonitor slvtimer(*_timerSolve);
615 #endif
616 
617  // tell bd_solver to iterate
618  while (1) {
619  try {
620  bd_solver->iterate();
621 
623  //
624  // check user-specified debug test; if it passed, return an exception
625  //
627  if (debugTest_ != Teuchos::null && debugTest_->getStatus() == Passed) {
628  throw AnasaziError("Anasazi::BlockDavidsonSolMgr::solve(): User-specified debug status test returned Passed.");
629  }
631  //
632  // check convergence next
633  //
635  else if (ordertest->getStatus() == Passed ) {
636  // we have convergence
637  // ordertest->whichVecs() tells us which vectors from lockvecs and solver state are the ones we want
638  // ordertest->howMany() will tell us how many
639  break;
640  }
642  //
643  // check for restarting before locking: if we need to lock, it will happen after the restart
644  //
646  else if ( bd_solver->getCurSubspaceDim() == bd_solver->getMaxSubspaceDim() ) {
647 
648 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
649  Teuchos::TimeMonitor restimer(*_timerRestarting);
650 #endif
651 
652  if ( numRestarts >= maxRestarts_ ) {
653  break; // break from while(1){bd_solver->iterate()}
654  }
655  numRestarts++;
656 
657  printer_->stream(IterationDetails) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << std::endl << std::endl;
658 
659  BlockDavidsonState<ScalarType,MV> state = bd_solver->getState();
660  int curdim = state.curDim;
661  int newdim = numRestartBlocks_*blockSize_;
662 
663 # ifdef TEUCHOS_DEBUG
664  {
665  std::vector<Value<ScalarType> > ritzvalues = bd_solver->getRitzValues();
666  *out << "Ritz values from solver:\n";
667  for (unsigned int i=0; i<ritzvalues.size(); i++) {
668  *out << ritzvalues[i].realpart << " ";
669  }
670  *out << "\n";
671  }
672 # endif
673 
674  //
675  // compute eigenvectors of the projected problem
677  std::vector<MagnitudeType> theta(curdim);
678  int rank = curdim;
679 # ifdef TEUCHOS_DEBUG
680  {
681  std::stringstream os;
682  os << "KK before HEEV...\n"
683  << printMat(*state.KK) << "\n";
684  *out << os.str();
685  }
686 # endif
687  int info = msutils::directSolver(curdim,*state.KK,Teuchos::null,S,theta,rank,10);
688  TEUCHOS_TEST_FOR_EXCEPTION(info != 0 ,std::logic_error,
689  "Anasazi::BlockDavidsonSolMgr::solve(): error calling SolverUtils::directSolver."); // this should never happen
690  TEUCHOS_TEST_FOR_EXCEPTION(rank != curdim,std::logic_error,
691  "Anasazi::BlockDavidsonSolMgr::solve(): direct solve did not compute all eigenvectors."); // this should never happen
692 
693 # ifdef TEUCHOS_DEBUG
694  {
695  std::stringstream os;
696  *out << "Ritz values from heev(KK):\n";
697  for (unsigned int i=0; i<theta.size(); i++) *out << theta[i] << " ";
698  os << "\nRitz vectors from heev(KK):\n"
699  << printMat(S) << "\n";
700  *out << os.str();
701  }
702 # endif
703 
704  //
705  // sort the eigenvalues (so that we can order the eigenvectors)
706  {
707  std::vector<int> order(curdim);
708  sorter->sort(theta,Teuchos::rcpFromRef(order),curdim);
709  //
710  // apply the same ordering to the primitive ritz vectors
711  msutils::permuteVectors(order,S);
712  }
713 # ifdef TEUCHOS_DEBUG
714  {
715  std::stringstream os;
716  *out << "Ritz values from heev(KK) after sorting:\n";
717  std::copy(theta.begin(), theta.end(), std::ostream_iterator<ScalarType>(*out, " "));
718  os << "\nRitz vectors from heev(KK) after sorting:\n"
719  << printMat(S) << "\n";
720  *out << os.str();
721  }
722 # endif
723  //
724  // select the significant primitive ritz vectors
726 # ifdef TEUCHOS_DEBUG
727  {
728  std::stringstream os;
729  os << "Significant primitive Ritz vectors:\n"
730  << printMat(Sr) << "\n";
731  *out << os.str();
732  }
733 # endif
734  //
735  // generate newKK = Sr'*KKold*Sr
736  Teuchos::SerialDenseMatrix<int,ScalarType> newKK(newdim,newdim);
737  {
738  Teuchos::SerialDenseMatrix<int,ScalarType> KKtmp(curdim,newdim),
739  KKold(Teuchos::View,*state.KK,curdim,curdim);
740  int teuchosRet;
741  // KKtmp = KKold*Sr
742  teuchosRet = KKtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,KKold,Sr,ZERO);
743  TEUCHOS_TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
744  "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
745  // newKK = Sr'*KKtmp = Sr'*KKold*Sr
746  teuchosRet = newKK.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,Sr,KKtmp,ZERO);
747  TEUCHOS_TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
748  "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
749  // make it Hermitian in memory
750  for (int j=0; j<newdim-1; ++j) {
751  for (int i=j+1; i<newdim; ++i) {
752  newKK(i,j) = SCT::conjugate(newKK(j,i));
753  }
754  }
755  }
756 # ifdef TEUCHOS_DEBUG
757  {
758  std::stringstream os;
759  os << "Sr'*KK*Sr:\n"
760  << printMat(newKK) << "\n";
761  *out << os.str();
762  }
763 # endif
764 
765  // prepare new state
767  rstate.curDim = newdim;
768  rstate.KK = Teuchos::rcpFromRef(newKK);
769  //
770  // we know that newX = newV*Sr(:,1:bS) = oldV*S(:1:bS) = oldX
771  // the restarting preserves the Ritz vectors and residual
772  // for the Ritz values, we want all of the values associated with newV.
773  // these have already been placed at the beginning of theta
774  rstate.X = state.X;
775  rstate.KX = state.KX;
776  rstate.MX = state.MX;
777  rstate.R = state.R;
778  rstate.T = Teuchos::rcp( new std::vector<MagnitudeType>(theta.begin(),theta.begin()+newdim) );
779 
780  if (inSituRestart_ == true) {
781 # ifdef TEUCHOS_DEBUG
782  *out << "Beginning in-situ restart...\n";
783 # endif
784  //
785  // get non-const pointer to solver's basis so we can work in situ
786  Teuchos::RCP<MV> solverbasis = Teuchos::rcp_const_cast<MV>(state.V);
787  //
788  // perform Householder QR of Sr = Q [D;0], where D is unit diag.
789  // WARNING: this will overwrite Sr; however, we do not need Sr anymore after this
790  std::vector<ScalarType> tau(newdim), work(newdim);
791  int geqrf_info;
792  lapack.GEQRF(curdim,newdim,Sr.values(),Sr.stride(),&tau[0],&work[0],work.size(),&geqrf_info);
793  TEUCHOS_TEST_FOR_EXCEPTION(geqrf_info != 0,std::logic_error,
794  "Anasazi::BlockDavidsonSolMgr::solve(): error calling GEQRF during restarting.");
795  if (printer_->isVerbosity(Debug)) {
797  for (int j=0; j<newdim; j++) {
798  R(j,j) = SCT::magnitude(R(j,j)) - 1.0;
799  for (int i=j+1; i<newdim; i++) {
800  R(i,j) = ZERO;
801  }
802  }
803  printer_->stream(Debug) << "||Triangular factor of Sr - I||: " << R.normFrobenius() << std::endl;
804  }
805  //
806  // perform implicit oldV*Sr
807  // this actually performs oldV*[Sr Su*M] = [newV truncV], for some unitary M
808  // we are actually interested in only the first newdim vectors of the result
809  {
810  std::vector<int> curind(curdim);
811  for (int i=0; i<curdim; i++) curind[i] = i;
812  Teuchos::RCP<MV> oldV = MVT::CloneViewNonConst(*solverbasis,curind);
813  msutils::applyHouse(newdim,*oldV,Sr,tau,workMV);
814  }
815  //
816  // put the new basis into the state for initialize()
817  // the new basis is contained in the the first newdim columns of solverbasis
818  // initialize() will recognize that pointer bd_solver.V_ == pointer rstate.V, and will neglect the copy.
819  rstate.V = solverbasis;
820  }
821  else { // inSituRestart == false)
822 # ifdef TEUCHOS_DEBUG
823  *out << "Beginning ex-situ restart...\n";
824 # endif
825  // newV = oldV*Sr, explicitly. workspace is in workMV
826  std::vector<int> curind(curdim), newind(newdim);
827  for (int i=0; i<curdim; i++) curind[i] = i;
828  for (int i=0; i<newdim; i++) newind[i] = i;
829  Teuchos::RCP<const MV> oldV = MVT::CloneView(*state.V,curind);
830  Teuchos::RCP<MV> newV = MVT::CloneViewNonConst(*workMV ,newind);
831 
832  MVT::MvTimesMatAddMv(ONE,*oldV,Sr,ZERO,*newV);
833  //
834  // put the new basis into the state for initialize()
835  rstate.V = newV;
836  }
837 
838  //
839  // send the new state to the solver
840  bd_solver->initialize(rstate);
841  } // end of restarting
843  //
844  // check locking if we didn't converge or restart
845  //
847  else if (locktest != Teuchos::null && locktest->getStatus() == Passed) {
848 
849 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
850  Teuchos::TimeMonitor lcktimer(*_timerLocking);
851 #endif
852 
853  //
854  // get current state
855  BlockDavidsonState<ScalarType,MV> state = bd_solver->getState();
856  const int curdim = state.curDim;
857 
858  //
859  // get number,indices of vectors to be locked
860  TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() <= 0,std::logic_error,
861  "Anasazi::BlockDavidsonSolMgr::solve(): status test mistake: howMany() non-positive.");
862  TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() != (int)locktest->whichVecs().size(),std::logic_error,
863  "Anasazi::BlockDavidsonSolMgr::solve(): status test mistake: howMany() not consistent with whichVecs().");
864  // we should have room to lock vectors; otherwise, locking should have been deactivated
865  TEUCHOS_TEST_FOR_EXCEPTION(curNumLocked == maxLocked_,std::logic_error,
866  "Anasazi::BlockDavidsonSolMgr::solve(): status test mistake: locking not deactivated.");
867  //
868  // don't lock more than maxLocked_; we didn't allocate enough space.
869  std::vector<int> tmp_vector_int;
870  if (curNumLocked + locktest->howMany() > maxLocked_) {
871  // just use the first of them
872  tmp_vector_int.insert(tmp_vector_int.begin(),locktest->whichVecs().begin(),locktest->whichVecs().begin()+maxLocked_-curNumLocked);
873  }
874  else {
875  tmp_vector_int = locktest->whichVecs();
876  }
877  const std::vector<int> lockind(tmp_vector_int);
878  const int numNewLocked = lockind.size();
879  //
880  // generate indices of vectors left unlocked
881  // curind = [0,...,curdim-1] = UNION( lockind, unlockind )
882  const int numUnlocked = curdim-numNewLocked;
883  tmp_vector_int.resize(curdim);
884  for (int i=0; i<curdim; i++) tmp_vector_int[i] = i;
885  const std::vector<int> curind(tmp_vector_int); // curind = [0 ... curdim-1]
886  tmp_vector_int.resize(numUnlocked);
887  std::set_difference(curind.begin(),curind.end(),lockind.begin(),lockind.end(),tmp_vector_int.begin());
888  const std::vector<int> unlockind(tmp_vector_int); // unlockind = [0 ... curdim-1] - lockind
889  tmp_vector_int.clear();
890 
891  //
892  // debug printing
893  if (printer_->isVerbosity(Debug)) {
894  printer_->print(Debug,"Locking vectors: ");
895  for (unsigned int i=0; i<lockind.size(); i++) {printer_->stream(Debug) << " " << lockind[i];}
896  printer_->print(Debug,"\n");
897  printer_->print(Debug,"Not locking vectors: ");
898  for (unsigned int i=0; i<unlockind.size(); i++) {printer_->stream(Debug) << " " << unlockind[i];}
899  printer_->print(Debug,"\n");
900  }
901 
902  //
903  // we need primitive ritz vectors/values:
904  // [S,L] = eig(oldKK)
905  //
906  // this will be partitioned as follows:
907  // locked: Sl = S(lockind) // we won't actually need Sl
908  // unlocked: Su = S(unlockind)
909  //
911  std::vector<MagnitudeType> theta(curdim);
912  {
913  int rank = curdim;
914  int info = msutils::directSolver(curdim,*state.KK,Teuchos::null,S,theta,rank,10);
915  TEUCHOS_TEST_FOR_EXCEPTION(info != 0 ,std::logic_error,
916  "Anasazi::BlockDavidsonSolMgr::solve(): error calling SolverUtils::directSolver."); // this should never happen
917  TEUCHOS_TEST_FOR_EXCEPTION(rank != curdim,std::logic_error,
918  "Anasazi::BlockDavidsonSolMgr::solve(): direct solve did not compute all eigenvectors."); // this should never happen
919  //
920  // sort the eigenvalues (so that we can order the eigenvectors)
921  std::vector<int> order(curdim);
922  sorter->sort(theta,Teuchos::rcpFromRef(order),curdim);
923  //
924  // apply the same ordering to the primitive ritz vectors
925  msutils::permuteVectors(order,S);
926  }
927  //
928  // select the unlocked ritz vectors
929  // the indexing in unlockind is relative to the ordered primitive ritz vectors
930  // (this is why we ordered theta,S above)
931  Teuchos::SerialDenseMatrix<int,ScalarType> Su(curdim,numUnlocked);
932  for (int i=0; i<numUnlocked; i++) {
933  blas.COPY(curdim, S[unlockind[i]], 1, Su[i], 1);
934  }
935 
936 
937  //
938  // newV has the following form:
939  // newV = [defV augV]
940  // - defV will be of size curdim - numNewLocked, and contain the generated basis: defV = oldV*Su
941  // - augV will be of size numNewLocked, and contain random directions to make up for the lost space
942  //
943  // we will need a pointer to defV below to generate the off-diagonal block of newKK
944  // go ahead and setup pointer to augV
945  //
946  Teuchos::RCP<MV> defV, augV;
947  if (inSituRestart_ == true) {
948  //
949  // get non-const pointer to solver's basis so we can work in situ
950  Teuchos::RCP<MV> solverbasis = Teuchos::rcp_const_cast<MV>(state.V);
951  //
952  // perform Householder QR of Su = Q [D;0], where D is unit diag.
953  // work on a copy of Su, since we need Su below to build newKK
955  std::vector<ScalarType> tau(numUnlocked), work(numUnlocked);
956  int info;
957  lapack.GEQRF(curdim,numUnlocked,copySu.values(),copySu.stride(),&tau[0],&work[0],work.size(),&info);
958  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,std::logic_error,
959  "Anasazi::BlockDavidsonSolMgr::solve(): error calling GEQRF during restarting.");
960  if (printer_->isVerbosity(Debug)) {
961  Teuchos::SerialDenseMatrix<int,ScalarType> R(Teuchos::Copy,copySu,numUnlocked,numUnlocked);
962  for (int j=0; j<numUnlocked; j++) {
963  R(j,j) = SCT::magnitude(R(j,j)) - 1.0;
964  for (int i=j+1; i<numUnlocked; i++) {
965  R(i,j) = ZERO;
966  }
967  }
968  printer_->stream(Debug) << "||Triangular factor of Su - I||: " << R.normFrobenius() << std::endl;
969  }
970  //
971  // perform implicit oldV*Su
972  // this actually performs oldV*[Su Sl*M] = [defV lockV], for some unitary M
973  // we are actually interested in only the first numUnlocked vectors of the result
974  {
975  Teuchos::RCP<MV> oldV = MVT::CloneViewNonConst(*solverbasis,curind);
976  msutils::applyHouse(numUnlocked,*oldV,copySu,tau,workMV);
977  }
978  std::vector<int> defind(numUnlocked), augind(numNewLocked);
979  for (int i=0; i<numUnlocked ; i++) defind[i] = i;
980  for (int i=0; i<numNewLocked; i++) augind[i] = numUnlocked+i;
981  defV = MVT::CloneViewNonConst(*solverbasis,defind);
982  augV = MVT::CloneViewNonConst(*solverbasis,augind);
983  }
984  else { // inSituRestart == false)
985  // defV = oldV*Su, explicitly. workspace is in workMV
986  std::vector<int> defind(numUnlocked), augind(numNewLocked);
987  for (int i=0; i<numUnlocked ; i++) defind[i] = i;
988  for (int i=0; i<numNewLocked; i++) augind[i] = numUnlocked+i;
989  Teuchos::RCP<const MV> oldV = MVT::CloneView(*state.V,curind);
990  defV = MVT::CloneViewNonConst(*workMV,defind);
991  augV = MVT::CloneViewNonConst(*workMV,augind);
992 
993  MVT::MvTimesMatAddMv(ONE,*oldV,Su,ZERO,*defV);
994  }
995 
996  //
997  // lockvecs will be partitioned as follows:
998  // lockvecs = [curlocked augTmp ...]
999  // - augTmp will be used for the storage of M*augV and K*augV
1000  // later, the locked vectors (stored in state.X and referenced via const MV view newLocked)
1001  // will be moved into lockvecs on top of augTmp when it is no longer needed as workspace.
1002  // - curlocked will be used in orthogonalization of augV
1003  //
1004  // newL is the new locked vectors; newL = oldV*Sl = RitzVectors(lockind)
1005  // we will not produce them, but instead retrieve them from RitzVectors
1006  //
1007  Teuchos::RCP<const MV> curlocked, newLocked;
1008  Teuchos::RCP<MV> augTmp;
1009  {
1010  // setup curlocked
1011  if (curNumLocked > 0) {
1012  std::vector<int> curlockind(curNumLocked);
1013  for (int i=0; i<curNumLocked; i++) curlockind[i] = i;
1014  curlocked = MVT::CloneView(*lockvecs,curlockind);
1015  }
1016  else {
1017  curlocked = Teuchos::null;
1018  }
1019  // setup augTmp
1020  std::vector<int> augtmpind(numNewLocked);
1021  for (int i=0; i<numNewLocked; i++) augtmpind[i] = curNumLocked+i;
1022  augTmp = MVT::CloneViewNonConst(*lockvecs,augtmpind);
1023  // setup newLocked
1024  newLocked = MVT::CloneView(*bd_solver->getRitzVectors(),lockind);
1025  }
1026 
1027  //
1028  // generate augV and perform orthogonalization
1029  //
1030  MVT::MvRandom(*augV);
1031  //
1032  // orthogonalize it against auxvecs, defV, and all locked vectors (new and current)
1033  // use augTmp as storage for M*augV, if hasM
1034  {
1037  if (probauxvecs != Teuchos::null) against.push_back(probauxvecs);
1038  if (curlocked != Teuchos::null) against.push_back(curlocked);
1039  against.push_back(newLocked);
1040  against.push_back(defV);
1041  if (problem_->getM() != Teuchos::null) {
1042  OPT::Apply(*problem_->getM(),*augV,*augTmp);
1043  }
1044  ortho->projectAndNormalizeMat(*augV,against,dummyC,Teuchos::null,augTmp);
1045  }
1046 
1047  //
1048  // form newKK
1049  //
1050  // newKK = newV'*K*newV = [Su'*KK*Su defV'*K*augV]
1051  // [augV'*K*defV augV'*K*augV]
1052  //
1053  // first, generate the principal submatrix, the projection of K onto the unlocked portion of oldV
1054  //
1055  Teuchos::SerialDenseMatrix<int,ScalarType> newKK(curdim,curdim);
1056  {
1057  Teuchos::SerialDenseMatrix<int,ScalarType> KKtmp(curdim,numUnlocked),
1058  KKold(Teuchos::View,*state.KK,curdim,curdim),
1059  KK11(Teuchos::View,newKK,numUnlocked,numUnlocked);
1060  int teuchosRet;
1061  // KKtmp = KKold*Su
1062  teuchosRet = KKtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,KKold,Su,ZERO);
1063  TEUCHOS_TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
1064  "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
1065  // KK11 = Su'*KKtmp = Su'*KKold*Su
1066  teuchosRet = KK11.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,Su,KKtmp,ZERO);
1067  TEUCHOS_TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
1068  "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
1069  }
1070  //
1071  // project the stiffness matrix on augV
1072  {
1073  OPT::Apply(*problem_->getOperator(),*augV,*augTmp);
1074  Teuchos::SerialDenseMatrix<int,ScalarType> KK12(Teuchos::View,newKK,numUnlocked,numNewLocked,0,numUnlocked),
1075  KK22(Teuchos::View,newKK,numNewLocked,numNewLocked,numUnlocked,numUnlocked);
1076  MVT::MvTransMv(ONE,*defV,*augTmp,KK12);
1077  MVT::MvTransMv(ONE,*augV,*augTmp,KK22);
1078  }
1079  //
1080  // done with defV,augV
1081  defV = Teuchos::null;
1082  augV = Teuchos::null;
1083  //
1084  // make it hermitian in memory (fill in KK21)
1085  for (int j=0; j<curdim; ++j) {
1086  for (int i=j+1; i<curdim; ++i) {
1087  newKK(i,j) = SCT::conjugate(newKK(j,i));
1088  }
1089  }
1090  //
1091  // we are done using augTmp as storage
1092  // put newLocked into lockvecs, new values into lockvals
1093  augTmp = Teuchos::null;
1094  {
1095  std::vector<Value<ScalarType> > allvals = bd_solver->getRitzValues();
1096  for (int i=0; i<numNewLocked; i++) {
1097  lockvals.push_back(allvals[lockind[i]].realpart);
1098  }
1099 
1100  std::vector<int> indlock(numNewLocked);
1101  for (int i=0; i<numNewLocked; i++) indlock[i] = curNumLocked+i;
1102  MVT::SetBlock(*newLocked,indlock,*lockvecs);
1103  newLocked = Teuchos::null;
1104 
1105  curNumLocked += numNewLocked;
1106  std::vector<int> curlockind(curNumLocked);
1107  for (int i=0; i<curNumLocked; i++) curlockind[i] = i;
1108  curlocked = MVT::CloneView(*lockvecs,curlockind);
1109  }
1110  // add locked vecs as aux vecs, along with aux vecs from problem
1111  // add lockvals to ordertest
1112  // disable locktest if curNumLocked == maxLocked
1113  {
1114  ordertest->setAuxVals(lockvals);
1115 
1117  if (probauxvecs != Teuchos::null) aux.push_back(probauxvecs);
1118  aux.push_back(curlocked);
1119  bd_solver->setAuxVecs(aux);
1120 
1121  if (curNumLocked == maxLocked_) {
1122  // disabled locking now
1123  combotest->removeTest(locktest);
1124  }
1125  }
1126 
1127  //
1128  // prepare new state
1130  rstate.curDim = curdim;
1131  if (inSituRestart_) {
1132  // data is already in the solver's memory
1133  rstate.V = state.V;
1134  }
1135  else {
1136  // data is in workspace and will be copied to solver memory
1137  rstate.V = workMV;
1138  }
1139  rstate.KK = Teuchos::rcpFromRef(newKK);
1140  //
1141  // pass new state to the solver
1142  bd_solver->initialize(rstate);
1143  } // end of locking
1145  //
1146  // we returned from iterate(), but none of our status tests Passed.
1147  // something is wrong, and it is probably our fault.
1148  //
1150  else {
1151  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): Invalid return from bd_solver::iterate().");
1152  }
1153  }
1154  catch (const AnasaziError &err) {
1155  printer_->stream(Errors)
1156  << "Anasazi::BlockDavidsonSolMgr::solve() caught unexpected exception from Anasazi::BlockDavidson::iterate() at iteration " << bd_solver->getNumIters() << std::endl
1157  << err.what() << std::endl
1158  << "Anasazi::BlockDavidsonSolMgr::solve() returning Unconverged with no solutions." << std::endl;
1159  return Unconverged;
1160  }
1161  }
1162 
1163  // clear temp space
1164  workMV = Teuchos::null;
1165 
1166  sol.numVecs = ordertest->howMany();
1167  if (sol.numVecs > 0) {
1168  sol.Evecs = MVT::Clone(*problem_->getInitVec(),sol.numVecs);
1169  sol.Espace = sol.Evecs;
1170  sol.Evals.resize(sol.numVecs);
1171  std::vector<MagnitudeType> vals(sol.numVecs);
1172 
1173  // copy them into the solution
1174  std::vector<int> which = ordertest->whichVecs();
1175  // indices between [0,blockSize) refer to vectors/values in the solver
1176  // indices between [-curNumLocked,-1] refer to locked vectors/values [0,curNumLocked)
1177  // everything has already been ordered by the solver; we just have to partition the two references
1178  std::vector<int> inlocked(0), insolver(0);
1179  for (unsigned int i=0; i<which.size(); i++) {
1180  if (which[i] >= 0) {
1181  TEUCHOS_TEST_FOR_EXCEPTION(which[i] >= blockSize_,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): positive indexing mistake from ordertest.");
1182  insolver.push_back(which[i]);
1183  }
1184  else {
1185  // sanity check
1186  TEUCHOS_TEST_FOR_EXCEPTION(which[i] < -curNumLocked,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): negative indexing mistake from ordertest.");
1187  inlocked.push_back(which[i] + curNumLocked);
1188  }
1189  }
1190 
1191  TEUCHOS_TEST_FOR_EXCEPTION(insolver.size() + inlocked.size() != (unsigned int)sol.numVecs,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): indexing mistake.");
1192 
1193  // set the vecs,vals in the solution
1194  if (insolver.size() > 0) {
1195  // set vecs
1196  int lclnum = insolver.size();
1197  std::vector<int> tosol(lclnum);
1198  for (int i=0; i<lclnum; i++) tosol[i] = i;
1199  Teuchos::RCP<const MV> v = MVT::CloneView(*bd_solver->getRitzVectors(),insolver);
1200  MVT::SetBlock(*v,tosol,*sol.Evecs);
1201  // set vals
1202  std::vector<Value<ScalarType> > fromsolver = bd_solver->getRitzValues();
1203  for (unsigned int i=0; i<insolver.size(); i++) {
1204  vals[i] = fromsolver[insolver[i]].realpart;
1205  }
1206  }
1207 
1208  // get the vecs,vals from locked storage
1209  if (inlocked.size() > 0) {
1210  int solnum = insolver.size();
1211  // set vecs
1212  int lclnum = inlocked.size();
1213  std::vector<int> tosol(lclnum);
1214  for (int i=0; i<lclnum; i++) tosol[i] = solnum + i;
1215  Teuchos::RCP<const MV> v = MVT::CloneView(*lockvecs,inlocked);
1216  MVT::SetBlock(*v,tosol,*sol.Evecs);
1217  // set vals
1218  for (unsigned int i=0; i<inlocked.size(); i++) {
1219  vals[i+solnum] = lockvals[inlocked[i]];
1220  }
1221  }
1222 
1223  // sort the eigenvalues and permute the eigenvectors appropriately
1224  {
1225  std::vector<int> order(sol.numVecs);
1226  sorter->sort(vals,Teuchos::rcpFromRef(order),sol.numVecs);
1227  // store the values in the Eigensolution
1228  for (int i=0; i<sol.numVecs; i++) {
1229  sol.Evals[i].realpart = vals[i];
1230  sol.Evals[i].imagpart = MT::zero();
1231  }
1232  // now permute the eigenvectors according to order
1233  msutils::permuteVectors(sol.numVecs,order,*sol.Evecs);
1234  }
1235 
1236  // setup sol.index, remembering that all eigenvalues are real so that index = {0,...,0}
1237  sol.index.resize(sol.numVecs,0);
1238  }
1239  }
1240 
1241  // print final summary
1242  bd_solver->currentStatus(printer_->stream(FinalSummary));
1243 
1244  // print timing information
1245 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1246  if ( printer_->isVerbosity( TimingDetails ) ) {
1247  Teuchos::TimeMonitor::summarize( printer_->stream( TimingDetails ) );
1248  }
1249 #endif
1250 
1251  problem_->setSolution(sol);
1252  printer_->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
1253 
1254  // get the number of iterations taken for this call to solve().
1255  numIters_ = bd_solver->getNumIters();
1256 
1257  if (sol.numVecs < nev) {
1258  return Unconverged; // return from BlockDavidsonSolMgr::solve()
1259  }
1260  return Converged; // return from BlockDavidsonSolMgr::solve()
1261 }
1262 
1263 
1264 template <class ScalarType, class MV, class OP>
1265 void
1267  const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global)
1268 {
1269  globalTest_ = global;
1270 }
1271 
1272 template <class ScalarType, class MV, class OP>
1275 {
1276  return globalTest_;
1277 }
1278 
1279 template <class ScalarType, class MV, class OP>
1280 void
1283 {
1284  debugTest_ = debug;
1285 }
1286 
1287 template <class ScalarType, class MV, class OP>
1290 {
1291  return debugTest_;
1292 }
1293 
1294 template <class ScalarType, class MV, class OP>
1295 void
1297  const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &locking)
1298 {
1299  lockingTest_ = locking;
1300 }
1301 
1302 template <class ScalarType, class MV, class OP>
1305 {
1306  return lockingTest_;
1307 }
1308 
1309 } // end Anasazi namespace
1310 
1311 #endif /* ANASAZI_BLOCKDAVIDSON_SOLMGR_HPP */
ScalarType * values() const
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. ...
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.
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...
SerialBandDenseMatrixPrinter< OrdinalType, ScalarType > printMat(const SerialBandDenseMatrix< OrdinalType, ScalarType > &obj)
basic_FancyOStream & setShowProcRank(const bool showProcRank)
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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.
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.
Output managers remove the need for the eigensolver to know any information about the required output...
This class implements a Block Davidson iteration, a preconditioned iteration for solving linear Hermi...
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Implementation of the block Davidson method.
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.
bool isParameter(const std::string &name) const
int numVecs
The number of computed eigenpairs.
int getNumIters() const
Get the iteration count for the most recent call to solve().
void COPY(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx, ScalarType *y, const OrdinalType &incy) const
Teuchos::RCP< MV > Espace
An orthonormal basis for the computed eigenspace.
Abstract class definition for Anasazi Output Managers.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Abstract base class which defines the interface required by an eigensolver and status test class to c...
static void summarize(Ptr< const Comm< int > > comm, std::ostream &out=std::cout, const bool alwaysWriteLocal=false, const bool writeGlobalStats=true, const bool writeZeroTimers=true, const ECounterSetOp setOp=Intersection, const std::string &filter="", const bool ignoreZeroTimers=false)
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.
Output managers remove the need for the eigensolver to know any information about the required output...
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.
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.
basic_FancyOStream & setShowAllFrontMatter(const bool showAllFrontMatter)
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
Orthogonalization manager based on the SVQB technique described in &quot;A Block Orthogonalization Procedu...
Struct for storing an eigenproblem solution.
int curDim
The current dimension of the solver.
void push_back(const value_type &x)
void GEQRF(const OrdinalType &m, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, ScalarType *TAU, ScalarType *WORK, const OrdinalType &lwork, OrdinalType *info) const
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.
Abstract class definition for Anasazi output stream.
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 ...
OrdinalType stride() const
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
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.