Anasazi  Version of the Day
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, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25 //
26 // ***********************************************************************
27 // @HEADER
28 
29 #ifndef ANASAZI_BLOCK_KRYLOV_SCHUR_SOLMGR_HPP
30 #define ANASAZI_BLOCK_KRYLOV_SCHUR_SOLMGR_HPP
31 
35 
36 #include "AnasaziConfigDefs.hpp"
37 #include "AnasaziTypes.hpp"
38 
39 #include "AnasaziEigenproblem.hpp"
40 #include "AnasaziSolverManager.hpp"
41 #include "AnasaziSolverUtils.hpp"
42 
44 #include "AnasaziBasicSort.hpp"
52 #include "Teuchos_BLAS.hpp"
53 #include "Teuchos_LAPACK.hpp"
54 #include "Teuchos_TimeMonitor.hpp"
55 
62 
110 namespace Anasazi {
111 
112 
139 template<class ScalarType, class MV, class OP>
140 class BlockKrylovSchurSolMgr : public SolverManager<ScalarType,MV,OP> {
141 
142  private:
146  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
148 
149  public:
150 
152 
153 
173 
177 
179 
180 
183  return *_problem;
184  }
185 
187  int getNumIters() const {
188  return _numIters;
189  }
190 
193  std::vector<Value<ScalarType> > getRitzValues() const {
194  std::vector<Value<ScalarType> > ret( _ritzValues );
195  return ret;
196  }
197 
205  return Teuchos::tuple(_timerSolve, _timerRestarting);
206  }
207 
209 
211 
212 
231  ReturnType solve();
232 
235 
238 
241 
244 
246 
247  private:
250 
251  std::string _whch, _ortho;
252  MagnitudeType _ortho_kappa;
253 
254  MagnitudeType _convtol;
255  int _maxRestarts;
256  bool _relconvtol,_conjSplit;
257  int _blockSize, _numBlocks, _stepSize, _nevBlocks, _xtra_nevBlocks;
258  int _numIters;
259  int _verbosity;
260  bool _inSituRestart, _dynXtraNev;
261 
262  std::vector<Value<ScalarType> > _ritzValues;
263 
264  int _printNum;
265  Teuchos::RCP<Teuchos::Time> _timerSolve, _timerRestarting;
266 
269 
270 };
271 
272 
273 // Constructor
274 template<class ScalarType, class MV, class OP>
277  Teuchos::ParameterList &pl ) :
278  _problem(problem),
279  _whch("LM"),
280  _ortho("SVQB"),
281  _ortho_kappa(-1.0),
282  _convtol(0),
283  _maxRestarts(20),
284  _relconvtol(true),
285  _conjSplit(false),
286  _blockSize(0),
287  _numBlocks(0),
288  _stepSize(0),
289  _nevBlocks(0),
290  _xtra_nevBlocks(0),
291  _numIters(0),
292  _verbosity(Anasazi::Errors),
293  _inSituRestart(false),
294  _dynXtraNev(false),
295  _printNum(-1)
296 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
297  ,_timerSolve(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchurSolMgr::solve()")),
298  _timerRestarting(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchurSolMgr restarting"))
299 #endif
300 {
301  TEUCHOS_TEST_FOR_EXCEPTION(_problem == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
302  TEUCHOS_TEST_FOR_EXCEPTION(!_problem->isProblemSet(), std::invalid_argument, "Problem not set.");
303  TEUCHOS_TEST_FOR_EXCEPTION(_problem->getInitVec() == Teuchos::null, std::invalid_argument, "Problem does not contain initial vectors to clone from.");
304 
305  const int nev = _problem->getNEV();
306 
307  // convergence tolerance
308  _convtol = pl.get("Convergence Tolerance",MT::prec());
309  _relconvtol = pl.get("Relative Convergence Tolerance",_relconvtol);
310 
311  // maximum number of restarts
312  _maxRestarts = pl.get("Maximum Restarts",_maxRestarts);
313 
314  // block size: default is 1
315  _blockSize = pl.get("Block Size",1);
316  TEUCHOS_TEST_FOR_EXCEPTION(_blockSize <= 0, std::invalid_argument,
317  "Anasazi::BlockKrylovSchurSolMgr: \"Block Size\" must be strictly positive.");
318 
319  // set the number of blocks we need to save to compute the nev eigenvalues of interest.
320  _xtra_nevBlocks = pl.get("Extra NEV Blocks",0);
321  if (nev%_blockSize) {
322  _nevBlocks = nev/_blockSize + 1;
323  } else {
324  _nevBlocks = nev/_blockSize;
325  }
326 
327  // determine if we should use the dynamic scheme for selecting the current number of retained eigenvalues.
328  // NOTE: This employs a technique similar to ARPACK in that it increases the number of retained eigenvalues
329  // by one for every converged eigenpair to accelerate convergence.
330  if (pl.isParameter("Dynamic Extra NEV")) {
331  if (Teuchos::isParameterType<bool>(pl,"Dynamic Extra NEV")) {
332  _dynXtraNev = pl.get("Dynamic Extra NEV",_dynXtraNev);
333  } else {
334  _dynXtraNev = ( Teuchos::getParameter<int>(pl,"Dynamic Extra NEV") != 0 );
335  }
336  }
337 
338  _numBlocks = pl.get("Num Blocks",3*_nevBlocks);
339  TEUCHOS_TEST_FOR_EXCEPTION(_numBlocks <= _nevBlocks, std::invalid_argument,
340  "Anasazi::BlockKrylovSchurSolMgr: \"Num Blocks\" must be strictly positive and large enough to compute the requested eigenvalues.");
341 
342  TEUCHOS_TEST_FOR_EXCEPTION(static_cast<ptrdiff_t>(_numBlocks)*_blockSize > MVT::GetGlobalLength(*_problem->getInitVec()),
343  std::invalid_argument,
344  "Anasazi::BlockKrylovSchurSolMgr: Potentially impossible orthogonality requests. Reduce basis size.");
345 
346  // step size: the default is _maxRestarts*_numBlocks, so that Ritz values are only computed every restart.
347  if (_maxRestarts) {
348  _stepSize = pl.get("Step Size", (_maxRestarts+1)*(_numBlocks+1));
349  } else {
350  _stepSize = pl.get("Step Size", _numBlocks+1);
351  }
352  TEUCHOS_TEST_FOR_EXCEPTION(_stepSize < 1, std::invalid_argument,
353  "Anasazi::BlockKrylovSchurSolMgr: \"Step Size\" must be strictly positive.");
354 
355  // get the sort manager
356  if (pl.isParameter("Sort Manager")) {
357  _sort = Teuchos::getParameter<Teuchos::RCP<Anasazi::SortManager<MagnitudeType> > >(pl,"Sort Manager");
358  } else {
359  // which values to solve for
360  _whch = pl.get("Which",_whch);
361  TEUCHOS_TEST_FOR_EXCEPTION(_whch != "SM" && _whch != "LM" && _whch != "SR" && _whch != "LR" && _whch != "SI" && _whch != "LI",
362  std::invalid_argument, "Invalid sorting string.");
363  _sort = Teuchos::rcp( new BasicSort<MagnitudeType>(_whch) );
364  }
365 
366  // which orthogonalization to use
367  _ortho = pl.get("Orthogonalization",_ortho);
368  if (_ortho != "DGKS" && _ortho != "SVQB") {
369  _ortho = "SVQB";
370  }
371 
372  // which orthogonalization constant to use
373  _ortho_kappa = pl.get("Orthogonalization Constant",_ortho_kappa);
374 
375  // verbosity level
376  if (pl.isParameter("Verbosity")) {
377  if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
378  _verbosity = pl.get("Verbosity", _verbosity);
379  } else {
380  _verbosity = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
381  }
382  }
383 
384  // restarting technique: V*Q or applyHouse(V,H,tau)
385  if (pl.isParameter("In Situ Restarting")) {
386  if (Teuchos::isParameterType<bool>(pl,"In Situ Restarting")) {
387  _inSituRestart = pl.get("In Situ Restarting",_inSituRestart);
388  } else {
389  _inSituRestart = ( Teuchos::getParameter<int>(pl,"In Situ Restarting") != 0 );
390  }
391  }
392 
393  _printNum = pl.get<int>("Print Number of Ritz Values",-1);
394 }
395 
396 
397 // solve()
398 template<class ScalarType, class MV, class OP>
401 
402  const int nev = _problem->getNEV();
403  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
404  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
405 
408  typedef SolverUtils<ScalarType,MV,OP> msutils;
409 
411  // Output manager
413 
415  // Status tests
416  //
417  // convergence
419  if (globalTest_ == Teuchos::null) {
420  convtest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(_convtol,nev,RITZRES_2NORM,_relconvtol) );
421  }
422  else {
423  convtest = globalTest_;
424  }
426  = Teuchos::rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,_sort,nev) );
427  // for a non-short-circuited OR test, the order doesn't matter
429  alltests.push_back(ordertest);
430 
431  if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
432 
435  // printing StatusTest
437  if ( printer->isVerbosity(Debug) ) {
438  outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,Passed+Failed+Undefined ) );
439  }
440  else {
441  outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,Passed ) );
442  }
443 
445  // Orthomanager
447  if (_ortho=="SVQB") {
448  ortho = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>(_problem->getM()) );
449  } else if (_ortho=="DGKS") {
450  if (_ortho_kappa <= 0) {
451  ortho = Teuchos::rcp( new BasicOrthoManager<ScalarType,MV,OP>(_problem->getM()) );
452  }
453  else {
454  ortho = Teuchos::rcp( new BasicOrthoManager<ScalarType,MV,OP>(_problem->getM(),_ortho_kappa) );
455  }
456  } else {
457  TEUCHOS_TEST_FOR_EXCEPTION(_ortho!="SVQB"&&_ortho!="DGKS",std::logic_error,"Anasazi::BlockKrylovSchurSolMgr::solve(): Invalid orthogonalization type.");
458  }
459 
461  // Parameter list
463  plist.set("Block Size",_blockSize);
464  plist.set("Num Blocks",_numBlocks);
465  plist.set("Step Size",_stepSize);
466  if (_printNum == -1) {
467  plist.set("Print Number of Ritz Values",_nevBlocks*_blockSize);
468  }
469  else {
470  plist.set("Print Number of Ritz Values",_printNum);
471  }
472 
474  // BlockKrylovSchur solver
476  = Teuchos::rcp( new BlockKrylovSchur<ScalarType,MV,OP>(_problem,_sort,printer,outputtest,ortho,plist) );
477  // set any auxiliary vectors defined in the problem
478  Teuchos::RCP< const MV > probauxvecs = _problem->getAuxVecs();
479  if (probauxvecs != Teuchos::null) {
480  bks_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
481  }
482 
483  // Create workspace for the Krylov basis generated during a restart
484  // Need at most (_nevBlocks*_blockSize+1) for the updated factorization and another block for the current factorization residual block (F).
485  // ---> (_nevBlocks*_blockSize+1) + _blockSize
486  // If Hermitian, this becomes _nevBlocks*_blockSize + _blockSize
487  // we only need this if there is the possibility of restarting, ex situ
488 
489  // Maximum allowable extra vectors that BKS can keep to accelerate convergence
490  int maxXtraBlocks = 0;
491  if ( _dynXtraNev ) maxXtraBlocks = ( ( bks_solver->getMaxSubspaceDim() - nev ) / _blockSize ) / 2;
492 
493  Teuchos::RCP<MV> workMV;
494  if (_maxRestarts > 0) {
495  if (_inSituRestart==true) {
496  // still need one work vector for applyHouse()
497  workMV = MVT::Clone( *_problem->getInitVec(), 1 );
498  }
499  else { // inSituRestart == false
500  if (_problem->isHermitian()) {
501  workMV = MVT::Clone( *_problem->getInitVec(), (_nevBlocks+maxXtraBlocks)*_blockSize + _blockSize );
502  } else {
503  // Increase workspace by 1 because of potential complex conjugate pairs on the Ritz vector boundary
504  workMV = MVT::Clone( *_problem->getInitVec(), (_nevBlocks+maxXtraBlocks)*_blockSize + 1 + _blockSize );
505  }
506  }
507  } else {
508  workMV = Teuchos::null;
509  }
510 
511  // go ahead and initialize the solution to nothing in case we throw an exception
513  sol.numVecs = 0;
514  _problem->setSolution(sol);
515 
516  int numRestarts = 0;
517  int cur_nevBlocks = 0;
518 
519  // enter solve() iterations
520  {
521 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
522  Teuchos::TimeMonitor slvtimer(*_timerSolve);
523 #endif
524 
525  // tell bks_solver to iterate
526  while (1) {
527  try {
528  bks_solver->iterate();
529 
531  //
532  // check convergence first
533  //
535  if ( ordertest->getStatus() == Passed ) {
536  // we have convergence
537  // ordertest->whichVecs() tells us which vectors from solver state are the ones we want
538  // ordertest->howMany() will tell us how many
539  break;
540  }
542  //
543  // check for restarting, i.e. the subspace is full
544  //
546  // this is for the Hermitian case, or non-Hermitian conjugate split situation.
547  // --> for the Hermitian case the current subspace dimension needs to match the maximum subspace dimension
548  // --> for the non-Hermitian case:
549  // --> if a conjugate pair was detected in the previous restart then the current subspace dimension needs to match the
550  // maximum subspace dimension (the BKS solver keeps one extra vector if the problem is non-Hermitian).
551  // --> if a conjugate pair was not detected in the previous restart then the current subspace dimension will be one less
552  // than the maximum subspace dimension.
553  else if ( (bks_solver->getCurSubspaceDim() == bks_solver->getMaxSubspaceDim()) ||
554  (!_problem->isHermitian() && !_conjSplit && (bks_solver->getCurSubspaceDim()+1 == bks_solver->getMaxSubspaceDim())) ) {
555 
556  // Update the Schur form of the projected eigenproblem, then sort it.
557  if (!bks_solver->isSchurCurrent()) {
558  bks_solver->computeSchurForm( true );
559 
560  // Check for convergence, just in case we wait for every restart to check
561  outputtest->checkStatus( &*bks_solver );
562  }
563 
564  // Don't bother to restart if we've converged or reached the maximum number of restarts
565  if ( numRestarts >= _maxRestarts || ordertest->getStatus() == Passed) {
566  break; // break from while(1){bks_solver->iterate()}
567  }
568 
569  // Start restarting timer and increment counter
570 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
571  Teuchos::TimeMonitor restimer(*_timerRestarting);
572 #endif
573  numRestarts++;
574 
575  int numConv = ordertest->howMany();
576  cur_nevBlocks = _nevBlocks*_blockSize;
577 
578  // Add in extra blocks for restarting if either static or dynamic boundaries are being used.
579  int moreNevBlocks = std::min( maxXtraBlocks, std::max( numConv/_blockSize, _xtra_nevBlocks) );
580  if ( _dynXtraNev )
581  cur_nevBlocks += moreNevBlocks * _blockSize;
582  else if ( _xtra_nevBlocks )
583  cur_nevBlocks += _xtra_nevBlocks * _blockSize;
584 /*
585  int cur_numConv = numConv;
586  while ( (cur_nevBlocks < (_nevBlocks + maxXtraVecs)) && cur_numConv > 0 ) {
587  cur_nevBlocks++;
588  cur_numConv--;
589 */
590 
591  printer->stream(Debug) << " Performing restart number " << numRestarts << " of " << _maxRestarts << std::endl << std::endl;
592  printer->stream(Debug) << " - Current NEV blocks is " << cur_nevBlocks << ", the minimum is " << _nevBlocks*_blockSize << std::endl;
593 
594  // Get the most current Ritz values before we continue.
595  _ritzValues = bks_solver->getRitzValues();
596 
597  // Get the state.
598  BlockKrylovSchurState<ScalarType,MV> oldState = bks_solver->getState();
599 
600  // Get the current dimension of the factorization
601  int curDim = oldState.curDim;
602 
603  // Determine if the storage for the nev eigenvalues of interest splits a complex conjugate pair.
604  std::vector<int> ritzIndex = bks_solver->getRitzIndex();
605  if (ritzIndex[cur_nevBlocks-1]==1) {
606  _conjSplit = true;
607  cur_nevBlocks++;
608  } else {
609  _conjSplit = false;
610  }
611 
612  // Print out a warning to the user if complex eigenvalues were found on the boundary of the restart subspace
613  // and the eigenproblem is Hermitian. This solver is not prepared to handle this situation.
614  if (_problem->isHermitian() && _conjSplit)
615  {
616  printer->stream(Warnings)
617  << " Eigenproblem is Hermitian, complex eigenvalues have been detected, and eigenvalues of interest split a conjugate pair!!!"
618  << std::endl
619  << " Block Krylov-Schur eigensolver cannot guarantee correct behavior in this situation, please turn Hermitian flag off!!!"
620  << std::endl;
621  }
622 
623  // Update the Krylov-Schur decomposition
624 
625  // Get a view of the Schur vectors of interest.
626  Teuchos::SerialDenseMatrix<int,ScalarType> Qnev(Teuchos::View, *(oldState.Q), curDim, cur_nevBlocks);
627 
628  // Get a view of the current Krylov basis.
629  std::vector<int> curind( curDim );
630  for (int i=0; i<curDim; i++) { curind[i] = i; }
631  Teuchos::RCP<const MV> basistemp = MVT::CloneView( *(oldState.V), curind );
632 
633  // Compute the new Krylov basis: Vnew = V*Qnev
634  //
635  // this will occur ex situ in workspace allocated for this purpose (tmpMV)
636  // or in situ in the solver's memory space.
637  //
638  // we will also set a pointer for the location that the current factorization residual block (F),
639  // currently located after the current basis in oldstate.V, will be moved to
640  //
641  Teuchos::RCP<MV> newF;
642  if (_inSituRestart) {
643  //
644  // get non-const pointer to solver's basis so we can work in situ
645  Teuchos::RCP<MV> solverbasis = Teuchos::rcp_const_cast<MV>(oldState.V);
647  //
648  // perform Householder QR of copyQnev = Q [D;0], where D is unit diag. We will want D below.
649  std::vector<ScalarType> tau(cur_nevBlocks), work(cur_nevBlocks);
650  int info;
651  lapack.GEQRF(curDim,cur_nevBlocks,copyQnev.values(),copyQnev.stride(),&tau[0],&work[0],work.size(),&info);
652  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,std::logic_error,
653  "Anasazi::BlockKrylovSchurSolMgr::solve(): error calling GEQRF during restarting.");
654  // we need to get the diagonal of D
655  std::vector<ScalarType> d(cur_nevBlocks);
656  for (int j=0; j<copyQnev.numCols(); j++) {
657  d[j] = copyQnev(j,j);
658  }
659  if (printer->isVerbosity(Debug)) {
660  Teuchos::SerialDenseMatrix<int,ScalarType> R(Teuchos::Copy,copyQnev,cur_nevBlocks,cur_nevBlocks);
661  for (int j=0; j<R.numCols(); j++) {
662  R(j,j) = SCT::magnitude(R(j,j)) - 1.0;
663  for (int i=j+1; i<R.numRows(); i++) {
664  R(i,j) = zero;
665  }
666  }
667  printer->stream(Debug) << "||Triangular factor of Su - I||: " << R.normFrobenius() << std::endl;
668  }
669  //
670  // perform implicit V*Qnev
671  // this actually performs V*[Qnev Qtrunc*M] = [newV truncV], for some unitary M
672  // we are interested in only the first cur_nevBlocks vectors of the result
673  curind.resize(curDim);
674  for (int i=0; i<curDim; i++) curind[i] = i;
675  {
676  Teuchos::RCP<MV> oldV = MVT::CloneViewNonConst(*solverbasis,curind);
677  msutils::applyHouse(cur_nevBlocks,*oldV,copyQnev,tau,workMV);
678  }
679  // multiply newV*D
680  // get pointer to new basis
681  curind.resize(cur_nevBlocks);
682  for (int i=0; i<cur_nevBlocks; i++) { curind[i] = i; }
683  {
684  Teuchos::RCP<MV> newV = MVT::CloneViewNonConst( *solverbasis, curind );
685  MVT::MvScale(*newV,d);
686  }
687  // get pointer to new location for F
688  curind.resize(_blockSize);
689  for (int i=0; i<_blockSize; i++) { curind[i] = cur_nevBlocks + i; }
690  newF = MVT::CloneViewNonConst( *solverbasis, curind );
691  }
692  else {
693  // get pointer to first part of work space
694  curind.resize(cur_nevBlocks);
695  for (int i=0; i<cur_nevBlocks; i++) { curind[i] = i; }
696  Teuchos::RCP<MV> tmp_newV = MVT::CloneViewNonConst(*workMV, curind );
697  // perform V*Qnev
698  MVT::MvTimesMatAddMv( one, *basistemp, Qnev, zero, *tmp_newV );
699  tmp_newV = Teuchos::null;
700  // get pointer to new location for F
701  curind.resize(_blockSize);
702  for (int i=0; i<_blockSize; i++) { curind[i] = cur_nevBlocks + i; }
703  newF = MVT::CloneViewNonConst( *workMV, curind );
704  }
705 
706  // Move the current factorization residual block (F) to the last block of newV.
707  curind.resize(_blockSize);
708  for (int i=0; i<_blockSize; i++) { curind[i] = curDim + i; }
709  Teuchos::RCP<const MV> oldF = MVT::CloneView( *(oldState.V), curind );
710  for (int i=0; i<_blockSize; i++) { curind[i] = i; }
711  MVT::SetBlock( *oldF, curind, *newF );
712  newF = Teuchos::null;
713 
714  // Update the Krylov-Schur quasi-triangular matrix.
715  //
716  // Create storage for the new Schur matrix of the Krylov-Schur factorization
717  // Copy over the current quasi-triangular factorization of oldState.H which is stored in oldState.S.
719  Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::Copy, *(oldState.S), cur_nevBlocks+_blockSize, cur_nevBlocks) );
720  //
721  // Get a view of the B block of the current factorization
722  Teuchos::SerialDenseMatrix<int,ScalarType> oldB(Teuchos::View, *(oldState.H), _blockSize, _blockSize, curDim, curDim-_blockSize);
723  //
724  // Get a view of the a block row of the Schur vectors.
725  Teuchos::SerialDenseMatrix<int,ScalarType> subQ(Teuchos::View, *(oldState.Q), _blockSize, cur_nevBlocks, curDim-_blockSize);
726  //
727  // Get a view of the new B block of the updated Krylov-Schur factorization
728  Teuchos::SerialDenseMatrix<int,ScalarType> newB(Teuchos::View, *newH, _blockSize, cur_nevBlocks, cur_nevBlocks);
729  //
730  // Compute the new B block.
731  blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, _blockSize, cur_nevBlocks, _blockSize, one,
732  oldB.values(), oldB.stride(), subQ.values(), subQ.stride(), zero, newB.values(), newB.stride() );
733 
734 
735  //
736  // Set the new state and initialize the solver.
738  if (_inSituRestart) {
739  newstate.V = oldState.V;
740  } else {
741  newstate.V = workMV;
742  }
743  newstate.H = newH;
744  newstate.curDim = cur_nevBlocks;
745  bks_solver->initialize(newstate);
746 
747  } // end of restarting
749  //
750  // we returned from iterate(), but none of our status tests Passed.
751  // something is wrong, and it is probably our fault.
752  //
754  else {
755  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::BlockKrylovSchurSolMgr::solve(): Invalid return from bks_solver::iterate().");
756  }
757  }
758  catch (const AnasaziError &err) {
759  printer->stream(Errors)
760  << "Anasazi::BlockKrylovSchurSolMgr::solve() caught unexpected exception from Anasazi::BlockKrylovSchur::iterate() at iteration " << bks_solver->getNumIters() << std::endl
761  << err.what() << std::endl
762  << "Anasazi::BlockKrylovSchurSolMgr::solve() returning Unconverged with no solutions." << std::endl;
763  return Unconverged;
764  }
765  }
766 
767  //
768  // free temporary space
769  workMV = Teuchos::null;
770 
771  // Get the most current Ritz values before we return
772  _ritzValues = bks_solver->getRitzValues();
773 
774  sol.numVecs = ordertest->howMany();
775  printer->stream(Debug) << "ordertest->howMany() : " << sol.numVecs << std::endl;
776  std::vector<int> whichVecs = ordertest->whichVecs();
777 
778  // Place any converged eigenpairs in the solution container.
779  if (sol.numVecs > 0) {
780 
781  // Next determine if there is a conjugate pair on the boundary and resize.
782  std::vector<int> tmpIndex = bks_solver->getRitzIndex();
783  for (int i=0; i<(int)_ritzValues.size(); ++i) {
784  printer->stream(Debug) << _ritzValues[i].realpart << " + i " << _ritzValues[i].imagpart << ", Index = " << tmpIndex[i] << std::endl;
785  }
786  printer->stream(Debug) << "Number of converged eigenpairs (before) = " << sol.numVecs << std::endl;
787  for (int i=0; i<sol.numVecs; ++i) {
788  printer->stream(Debug) << "whichVecs[" << i << "] = " << whichVecs[i] << ", tmpIndex[" << whichVecs[i] << "] = " << tmpIndex[whichVecs[i]] << std::endl;
789  }
790  if (tmpIndex[whichVecs[sol.numVecs-1]]==1) {
791  printer->stream(Debug) << "There is a conjugate pair on the boundary, resizing sol.numVecs" << std::endl;
792  whichVecs.push_back(whichVecs[sol.numVecs-1]+1);
793  sol.numVecs++;
794  for (int i=0; i<sol.numVecs; ++i) {
795  printer->stream(Debug) << "whichVecs[" << i << "] = " << whichVecs[i] << ", tmpIndex[" << whichVecs[i] << "] = " << tmpIndex[whichVecs[i]] << std::endl;
796  }
797  }
798 
799  bool keepMore = false;
800  int numEvecs = sol.numVecs;
801  printer->stream(Debug) << "Number of converged eigenpairs (after) = " << sol.numVecs << std::endl;
802  printer->stream(Debug) << "whichVecs[sol.numVecs-1] > sol.numVecs-1 : " << whichVecs[sol.numVecs-1] << " > " << sol.numVecs-1 << std::endl;
803  if (whichVecs[sol.numVecs-1] > (sol.numVecs-1)) {
804  keepMore = true;
805  numEvecs = whichVecs[sol.numVecs-1]+1; // Add 1 to fix zero-based indexing
806  printer->stream(Debug) << "keepMore = true; numEvecs = " << numEvecs << std::endl;
807  }
808 
809  // Next set the number of Ritz vectors that the iteration must compute and compute them.
810  bks_solver->setNumRitzVectors(numEvecs);
811  bks_solver->computeRitzVectors();
812 
813  // If the leading Ritz pairs are the converged ones, get the information
814  // from the iteration to the solution container. Otherwise copy the necessary
815  // information using 'whichVecs'.
816  if (!keepMore) {
817  sol.index = bks_solver->getRitzIndex();
818  sol.Evals = bks_solver->getRitzValues();
819  sol.Evecs = MVT::CloneCopy( *(bks_solver->getRitzVectors()) );
820  }
821 
822  // Resize based on the number of solutions being returned and set the number of Ritz
823  // vectors for the iteration to compute.
824  sol.Evals.resize(sol.numVecs);
825  sol.index.resize(sol.numVecs);
826 
827  // If the converged Ritz pairs are not the leading ones, copy over the information directly.
828  if (keepMore) {
829  std::vector<Anasazi::Value<ScalarType> > tmpEvals = bks_solver->getRitzValues();
830  for (int vec_i=0; vec_i<sol.numVecs; ++vec_i) {
831  sol.index[vec_i] = tmpIndex[whichVecs[vec_i]];
832  sol.Evals[vec_i] = tmpEvals[whichVecs[vec_i]];
833  }
834  sol.Evecs = MVT::CloneCopy( *(bks_solver->getRitzVectors()), whichVecs );
835  }
836 
837  // Set the solution space to be the Ritz vectors at this time.
838  sol.Espace = sol.Evecs;
839  }
840  }
841 
842  // print final summary
843  bks_solver->currentStatus(printer->stream(FinalSummary));
844 
845  // print timing information
846 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
847  if ( printer->isVerbosity( TimingDetails ) ) {
848  Teuchos::TimeMonitor::summarize( printer->stream( TimingDetails ) );
849  }
850 #endif
851 
852  _problem->setSolution(sol);
853  printer->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
854 
855  // get the number of iterations performed during this solve.
856  _numIters = bks_solver->getNumIters();
857 
858  if (sol.numVecs < nev) {
859  return Unconverged; // return from BlockKrylovSchurSolMgr::solve()
860  }
861  return Converged; // return from BlockKrylovSchurSolMgr::solve()
862 }
863 
864 
865 template <class ScalarType, class MV, class OP>
866 void
869 {
870  globalTest_ = global;
871 }
872 
873 template <class ScalarType, class MV, class OP>
876 {
877  return globalTest_;
878 }
879 
880 template <class ScalarType, class MV, class OP>
881 void
884 {
885  debugTest_ = debug;
886 }
887 
888 template <class ScalarType, class MV, class OP>
891 {
892  return debugTest_;
893 }
894 
895 } // end Anasazi namespace
896 
897 #endif /* ANASAZI_BLOCK_KRYLOV_SCHUR_SOLMGR_HPP */
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
Pure virtual base class which describes the basic interface for a solver manager. ...
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.
ScalarType * values() const
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.
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 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
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.
An exception class parent to all Anasazi exceptions.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getDebugStatusTest() const
Get the status test for debugging.
The Anasazi::BlockKrylovSchurSolMgr provides a flexible solver manager over the BlockKrylovSchur eige...
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
Anasazi&#39;s templated, static class providing utilities for the solvers.
void GEQRF(const OrdinalType m, const OrdinalType n, ScalarType *A, const OrdinalType lda, ScalarType *TAU, ScalarType *WORK, const OrdinalType lwork, OrdinalType *info) const
int numVecs
The number of computed eigenpairs.
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.
Basic output manager for sending information of select verbosity levels to the appropriate output str...
Teuchos::RCP< MV > Espace
An orthonormal basis for the computed eigenspace.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
OrdinalType numRows() const
Anasazi&#39;s basic output manager for sending information of select verbosity levels to the appropriate ...
Abstract base class which defines the interface required by an eigensolver and status test class to c...
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()
ReturnType
Enumerated type used to pass back information from a solver manager.
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...
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
Orthogonalization manager based on the SVQB technique described in "A Block Orthogonalization Procedu...
Struct for storing an eigenproblem solution.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
std::vector< Value< ScalarType > > getRitzValues() const
Return the Ritz values from the most recent solve.
void push_back(const value_type &x)
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.
int getNumIters() const
Get the iteration count for the most recent call to solve().
bool isParameter(const std::string &name) const
OrdinalType stride() const
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...
OrdinalType numCols() const
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.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > S
The current Schur form reduction of the valid part of H.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getGlobalStatusTest() const
Get the status test defining global convergence.
Class which provides internal utilities for the Anasazi solvers.