Anasazi Version of the Day
Loading...
Searching...
No Matches
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
55
57#include "AnasaziBasicSort.hpp"
67
68#include "Teuchos_BLAS.hpp"
69#include "Teuchos_LAPACK.hpp"
70#include "Teuchos_TimeMonitor.hpp"
71#include "Teuchos_FancyOStream.hpp"
72
109
127namespace Anasazi {
128
129
156template<class ScalarType, class MV, class OP>
157class BlockKrylovSchurSolMgr : public SolverManager<ScalarType,MV,OP> {
158
159 private:
162 typedef Teuchos::ScalarTraits<ScalarType> SCT;
163 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
164 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
165
166 public:
167
169
170
191 BlockKrylovSchurSolMgr( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
192 Teuchos::ParameterList &pl );
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
224 Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
225 return Teuchos::tuple(timerSolve_, timerRestarting_);
226 }
227
229
231
232
252
254 void setGlobalStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global);
255
257 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getGlobalStatusTest() const;
258
260 void setDebugStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &debug);
261
263 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getDebugStatusTest() const;
264
266
267 private:
268 Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
269 Teuchos::RCP<SortManager<MagnitudeType> > sort_;
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
288 Teuchos::RCP<Teuchos::FancyOStream> osp_;
289
290 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > globalTest_;
291 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > debugTest_;
292
293};
294
295
296// Constructor
297template<class ScalarType, class MV, class OP>
299 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
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()
434template<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
442 Teuchos::BLAS<int,ScalarType> blas;
443 Teuchos::LAPACK<int,ScalarType> lapack;
444
446 // Output manager
447 Teuchos::RCP<OutputManager<ScalarType> > printer = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_,osp_) );
448
450 // Status tests
451 //
452 // convergence
453 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convtest;
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 }
460 Teuchos::RCP<StatusTestWithOrdering<ScalarType,MV,OP> > ordertest
461 = Teuchos::rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sort_,nev) );
462 // for a non-short-circuited OR test, the order doesn't matter
463 Teuchos::Array<Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
464 alltests.push_back(ordertest);
465
466 if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
467
468 Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > combotest
470 // printing StatusTest
471 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
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
481 Teuchos::RCP<OrthoManager<ScalarType,MV> > ortho;
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
498 Teuchos::ParameterList plist;
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
511 Teuchos::RCP<BlockKrylovSchur<ScalarType,MV,OP> > bks_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);
682 Teuchos::SerialDenseMatrix<int,ScalarType> copyQnev(Teuchos::Copy, Qnev);
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.
754 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > newH =
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
901template <class ScalarType, class MV, class OP>
902void
904 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global)
905{
906 globalTest_ = global;
907}
908
909template <class ScalarType, class MV, class OP>
910const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
912{
913 return globalTest_;
914}
915
916template <class ScalarType, class MV, class OP>
917void
919 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &debug)
920{
921 debugTest_ = debug;
922}
923
924template <class ScalarType, class MV, class OP>
925const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
927{
928 return debugTest_;
929}
930
931} // end Anasazi namespace
932
933#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.
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
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...