Anasazi  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
AnasaziBlockKrylovSchurSolMgr.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Anasazi: Block Eigensolvers Package
4 //
5 // Copyright 2004 NTESS and the Anasazi contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef ANASAZI_BLOCK_KRYLOV_SCHUR_SOLMGR_HPP
11 #define ANASAZI_BLOCK_KRYLOV_SCHUR_SOLMGR_HPP
12 
16 
17 #include "AnasaziConfigDefs.hpp"
18 #include "AnasaziTypes.hpp"
19 
20 #include "AnasaziEigenproblem.hpp"
21 #include "AnasaziSolverManager.hpp"
22 #include "AnasaziSolverUtils.hpp"
23 
25 #include "AnasaziBasicSort.hpp"
33 #include "AnasaziOutputManager.hpp"
35 
36 #include "Teuchos_BLAS.hpp"
37 #include "Teuchos_LAPACK.hpp"
38 #include "Teuchos_TimeMonitor.hpp"
39 #include "Teuchos_FancyOStream.hpp"
40 
47 
95 namespace Anasazi {
96 
97 
124 template<class ScalarType, class MV, class OP>
125 class BlockKrylovSchurSolMgr : public SolverManager<ScalarType,MV,OP> {
126 
127  private:
131  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
133 
134  public:
135 
137 
138 
161 
165 
167 
168 
171  return *problem_;
172  }
173 
175  int getNumIters() const {
176  return numIters_;
177  }
178 
181  std::vector<Value<ScalarType> > getRitzValues() const {
182  std::vector<Value<ScalarType> > ret( ritzValues_ );
183  return ret;
184  }
185 
193  return Teuchos::tuple(timerSolve_, timerRestarting_);
194  }
195 
197 
199 
200 
219  ReturnType solve();
220 
223 
226 
229 
232 
234 
235  private:
238 
239  std::string whch_, ortho_;
240  MagnitudeType ortho_kappa_;
241 
242  MagnitudeType convtol_;
243  int maxRestarts_;
244  bool relconvtol_,conjSplit_;
245  int blockSize_, numBlocks_, stepSize_, nevBlocks_, xtra_nevBlocks_;
246  int numIters_;
247  int verbosity_;
248  bool inSituRestart_, dynXtraNev_;
249  int osProc_;
250 
251  std::vector<Value<ScalarType> > ritzValues_;
252 
253  int printNum_;
254  Teuchos::RCP<Teuchos::Time> timerSolve_, timerRestarting_;
255 
257 
260 
261 };
262 
263 
264 // Constructor
265 template<class ScalarType, class MV, class OP>
268  Teuchos::ParameterList &pl ) :
269  problem_(problem),
270  whch_("LM"),
271  ortho_("SVQB"),
272  ortho_kappa_(-1.0),
273  convtol_(0),
274  maxRestarts_(20),
275  relconvtol_(true),
276  conjSplit_(false),
277  blockSize_(0),
278  numBlocks_(0),
279  stepSize_(0),
280  nevBlocks_(0),
281  xtra_nevBlocks_(0),
282  numIters_(0),
283  verbosity_(Anasazi::Errors),
284  inSituRestart_(false),
285  dynXtraNev_(false),
286  osProc_(0),
287  printNum_(-1)
288 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
289  ,timerSolve_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchurSolMgr::solve()")),
290  timerRestarting_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchurSolMgr restarting"))
291 #endif
292 {
293  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
294  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument, "Problem not set.");
295  TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null, std::invalid_argument, "Problem does not contain initial vectors to clone from.");
296 
297  const int nev = problem_->getNEV();
298 
299  // convergence tolerance
300  convtol_ = pl.get("Convergence Tolerance",MT::prec());
301  relconvtol_ = pl.get("Relative Convergence Tolerance",relconvtol_);
302 
303  // maximum number of restarts
304  maxRestarts_ = pl.get("Maximum Restarts",maxRestarts_);
305 
306  // block size: default is 1
307  blockSize_ = pl.get("Block Size",1);
308  TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
309  "Anasazi::BlockKrylovSchurSolMgr: \"Block Size\" must be strictly positive.");
310 
311  // set the number of blocks we need to save to compute the nev eigenvalues of interest.
312  xtra_nevBlocks_ = pl.get("Extra NEV Blocks",0);
313  if (nev%blockSize_) {
314  nevBlocks_ = nev/blockSize_ + 1;
315  } else {
316  nevBlocks_ = nev/blockSize_;
317  }
318 
319  // determine if we should use the dynamic scheme for selecting the current number of retained eigenvalues.
320  // NOTE: This employs a technique similar to ARPACK in that it increases the number of retained eigenvalues
321  // by one for every converged eigenpair to accelerate convergence.
322  if (pl.isParameter("Dynamic Extra NEV")) {
323  if (Teuchos::isParameterType<bool>(pl,"Dynamic Extra NEV")) {
324  dynXtraNev_ = pl.get("Dynamic Extra NEV",dynXtraNev_);
325  } else {
326  dynXtraNev_ = ( Teuchos::getParameter<int>(pl,"Dynamic Extra NEV") != 0 );
327  }
328  }
329 
330  numBlocks_ = pl.get("Num Blocks",3*nevBlocks_);
331  TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= nevBlocks_, std::invalid_argument,
332  "Anasazi::BlockKrylovSchurSolMgr: \"Num Blocks\" must be strictly positive and large enough to compute the requested eigenvalues.");
333 
334  TEUCHOS_TEST_FOR_EXCEPTION(static_cast<ptrdiff_t>(numBlocks_)*blockSize_ > MVT::GetGlobalLength(*problem_->getInitVec()),
335  std::invalid_argument,
336  "Anasazi::BlockKrylovSchurSolMgr: Potentially impossible orthogonality requests. Reduce basis size.");
337 
338  // step size: the default is maxRestarts_*numBlocks_, so that Ritz values are only computed every restart.
339  if (maxRestarts_) {
340  stepSize_ = pl.get("Step Size", (maxRestarts_+1)*(numBlocks_+1));
341  } else {
342  stepSize_ = pl.get("Step Size", numBlocks_+1);
343  }
344  TEUCHOS_TEST_FOR_EXCEPTION(stepSize_ < 1, std::invalid_argument,
345  "Anasazi::BlockKrylovSchurSolMgr: \"Step Size\" must be strictly positive.");
346 
347  // get the sort manager
348  if (pl.isParameter("Sort Manager")) {
349  sort_ = Teuchos::getParameter<Teuchos::RCP<Anasazi::SortManager<MagnitudeType> > >(pl,"Sort Manager");
350  } else {
351  // which values to solve for
352  whch_ = pl.get("Which",whch_);
353  TEUCHOS_TEST_FOR_EXCEPTION(whch_ != "SM" && whch_ != "LM" && whch_ != "SR" && whch_ != "LR" && whch_ != "SI" && whch_ != "LI",
354  std::invalid_argument, "Invalid sorting string.");
355  sort_ = Teuchos::rcp( new BasicSort<MagnitudeType>(whch_) );
356  }
357 
358  // which orthogonalization to use
359  ortho_ = pl.get("Orthogonalization",ortho_);
360  if (ortho_ != "DGKS" && ortho_ != "SVQB" && ortho_ != "ICGS") {
361  ortho_ = "SVQB";
362  }
363 
364  // which orthogonalization constant to use
365  ortho_kappa_ = pl.get("Orthogonalization Constant",ortho_kappa_);
366 
367  // Create a formatted output stream to print to.
368  // See if user requests output processor.
369  osProc_ = pl.get("Output Processor", osProc_);
370 
371  // If not passed in by user, it will be chosen based upon operator type.
372  if (pl.isParameter("Output Stream")) {
373  osp_ = Teuchos::getParameter<Teuchos::RCP<Teuchos::FancyOStream> >(pl,"Output Stream");
374  }
375  else {
376  osp_ = OutputStreamTraits<OP>::getOutputStream (*problem_->getOperator(),osProc_);
377  }
378 
379  // verbosity level
380  if (pl.isParameter("Verbosity")) {
381  if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
382  verbosity_ = pl.get("Verbosity", verbosity_);
383  } else {
384  verbosity_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
385  }
386  }
387 
388  // restarting technique: V*Q or applyHouse(V,H,tau)
389  if (pl.isParameter("In Situ Restarting")) {
390  if (Teuchos::isParameterType<bool>(pl,"In Situ Restarting")) {
391  inSituRestart_ = pl.get("In Situ Restarting",inSituRestart_);
392  } else {
393  inSituRestart_ = ( Teuchos::getParameter<int>(pl,"In Situ Restarting") != 0 );
394  }
395  }
396 
397  printNum_ = pl.get<int>("Print Number of Ritz Values",printNum_);
398 }
399 
400 
401 // solve()
402 template<class ScalarType, class MV, class OP>
405 
406  const int nev = problem_->getNEV();
407  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
408  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
409 
412 
414  // Output manager
416 
418  // Status tests
419  //
420  // convergence
422  if (globalTest_ == Teuchos::null) {
423  convtest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(convtol_,nev,RITZRES_2NORM,relconvtol_) );
424  }
425  else {
426  convtest = globalTest_;
427  }
429  = Teuchos::rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sort_,nev) );
430  // for a non-short-circuited OR test, the order doesn't matter
432  alltests.push_back(ordertest);
433 
434  if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
435 
438  // printing StatusTest
440  if ( printer->isVerbosity(Debug) ) {
441  outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,Passed+Failed+Undefined ) );
442  }
443  else {
444  outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,Passed ) );
445  }
446 
448  // Orthomanager
450  if (ortho_=="SVQB") {
451  ortho = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
452  } else if (ortho_=="DGKS") {
453  if (ortho_kappa_ <= 0) {
454  ortho = Teuchos::rcp( new BasicOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
455  } else {
456  ortho = Teuchos::rcp( new BasicOrthoManager<ScalarType,MV,OP>(problem_->getM(),ortho_kappa_) );
457  }
458  } else if (ortho_=="ICGS") {
459  ortho = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
460  } else {
461  TEUCHOS_TEST_FOR_EXCEPTION(ortho_!="SVQB"&&ortho_!="DGKS"&&ortho_!="ICGS",std::logic_error,"Anasazi::BlockKrylovSchurSolMgr::solve(): Invalid orthogonalization type.");
462  }
463 
465  // Parameter list
467  plist.set("Block Size",blockSize_);
468  plist.set("Num Blocks",numBlocks_);
469  plist.set("Step Size",stepSize_);
470  if (printNum_ == -1) {
471  plist.set("Print Number of Ritz Values",nevBlocks_*blockSize_);
472  }
473  else {
474  plist.set("Print Number of Ritz Values",printNum_);
475  }
476 
478  // BlockKrylovSchur solver
480  = Teuchos::rcp( new BlockKrylovSchur<ScalarType,MV,OP>(problem_,sort_,printer,outputtest,ortho,plist) );
481  // set any auxiliary vectors defined in the problem
482  Teuchos::RCP< const MV > probauxvecs = problem_->getAuxVecs();
483  if (probauxvecs != Teuchos::null) {
484  bks_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
485  }
486 
487  // Create workspace for the Krylov basis generated during a restart
488  // Need at most (nevBlocks_*blockSize_+1) for the updated factorization and another block for the current factorization residual block (F).
489  // ---> (nevBlocks_*blockSize_+1) + blockSize_
490  // If Hermitian, this becomes nevBlocks_*blockSize_ + blockSize_
491  // we only need this if there is the possibility of restarting, ex situ
492 
493  // Maximum allowable extra vectors that BKS can keep to accelerate convergence
494  int maxXtraBlocks = 0;
495  if ( dynXtraNev_ ) maxXtraBlocks = ( ( bks_solver->getMaxSubspaceDim() - nev ) / blockSize_ ) / 2;
496 
497  Teuchos::RCP<MV> workMV;
498  if (maxRestarts_ > 0) {
499  if (inSituRestart_==true) {
500  // still need one work vector for applyHouse()
501  workMV = MVT::Clone( *problem_->getInitVec(), 1 );
502  }
503  else { // inSituRestart == false
504  if (problem_->isHermitian()) {
505  workMV = MVT::Clone( *problem_->getInitVec(), (nevBlocks_+maxXtraBlocks)*blockSize_ + blockSize_ );
506  } else {
507  // Increase workspace by 1 because of potential complex conjugate pairs on the Ritz vector boundary
508  workMV = MVT::Clone( *problem_->getInitVec(), (nevBlocks_+maxXtraBlocks)*blockSize_ + 1 + blockSize_ );
509  }
510  }
511  } else {
512  workMV = Teuchos::null;
513  }
514 
515  // go ahead and initialize the solution to nothing in case we throw an exception
517  sol.numVecs = 0;
518  problem_->setSolution(sol);
519 
520  int numRestarts = 0;
521  int cur_nevBlocks = 0;
522 
523  // enter solve() iterations
524  {
525 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
526  Teuchos::TimeMonitor slvtimer(*timerSolve_);
527 #endif
528 
529  // tell bks_solver to iterate
530  while (1) {
531  try {
532  bks_solver->iterate();
533 
535  //
536  // check convergence first
537  //
539  if ( ordertest->getStatus() == Passed ) {
540  // we have convergence
541  // ordertest->whichVecs() tells us which vectors from solver state are the ones we want
542  // ordertest->howMany() will tell us how many
543  break;
544  }
546  //
547  // check for restarting, i.e. the subspace is full
548  //
550  // this is for the Hermitian case, or non-Hermitian conjugate split situation.
551  // --> for the Hermitian case the current subspace dimension needs to match the maximum subspace dimension
552  // --> for the non-Hermitian case:
553  // --> if a conjugate pair was detected in the previous restart then the current subspace dimension needs to match the
554  // maximum subspace dimension (the BKS solver keeps one extra vector if the problem is non-Hermitian).
555  // --> if a conjugate pair was not detected in the previous restart then the current subspace dimension will be one less
556  // than the maximum subspace dimension.
557  else if ( (bks_solver->getCurSubspaceDim() == bks_solver->getMaxSubspaceDim()) ||
558  (!problem_->isHermitian() && !conjSplit_ && (bks_solver->getCurSubspaceDim()+1 == bks_solver->getMaxSubspaceDim())) ) {
559 
560  // Update the Schur form of the projected eigenproblem, then sort it.
561  if (!bks_solver->isSchurCurrent()) {
562  bks_solver->computeSchurForm( true );
563 
564  // Check for convergence, just in case we wait for every restart to check
565  outputtest->checkStatus( &*bks_solver );
566  }
567 
568  // Don't bother to restart if we've converged or reached the maximum number of restarts
569  if ( numRestarts >= maxRestarts_ || ordertest->getStatus() == Passed) {
570  break; // break from while(1){bks_solver->iterate()}
571  }
572 
573  // Start restarting timer and increment counter
574 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
575  Teuchos::TimeMonitor restimer(*timerRestarting_);
576 #endif
577  numRestarts++;
578 
579  int numConv = ordertest->howMany();
580  cur_nevBlocks = nevBlocks_*blockSize_;
581 
582  // Add in extra blocks for restarting if either static or dynamic boundaries are being used.
583  int moreNevBlocks = std::min( maxXtraBlocks, std::max( numConv/blockSize_, xtra_nevBlocks_) );
584  if ( dynXtraNev_ )
585  cur_nevBlocks += moreNevBlocks * blockSize_;
586  else if ( xtra_nevBlocks_ )
587  cur_nevBlocks += xtra_nevBlocks_ * blockSize_;
588 /*
589  int cur_numConv = numConv;
590  while ( (cur_nevBlocks < (_nevBlocks + maxXtraVecs)) && cur_numConv > 0 ) {
591  cur_nevBlocks++;
592  cur_numConv--;
593 */
594 
595  printer->stream(Debug) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << std::endl << std::endl;
596  printer->stream(Debug) << " - Current NEV blocks is " << cur_nevBlocks << ", the minimum is " << nevBlocks_*blockSize_ << std::endl;
597 
598  // Get the most current Ritz values before we continue.
599  ritzValues_ = bks_solver->getRitzValues();
600 
601  // Get the state.
602  BlockKrylovSchurState<ScalarType,MV> oldState = bks_solver->getState();
603 
604  // Get the current dimension of the factorization
605  int curDim = oldState.curDim;
606 
607  // Determine if the storage for the nev eigenvalues of interest splits a complex conjugate pair.
608  std::vector<int> ritzIndex = bks_solver->getRitzIndex();
609  if (ritzIndex[cur_nevBlocks-1]==1) {
610  conjSplit_ = true;
611  cur_nevBlocks++;
612  } else {
613  conjSplit_ = false;
614  }
615 
616  // Print out a warning to the user if complex eigenvalues were found on the boundary of the restart subspace
617  // and the eigenproblem is Hermitian. This solver is not prepared to handle this situation.
618  if (problem_->isHermitian() && conjSplit_)
619  {
620  printer->stream(Warnings)
621  << " Eigenproblem is Hermitian, complex eigenvalues have been detected, and eigenvalues of interest split a conjugate pair!!!"
622  << std::endl
623  << " Block Krylov-Schur eigensolver cannot guarantee correct behavior in this situation, please turn Hermitian flag off!!!"
624  << std::endl;
625  }
626 
627  // Update the Krylov-Schur decomposition
628 
629  // Get a view of the Schur vectors of interest.
630  Teuchos::SerialDenseMatrix<int,ScalarType> Qnev(Teuchos::View, *(oldState.Q), curDim, cur_nevBlocks);
631 
632  // Get a view of the current Krylov basis.
633  std::vector<int> curind( curDim );
634  for (int i=0; i<curDim; i++) { curind[i] = i; }
635  Teuchos::RCP<const MV> basistemp = MVT::CloneView( *(oldState.V), curind );
636 
637  // Compute the new Krylov basis: Vnew = V*Qnev
638  //
639  // this will occur ex situ in workspace allocated for this purpose (tmpMV)
640  // or in situ in the solver's memory space.
641  //
642  // we will also set a pointer for the location that the current factorization residual block (F),
643  // currently located after the current basis in oldstate.V, will be moved to
644  //
645  Teuchos::RCP<MV> newF;
646  if (inSituRestart_) {
647  //
648  // get non-const pointer to solver's basis so we can work in situ
649  Teuchos::RCP<MV> solverbasis = Teuchos::rcp_const_cast<MV>(oldState.V);
651  //
652  // perform Householder QR of copyQnev = Q [D;0], where D is unit diag. We will want D below.
653  std::vector<ScalarType> tau(cur_nevBlocks), work(cur_nevBlocks);
654  int info;
655  lapack.GEQRF(curDim,cur_nevBlocks,copyQnev.values(),copyQnev.stride(),&tau[0],&work[0],work.size(),&info);
656  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,std::logic_error,
657  "Anasazi::BlockKrylovSchurSolMgr::solve(): error calling GEQRF during restarting.");
658  // we need to get the diagonal of D
659  std::vector<ScalarType> d(cur_nevBlocks);
660  for (int j=0; j<copyQnev.numCols(); j++) {
661  d[j] = copyQnev(j,j);
662  }
663  if (printer->isVerbosity(Debug)) {
664  Teuchos::SerialDenseMatrix<int,ScalarType> R(Teuchos::Copy,copyQnev,cur_nevBlocks,cur_nevBlocks);
665  for (int j=0; j<R.numCols(); j++) {
666  R(j,j) = SCT::magnitude(R(j,j)) - 1.0;
667  for (int i=j+1; i<R.numRows(); i++) {
668  R(i,j) = zero;
669  }
670  }
671  printer->stream(Debug) << "||Triangular factor of Su - I||: " << R.normFrobenius() << std::endl;
672  }
673  //
674  // perform implicit V*Qnev
675  // this actually performs V*[Qnev Qtrunc*M] = [newV truncV], for some unitary M
676  // we are interested in only the first cur_nevBlocks vectors of the result
677  curind.resize(curDim);
678  for (int i=0; i<curDim; i++) curind[i] = i;
679  {
680  Teuchos::RCP<MV> oldV = MVT::CloneViewNonConst(*solverbasis,curind);
681  SolverUtils<ScalarType,MV,OP>::applyHouse(cur_nevBlocks,*oldV,copyQnev,tau,workMV);
682  }
683  // multiply newV*D
684  // get pointer to new basis
685  curind.resize(cur_nevBlocks);
686  for (int i=0; i<cur_nevBlocks; i++) { curind[i] = i; }
687  {
688  Teuchos::RCP<MV> newV = MVT::CloneViewNonConst( *solverbasis, curind );
689  MVT::MvScale(*newV,d);
690  }
691  // get pointer to new location for F
692  curind.resize(blockSize_);
693  for (int i=0; i<blockSize_; i++) { curind[i] = cur_nevBlocks + i; }
694  newF = MVT::CloneViewNonConst( *solverbasis, curind );
695  }
696  else {
697  // get pointer to first part of work space
698  curind.resize(cur_nevBlocks);
699  for (int i=0; i<cur_nevBlocks; i++) { curind[i] = i; }
700  Teuchos::RCP<MV> tmp_newV = MVT::CloneViewNonConst(*workMV, curind );
701  // perform V*Qnev
702  MVT::MvTimesMatAddMv( one, *basistemp, Qnev, zero, *tmp_newV );
703  tmp_newV = Teuchos::null;
704  // get pointer to new location for F
705  curind.resize(blockSize_);
706  for (int i=0; i<blockSize_; i++) { curind[i] = cur_nevBlocks + i; }
707  newF = MVT::CloneViewNonConst( *workMV, curind );
708  }
709 
710  // Move the current factorization residual block (F) to the last block of newV.
711  curind.resize(blockSize_);
712  for (int i=0; i<blockSize_; i++) { curind[i] = curDim + i; }
713  Teuchos::RCP<const MV> oldF = MVT::CloneView( *(oldState.V), curind );
714  for (int i=0; i<blockSize_; i++) { curind[i] = i; }
715  MVT::SetBlock( *oldF, curind, *newF );
716  newF = Teuchos::null;
717 
718  // Update the Krylov-Schur quasi-triangular matrix.
719  //
720  // Create storage for the new Schur matrix of the Krylov-Schur factorization
721  // Copy over the current quasi-triangular factorization of oldState.H which is stored in oldState.S.
723  Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::Copy, *(oldState.S), cur_nevBlocks+blockSize_, cur_nevBlocks) );
724  //
725  // Get a view of the B block of the current factorization
726  Teuchos::SerialDenseMatrix<int,ScalarType> oldB(Teuchos::View, *(oldState.H), blockSize_, blockSize_, curDim, curDim-blockSize_);
727  //
728  // Get a view of the a block row of the Schur vectors.
729  Teuchos::SerialDenseMatrix<int,ScalarType> subQ(Teuchos::View, *(oldState.Q), blockSize_, cur_nevBlocks, curDim-blockSize_);
730  //
731  // Get a view of the new B block of the updated Krylov-Schur factorization
732  Teuchos::SerialDenseMatrix<int,ScalarType> newB(Teuchos::View, *newH, blockSize_, cur_nevBlocks, cur_nevBlocks);
733  //
734  // Compute the new B block.
735  blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, blockSize_, cur_nevBlocks, blockSize_, one,
736  oldB.values(), oldB.stride(), subQ.values(), subQ.stride(), zero, newB.values(), newB.stride() );
737 
738 
739  //
740  // Set the new state and initialize the solver.
742  if (inSituRestart_) {
743  newstate.V = oldState.V;
744  } else {
745  newstate.V = workMV;
746  }
747  newstate.H = newH;
748  newstate.curDim = cur_nevBlocks;
749  bks_solver->initialize(newstate);
750 
751  } // end of restarting
753  //
754  // we returned from iterate(), but none of our status tests Passed.
755  // something is wrong, and it is probably our fault.
756  //
758  else {
759  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::BlockKrylovSchurSolMgr::solve(): Invalid return from bks_solver::iterate().");
760  }
761  }
762  catch (const AnasaziError &err) {
763  printer->stream(Errors)
764  << "Anasazi::BlockKrylovSchurSolMgr::solve() caught unexpected exception from Anasazi::BlockKrylovSchur::iterate() at iteration " << bks_solver->getNumIters() << std::endl
765  << err.what() << std::endl
766  << "Anasazi::BlockKrylovSchurSolMgr::solve() returning Unconverged with no solutions." << std::endl;
767  return Unconverged;
768  }
769  }
770 
771  //
772  // free temporary space
773  workMV = Teuchos::null;
774 
775  // Get the most current Ritz values before we return
776  ritzValues_ = bks_solver->getRitzValues();
777 
778  sol.numVecs = ordertest->howMany();
779  printer->stream(Debug) << "ordertest->howMany() : " << sol.numVecs << std::endl;
780  std::vector<int> whichVecs = ordertest->whichVecs();
781 
782  // Place any converged eigenpairs in the solution container.
783  if (sol.numVecs > 0) {
784 
785  // Next determine if there is a conjugate pair on the boundary and resize.
786  std::vector<int> tmpIndex = bks_solver->getRitzIndex();
787  for (int i=0; i<(int)ritzValues_.size(); ++i) {
788  printer->stream(Debug) << ritzValues_[i].realpart << " + i " << ritzValues_[i].imagpart << ", Index = " << tmpIndex[i] << std::endl;
789  }
790  printer->stream(Debug) << "Number of converged eigenpairs (before) = " << sol.numVecs << std::endl;
791  for (int i=0; i<sol.numVecs; ++i) {
792  printer->stream(Debug) << "whichVecs[" << i << "] = " << whichVecs[i] << ", tmpIndex[" << whichVecs[i] << "] = " << tmpIndex[whichVecs[i]] << std::endl;
793  }
794  if (tmpIndex[whichVecs[sol.numVecs-1]]==1) {
795  printer->stream(Debug) << "There is a conjugate pair on the boundary, resizing sol.numVecs" << std::endl;
796  whichVecs.push_back(whichVecs[sol.numVecs-1]+1);
797  sol.numVecs++;
798  for (int i=0; i<sol.numVecs; ++i) {
799  printer->stream(Debug) << "whichVecs[" << i << "] = " << whichVecs[i] << ", tmpIndex[" << whichVecs[i] << "] = " << tmpIndex[whichVecs[i]] << std::endl;
800  }
801  }
802 
803  bool keepMore = false;
804  int numEvecs = sol.numVecs;
805  printer->stream(Debug) << "Number of converged eigenpairs (after) = " << sol.numVecs << std::endl;
806  printer->stream(Debug) << "whichVecs[sol.numVecs-1] > sol.numVecs-1 : " << whichVecs[sol.numVecs-1] << " > " << sol.numVecs-1 << std::endl;
807  if (whichVecs[sol.numVecs-1] > (sol.numVecs-1)) {
808  keepMore = true;
809  numEvecs = whichVecs[sol.numVecs-1]+1; // Add 1 to fix zero-based indexing
810  printer->stream(Debug) << "keepMore = true; numEvecs = " << numEvecs << std::endl;
811  }
812 
813  // Next set the number of Ritz vectors that the iteration must compute and compute them.
814  bks_solver->setNumRitzVectors(numEvecs);
815  bks_solver->computeRitzVectors();
816 
817  // If the leading Ritz pairs are the converged ones, get the information
818  // from the iteration to the solution container. Otherwise copy the necessary
819  // information using 'whichVecs'.
820  if (!keepMore) {
821  sol.index = bks_solver->getRitzIndex();
822  sol.Evals = bks_solver->getRitzValues();
823  sol.Evecs = MVT::CloneCopy( *(bks_solver->getRitzVectors()) );
824  }
825 
826  // Resize based on the number of solutions being returned and set the number of Ritz
827  // vectors for the iteration to compute.
828  sol.Evals.resize(sol.numVecs);
829  sol.index.resize(sol.numVecs);
830 
831  // If the converged Ritz pairs are not the leading ones, copy over the information directly.
832  if (keepMore) {
833  std::vector<Anasazi::Value<ScalarType> > tmpEvals = bks_solver->getRitzValues();
834  for (int vec_i=0; vec_i<sol.numVecs; ++vec_i) {
835  sol.index[vec_i] = tmpIndex[whichVecs[vec_i]];
836  sol.Evals[vec_i] = tmpEvals[whichVecs[vec_i]];
837  }
838  sol.Evecs = MVT::CloneCopy( *(bks_solver->getRitzVectors()), whichVecs );
839  }
840 
841  // Set the solution space to be the Ritz vectors at this time.
842  sol.Espace = sol.Evecs;
843  }
844  }
845 
846  // print final summary
847  bks_solver->currentStatus(printer->stream(FinalSummary));
848 
849  // print timing information
850 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
851  if ( printer->isVerbosity( TimingDetails ) ) {
852  Teuchos::TimeMonitor::summarize( printer->stream( TimingDetails ) );
853  }
854 #endif
855 
856  problem_->setSolution(sol);
857  printer->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
858 
859  // get the number of iterations performed during this solve.
860  numIters_ = bks_solver->getNumIters();
861 
862  if (sol.numVecs < nev) {
863  return Unconverged; // return from BlockKrylovSchurSolMgr::solve()
864  }
865  return Converged; // return from BlockKrylovSchurSolMgr::solve()
866 }
867 
868 
869 template <class ScalarType, class MV, class OP>
870 void
873 {
874  globalTest_ = global;
875 }
876 
877 template <class ScalarType, class MV, class OP>
880 {
881  return globalTest_;
882 }
883 
884 template <class ScalarType, class MV, class OP>
885 void
888 {
889  debugTest_ = debug;
890 }
891 
892 template <class ScalarType, class MV, class OP>
895 {
896  return debugTest_;
897 }
898 
899 } // end Anasazi namespace
900 
901 #endif /* ANASAZI_BLOCK_KRYLOV_SCHUR_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. ...
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.
A special StatusTest for printing other status tests.
This class defines the interface required by an eigensolver and status test class to compute solution...
An implementation of the Anasazi::SortManager that performs a collection of common sorting techniques...
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)
#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.
An implementation of the Anasazi::GenOrthoManager that performs orthogonalization using iterated clas...
ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver&#39;s iterate() routine until ...
Status test for forming logical combinations of other status tests.
The Anasazi::SolverManager is a templated virtual base class that defines the basic interface that an...
void setDebugStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &debug)
Set 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...
Structure to contain pointers to BlockKrylovSchur state variables.
void GEMM(ETransp transa, ETransp transb, const OrdinalType &m, const OrdinalType &n, const OrdinalType &k, const alpha_type alpha, const A_type *A, const OrdinalType &lda, const B_type *B, const OrdinalType &ldb, const beta_type beta, ScalarType *C, const OrdinalType &ldc) const
An exception class parent to all Anasazi exceptions.
Output managers remove the need for the eigensolver to know any information about the required output...
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
The Anasazi::BlockKrylovSchurSolMgr provides a flexible solver manager over the BlockKrylovSchur eige...
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
bool isParameter(const std::string &name) const
int numVecs
The number of computed eigenpairs.
std::vector< Value< ScalarType > > getRitzValues() const
Return the Ritz values from the most recent solve.
void setGlobalStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &global)
Set the status test defining global convergence.
BlockKrylovSchurSolMgr(const Teuchos::RCP< Eigenproblem< ScalarType, MV, OP > > &problem, Teuchos::ParameterList &pl)
Basic constructor for BlockKrylovSchurSolMgr.
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)
static magnitudeType prec()
static void applyHouse(int k, MV &V, const Teuchos::SerialDenseMatrix< int, ScalarType > &H, const std::vector< ScalarType > &tau, Teuchos::RCP< MV > workMV=Teuchos::null)
Apply a sequence of Householder reflectors (from GEQRF) to a multivector, using minimal workspace...
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...
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 Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
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.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getGlobalStatusTest() const
Get the status test defining global convergence.
OrdinalType numCols() const
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
Teuchos::RCP< const MulVec > V
The current Krylov basis.
Special StatusTest for printing status tests.
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
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.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > Q
The current Schur vectors of the valid part of H.
Implementation of a block Krylov-Schur eigensolver.
This class implements the block Krylov-Schur iteration, for solving linear eigenvalue problems...
int curDim
The current dimension of the reduction.
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.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getDebugStatusTest() const
Get the status test for debugging.
Basic implementation of the Anasazi::OrthoManager class.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > S
The current Schur form reduction of the valid part of H.
OrdinalType stride() const
int getNumIters() const
Get the iteration count for the most recent call to solve().
OrdinalType numRows() const
Class which provides internal utilities for the Anasazi solvers.