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 //
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_BLOCK_KRYLOV_SCHUR_SOLMGR_HPP
43 #define ANASAZI_BLOCK_KRYLOV_SCHUR_SOLMGR_HPP
44 
48 
49 #include "AnasaziConfigDefs.hpp"
50 #include "AnasaziTypes.hpp"
51 
52 #include "AnasaziEigenproblem.hpp"
53 #include "AnasaziSolverManager.hpp"
54 #include "AnasaziSolverUtils.hpp"
55 
57 #include "AnasaziBasicSort.hpp"
65 #include "AnasaziOutputManager.hpp"
67 
68 #include "Teuchos_BLAS.hpp"
69 #include "Teuchos_LAPACK.hpp"
70 #include "Teuchos_TimeMonitor.hpp"
71 #include "Teuchos_FancyOStream.hpp"
72 
79 
127 namespace Anasazi {
128 
129 
156 template<class ScalarType, class MV, class OP>
157 class BlockKrylovSchurSolMgr : public SolverManager<ScalarType,MV,OP> {
158 
159  private:
163  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
165 
166  public:
167 
169 
170 
193 
197 
199 
200 
203  return *problem_;
204  }
205 
207  int getNumIters() const {
208  return numIters_;
209  }
210 
213  std::vector<Value<ScalarType> > getRitzValues() const {
214  std::vector<Value<ScalarType> > ret( ritzValues_ );
215  return ret;
216  }
217 
225  return Teuchos::tuple(timerSolve_, timerRestarting_);
226  }
227 
229 
231 
232 
251  ReturnType solve();
252 
255 
258 
261 
264 
266 
267  private:
270 
271  std::string whch_, ortho_;
272  MagnitudeType ortho_kappa_;
273 
274  MagnitudeType convtol_;
275  int maxRestarts_;
276  bool relconvtol_,conjSplit_;
277  int blockSize_, numBlocks_, stepSize_, nevBlocks_, xtra_nevBlocks_;
278  int numIters_;
279  int verbosity_;
280  bool inSituRestart_, dynXtraNev_;
281  int osProc_;
282 
283  std::vector<Value<ScalarType> > ritzValues_;
284 
285  int printNum_;
286  Teuchos::RCP<Teuchos::Time> timerSolve_, timerRestarting_;
287 
289 
292 
293 };
294 
295 
296 // Constructor
297 template<class ScalarType, class MV, class OP>
300  Teuchos::ParameterList &pl ) :
301  problem_(problem),
302  whch_("LM"),
303  ortho_("SVQB"),
304  ortho_kappa_(-1.0),
305  convtol_(0),
306  maxRestarts_(20),
307  relconvtol_(true),
308  conjSplit_(false),
309  blockSize_(0),
310  numBlocks_(0),
311  stepSize_(0),
312  nevBlocks_(0),
313  xtra_nevBlocks_(0),
314  numIters_(0),
315  verbosity_(Anasazi::Errors),
316  inSituRestart_(false),
317  dynXtraNev_(false),
318  osProc_(0),
319  printNum_(-1)
320 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
321  ,timerSolve_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchurSolMgr::solve()")),
322  timerRestarting_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchurSolMgr restarting"))
323 #endif
324 {
325  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
326  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument, "Problem not set.");
327  TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null, std::invalid_argument, "Problem does not contain initial vectors to clone from.");
328 
329  const int nev = problem_->getNEV();
330 
331  // convergence tolerance
332  convtol_ = pl.get("Convergence Tolerance",MT::prec());
333  relconvtol_ = pl.get("Relative Convergence Tolerance",relconvtol_);
334 
335  // maximum number of restarts
336  maxRestarts_ = pl.get("Maximum Restarts",maxRestarts_);
337 
338  // block size: default is 1
339  blockSize_ = pl.get("Block Size",1);
340  TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
341  "Anasazi::BlockKrylovSchurSolMgr: \"Block Size\" must be strictly positive.");
342 
343  // set the number of blocks we need to save to compute the nev eigenvalues of interest.
344  xtra_nevBlocks_ = pl.get("Extra NEV Blocks",0);
345  if (nev%blockSize_) {
346  nevBlocks_ = nev/blockSize_ + 1;
347  } else {
348  nevBlocks_ = nev/blockSize_;
349  }
350 
351  // determine if we should use the dynamic scheme for selecting the current number of retained eigenvalues.
352  // NOTE: This employs a technique similar to ARPACK in that it increases the number of retained eigenvalues
353  // by one for every converged eigenpair to accelerate convergence.
354  if (pl.isParameter("Dynamic Extra NEV")) {
355  if (Teuchos::isParameterType<bool>(pl,"Dynamic Extra NEV")) {
356  dynXtraNev_ = pl.get("Dynamic Extra NEV",dynXtraNev_);
357  } else {
358  dynXtraNev_ = ( Teuchos::getParameter<int>(pl,"Dynamic Extra NEV") != 0 );
359  }
360  }
361 
362  numBlocks_ = pl.get("Num Blocks",3*nevBlocks_);
363  TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= nevBlocks_, std::invalid_argument,
364  "Anasazi::BlockKrylovSchurSolMgr: \"Num Blocks\" must be strictly positive and large enough to compute the requested eigenvalues.");
365 
366  TEUCHOS_TEST_FOR_EXCEPTION(static_cast<ptrdiff_t>(numBlocks_)*blockSize_ > MVT::GetGlobalLength(*problem_->getInitVec()),
367  std::invalid_argument,
368  "Anasazi::BlockKrylovSchurSolMgr: Potentially impossible orthogonality requests. Reduce basis size.");
369 
370  // step size: the default is maxRestarts_*numBlocks_, so that Ritz values are only computed every restart.
371  if (maxRestarts_) {
372  stepSize_ = pl.get("Step Size", (maxRestarts_+1)*(numBlocks_+1));
373  } else {
374  stepSize_ = pl.get("Step Size", numBlocks_+1);
375  }
376  TEUCHOS_TEST_FOR_EXCEPTION(stepSize_ < 1, std::invalid_argument,
377  "Anasazi::BlockKrylovSchurSolMgr: \"Step Size\" must be strictly positive.");
378 
379  // get the sort manager
380  if (pl.isParameter("Sort Manager")) {
381  sort_ = Teuchos::getParameter<Teuchos::RCP<Anasazi::SortManager<MagnitudeType> > >(pl,"Sort Manager");
382  } else {
383  // which values to solve for
384  whch_ = pl.get("Which",whch_);
385  TEUCHOS_TEST_FOR_EXCEPTION(whch_ != "SM" && whch_ != "LM" && whch_ != "SR" && whch_ != "LR" && whch_ != "SI" && whch_ != "LI",
386  std::invalid_argument, "Invalid sorting string.");
387  sort_ = Teuchos::rcp( new BasicSort<MagnitudeType>(whch_) );
388  }
389 
390  // which orthogonalization to use
391  ortho_ = pl.get("Orthogonalization",ortho_);
392  if (ortho_ != "DGKS" && ortho_ != "SVQB" && ortho_ != "ICGS") {
393  ortho_ = "SVQB";
394  }
395 
396  // which orthogonalization constant to use
397  ortho_kappa_ = pl.get("Orthogonalization Constant",ortho_kappa_);
398 
399  // Create a formatted output stream to print to.
400  // See if user requests output processor.
401  osProc_ = pl.get("Output Processor", osProc_);
402 
403  // If not passed in by user, it will be chosen based upon operator type.
404  if (pl.isParameter("Output Stream")) {
405  osp_ = Teuchos::getParameter<Teuchos::RCP<Teuchos::FancyOStream> >(pl,"Output Stream");
406  }
407  else {
408  osp_ = OutputStreamTraits<OP>::getOutputStream (*problem_->getOperator(),osProc_);
409  }
410 
411  // verbosity level
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 
420  // restarting technique: V*Q or applyHouse(V,H,tau)
421  if (pl.isParameter("In Situ Restarting")) {
422  if (Teuchos::isParameterType<bool>(pl,"In Situ Restarting")) {
423  inSituRestart_ = pl.get("In Situ Restarting",inSituRestart_);
424  } else {
425  inSituRestart_ = ( Teuchos::getParameter<int>(pl,"In Situ Restarting") != 0 );
426  }
427  }
428 
429  printNum_ = pl.get<int>("Print Number of Ritz Values",printNum_);
430 }
431 
432 
433 // solve()
434 template<class ScalarType, class MV, class OP>
437 
438  const int nev = problem_->getNEV();
439  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
440  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
441 
444 
446  // Output manager
448 
450  // Status tests
451  //
452  // convergence
454  if (globalTest_ == Teuchos::null) {
455  convtest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(convtol_,nev,RITZRES_2NORM,relconvtol_) );
456  }
457  else {
458  convtest = globalTest_;
459  }
461  = Teuchos::rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sort_,nev) );
462  // for a non-short-circuited OR test, the order doesn't matter
464  alltests.push_back(ordertest);
465 
466  if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
467 
470  // printing StatusTest
472  if ( printer->isVerbosity(Debug) ) {
473  outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,Passed+Failed+Undefined ) );
474  }
475  else {
476  outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,Passed ) );
477  }
478 
480  // Orthomanager
482  if (ortho_=="SVQB") {
483  ortho = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
484  } else if (ortho_=="DGKS") {
485  if (ortho_kappa_ <= 0) {
486  ortho = Teuchos::rcp( new BasicOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
487  } else {
488  ortho = Teuchos::rcp( new BasicOrthoManager<ScalarType,MV,OP>(problem_->getM(),ortho_kappa_) );
489  }
490  } else if (ortho_=="ICGS") {
491  ortho = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
492  } else {
493  TEUCHOS_TEST_FOR_EXCEPTION(ortho_!="SVQB"&&ortho_!="DGKS"&&ortho_!="ICGS",std::logic_error,"Anasazi::BlockKrylovSchurSolMgr::solve(): Invalid orthogonalization type.");
494  }
495 
497  // Parameter list
499  plist.set("Block Size",blockSize_);
500  plist.set("Num Blocks",numBlocks_);
501  plist.set("Step Size",stepSize_);
502  if (printNum_ == -1) {
503  plist.set("Print Number of Ritz Values",nevBlocks_*blockSize_);
504  }
505  else {
506  plist.set("Print Number of Ritz Values",printNum_);
507  }
508 
510  // BlockKrylovSchur solver
512  = Teuchos::rcp( new BlockKrylovSchur<ScalarType,MV,OP>(problem_,sort_,printer,outputtest,ortho,plist) );
513  // set any auxiliary vectors defined in the problem
514  Teuchos::RCP< const MV > probauxvecs = problem_->getAuxVecs();
515  if (probauxvecs != Teuchos::null) {
516  bks_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
517  }
518 
519  // Create workspace for the Krylov basis generated during a restart
520  // Need at most (nevBlocks_*blockSize_+1) for the updated factorization and another block for the current factorization residual block (F).
521  // ---> (nevBlocks_*blockSize_+1) + blockSize_
522  // If Hermitian, this becomes nevBlocks_*blockSize_ + blockSize_
523  // we only need this if there is the possibility of restarting, ex situ
524 
525  // Maximum allowable extra vectors that BKS can keep to accelerate convergence
526  int maxXtraBlocks = 0;
527  if ( dynXtraNev_ ) maxXtraBlocks = ( ( bks_solver->getMaxSubspaceDim() - nev ) / blockSize_ ) / 2;
528 
529  Teuchos::RCP<MV> workMV;
530  if (maxRestarts_ > 0) {
531  if (inSituRestart_==true) {
532  // still need one work vector for applyHouse()
533  workMV = MVT::Clone( *problem_->getInitVec(), 1 );
534  }
535  else { // inSituRestart == false
536  if (problem_->isHermitian()) {
537  workMV = MVT::Clone( *problem_->getInitVec(), (nevBlocks_+maxXtraBlocks)*blockSize_ + blockSize_ );
538  } else {
539  // Increase workspace by 1 because of potential complex conjugate pairs on the Ritz vector boundary
540  workMV = MVT::Clone( *problem_->getInitVec(), (nevBlocks_+maxXtraBlocks)*blockSize_ + 1 + blockSize_ );
541  }
542  }
543  } else {
544  workMV = Teuchos::null;
545  }
546 
547  // go ahead and initialize the solution to nothing in case we throw an exception
549  sol.numVecs = 0;
550  problem_->setSolution(sol);
551 
552  int numRestarts = 0;
553  int cur_nevBlocks = 0;
554 
555  // enter solve() iterations
556  {
557 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
558  Teuchos::TimeMonitor slvtimer(*timerSolve_);
559 #endif
560 
561  // tell bks_solver to iterate
562  while (1) {
563  try {
564  bks_solver->iterate();
565 
567  //
568  // check convergence first
569  //
571  if ( ordertest->getStatus() == Passed ) {
572  // we have convergence
573  // ordertest->whichVecs() tells us which vectors from solver state are the ones we want
574  // ordertest->howMany() will tell us how many
575  break;
576  }
578  //
579  // check for restarting, i.e. the subspace is full
580  //
582  // this is for the Hermitian case, or non-Hermitian conjugate split situation.
583  // --> for the Hermitian case the current subspace dimension needs to match the maximum subspace dimension
584  // --> for the non-Hermitian case:
585  // --> if a conjugate pair was detected in the previous restart then the current subspace dimension needs to match the
586  // maximum subspace dimension (the BKS solver keeps one extra vector if the problem is non-Hermitian).
587  // --> if a conjugate pair was not detected in the previous restart then the current subspace dimension will be one less
588  // than the maximum subspace dimension.
589  else if ( (bks_solver->getCurSubspaceDim() == bks_solver->getMaxSubspaceDim()) ||
590  (!problem_->isHermitian() && !conjSplit_ && (bks_solver->getCurSubspaceDim()+1 == bks_solver->getMaxSubspaceDim())) ) {
591 
592  // Update the Schur form of the projected eigenproblem, then sort it.
593  if (!bks_solver->isSchurCurrent()) {
594  bks_solver->computeSchurForm( true );
595 
596  // Check for convergence, just in case we wait for every restart to check
597  outputtest->checkStatus( &*bks_solver );
598  }
599 
600  // Don't bother to restart if we've converged or reached the maximum number of restarts
601  if ( numRestarts >= maxRestarts_ || ordertest->getStatus() == Passed) {
602  break; // break from while(1){bks_solver->iterate()}
603  }
604 
605  // Start restarting timer and increment counter
606 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
607  Teuchos::TimeMonitor restimer(*timerRestarting_);
608 #endif
609  numRestarts++;
610 
611  int numConv = ordertest->howMany();
612  cur_nevBlocks = nevBlocks_*blockSize_;
613 
614  // Add in extra blocks for restarting if either static or dynamic boundaries are being used.
615  int moreNevBlocks = std::min( maxXtraBlocks, std::max( numConv/blockSize_, xtra_nevBlocks_) );
616  if ( dynXtraNev_ )
617  cur_nevBlocks += moreNevBlocks * blockSize_;
618  else if ( xtra_nevBlocks_ )
619  cur_nevBlocks += xtra_nevBlocks_ * blockSize_;
620 /*
621  int cur_numConv = numConv;
622  while ( (cur_nevBlocks < (_nevBlocks + maxXtraVecs)) && cur_numConv > 0 ) {
623  cur_nevBlocks++;
624  cur_numConv--;
625 */
626 
627  printer->stream(Debug) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << std::endl << std::endl;
628  printer->stream(Debug) << " - Current NEV blocks is " << cur_nevBlocks << ", the minimum is " << nevBlocks_*blockSize_ << std::endl;
629 
630  // Get the most current Ritz values before we continue.
631  ritzValues_ = bks_solver->getRitzValues();
632 
633  // Get the state.
634  BlockKrylovSchurState<ScalarType,MV> oldState = bks_solver->getState();
635 
636  // Get the current dimension of the factorization
637  int curDim = oldState.curDim;
638 
639  // Determine if the storage for the nev eigenvalues of interest splits a complex conjugate pair.
640  std::vector<int> ritzIndex = bks_solver->getRitzIndex();
641  if (ritzIndex[cur_nevBlocks-1]==1) {
642  conjSplit_ = true;
643  cur_nevBlocks++;
644  } else {
645  conjSplit_ = false;
646  }
647 
648  // Print out a warning to the user if complex eigenvalues were found on the boundary of the restart subspace
649  // and the eigenproblem is Hermitian. This solver is not prepared to handle this situation.
650  if (problem_->isHermitian() && conjSplit_)
651  {
652  printer->stream(Warnings)
653  << " Eigenproblem is Hermitian, complex eigenvalues have been detected, and eigenvalues of interest split a conjugate pair!!!"
654  << std::endl
655  << " Block Krylov-Schur eigensolver cannot guarantee correct behavior in this situation, please turn Hermitian flag off!!!"
656  << std::endl;
657  }
658 
659  // Update the Krylov-Schur decomposition
660 
661  // Get a view of the Schur vectors of interest.
662  Teuchos::SerialDenseMatrix<int,ScalarType> Qnev(Teuchos::View, *(oldState.Q), curDim, cur_nevBlocks);
663 
664  // Get a view of the current Krylov basis.
665  std::vector<int> curind( curDim );
666  for (int i=0; i<curDim; i++) { curind[i] = i; }
667  Teuchos::RCP<const MV> basistemp = MVT::CloneView( *(oldState.V), curind );
668 
669  // Compute the new Krylov basis: Vnew = V*Qnev
670  //
671  // this will occur ex situ in workspace allocated for this purpose (tmpMV)
672  // or in situ in the solver's memory space.
673  //
674  // we will also set a pointer for the location that the current factorization residual block (F),
675  // currently located after the current basis in oldstate.V, will be moved to
676  //
677  Teuchos::RCP<MV> newF;
678  if (inSituRestart_) {
679  //
680  // get non-const pointer to solver's basis so we can work in situ
681  Teuchos::RCP<MV> solverbasis = Teuchos::rcp_const_cast<MV>(oldState.V);
683  //
684  // perform Householder QR of copyQnev = Q [D;0], where D is unit diag. We will want D below.
685  std::vector<ScalarType> tau(cur_nevBlocks), work(cur_nevBlocks);
686  int info;
687  lapack.GEQRF(curDim,cur_nevBlocks,copyQnev.values(),copyQnev.stride(),&tau[0],&work[0],work.size(),&info);
688  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,std::logic_error,
689  "Anasazi::BlockKrylovSchurSolMgr::solve(): error calling GEQRF during restarting.");
690  // we need to get the diagonal of D
691  std::vector<ScalarType> d(cur_nevBlocks);
692  for (int j=0; j<copyQnev.numCols(); j++) {
693  d[j] = copyQnev(j,j);
694  }
695  if (printer->isVerbosity(Debug)) {
696  Teuchos::SerialDenseMatrix<int,ScalarType> R(Teuchos::Copy,copyQnev,cur_nevBlocks,cur_nevBlocks);
697  for (int j=0; j<R.numCols(); j++) {
698  R(j,j) = SCT::magnitude(R(j,j)) - 1.0;
699  for (int i=j+1; i<R.numRows(); i++) {
700  R(i,j) = zero;
701  }
702  }
703  printer->stream(Debug) << "||Triangular factor of Su - I||: " << R.normFrobenius() << std::endl;
704  }
705  //
706  // perform implicit V*Qnev
707  // this actually performs V*[Qnev Qtrunc*M] = [newV truncV], for some unitary M
708  // we are interested in only the first cur_nevBlocks vectors of the result
709  curind.resize(curDim);
710  for (int i=0; i<curDim; i++) curind[i] = i;
711  {
712  Teuchos::RCP<MV> oldV = MVT::CloneViewNonConst(*solverbasis,curind);
713  SolverUtils<ScalarType,MV,OP>::applyHouse(cur_nevBlocks,*oldV,copyQnev,tau,workMV);
714  }
715  // multiply newV*D
716  // get pointer to new basis
717  curind.resize(cur_nevBlocks);
718  for (int i=0; i<cur_nevBlocks; i++) { curind[i] = i; }
719  {
720  Teuchos::RCP<MV> newV = MVT::CloneViewNonConst( *solverbasis, curind );
721  MVT::MvScale(*newV,d);
722  }
723  // get pointer to new location for F
724  curind.resize(blockSize_);
725  for (int i=0; i<blockSize_; i++) { curind[i] = cur_nevBlocks + i; }
726  newF = MVT::CloneViewNonConst( *solverbasis, curind );
727  }
728  else {
729  // get pointer to first part of work space
730  curind.resize(cur_nevBlocks);
731  for (int i=0; i<cur_nevBlocks; i++) { curind[i] = i; }
732  Teuchos::RCP<MV> tmp_newV = MVT::CloneViewNonConst(*workMV, curind );
733  // perform V*Qnev
734  MVT::MvTimesMatAddMv( one, *basistemp, Qnev, zero, *tmp_newV );
735  tmp_newV = Teuchos::null;
736  // get pointer to new location for F
737  curind.resize(blockSize_);
738  for (int i=0; i<blockSize_; i++) { curind[i] = cur_nevBlocks + i; }
739  newF = MVT::CloneViewNonConst( *workMV, curind );
740  }
741 
742  // Move the current factorization residual block (F) to the last block of newV.
743  curind.resize(blockSize_);
744  for (int i=0; i<blockSize_; i++) { curind[i] = curDim + i; }
745  Teuchos::RCP<const MV> oldF = MVT::CloneView( *(oldState.V), curind );
746  for (int i=0; i<blockSize_; i++) { curind[i] = i; }
747  MVT::SetBlock( *oldF, curind, *newF );
748  newF = Teuchos::null;
749 
750  // Update the Krylov-Schur quasi-triangular matrix.
751  //
752  // Create storage for the new Schur matrix of the Krylov-Schur factorization
753  // Copy over the current quasi-triangular factorization of oldState.H which is stored in oldState.S.
755  Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::Copy, *(oldState.S), cur_nevBlocks+blockSize_, cur_nevBlocks) );
756  //
757  // Get a view of the B block of the current factorization
758  Teuchos::SerialDenseMatrix<int,ScalarType> oldB(Teuchos::View, *(oldState.H), blockSize_, blockSize_, curDim, curDim-blockSize_);
759  //
760  // Get a view of the a block row of the Schur vectors.
761  Teuchos::SerialDenseMatrix<int,ScalarType> subQ(Teuchos::View, *(oldState.Q), blockSize_, cur_nevBlocks, curDim-blockSize_);
762  //
763  // Get a view of the new B block of the updated Krylov-Schur factorization
764  Teuchos::SerialDenseMatrix<int,ScalarType> newB(Teuchos::View, *newH, blockSize_, cur_nevBlocks, cur_nevBlocks);
765  //
766  // Compute the new B block.
767  blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, blockSize_, cur_nevBlocks, blockSize_, one,
768  oldB.values(), oldB.stride(), subQ.values(), subQ.stride(), zero, newB.values(), newB.stride() );
769 
770 
771  //
772  // Set the new state and initialize the solver.
774  if (inSituRestart_) {
775  newstate.V = oldState.V;
776  } else {
777  newstate.V = workMV;
778  }
779  newstate.H = newH;
780  newstate.curDim = cur_nevBlocks;
781  bks_solver->initialize(newstate);
782 
783  } // end of restarting
785  //
786  // we returned from iterate(), but none of our status tests Passed.
787  // something is wrong, and it is probably our fault.
788  //
790  else {
791  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::BlockKrylovSchurSolMgr::solve(): Invalid return from bks_solver::iterate().");
792  }
793  }
794  catch (const AnasaziError &err) {
795  printer->stream(Errors)
796  << "Anasazi::BlockKrylovSchurSolMgr::solve() caught unexpected exception from Anasazi::BlockKrylovSchur::iterate() at iteration " << bks_solver->getNumIters() << std::endl
797  << err.what() << std::endl
798  << "Anasazi::BlockKrylovSchurSolMgr::solve() returning Unconverged with no solutions." << std::endl;
799  return Unconverged;
800  }
801  }
802 
803  //
804  // free temporary space
805  workMV = Teuchos::null;
806 
807  // Get the most current Ritz values before we return
808  ritzValues_ = bks_solver->getRitzValues();
809 
810  sol.numVecs = ordertest->howMany();
811  printer->stream(Debug) << "ordertest->howMany() : " << sol.numVecs << std::endl;
812  std::vector<int> whichVecs = ordertest->whichVecs();
813 
814  // Place any converged eigenpairs in the solution container.
815  if (sol.numVecs > 0) {
816 
817  // Next determine if there is a conjugate pair on the boundary and resize.
818  std::vector<int> tmpIndex = bks_solver->getRitzIndex();
819  for (int i=0; i<(int)ritzValues_.size(); ++i) {
820  printer->stream(Debug) << ritzValues_[i].realpart << " + i " << ritzValues_[i].imagpart << ", Index = " << tmpIndex[i] << std::endl;
821  }
822  printer->stream(Debug) << "Number of converged eigenpairs (before) = " << sol.numVecs << std::endl;
823  for (int i=0; i<sol.numVecs; ++i) {
824  printer->stream(Debug) << "whichVecs[" << i << "] = " << whichVecs[i] << ", tmpIndex[" << whichVecs[i] << "] = " << tmpIndex[whichVecs[i]] << std::endl;
825  }
826  if (tmpIndex[whichVecs[sol.numVecs-1]]==1) {
827  printer->stream(Debug) << "There is a conjugate pair on the boundary, resizing sol.numVecs" << std::endl;
828  whichVecs.push_back(whichVecs[sol.numVecs-1]+1);
829  sol.numVecs++;
830  for (int i=0; i<sol.numVecs; ++i) {
831  printer->stream(Debug) << "whichVecs[" << i << "] = " << whichVecs[i] << ", tmpIndex[" << whichVecs[i] << "] = " << tmpIndex[whichVecs[i]] << std::endl;
832  }
833  }
834 
835  bool keepMore = false;
836  int numEvecs = sol.numVecs;
837  printer->stream(Debug) << "Number of converged eigenpairs (after) = " << sol.numVecs << std::endl;
838  printer->stream(Debug) << "whichVecs[sol.numVecs-1] > sol.numVecs-1 : " << whichVecs[sol.numVecs-1] << " > " << sol.numVecs-1 << std::endl;
839  if (whichVecs[sol.numVecs-1] > (sol.numVecs-1)) {
840  keepMore = true;
841  numEvecs = whichVecs[sol.numVecs-1]+1; // Add 1 to fix zero-based indexing
842  printer->stream(Debug) << "keepMore = true; numEvecs = " << numEvecs << std::endl;
843  }
844 
845  // Next set the number of Ritz vectors that the iteration must compute and compute them.
846  bks_solver->setNumRitzVectors(numEvecs);
847  bks_solver->computeRitzVectors();
848 
849  // If the leading Ritz pairs are the converged ones, get the information
850  // from the iteration to the solution container. Otherwise copy the necessary
851  // information using 'whichVecs'.
852  if (!keepMore) {
853  sol.index = bks_solver->getRitzIndex();
854  sol.Evals = bks_solver->getRitzValues();
855  sol.Evecs = MVT::CloneCopy( *(bks_solver->getRitzVectors()) );
856  }
857 
858  // Resize based on the number of solutions being returned and set the number of Ritz
859  // vectors for the iteration to compute.
860  sol.Evals.resize(sol.numVecs);
861  sol.index.resize(sol.numVecs);
862 
863  // If the converged Ritz pairs are not the leading ones, copy over the information directly.
864  if (keepMore) {
865  std::vector<Anasazi::Value<ScalarType> > tmpEvals = bks_solver->getRitzValues();
866  for (int vec_i=0; vec_i<sol.numVecs; ++vec_i) {
867  sol.index[vec_i] = tmpIndex[whichVecs[vec_i]];
868  sol.Evals[vec_i] = tmpEvals[whichVecs[vec_i]];
869  }
870  sol.Evecs = MVT::CloneCopy( *(bks_solver->getRitzVectors()), whichVecs );
871  }
872 
873  // Set the solution space to be the Ritz vectors at this time.
874  sol.Espace = sol.Evecs;
875  }
876  }
877 
878  // print final summary
879  bks_solver->currentStatus(printer->stream(FinalSummary));
880 
881  // print timing information
882 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
883  if ( printer->isVerbosity( TimingDetails ) ) {
884  Teuchos::TimeMonitor::summarize( printer->stream( TimingDetails ) );
885  }
886 #endif
887 
888  problem_->setSolution(sol);
889  printer->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
890 
891  // get the number of iterations performed during this solve.
892  numIters_ = bks_solver->getNumIters();
893 
894  if (sol.numVecs < nev) {
895  return Unconverged; // return from BlockKrylovSchurSolMgr::solve()
896  }
897  return Converged; // return from BlockKrylovSchurSolMgr::solve()
898 }
899 
900 
901 template <class ScalarType, class MV, class OP>
902 void
905 {
906  globalTest_ = global;
907 }
908 
909 template <class ScalarType, class MV, class OP>
912 {
913  return globalTest_;
914 }
915 
916 template <class ScalarType, class MV, class OP>
917 void
920 {
921  debugTest_ = debug;
922 }
923 
924 template <class ScalarType, class MV, class OP>
927 {
928  return debugTest_;
929 }
930 
931 } // end Anasazi namespace
932 
933 #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.