Anasazi Version of the Day
Loading...
Searching...
No Matches
AnasaziLOBPCGSolMgr.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_LOBPCG_SOLMGR_HPP
43#define ANASAZI_LOBPCG_SOLMGR_HPP
44
49#include "AnasaziConfigDefs.hpp"
50#include "AnasaziTypes.hpp"
51
55
56#include "AnasaziLOBPCG.hpp"
57#include "AnasaziBasicSort.hpp"
68
69#include "Teuchos_BLAS.hpp"
70#include "Teuchos_TimeMonitor.hpp"
71#include "Teuchos_FancyOStream.hpp"
72
123
124namespace Anasazi {
125
173template<class ScalarType, class MV, class OP>
174class LOBPCGSolMgr : public SolverManager<ScalarType,MV,OP> {
175
176 private:
179 typedef Teuchos::ScalarTraits<ScalarType> SCT;
180 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
181 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
182
183 public:
184
186
187
215 LOBPCGSolMgr( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
216 Teuchos::ParameterList &pl );
217
219 virtual ~LOBPCGSolMgr() {};
221
223
224
227 return *problem_;
228 }
229
231 int getNumIters() const {
232 return numIters_;
233 }
234
241 Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
242 return Teuchos::tuple(_timerSolve, _timerLocking);
243 }
244
245
247
249
250
278
280 void setGlobalStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global);
281
283 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getGlobalStatusTest() const;
284
286 void setLockingStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &locking);
287
289 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getLockingStatusTest() const;
290
292 void setDebugStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &debug);
293
295 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getDebugStatusTest() const;
296
298
299 private:
300 Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
301
302 std::string whch_, ortho_;
303
304 MagnitudeType convtol_, locktol_;
305 int maxIters_, numIters_;
306 bool useLocking_;
307 bool relconvtol_, rellocktol_;
308 int blockSize_;
309 bool fullOrtho_;
310 int maxLocked_;
311 int verbosity_;
312 int lockQuorum_;
313 bool recover_;
314 Teuchos::RCP<LOBPCGState<ScalarType,MV> > state_;
315 enum ResType convNorm_, lockNorm_;
316 int osProc_;
317
318 Teuchos::RCP<Teuchos::Time> _timerSolve, _timerLocking;
319
320 Teuchos::RCP<Teuchos::FancyOStream> osp_;
321
322 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > globalTest_;
323 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > lockingTest_;
324 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > debugTest_;
325};
326
327
328// Constructor
329template<class ScalarType, class MV, class OP>
331 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
332 Teuchos::ParameterList &pl ) :
333 problem_(problem),
334 whch_("SR"),
335 ortho_("SVQB"),
336 convtol_(MT::prec()),
337 maxIters_(100),
338 numIters_(0),
339 useLocking_(false),
340 relconvtol_(true),
341 rellocktol_(true),
342 blockSize_(0),
343 fullOrtho_(true),
344 maxLocked_(0),
345 verbosity_(Anasazi::Errors),
346 lockQuorum_(1),
347 recover_(true),
348 osProc_(0)
349#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
350 , _timerSolve(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCGSolMgr::solve()")),
351 _timerLocking(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCGSolMgr locking"))
352#endif
353{
354 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
355 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument, "Problem not set.");
356 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isHermitian(), std::invalid_argument, "Problem not symmetric.");
357 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument, "Problem does not contain initial vectors to clone from.");
358
359
360 std::string strtmp;
361
362 // which values to solve for
363 whch_ = pl.get("Which",whch_);
364 TEUCHOS_TEST_FOR_EXCEPTION(whch_ != "SM" && whch_ != "LM" && whch_ != "SR" && whch_ != "LR",
365 std::invalid_argument, "Anasazi::LOBPCGSolMgr: Invalid sorting string.");
366
367 // which orthogonalization to use
368 ortho_ = pl.get("Orthogonalization",ortho_);
369 if (ortho_ != "DGKS" && ortho_ != "SVQB" && ortho_ != "ICGS") {
370 ortho_ = "SVQB";
371 }
372
373 // convergence tolerance
374 convtol_ = pl.get("Convergence Tolerance",convtol_);
375 relconvtol_ = pl.get("Relative Convergence Tolerance",relconvtol_);
376 strtmp = pl.get("Convergence Norm",std::string("2"));
377 if (strtmp == "2") {
378 convNorm_ = RES_2NORM;
379 }
380 else if (strtmp == "M") {
381 convNorm_ = RES_ORTH;
382 }
383 else {
384 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
385 "Anasazi::LOBPCGSolMgr: Invalid Convergence Norm.");
386 }
387
388
389 // locking tolerance
390 useLocking_ = pl.get("Use Locking",useLocking_);
391 rellocktol_ = pl.get("Relative Locking Tolerance",rellocktol_);
392 // default: should be less than convtol_
393 locktol_ = convtol_/10;
394 locktol_ = pl.get("Locking Tolerance",locktol_);
395 strtmp = pl.get("Locking Norm",std::string("2"));
396 if (strtmp == "2") {
397 lockNorm_ = RES_2NORM;
398 }
399 else if (strtmp == "M") {
400 lockNorm_ = RES_ORTH;
401 }
402 else {
403 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
404 "Anasazi::LOBPCGSolMgr: Invalid Locking Norm.");
405 }
406
407 // maximum number of iterations
408 maxIters_ = pl.get("Maximum Iterations",maxIters_);
409
410 // block size: default is nev()
411 blockSize_ = pl.get("Block Size",problem_->getNEV());
412 TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
413 "Anasazi::LOBPCGSolMgr: \"Block Size\" must be strictly positive.");
414
415 // max locked: default is nev(), must satisfy maxLocked_ + blockSize_ >= nev
416 if (useLocking_) {
417 maxLocked_ = pl.get("Max Locked",problem_->getNEV());
418 }
419 else {
420 maxLocked_ = 0;
421 }
422 if (maxLocked_ == 0) {
423 useLocking_ = false;
424 }
425 TEUCHOS_TEST_FOR_EXCEPTION(maxLocked_ < 0, std::invalid_argument,
426 "Anasazi::LOBPCGSolMgr: \"Max Locked\" must be positive.");
427 TEUCHOS_TEST_FOR_EXCEPTION(maxLocked_ + blockSize_ < problem_->getNEV(),
428 std::invalid_argument,
429 "Anasazi::LOBPCGSolMgr: Not enough storage space for requested number of eigenpairs.");
430
431 if (useLocking_) {
432 lockQuorum_ = pl.get("Locking Quorum",lockQuorum_);
433 TEUCHOS_TEST_FOR_EXCEPTION(lockQuorum_ <= 0,
434 std::invalid_argument,
435 "Anasazi::LOBPCGSolMgr: \"Locking Quorum\" must be strictly positive.");
436 }
437
438 // full orthogonalization: default true
439 fullOrtho_ = pl.get("Full Ortho",fullOrtho_);
440
441 // Create a formatted output stream to print to.
442 // See if user requests output processor.
443 osProc_ = pl.get("Output Processor", osProc_);
444
445 // If not passed in by user, it will be chosen based upon operator type.
446 if (pl.isParameter("Output Stream")) {
447 osp_ = Teuchos::getParameter<Teuchos::RCP<Teuchos::FancyOStream> >(pl,"Output Stream");
448 }
449 else {
450 osp_ = OutputStreamTraits<OP>::getOutputStream (*problem_->getOperator(), osProc_);
451 }
452
453 // verbosity level
454 if (pl.isParameter("Verbosity")) {
455 if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
456 verbosity_ = pl.get("Verbosity", verbosity_);
457 } else {
458 verbosity_ = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
459 }
460 }
461
462 // recover from LOBPCGRitzFailure
463 recover_ = pl.get("Recover",recover_);
464
465 // get (optionally) an initial state
466 if (pl.isParameter("Init")) {
467 state_ = Teuchos::getParameter<Teuchos::RCP<Anasazi::LOBPCGState<ScalarType,MV> > >(pl,"Init");
468 }
469}
470
471
472// solve()
473template<class ScalarType, class MV, class OP>
476
477 typedef SolverUtils<ScalarType,MV,OP> msutils;
478
479 const int nev = problem_->getNEV();
480
481
482
484 // Sort manager
485 Teuchos::RCP<BasicSort<MagnitudeType> > sorter = Teuchos::rcp( new BasicSort<MagnitudeType>(whch_) );
486
488 // Output manager
489 Teuchos::RCP<OutputManager<ScalarType> > printer = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_,osp_) );
490
492 // Status tests
493 //
494 // maximum number of iterations: optional test
495 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxtest;
496 if (maxIters_ > 0) {
497 maxtest = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>(maxIters_) );
498 }
499 // convergence
500 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convtest;
501 if (globalTest_ == Teuchos::null) {
502 convtest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(convtol_,nev,convNorm_,relconvtol_) );
503 }
504 else {
505 convtest = globalTest_;
506 }
507 Teuchos::RCP<StatusTestWithOrdering<ScalarType,MV,OP> > ordertest
508 = Teuchos::rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sorter,nev) );
509 // locking
510 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > locktest;
511 if (useLocking_) {
512 if (lockingTest_ == Teuchos::null) {
513 locktest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(locktol_,lockQuorum_,lockNorm_,rellocktol_) );
514 }
515 else {
516 locktest = lockingTest_;
517 }
518 }
519 // for a non-short-circuited OR test, the order doesn't matter
520 Teuchos::Array<Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
521 alltests.push_back(ordertest);
522 if (locktest != Teuchos::null) alltests.push_back(locktest);
523 if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
524 if (maxtest != Teuchos::null) alltests.push_back(maxtest);
525 Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > combotest
527 // printing StatusTest
528 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
529 if ( printer->isVerbosity(Debug) ) {
530 outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,Passed+Failed+Undefined ) );
531 }
532 else {
533 outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer,combotest,1,Passed ) );
534 }
535
537 // Orthomanager
538 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho;
539 if (ortho_=="SVQB") {
540 ortho = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
541 } else if (ortho_=="DGKS") {
542 ortho = Teuchos::rcp( new BasicOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
543 } else if (ortho_=="ICGS") {
544 ortho = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
545 } else {
546 TEUCHOS_TEST_FOR_EXCEPTION(ortho_!="SVQB"&&ortho_!="DGKS"&&ortho_!="ICGS",std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): Invalid orthogonalization type.");
547 }
548
550 // Parameter list
551 Teuchos::ParameterList plist;
552 plist.set("Block Size",blockSize_);
553 plist.set("Full Ortho",fullOrtho_);
554
556 // LOBPCG solver
557 Teuchos::RCP<LOBPCG<ScalarType,MV,OP> > lobpcg_solver
558 = Teuchos::rcp( new LOBPCG<ScalarType,MV,OP>(problem_,sorter,printer,outputtest,ortho,plist) );
559 // set any auxiliary vectors defined in the problem
560 Teuchos::RCP< const MV > probauxvecs = problem_->getAuxVecs();
561 if (probauxvecs != Teuchos::null) {
562 lobpcg_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
563 }
564
566 // Storage
567 //
568 // lockvecs will contain eigenvectors that have been determined "locked" by the status test
569 int curNumLocked = 0;
570 Teuchos::RCP<MV> lockvecs;
571 if (useLocking_) {
572 lockvecs = MVT::Clone(*problem_->getInitVec(),maxLocked_);
573 }
574 std::vector<MagnitudeType> lockvals;
575 // workMV will be used as work space for LOBPCGRitzFailure recovery and locking
576 // it will be partitioned in these cases as follows:
577 // for LOBPCGRitzFailure recovery:
578 // workMV = [X H P OpX OpH OpP], where OpX OpH OpP will be used for K and M
579 // total size: 2*3*blocksize
580 // for locking
581 // workMV = [X P MX MP], with MX,MP needing storage only if hasM==true
582 // total size: 2*blocksize or 4*blocksize
583 Teuchos::RCP<MV> workMV;
584 if (fullOrtho_ == false && recover_ == true) {
585 workMV = MVT::Clone(*problem_->getInitVec(),2*3*blockSize_);
586 }
587 else if (useLocking_) {
588 if (problem_->getM() != Teuchos::null) {
589 workMV = MVT::Clone(*problem_->getInitVec(),4*blockSize_);
590 }
591 else {
592 workMV = MVT::Clone(*problem_->getInitVec(),2*blockSize_);
593 }
594 }
595
596 // initialize the solution to nothing in case we throw an exception
598 sol.numVecs = 0;
599 problem_->setSolution(sol);
600
601 // initialize the solver if the user specified a state
602 if (state_ != Teuchos::null) {
603 lobpcg_solver->initialize(*state_);
604 }
605
606 // enter solve() iterations
607 {
608#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
609 Teuchos::TimeMonitor slvtimer(*_timerSolve);
610#endif
611
612 // tell the lobpcg_solver to iterate
613 while (1) {
614 try {
615 lobpcg_solver->iterate();
616
618 //
619 // check user-specified debug test; if it passed, return an exception
620 //
622 if (debugTest_ != Teuchos::null && debugTest_->getStatus() == Passed) {
623 throw AnasaziError("Anasazi::LOBPCGSolMgr::solve(): User-specified debug status test returned Passed.");
624 }
626 //
627 // check convergence first
628 //
630 else if (ordertest->getStatus() == Passed || (maxtest != Teuchos::null && maxtest->getStatus() == Passed) ) {
631 // we have convergence or not
632 // ordertest->whichVecs() tells us which vectors from lockvecs and solver state are the ones we want
633 // ordertest->howMany() will tell us how many
634 break;
635 }
637 //
638 // check locking if we didn't converge
639 //
641 else if (locktest != Teuchos::null && locktest->getStatus() == Passed) {
642
643#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
644 Teuchos::TimeMonitor lcktimer(*_timerLocking);
645#endif
646
647 // remove the locked vectors,values from lobpcg_solver: put them in newvecs, newvals
648 TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() <= 0,std::logic_error,
649 "Anasazi::LOBPCGSolMgr::solve(): status test mistake: howMany() non-positive.");
650 TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() != (int)locktest->whichVecs().size(),std::logic_error,
651 "Anasazi::LOBPCGSolMgr::solve(): status test mistake: howMany() not consistent with whichVecs().");
652 TEUCHOS_TEST_FOR_EXCEPTION(curNumLocked == maxLocked_,std::logic_error,
653 "Anasazi::LOBPCGSolMgr::solve(): status test mistake: locking not deactivated.");
654 // get the indices
655 int numnew = locktest->howMany();
656 std::vector<int> indnew = locktest->whichVecs();
657
658 // don't lock more than maxLocked_; we didn't allocate enough space.
659 if (curNumLocked + numnew > maxLocked_) {
660 numnew = maxLocked_ - curNumLocked;
661 indnew.resize(numnew);
662 }
663
664 // the call below to lobpcg_solver->setAuxVecs() will reset the solver to unitialized with hasP() == false
665 // store the hasP() state for use below
666 bool hadP = lobpcg_solver->hasP();
667
668 {
669 // debug printing
670 printer->print(Debug,"Locking vectors: ");
671 for (unsigned int i=0; i<indnew.size(); i++) {printer->stream(Debug) << " " << indnew[i];}
672 printer->print(Debug,"\n");
673 }
674 std::vector<MagnitudeType> newvals(numnew);
675 Teuchos::RCP<const MV> newvecs;
676 {
677 // work in a local scope, to hide the variabes needed for extracting this info
678 // get the vectors
679 newvecs = MVT::CloneView(*lobpcg_solver->getRitzVectors(),indnew);
680 // get the values
681 std::vector<Value<ScalarType> > allvals = lobpcg_solver->getRitzValues();
682 for (int i=0; i<numnew; i++) {
683 newvals[i] = allvals[indnew[i]].realpart;
684 }
685 }
686 // put newvecs into lockvecs
687 {
688 std::vector<int> indlock(numnew);
689 for (int i=0; i<numnew; i++) indlock[i] = curNumLocked+i;
690 MVT::SetBlock(*newvecs,indlock,*lockvecs);
691 newvecs = Teuchos::null;
692 }
693 // put newvals into lockvals
694 lockvals.insert(lockvals.end(),newvals.begin(),newvals.end());
695 curNumLocked += numnew;
696 // add locked vecs as aux vecs, along with aux vecs from problem
697 {
698 std::vector<int> indlock(curNumLocked);
699 for (int i=0; i<curNumLocked; i++) indlock[i] = i;
700 Teuchos::RCP<const MV> curlocked = MVT::CloneView(*lockvecs,indlock);
701 if (probauxvecs != Teuchos::null) {
702 lobpcg_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs,curlocked) );
703 }
704 else {
705 lobpcg_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(curlocked) );
706 }
707 }
708 // add locked vals to ordertest
709 ordertest->setAuxVals(lockvals);
710 // fill out the empty state in the solver
711 {
712 LOBPCGState<ScalarType,MV> state = lobpcg_solver->getState();
713 Teuchos::RCP<MV> newstateX, newstateMX, newstateP, newstateMP;
714 //
715 // workMV will be partitioned as follows: workMV = [X P MX MP],
716 //
717 // make a copy of the current X,MX state
718 std::vector<int> bsind(blockSize_);
719 for (int i=0; i<blockSize_; i++) bsind[i] = i;
720 newstateX = MVT::CloneViewNonConst(*workMV,bsind);
721 MVT::SetBlock(*state.X,bsind,*newstateX);
722
723 if (state.MX != Teuchos::null) {
724 std::vector<int> block3(blockSize_);
725 for (int i=0; i<blockSize_; i++) block3[i] = 2*blockSize_+i;
726 newstateMX = MVT::CloneViewNonConst(*workMV,block3);
727 MVT::SetBlock(*state.MX,bsind,*newstateMX);
728 }
729 //
730 // get select part, set to random, apply M
731 {
732 Teuchos::RCP<MV> newX = MVT::CloneViewNonConst(*newstateX,indnew);
733 MVT::MvRandom(*newX);
734
735 if (newstateMX != Teuchos::null) {
736 Teuchos::RCP<MV> newMX = MVT::CloneViewNonConst(*newstateMX,indnew);
737 OPT::Apply(*problem_->getM(),*newX,*newMX);
738 }
739 }
740
741 Teuchos::Array<Teuchos::RCP<const MV> > curauxvecs = lobpcg_solver->getAuxVecs();
742 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
743 // ortho X against the aux vectors
744 ortho->projectAndNormalizeMat(*newstateX,curauxvecs,dummyC,Teuchos::null,newstateMX);
745
746 if (hadP) {
747 //
748 // get P and optionally MP, orthogonalize against X and auxiliary vectors
749 std::vector<int> block2(blockSize_);
750 for (int i=0; i<blockSize_; i++) block2[i] = blockSize_+i;
751 newstateP = MVT::CloneViewNonConst(*workMV,block2);
752 MVT::SetBlock(*state.P,bsind,*newstateP);
753
754 if (state.MP != Teuchos::null) {
755 std::vector<int> block4(blockSize_);
756 for (int i=0; i<blockSize_; i++) block4[i] = 3*blockSize_+i;
757 newstateMP = MVT::CloneViewNonConst(*workMV,block4);
758 MVT::SetBlock(*state.MP,bsind,*newstateMP);
759 }
760
761 if (fullOrtho_) {
762 // ortho P against the new aux vectors and new X
763 curauxvecs.push_back(newstateX);
764 ortho->projectAndNormalizeMat(*newstateP,curauxvecs,dummyC,Teuchos::null,newstateMP);
765 }
766 else {
767 // ortho P against the new aux vectors
768 ortho->projectAndNormalizeMat(*newstateP,curauxvecs,dummyC,Teuchos::null,newstateMP);
769 }
770 }
771 // set the new state
773 newstate.X = newstateX;
774 newstate.MX = newstateMX;
775 newstate.P = newstateP;
776 newstate.MP = newstateMP;
777 lobpcg_solver->initialize(newstate);
778 }
779
780 if (curNumLocked == maxLocked_) {
781 // disable locking now; remove locking test from combo test
782 combotest->removeTest(locktest);
783 }
784 }
785 else {
786 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): Invalid return from lobpcg_solver::iterate().");
787 }
788 }
790 //
791 // check Ritz Failure
792 //
794 catch (const LOBPCGRitzFailure &re) {
795 if (fullOrtho_==true || recover_==false) {
796 // if we are already using full orthogonalization, there isn't much we can do here.
797 // the most recent information in the status tests is still valid, and can be used to extract/return the
798 // eigenpairs that have converged.
799 printer->stream(Warnings) << "Error! Caught LOBPCGRitzFailure at iteration " << lobpcg_solver->getNumIters() << std::endl
800 << "Will not try to recover." << std::endl;
801 break; // while(1)
802 }
803 printer->stream(Warnings) << "Error! Caught LOBPCGRitzFailure at iteration " << lobpcg_solver->getNumIters() << std::endl
804 << "Full orthogonalization is off; will try to recover." << std::endl;
805 // get the current "basis" from the solver, orthonormalize it, do a rayleigh-ritz, and restart with the ritz vectors
806 // if there aren't enough, break and quit with what we have
807 //
808 // workMV = [X H P OpX OpH OpP], where OpX OpH OpP will be used for K and M
809 LOBPCGState<ScalarType,MV> curstate = lobpcg_solver->getState();
810 Teuchos::RCP<MV> restart, Krestart, Mrestart;
811 int localsize = lobpcg_solver->hasP() ? 3*blockSize_ : 2*blockSize_;
812 bool hasM = problem_->getM() != Teuchos::null;
813 {
814 std::vector<int> recind(localsize);
815 for (int i=0; i<localsize; i++) recind[i] = i;
816 restart = MVT::CloneViewNonConst(*workMV,recind);
817 }
818 {
819 std::vector<int> recind(localsize);
820 for (int i=0; i<localsize; i++) recind[i] = localsize+i;
821 Krestart = MVT::CloneViewNonConst(*workMV,recind);
822 }
823 if (hasM) {
824 Mrestart = Krestart;
825 }
826 else {
827 Mrestart = restart;
828 }
829 //
830 // set restart = [X H P] and Mrestart = M*[X H P]
831 //
832 // put X into [0 , blockSize)
833 {
834 std::vector<int> blk1(blockSize_);
835 for (int i=0; i < blockSize_; i++) blk1[i] = i;
836 MVT::SetBlock(*curstate.X,blk1,*restart);
837
838 // put MX into [0 , blockSize)
839 if (hasM) {
840 MVT::SetBlock(*curstate.MX,blk1,*Mrestart);
841 }
842 }
843 //
844 // put H into [blockSize_ , 2*blockSize)
845 {
846 std::vector<int> blk2(blockSize_);
847 for (int i=0; i < blockSize_; i++) blk2[i] = blockSize_+i;
848 MVT::SetBlock(*curstate.H,blk2,*restart);
849
850 // put MH into [blockSize_ , 2*blockSize)
851 if (hasM) {
852 MVT::SetBlock(*curstate.MH,blk2,*Mrestart);
853 }
854 }
855 // optionally, put P into [2*blockSize,3*blockSize)
856 if (localsize == 3*blockSize_) {
857 std::vector<int> blk3(blockSize_);
858 for (int i=0; i < blockSize_; i++) blk3[i] = 2*blockSize_+i;
859 MVT::SetBlock(*curstate.P,blk3,*restart);
860
861 // put MP into [2*blockSize,3*blockSize)
862 if (hasM) {
863 MVT::SetBlock(*curstate.MP,blk3,*Mrestart);
864 }
865 }
866 // project against auxvecs and locked vecs, and orthonormalize the basis
867 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
868 Teuchos::Array<Teuchos::RCP<const MV> > Q;
869 {
870 if (curNumLocked > 0) {
871 std::vector<int> indlock(curNumLocked);
872 for (int i=0; i<curNumLocked; i++) indlock[i] = i;
873 Teuchos::RCP<const MV> curlocked = MVT::CloneView(*lockvecs,indlock);
874 Q.push_back(curlocked);
875 }
876 if (probauxvecs != Teuchos::null) {
877 Q.push_back(probauxvecs);
878 }
879 }
880 int rank = ortho->projectAndNormalizeMat(*restart,Q,dummyC,Teuchos::null,Mrestart);
881 if (rank < blockSize_) {
882 // quit
883 printer->stream(Errors) << "Error! Recovered basis only rank " << rank << ". Block size is " << blockSize_ << ".\n"
884 << "Recovery failed." << std::endl;
885 break;
886 }
887 // reduce multivec size if necessary
888 if (rank < localsize) {
889 localsize = rank;
890 std::vector<int> redind(localsize);
891 for (int i=0; i<localsize; i++) redind[i] = i;
892 // grab the first part of restart,Krestart
893 restart = MVT::CloneViewNonConst(*restart,redind);
894 Krestart = MVT::CloneViewNonConst(*Krestart,redind);
895 if (hasM) {
896 Mrestart = Krestart;
897 }
898 else {
899 Mrestart = restart;
900 }
901 }
902 Teuchos::SerialDenseMatrix<int,ScalarType> KK(localsize,localsize), MM(localsize,localsize), S(localsize,localsize);
903 std::vector<MagnitudeType> theta(localsize);
904 // project the matrices
905 //
906 // MM = restart^H M restart
907 MVT::MvTransMv(1.0,*restart,*Mrestart,MM);
908 //
909 // compute Krestart = K*restart
910 OPT::Apply(*problem_->getOperator(),*restart,*Krestart);
911 //
912 // KK = restart^H K restart
913 MVT::MvTransMv(1.0,*restart,*Krestart,KK);
914 rank = localsize;
915 msutils::directSolver(localsize,KK,Teuchos::rcpFromRef(MM),S,theta,rank,1);
916 if (rank < blockSize_) {
917 printer->stream(Errors) << "Error! Recovered basis of rank " << rank << " produced only " << rank << "ritz vectors.\n"
918 << "Block size is " << blockSize_ << ".\n"
919 << "Recovery failed." << std::endl;
920 break;
921 }
922 theta.resize(rank);
923 //
924 // sort the ritz values using the sort manager
925 {
926 Teuchos::BLAS<int,ScalarType> blas;
927 std::vector<int> order(rank);
928 // sort
929 sorter->sort( theta, Teuchos::rcpFromRef(order),rank ); // don't catch exception
930 // Sort the primitive ritz vectors
931 Teuchos::SerialDenseMatrix<int,ScalarType> curS(Teuchos::View,S,rank,rank);
932 msutils::permuteVectors(order,curS);
933 }
934 //
935 Teuchos::SerialDenseMatrix<int,ScalarType> S1(Teuchos::View,S,localsize,blockSize_);
936 //
937 // compute the ritz vectors: store them in Krestart
939 Teuchos::RCP<MV> newX;
940 {
941 std::vector<int> bsind(blockSize_);
942 for (int i=0; i<blockSize_; i++) bsind[i] = i;
943 newX = MVT::CloneViewNonConst(*Krestart,bsind);
944 }
945 MVT::MvTimesMatAddMv(1.0,*restart,S1,0.0,*newX);
946 // send X and theta into the solver
947 newstate.X = newX;
948 theta.resize(blockSize_);
949 newstate.T = Teuchos::rcpFromRef(theta);
950 // initialize
951 lobpcg_solver->initialize(newstate);
952 }
953 catch (const AnasaziError &err) {
954 printer->stream(Errors)
955 << "Anasazi::LOBPCGSolMgr::solve() caught unexpected exception from Anasazi::LOBPCG::iterate() at iteration " << lobpcg_solver->getNumIters() << std::endl
956 << err.what() << std::endl
957 << "Anasazi::LOBPCGSolMgr::solve() returning Unconverged with no solutions." << std::endl;
958 return Unconverged;
959 }
960 // don't catch any other exceptions
961 }
962
963 sol.numVecs = ordertest->howMany();
964 if (sol.numVecs > 0) {
965 sol.Evecs = MVT::Clone(*problem_->getInitVec(),sol.numVecs);
966 sol.Espace = sol.Evecs;
967 sol.Evals.resize(sol.numVecs);
968 std::vector<MagnitudeType> vals(sol.numVecs);
969
970 // copy them into the solution
971 std::vector<int> which = ordertest->whichVecs();
972 // indices between [0,blockSize) refer to vectors/values in the solver
973 // indices between [-curNumLocked,-1] refer to locked vectors/values [0,curNumLocked)
974 // everything has already been ordered by the solver; we just have to partition the two references
975 std::vector<int> inlocked(0), insolver(0);
976 for (unsigned int i=0; i<which.size(); i++) {
977 if (which[i] >= 0) {
978 TEUCHOS_TEST_FOR_EXCEPTION(which[i] >= blockSize_,std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): positive indexing mistake from ordertest.");
979 insolver.push_back(which[i]);
980 }
981 else {
982 // sanity check
983 TEUCHOS_TEST_FOR_EXCEPTION(which[i] < -curNumLocked,std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): negative indexing mistake from ordertest.");
984 inlocked.push_back(which[i] + curNumLocked);
985 }
986 }
987
988 TEUCHOS_TEST_FOR_EXCEPTION(insolver.size() + inlocked.size() != (unsigned int)sol.numVecs,std::logic_error,"Anasazi::LOBPCGSolMgr::solve(): indexing mistake.");
989
990 // set the vecs,vals in the solution
991 if (insolver.size() > 0) {
992 // set vecs
993 int lclnum = insolver.size();
994 std::vector<int> tosol(lclnum);
995 for (int i=0; i<lclnum; i++) tosol[i] = i;
996 Teuchos::RCP<const MV> v = MVT::CloneView(*lobpcg_solver->getRitzVectors(),insolver);
997 MVT::SetBlock(*v,tosol,*sol.Evecs);
998 // set vals
999 std::vector<Value<ScalarType> > fromsolver = lobpcg_solver->getRitzValues();
1000 for (unsigned int i=0; i<insolver.size(); i++) {
1001 vals[i] = fromsolver[insolver[i]].realpart;
1002 }
1003 }
1004
1005 // get the vecs,vals from locked storage
1006 if (inlocked.size() > 0) {
1007 int solnum = insolver.size();
1008 // set vecs
1009 int lclnum = inlocked.size();
1010 std::vector<int> tosol(lclnum);
1011 for (int i=0; i<lclnum; i++) tosol[i] = solnum + i;
1012 Teuchos::RCP<const MV> v = MVT::CloneView(*lockvecs,inlocked);
1013 MVT::SetBlock(*v,tosol,*sol.Evecs);
1014 // set vals
1015 for (unsigned int i=0; i<inlocked.size(); i++) {
1016 vals[i+solnum] = lockvals[inlocked[i]];
1017 }
1018 }
1019
1020 // sort the eigenvalues and permute the eigenvectors appropriately
1021 {
1022 std::vector<int> order(sol.numVecs);
1023 sorter->sort( vals, Teuchos::rcpFromRef(order), sol.numVecs);
1024 // store the values in the Eigensolution
1025 for (int i=0; i<sol.numVecs; i++) {
1026 sol.Evals[i].realpart = vals[i];
1027 sol.Evals[i].imagpart = MT::zero();
1028 }
1029 // now permute the eigenvectors according to order
1030 msutils::permuteVectors(sol.numVecs,order,*sol.Evecs);
1031 }
1032
1033 // setup sol.index, remembering that all eigenvalues are real so that index = {0,...,0}
1034 sol.index.resize(sol.numVecs,0);
1035 }
1036 }
1037
1038 // print final summary
1039 lobpcg_solver->currentStatus(printer->stream(FinalSummary));
1040
1041 // print timing information
1042#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1043 if ( printer->isVerbosity( TimingDetails ) ) {
1044 Teuchos::TimeMonitor::summarize( printer->stream( TimingDetails ) );
1045 }
1046#endif
1047
1048 problem_->setSolution(sol);
1049 printer->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
1050
1051 // get the number of iterations performed in this call to solve.
1052 numIters_ = lobpcg_solver->getNumIters();
1053
1054 if (sol.numVecs < nev) {
1055 return Unconverged; // return from LOBPCGSolMgr::solve()
1056 }
1057 return Converged; // return from LOBPCGSolMgr::solve()
1058}
1059
1060
1061template <class ScalarType, class MV, class OP>
1062void
1064 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global)
1065{
1066 globalTest_ = global;
1067}
1068
1069template <class ScalarType, class MV, class OP>
1070const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
1072{
1073 return globalTest_;
1074}
1075
1076template <class ScalarType, class MV, class OP>
1077void
1079 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &debug)
1080{
1081 debugTest_ = debug;
1082}
1083
1084template <class ScalarType, class MV, class OP>
1085const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
1087{
1088 return debugTest_;
1089}
1090
1091template <class ScalarType, class MV, class OP>
1092void
1094 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &locking)
1095{
1096 lockingTest_ = locking;
1097}
1098
1099template <class ScalarType, class MV, class OP>
1100const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
1102{
1103 return lockingTest_;
1104}
1105
1106} // end Anasazi namespace
1107
1108#endif /* ANASAZI_LOBPCG_SOLMGR_HPP */
Basic implementation of the Anasazi::OrthoManager class.
Basic implementation of the Anasazi::SortManager class.
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.
Implementation of the locally-optimal block preconditioned conjugate gradient (LOBPCG) method.
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.
Status test for testing the number of iterations.
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...
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...
LOBPCGRitzFailure is thrown when the LOBPCG solver is unable to continue a call to LOBPCG::iterate() ...
User interface for the LOBPCG eigensolver.
virtual ~LOBPCGSolMgr()
Destructor.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
void setLockingStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &locking)
Set the status test defining locking.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getGlobalStatusTest() const
Get the status test defining global convergence.
void setGlobalStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &global)
Set the status test defining global convergence.
ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver's iterate() routine until ...
int getNumIters() const
Get the iteration count for the most recent call to solve().
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getDebugStatusTest() const
Get the status test for debugging.
void setDebugStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &debug)
Set the status test for debugging.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getLockingStatusTest() const
Get the status test defining locking.
LOBPCGSolMgr(const Teuchos::RCP< Eigenproblem< ScalarType, MV, OP > > &problem, Teuchos::ParameterList &pl)
Basic constructor for LOBPCGSolMgr.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
This class provides the Locally Optimal Block Preconditioned Conjugate Gradient (LOBPCG) iteration,...
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...
Anasazi's templated, static class providing utilities for the solvers.
Status test for forming logical combinations of other status tests.
A status test for testing the number of iterations.
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.
ResType
Enumerated type used to specify which residual norm used by residual norm status tests.
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.
Structure to contain pointers to Anasazi state variables.
Teuchos::RCP< const MultiVector > X
The current eigenvectors.
Teuchos::RCP< const MultiVector > MX
The image of the current eigenvectors under M, or Teuchos::null if M was not specified.
Teuchos::RCP< const MultiVector > MP
The image of the current search direction under M, or Teuchos::null if M was not specified.
Teuchos::RCP< const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > T
The current Ritz values.
Teuchos::RCP< const MultiVector > H
The current preconditioned residual vectors.
Teuchos::RCP< const MultiVector > P
The current search direction.
Teuchos::RCP< const MultiVector > MH
The image of the current preconditioned residual vectors under M, or Teuchos::null if M was not speci...
Output managers remove the need for the eigensolver to know any information about the required output...