Anasazi Version of the Day
Loading...
Searching...
No Matches
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
23
25#include "AnasaziBasicSort.hpp"
35
36#include "Teuchos_BLAS.hpp"
37#include "Teuchos_LAPACK.hpp"
38#include "Teuchos_TimeMonitor.hpp"
39#include "Teuchos_FancyOStream.hpp"
40
77
95namespace Anasazi {
96
97
124template<class ScalarType, class MV, class OP>
125class BlockKrylovSchurSolMgr : public SolverManager<ScalarType,MV,OP> {
126
127 private:
130 typedef Teuchos::ScalarTraits<ScalarType> SCT;
131 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
132 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
133
134 public:
135
137
138
159 BlockKrylovSchurSolMgr( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
160 Teuchos::ParameterList &pl );
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
192 Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
193 return Teuchos::tuple(timerSolve_, timerRestarting_);
194 }
195
197
199
200
220
222 void setGlobalStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global);
223
225 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getGlobalStatusTest() const;
226
228 void setDebugStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &debug);
229
231 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getDebugStatusTest() const;
232
234
235 private:
236 Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
237 Teuchos::RCP<SortManager<MagnitudeType> > sort_;
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
256 Teuchos::RCP<Teuchos::FancyOStream> osp_;
257
258 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > globalTest_;
259 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > debugTest_;
260
261};
262
263
264// Constructor
265template<class ScalarType, class MV, class OP>
267 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
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()
402template<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
410 Teuchos::BLAS<int,ScalarType> blas;
411 Teuchos::LAPACK<int,ScalarType> lapack;
412
414 // Output manager
415 Teuchos::RCP<OutputManager<ScalarType> > printer = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_,osp_) );
416
418 // Status tests
419 //
420 // convergence
421 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convtest;
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 }
428 Teuchos::RCP<StatusTestWithOrdering<ScalarType,MV,OP> > ordertest
429 = Teuchos::rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sort_,nev) );
430 // for a non-short-circuited OR test, the order doesn't matter
431 Teuchos::Array<Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
432 alltests.push_back(ordertest);
433
434 if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
435
436 Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > combotest
438 // printing StatusTest
439 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
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
449 Teuchos::RCP<OrthoManager<ScalarType,MV> > ortho;
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
466 Teuchos::ParameterList plist;
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
479 Teuchos::RCP<BlockKrylovSchur<ScalarType,MV,OP> > bks_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);
650 Teuchos::SerialDenseMatrix<int,ScalarType> copyQnev(Teuchos::Copy, Qnev);
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.
722 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > newH =
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
869template <class ScalarType, class MV, class OP>
870void
872 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global)
873{
874 globalTest_ = global;
875}
876
877template <class ScalarType, class MV, class OP>
878const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
883
884template <class ScalarType, class MV, class OP>
885void
887 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &debug)
888{
889 debugTest_ = debug;
890}
891
892template <class ScalarType, class MV, class OP>
893const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
898
899} // end Anasazi namespace
900
901#endif /* ANASAZI_BLOCK_KRYLOV_SCHUR_SOLMGR_HPP */
Basic implementation of the Anasazi::OrthoManager class.
Basic implementation of the Anasazi::SortManager class.
Implementation of a block Krylov-Schur eigensolver.
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
Abstract base class which defines the interface required by an eigensolver and status test class to c...
Basic implementation of the Anasazi::OrthoManager class.
Abstract class definition for Anasazi Output Managers.
Abstract class definition for Anasazi output stream.
Orthogonalization manager based on the SVQB technique described in "A Block Orthogonalization Procedu...
Pure virtual base class which describes the basic interface for a solver manager.
Class which provides internal utilities for the Anasazi solvers.
Status test for forming logical combinations of other status tests.
Special StatusTest for printing status tests.
A status test for testing the norm of the eigenvectors residuals.
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
Types and exceptions used within Anasazi solvers and interfaces.
An exception class parent to all Anasazi exceptions.
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using (potentially)...
An implementation of the Anasazi::SortManager that performs a collection of common sorting techniques...
The Anasazi::BlockKrylovSchurSolMgr provides a flexible solver manager over the BlockKrylovSchur eige...
int getNumIters() const
Get the iteration count for the most recent call to solve().
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getGlobalStatusTest() const
Get the status test defining global convergence.
void setDebugStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &debug)
Set the status test for debugging.
BlockKrylovSchurSolMgr(const Teuchos::RCP< Eigenproblem< ScalarType, MV, OP > > &problem, Teuchos::ParameterList &pl)
Basic constructor for BlockKrylovSchurSolMgr.
std::vector< Value< ScalarType > > getRitzValues() const
Return the Ritz values from the most recent solve.
ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver's iterate() routine until ...
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getDebugStatusTest() const
Get the status test for debugging.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
void setGlobalStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &global)
Set the status test defining global convergence.
This class implements the block Krylov-Schur iteration, for solving linear eigenvalue problems.
This class defines the interface required by an eigensolver and status test class to compute solution...
An implementation of the Anasazi::GenOrthoManager that performs orthogonalization using iterated clas...
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
Output managers remove the need for the eigensolver to know any information about the required output...
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using the SVQB iter...
The Anasazi::SolverManager is a templated virtual base class that defines the basic interface that an...
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.
Status test for forming logical combinations of other status tests.
A special StatusTest for printing other status tests.
A status test for testing the norm of the eigenvectors residuals.
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
Common interface of stopping criteria for Anasazi's solvers.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
ReturnType
Enumerated type used to pass back information from a solver manager.
Structure to contain pointers to BlockKrylovSchur state variables.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > S
The current Schur form reduction of the valid part of H.
int curDim
The current dimension of the reduction.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > Q
The current Schur vectors of the valid part of H.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
Teuchos::RCP< const MulVec > V
The current Krylov basis.
Struct for storing an eigenproblem solution.
Teuchos::RCP< MV > Evecs
The computed eigenvectors.
std::vector< int > index
An index into Evecs to allow compressed storage of eigenvectors for real, non-Hermitian problems.
int numVecs
The number of computed eigenpairs.
Teuchos::RCP< MV > Espace
An orthonormal basis for the computed eigenspace.
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.
Output managers remove the need for the eigensolver to know any information about the required output...