Anasazi Version of the Day
Loading...
Searching...
No Matches
AnasaziRTRBase.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
46#ifndef ANASAZI_RTRBASE_HPP
47#define ANASAZI_RTRBASE_HPP
48
49#include "AnasaziTypes.hpp"
50
54#include "Teuchos_ScalarTraits.hpp"
55
58
59#include "Teuchos_LAPACK.hpp"
60#include "Teuchos_BLAS.hpp"
61#include "Teuchos_SerialDenseMatrix.hpp"
62#include "Teuchos_ParameterList.hpp"
63#include "Teuchos_TimeMonitor.hpp"
64
114namespace Anasazi {
115
117
118
123 template <class ScalarType, class MV>
124 struct RTRState {
126 Teuchos::RCP<const MV> X;
128 Teuchos::RCP<const MV> AX;
130 Teuchos::RCP<const MV> BX;
132 Teuchos::RCP<const MV> R;
134 Teuchos::RCP<const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > T;
138 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType rho;
139 RTRState() : X(Teuchos::null),AX(Teuchos::null),BX(Teuchos::null),
140 R(Teuchos::null),T(Teuchos::null),rho(0) {};
141 };
142
144
146
147
161 class RTRRitzFailure : public AnasaziError {public:
162 RTRRitzFailure(const std::string& what_arg) : AnasaziError(what_arg)
163 {}};
164
173 class RTRInitFailure : public AnasaziError {public:
174 RTRInitFailure(const std::string& what_arg) : AnasaziError(what_arg)
175 {}};
176
193 class RTROrthoFailure : public AnasaziError {public:
194 RTROrthoFailure(const std::string& what_arg) : AnasaziError(what_arg)
195 {}};
196
197
199
200
201 template <class ScalarType, class MV, class OP>
202 class RTRBase : public Eigensolver<ScalarType,MV,OP> {
203 public:
204
206
207
213 RTRBase(const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
214 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
215 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
216 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
217 const Teuchos::RCP<GenOrthoManager<ScalarType,MV,OP> > &ortho,
218 Teuchos::ParameterList &params,
219 const std::string &solverLabel, bool skinnySolver
220 );
221
223 virtual ~RTRBase() {};
224
226
228
229
251 virtual void iterate() = 0;
252
277 void initialize(RTRState<ScalarType,MV>& newstate);
278
282 void initialize();
283
296 bool isInitialized() const;
297
306
308
310
311
313 int getNumIters() const;
314
316 void resetNumIters();
317
325 Teuchos::RCP<const MV> getRitzVectors();
326
332 std::vector<Value<ScalarType> > getRitzValues();
333
341 std::vector<int> getRitzIndex();
342
348 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getResNorms();
349
350
356 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRes2Norms();
357
358
363 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRitzRes2Norms();
364
365
374 int getCurSubspaceDim() const;
375
378 int getMaxSubspaceDim() const;
379
381
383
384
386 void setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test);
387
389 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > getStatusTest() const;
390
393
394
403 void setBlockSize(int blockSize);
404
405
407 int getBlockSize() const;
408
409
430 void setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs);
431
433 Teuchos::Array<Teuchos::RCP<const MV> > getAuxVecs() const;
434
436
438
439
441 virtual void currentStatus(std::ostream &os);
442
444
445 protected:
446 //
447 // Convenience typedefs
448 //
452 typedef Teuchos::ScalarTraits<ScalarType> SCT;
453 typedef typename SCT::magnitudeType MagnitudeType;
454 typedef Teuchos::ScalarTraits<MagnitudeType> MAT;
455 const MagnitudeType ONE;
456 const MagnitudeType ZERO;
457 const MagnitudeType NANVAL;
458 typedef typename std::vector<MagnitudeType>::iterator vecMTiter;
459 typedef typename std::vector<ScalarType>::iterator vecSTiter;
460 //
461 // Internal structs
462 //
463 struct CheckList {
464 bool checkX, checkAX, checkBX;
465 bool checkEta, checkAEta, checkBEta;
466 bool checkR, checkQ, checkBR;
467 bool checkZ, checkPBX;
468 CheckList() : checkX(false),checkAX(false),checkBX(false),
469 checkEta(false),checkAEta(false),checkBEta(false),
470 checkR(false),checkQ(false),checkBR(false),
471 checkZ(false), checkPBX(false) {};
472 };
473 //
474 // Internal methods
475 //
476 std::string accuracyCheck(const CheckList &chk, const std::string &where) const;
477 // solves the model minimization
478 virtual void solveTRSubproblem() = 0;
479 // Riemannian metric
480 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType ginner(const MV &xi) const;
481 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType ginner(const MV &xi, const MV &zeta) const;
482 void ginnersep(const MV &xi, std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &d) const;
483 void ginnersep(const MV &xi, const MV &zeta, std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &d) const;
484 //
485 // Classes input through constructor that define the eigenproblem to be solved.
486 //
487 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
488 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sm_;
489 const Teuchos::RCP<OutputManager<ScalarType> > om_;
490 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > tester_;
491 const Teuchos::RCP<GenOrthoManager<ScalarType,MV,OP> > orthman_;
492 //
493 // Information obtained from the eigenproblem
494 //
495 Teuchos::RCP<const OP> AOp_;
496 Teuchos::RCP<const OP> BOp_;
497 Teuchos::RCP<const OP> Prec_;
498 bool hasBOp_, hasPrec_, olsenPrec_;
499 //
500 // Internal timers
501 //
502 Teuchos::RCP<Teuchos::Time> timerAOp_, timerBOp_, timerPrec_,
503 timerSort_,
504 timerLocalProj_, timerDS_,
505 timerLocalUpdate_, timerCompRes_,
506 timerOrtho_, timerInit_;
507 //
508 // Counters
509 //
510 // Number of operator applications
511 int counterAOp_, counterBOp_, counterPrec_;
512
513 //
514 // Algorithmic parameters.
515 //
516 // blockSize_ is the solver block size
517 int blockSize_;
518 //
519 // Current solver state
520 //
521 // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
522 // is capable of running; _initialize is controlled by the initialize() member method
523 // For the implications of the state of initialized_, please see documentation for initialize()
524 bool initialized_;
525 //
526 // nevLocal_ reflects how much of the current basis is valid (0 <= nevLocal_ <= blockSize_)
527 // this tells us how many of the values in theta_ are valid Ritz values
528 int nevLocal_;
529 //
530 // are we implementing a skinny solver? (SkinnyIRTR)
531 bool isSkinny_;
532 //
533 // are we computing leftmost or rightmost eigenvalue?
534 bool leftMost_;
535 //
536 // State Multivecs
537 //
538 // if we are implementing a skinny solver (SkinnyIRTR),
539 // then some of these will never be allocated
540 //
541 // In order to handle auxiliary vectors, we need to handle the projector
542 // P_{[BQ BX],[BQ BX]}
543 // Using an orthomanager with B-inner product, this requires calling with multivectors
544 // [BQ,BX] and [Q,X].
545 // These multivectors must be combined because <[BQ,BX],[Q,X]>_B != I
546 // Hence, we will create two multivectors V and BV, which store
547 // V = [Q,X]
548 // BV = [BQ,BX]
549 //
550 // In the context of preconditioning, we may need to apply the projector
551 // P_{prec*[BQ,BX],[BQ,BX]}
552 // Because [BQ,BX] do not change during the supproblem solver, we will apply
553 // the preconditioner to [BQ,BX] only once. This result is stored in PBV.
554 //
555 // X,BX are views into V,BV
556 // We don't need views for Q,BQ
557 // Inside the subproblem solver, X,BX are constant, so we can allow these
558 // views to exist alongside the full view of V,BV without violating
559 // view semantics.
560 //
561 // Skinny solver allocates 6/7/8 multivectors:
562 // V_, BV_ (if hasB)
563 // PBV_ (if hasPrec and olsenPrec)
564 // R_, Z_ (regardless of hasPrec)
565 // eta_, delta_, Hdelta_
566 //
567 // Hefty solver allocates 10/11/12/13 multivectors:
568 // V_, AX_, BV_ (if hasB)
569 // PBV_ (if hasPrec and olsenPrec)
570 // R_, Z_ (if hasPrec)
571 // eta_, Aeta_, Beta_
572 // delta_, Adelta_, Bdelta_, Hdelta_
573 //
574 Teuchos::RCP<MV> V_, BV_, PBV_, // V = [Q,X]; B*V; Prec*B*V
575 AX_, R_, // A*X_; residual, gradient, and residual of model minimization
576 eta_, Aeta_, Beta_, // update vector from model minimization
577 delta_, Adelta_, Bdelta_, Hdelta_, // search direction in model minimization
578 Z_; // preconditioned residual
579 Teuchos::RCP<const MV> X_, BX_;
580 //
581 // auxiliary vectors
582 Teuchos::Array<Teuchos::RCP<const MV> > auxVecs_;
583 int numAuxVecs_;
584 //
585 // Number of iterations that have been performed.
586 int iter_;
587 //
588 // Current eigenvalues, residual norms
589 std::vector<MagnitudeType> theta_, Rnorms_, R2norms_, ritz2norms_;
590 //
591 // are the residual norms current with the residual?
592 bool Rnorms_current_, R2norms_current_;
593 //
594 // parameters solver and inner solver
595 MagnitudeType conv_kappa_, conv_theta_;
596 MagnitudeType rho_;
597 //
598 // current objective function value
599 MagnitudeType fx_;
600 };
601
602
603
604
606 // Constructor
607 template <class ScalarType, class MV, class OP>
609 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
610 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
611 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
612 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
613 const Teuchos::RCP<GenOrthoManager<ScalarType,MV,OP> > &ortho,
614 Teuchos::ParameterList &params,
615 const std::string &solverLabel, bool skinnySolver
616 ) :
617 ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
618 ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
619 NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan()),
620 // problem, tools
621 problem_(problem),
622 sm_(sorter),
623 om_(printer),
624 tester_(tester),
625 orthman_(ortho),
626#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
627 timerAOp_(Teuchos::TimeMonitor::getNewTimer("Anasazi: "+solverLabel+"::Operation A*x")),
628 timerBOp_(Teuchos::TimeMonitor::getNewTimer("Anasazi: "+solverLabel+"::Operation B*x")),
629 timerPrec_(Teuchos::TimeMonitor::getNewTimer("Anasazi: "+solverLabel+"::Operation Prec*x")),
630 timerSort_(Teuchos::TimeMonitor::getNewTimer("Anasazi: "+solverLabel+"::Sorting eigenvalues")),
631 timerLocalProj_(Teuchos::TimeMonitor::getNewTimer("Anasazi: "+solverLabel+"::Local projection")),
632 timerDS_(Teuchos::TimeMonitor::getNewTimer("Anasazi: "+solverLabel+"::Direct solve")),
633 timerLocalUpdate_(Teuchos::TimeMonitor::getNewTimer("Anasazi: "+solverLabel+"::Local update")),
634 timerCompRes_(Teuchos::TimeMonitor::getNewTimer("Anasazi: "+solverLabel+"::Computing residuals")),
635 timerOrtho_(Teuchos::TimeMonitor::getNewTimer("Anasazi: "+solverLabel+"::Orthogonalization")),
636 timerInit_(Teuchos::TimeMonitor::getNewTimer("Anasazi: "+solverLabel+"::Initialization")),
637#endif
638 counterAOp_(0),
639 counterBOp_(0),
640 counterPrec_(0),
641 // internal data
642 blockSize_(0),
643 initialized_(false),
644 nevLocal_(0),
645 isSkinny_(skinnySolver),
646 leftMost_(true),
647 auxVecs_( Teuchos::Array<Teuchos::RCP<const MV> >(0) ),
648 numAuxVecs_(0),
649 iter_(0),
650 Rnorms_current_(false),
651 R2norms_current_(false),
652 conv_kappa_(.1),
653 conv_theta_(1),
654 rho_( MAT::nan() ),
655 fx_( MAT::nan() )
656 {
657 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument,
658 "Anasazi::RTRBase::constructor: user passed null problem pointer.");
659 TEUCHOS_TEST_FOR_EXCEPTION(sm_ == Teuchos::null,std::invalid_argument,
660 "Anasazi::RTRBase::constructor: user passed null sort manager pointer.");
661 TEUCHOS_TEST_FOR_EXCEPTION(om_ == Teuchos::null,std::invalid_argument,
662 "Anasazi::RTRBase::constructor: user passed null output manager pointer.");
663 TEUCHOS_TEST_FOR_EXCEPTION(tester_ == Teuchos::null,std::invalid_argument,
664 "Anasazi::RTRBase::constructor: user passed null status test pointer.");
665 TEUCHOS_TEST_FOR_EXCEPTION(orthman_ == Teuchos::null,std::invalid_argument,
666 "Anasazi::RTRBase::constructor: user passed null orthogonalization manager pointer.");
667 TEUCHOS_TEST_FOR_EXCEPTION(problem_->isProblemSet() == false, std::invalid_argument,
668 "Anasazi::RTRBase::constructor: problem is not set.");
669 TEUCHOS_TEST_FOR_EXCEPTION(problem_->isHermitian() == false, std::invalid_argument,
670 "Anasazi::RTRBase::constructor: problem is not Hermitian.");
671
672 // get the problem operators
673 AOp_ = problem_->getOperator();
674 TEUCHOS_TEST_FOR_EXCEPTION(AOp_ == Teuchos::null, std::invalid_argument,
675 "Anasazi::RTRBase::constructor: problem provides no A matrix.");
676 BOp_ = problem_->getM();
677 Prec_ = problem_->getPrec();
678 hasBOp_ = (BOp_ != Teuchos::null);
679 hasPrec_ = (Prec_ != Teuchos::null);
680 olsenPrec_ = params.get<bool>("Olsen Prec", true);
681
682 TEUCHOS_TEST_FOR_EXCEPTION(orthman_->getOp() != BOp_,std::invalid_argument,
683 "Anasazi::RTRBase::constructor: orthogonalization manager must use mass matrix for inner product.");
684
685 // set the block size and allocate data
686 int bs = params.get("Block Size", problem_->getNEV());
687 setBlockSize(bs);
688
689 // leftmost or rightmost?
690 leftMost_ = params.get("Leftmost",leftMost_);
691
692 conv_kappa_ = params.get("Kappa Convergence",conv_kappa_);
693 TEUCHOS_TEST_FOR_EXCEPTION(conv_kappa_ <= 0 || conv_kappa_ >= 1,std::invalid_argument,
694 "Anasazi::RTRBase::constructor: kappa must be in (0,1).");
695 conv_theta_ = params.get("Theta Convergence",conv_theta_);
696 TEUCHOS_TEST_FOR_EXCEPTION(conv_theta_ <= 0,std::invalid_argument,
697 "Anasazi::RTRBase::constructor: theta must be strictly postitive.");
698 }
699
700
702 // Set the block size and make necessary adjustments.
703 template <class ScalarType, class MV, class OP>
705 {
706 // time spent here counts towards timerInit_
707#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
708 Teuchos::TimeMonitor lcltimer( *timerInit_ );
709#endif
710
711 // This routine only allocates space; it doesn't not perform any computation
712 // if solver is initialized and size is to be decreased, take the first blockSize vectors of all to preserve state
713 // otherwise, shrink/grow/allocate space and set solver to unitialized
714
715 Teuchos::RCP<const MV> tmp;
716 // grab some Multivector to Clone
717 // in practice, getInitVec() should always provide this, but it is possible to use a
718 // Eigenproblem with nothing in getInitVec() by manually initializing with initialize();
719 // in case of that strange scenario, we will try to Clone from R_
720 // we like R_ for this, because it has minimal size (blockSize_), as opposed to V_ (blockSize_+numAuxVecs_)
721 if (blockSize_ > 0) {
722 tmp = R_;
723 }
724 else {
725 tmp = problem_->getInitVec();
726 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::logic_error,
727 "Anasazi::RTRBase::setBlockSize(): Eigenproblem did not specify initial vectors to clone from");
728 }
729
730 TEUCHOS_TEST_FOR_EXCEPTION(blockSize <= 0 || blockSize > MVT::GetGlobalLength(*tmp), std::invalid_argument,
731 "Anasazi::RTRBase::setBlockSize was passed a non-positive block size");
732
733 // last chance to quit before causing side-effects
734 if (blockSize == blockSize_) {
735 // do nothing
736 return;
737 }
738
739 // clear views
740 X_ = Teuchos::null;
741 BX_ = Teuchos::null;
742
743 // regardless of whether we preserve any data, any content in V, BV and PBV corresponding to the
744 // auxiliary vectors must be retained
745 // go ahead and do these first
746 //
747 // two cases here:
748 // * we are growing (possibly, from empty)
749 // any aux data must be copied over, nothing from X need copying
750 // * we are shrinking
751 // any aux data must be copied over, go ahead and copy X material (initialized or not)
752 //
753 if (blockSize > blockSize_)
754 {
755 // GROWING
756 // get a pointer for Q-related material, and an index for extracting/setting it
757 Teuchos::RCP<const MV> Q;
758 std::vector<int> indQ(numAuxVecs_);
759 for (int i=0; i<numAuxVecs_; i++) indQ[i] = i;
760 // if numAuxVecs_ > 0, then necessarily blockSize_ > 0 (we have already been allocated once)
761 TEUCHOS_TEST_FOR_EXCEPTION(numAuxVecs_ > 0 && blockSize_ == 0, std::logic_error,
762 "Anasazi::RTRBase::setSize(): logic error. Please contact Anasazi team.");
763 // V
764 if (numAuxVecs_ > 0) Q = MVT::CloneView(*V_,indQ);
765 V_ = MVT::Clone(*tmp,numAuxVecs_ + blockSize);
766 if (numAuxVecs_ > 0) MVT::SetBlock(*Q,indQ,*V_);
767 // BV
768 if (hasBOp_) {
769 if (numAuxVecs_ > 0) Q = MVT::CloneView(*BV_,indQ);
770 BV_ = MVT::Clone(*tmp,numAuxVecs_ + blockSize);
771 if (numAuxVecs_ > 0) MVT::SetBlock(*Q,indQ,*BV_);
772 }
773 else {
774 BV_ = V_;
775 }
776 // PBV
777 if (hasPrec_ && olsenPrec_) {
778 if (numAuxVecs_ > 0) Q = MVT::CloneView(*PBV_,indQ);
779 PBV_ = MVT::Clone(*tmp,numAuxVecs_ + blockSize);
780 if (numAuxVecs_ > 0) MVT::SetBlock(*Q,indQ,*PBV_);
781 }
782 }
783 else
784 {
785 // SHRINKING
786 std::vector<int> indV(numAuxVecs_+blockSize);
787 for (int i=0; i<numAuxVecs_+blockSize; i++) indV[i] = i;
788 // V
789 V_ = MVT::CloneCopy(*V_,indV);
790 // BV
791 if (hasBOp_) {
792 BV_ = MVT::CloneCopy(*BV_,indV);
793 }
794 else {
795 BV_ = V_;
796 }
797 // PBV
798 if (hasPrec_ && olsenPrec_) {
799 PBV_ = MVT::CloneCopy(*PBV_,indV);
800 }
801 }
802
803 if (blockSize < blockSize_) {
804 // shrink vectors
805 blockSize_ = blockSize;
806
807 theta_.resize(blockSize_);
808 ritz2norms_.resize(blockSize_);
809 Rnorms_.resize(blockSize_);
810 R2norms_.resize(blockSize_);
811
812 if (initialized_) {
813 // shrink multivectors with copy
814 std::vector<int> ind(blockSize_);
815 for (int i=0; i<blockSize_; i++) ind[i] = i;
816
817 // Z can be deleted, no important info there
818 Z_ = Teuchos::null;
819
820 // we will not use tmp for cloning, clear it and free some space
821 tmp = Teuchos::null;
822
823 R_ = MVT::CloneCopy(*R_ ,ind);
824 eta_ = MVT::CloneCopy(*eta_ ,ind);
825 delta_ = MVT::CloneCopy(*delta_ ,ind);
826 Hdelta_ = MVT::CloneCopy(*Hdelta_,ind);
827 if (!isSkinny_) {
828 AX_ = MVT::CloneCopy(*AX_ ,ind);
829 Aeta_ = MVT::CloneCopy(*Aeta_ ,ind);
830 Adelta_ = MVT::CloneCopy(*Adelta_,ind);
831 }
832 else {
833 AX_ = Teuchos::null;
834 Aeta_ = Teuchos::null;
835 Adelta_ = Teuchos::null;
836 }
837
838 if (hasBOp_) {
839 if (!isSkinny_) {
840 Beta_ = MVT::CloneCopy(*Beta_,ind);
841 Bdelta_ = MVT::CloneCopy(*Bdelta_,ind);
842 }
843 else {
844 Beta_ = Teuchos::null;
845 Bdelta_ = Teuchos::null;
846 }
847 }
848 else {
849 Beta_ = eta_;
850 Bdelta_ = delta_;
851 }
852
853 // we need Z if we have a preconditioner
854 // we also use Z for temp storage in the skinny solvers
855 if (hasPrec_ || isSkinny_) {
856 Z_ = MVT::Clone(*V_,blockSize_);
857 }
858 else {
859 Z_ = R_;
860 }
861
862 }
863 else {
864 // release previous multivectors
865 // shrink multivectors without copying
866 AX_ = Teuchos::null;
867 R_ = Teuchos::null;
868 eta_ = Teuchos::null;
869 Aeta_ = Teuchos::null;
870 delta_ = Teuchos::null;
871 Adelta_ = Teuchos::null;
872 Hdelta_ = Teuchos::null;
873 Beta_ = Teuchos::null;
874 Bdelta_ = Teuchos::null;
875 Z_ = Teuchos::null;
876
877 R_ = MVT::Clone(*tmp,blockSize_);
878 eta_ = MVT::Clone(*tmp,blockSize_);
879 delta_ = MVT::Clone(*tmp,blockSize_);
880 Hdelta_ = MVT::Clone(*tmp,blockSize_);
881 if (!isSkinny_) {
882 AX_ = MVT::Clone(*tmp,blockSize_);
883 Aeta_ = MVT::Clone(*tmp,blockSize_);
884 Adelta_ = MVT::Clone(*tmp,blockSize_);
885 }
886
887 if (hasBOp_) {
888 if (!isSkinny_) {
889 Beta_ = MVT::Clone(*tmp,blockSize_);
890 Bdelta_ = MVT::Clone(*tmp,blockSize_);
891 }
892 }
893 else {
894 Beta_ = eta_;
895 Bdelta_ = delta_;
896 }
897
898 // we need Z if we have a preconditioner
899 // we also use Z for temp storage in the skinny solvers
900 if (hasPrec_ || isSkinny_) {
901 Z_ = MVT::Clone(*tmp,blockSize_);
902 }
903 else {
904 Z_ = R_;
905 }
906 } // if initialized_
907 } // if blockSize is shrinking
908 else { // if blockSize > blockSize_
909 // this is also the scenario for our initial call to setBlockSize(), in the constructor
910 initialized_ = false;
911
912 // grow/allocate vectors
913 theta_.resize(blockSize,NANVAL);
914 ritz2norms_.resize(blockSize,NANVAL);
915 Rnorms_.resize(blockSize,NANVAL);
916 R2norms_.resize(blockSize,NANVAL);
917
918 // deallocate old multivectors
919 AX_ = Teuchos::null;
920 R_ = Teuchos::null;
921 eta_ = Teuchos::null;
922 Aeta_ = Teuchos::null;
923 delta_ = Teuchos::null;
924 Adelta_ = Teuchos::null;
925 Hdelta_ = Teuchos::null;
926 Beta_ = Teuchos::null;
927 Bdelta_ = Teuchos::null;
928 Z_ = Teuchos::null;
929
930 // clone multivectors off of tmp
931 R_ = MVT::Clone(*tmp,blockSize);
932 eta_ = MVT::Clone(*tmp,blockSize);
933 delta_ = MVT::Clone(*tmp,blockSize);
934 Hdelta_ = MVT::Clone(*tmp,blockSize);
935 if (!isSkinny_) {
936 AX_ = MVT::Clone(*tmp,blockSize);
937 Aeta_ = MVT::Clone(*tmp,blockSize);
938 Adelta_ = MVT::Clone(*tmp,blockSize);
939 }
940
941 if (hasBOp_) {
942 if (!isSkinny_) {
943 Beta_ = MVT::Clone(*tmp,blockSize);
944 Bdelta_ = MVT::Clone(*tmp,blockSize);
945 }
946 }
947 else {
948 Beta_ = eta_;
949 Bdelta_ = delta_;
950 }
951 if (hasPrec_ || isSkinny_) {
952 Z_ = MVT::Clone(*tmp,blockSize);
953 }
954 else {
955 Z_ = R_;
956 }
957 blockSize_ = blockSize;
958 }
959
960 // get view of X from V, BX from BV
961 // these are located after the first numAuxVecs columns
962 {
963 std::vector<int> indX(blockSize_);
964 for (int i=0; i<blockSize_; i++) indX[i] = numAuxVecs_+i;
965 X_ = MVT::CloneView(*V_,indX);
966 if (hasBOp_) {
967 BX_ = MVT::CloneView(*BV_,indX);
968 }
969 else {
970 BX_ = X_;
971 }
972 }
973 }
974
975
977 // Set a new StatusTest for the solver.
978 template <class ScalarType, class MV, class OP>
980 TEUCHOS_TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument,
981 "Anasazi::RTRBase::setStatusTest() was passed a null StatusTest.");
982 tester_ = test;
983 }
984
985
987 // Get the current StatusTest used by the solver.
988 template <class ScalarType, class MV, class OP>
989 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > RTRBase<ScalarType,MV,OP>::getStatusTest() const {
990 return tester_;
991 }
992
993
995 // Set the auxiliary vectors
996 template <class ScalarType, class MV, class OP>
997 void RTRBase<ScalarType,MV,OP>::setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs) {
998 typedef typename Teuchos::Array<Teuchos::RCP<const MV> >::const_iterator tarcpmv;
999
1000 // set new auxiliary vectors
1001 auxVecs_.resize(0);
1002 auxVecs_.reserve(auxvecs.size());
1003
1004 numAuxVecs_ = 0;
1005 for (tarcpmv v=auxvecs.begin(); v != auxvecs.end(); ++v) {
1006 numAuxVecs_ += MVT::GetNumberVecs(**v);
1007 }
1008
1009 // If the solver has been initialized, X is not necessarily orthogonal to new auxiliary vectors
1010 if (numAuxVecs_ > 0 && initialized_) {
1011 initialized_ = false;
1012 }
1013
1014 // clear X,BX views
1015 X_ = Teuchos::null;
1016 BX_ = Teuchos::null;
1017 // deallocate, we'll clone off R_ below
1018 V_ = Teuchos::null;
1019 BV_ = Teuchos::null;
1020 PBV_ = Teuchos::null;
1021
1022 // put auxvecs into V, update BV and PBV if necessary
1023 if (numAuxVecs_ > 0) {
1024 V_ = MVT::Clone(*R_,numAuxVecs_ + blockSize_);
1025 int numsofar = 0;
1026 for (tarcpmv v=auxvecs.begin(); v != auxvecs.end(); ++v) {
1027 std::vector<int> ind(MVT::GetNumberVecs(**v));
1028 for (size_t j=0; j<ind.size(); j++) ind[j] = numsofar++;
1029 MVT::SetBlock(**v,ind,*V_);
1030 auxVecs_.push_back(MVT::CloneView(*Teuchos::rcp_static_cast<const MV>(V_),ind));
1031 }
1032 TEUCHOS_TEST_FOR_EXCEPTION(numsofar != numAuxVecs_, std::logic_error,
1033 "Anasazi::RTRBase::setAuxVecs(): Logic error. Please contact Anasazi team.");
1034 // compute B*V, Prec*B*V
1035 if (hasBOp_) {
1036 BV_ = MVT::Clone(*R_,numAuxVecs_ + blockSize_);
1037 OPT::Apply(*BOp_,*V_,*BV_);
1038 }
1039 else {
1040 BV_ = V_;
1041 }
1042 if (hasPrec_ && olsenPrec_) {
1043 PBV_ = MVT::Clone(*R_,numAuxVecs_ + blockSize_);
1044 OPT::Apply(*Prec_,*BV_,*V_);
1045 }
1046 }
1047 //
1048
1049 if (om_->isVerbosity( Debug ) ) {
1050 // Check almost everything here
1051 CheckList chk;
1052 chk.checkQ = true;
1053 om_->print( Debug, accuracyCheck(chk, "in setAuxVecs()") );
1054 }
1055 }
1056
1057
1059 /* Initialize the state of the solver
1060 *
1061 * POST-CONDITIONS:
1062 *
1063 * initialized_ == true
1064 * X is orthonormal, orthogonal to auxVecs_
1065 * AX = A*X if not skinnny
1066 * BX = B*X if hasBOp_
1067 * theta_ contains Ritz values of X
1068 * R = AX - BX*diag(theta_)
1069 */
1070 template <class ScalarType, class MV, class OP>
1072 {
1073 // NOTE: memory has been allocated by setBlockSize(). Use SetBlock below; do not Clone
1074 // NOTE: Overall time spent in this routine is counted to timerInit_; portions will also be counted towards other primitives
1075
1076 // clear const views to X,BX in V,BV
1077 // work with temporary non-const views
1078 X_ = Teuchos::null;
1079 BX_ = Teuchos::null;
1080 Teuchos::RCP<MV> X, BX;
1081 {
1082 std::vector<int> ind(blockSize_);
1083 for (int i=0; i<blockSize_; ++i) ind[i] = numAuxVecs_+i;
1084 X = MVT::CloneViewNonConst(*V_,ind);
1085 if (hasBOp_) {
1086 BX = MVT::CloneViewNonConst(*BV_,ind);
1087 }
1088 else {
1089 BX = X;
1090 }
1091 }
1092
1093#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1094 Teuchos::TimeMonitor inittimer( *timerInit_ );
1095#endif
1096
1097 std::vector<int> bsind(blockSize_);
1098 for (int i=0; i<blockSize_; i++) bsind[i] = i;
1099
1100 // in RTR, X (the subspace iterate) is primary
1101 // the order of dependence follows like so.
1102 // --init-> X
1103 // --op apply-> AX,BX
1104 // --ritz analysis-> theta
1105 //
1106 // if the user specifies all data for a level, we will accept it.
1107 // otherwise, we will generate the whole level, and all subsequent levels.
1108 //
1109 // the data members are ordered based on dependence, and the levels are
1110 // partitioned according to the amount of work required to produce the
1111 // items in a level.
1112 //
1113 // inconsitent multivectors widths and lengths will not be tolerated, and
1114 // will be treated with exceptions.
1115
1116 // set up X, AX, BX: get them from "state" if user specified them
1117 if (newstate.X != Teuchos::null) {
1118 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.X) != MVT::GetGlobalLength(*X),
1119 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): vector length of newstate.X not correct." );
1120 // newstate.X must have blockSize_ vectors; any more will be ignored
1121 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.X) < blockSize_,
1122 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.X must have at least block size vectors.");
1123
1124 // put data in X
1125 MVT::SetBlock(*newstate.X,bsind,*X);
1126
1127 // put data in AX
1128 // if we are implementing a skinny solver, then we don't have memory allocated for AX
1129 // in this case, point AX at Z (skinny solvers allocate Z) and use it for temporary storage
1130 // we will clear the pointer later
1131 if (isSkinny_) {
1132 AX_ = Z_;
1133 }
1134 if (newstate.AX != Teuchos::null) {
1135 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.AX) != MVT::GetGlobalLength(*X),
1136 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): vector length of newstate.AX not correct." );
1137 // newstate.AX must have blockSize_ vectors; any more will be ignored
1138 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.AX) < blockSize_,
1139 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.AX must have at least block size vectors.");
1140 MVT::SetBlock(*newstate.AX,bsind,*AX_);
1141 }
1142 else {
1143 {
1144#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1145 Teuchos::TimeMonitor lcltimer( *timerAOp_ );
1146#endif
1147 OPT::Apply(*AOp_,*X,*AX_);
1148 counterAOp_ += blockSize_;
1149 }
1150 // we generated AX; we will generate R as well
1151 newstate.R = Teuchos::null;
1152 }
1153
1154 // put data in BX
1155 // skinny solvers always allocate BX if hasB, so this is unconditionally appropriate
1156 if (hasBOp_) {
1157 if (newstate.BX != Teuchos::null) {
1158 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.BX) != MVT::GetGlobalLength(*X),
1159 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): vector length of newstate.BX not correct." );
1160 // newstate.BX must have blockSize_ vectors; any more will be ignored
1161 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.BX) < blockSize_,
1162 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.BX must have at least block size vectors.");
1163 MVT::SetBlock(*newstate.BX,bsind,*BX);
1164 }
1165 else {
1166 {
1167#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1168 Teuchos::TimeMonitor lcltimer( *timerBOp_ );
1169#endif
1170 OPT::Apply(*BOp_,*X,*BX);
1171 counterBOp_ += blockSize_;
1172 }
1173 // we generated BX; we will generate R as well
1174 newstate.R = Teuchos::null;
1175 }
1176 }
1177 else {
1178 // the assignment BX_==X_ would be redundant; take advantage of this opportunity to debug a little
1179 TEUCHOS_TEST_FOR_EXCEPTION(BX != X, std::logic_error, "Anasazi::RTRBase::initialize(): solver invariant not satisfied (BX==X).");
1180 }
1181
1182 }
1183 else {
1184 // user did not specify X
1185
1186 // clear state so we won't use any data from it below
1187 newstate.R = Teuchos::null;
1188 newstate.T = Teuchos::null;
1189
1190 // generate something and projectAndNormalize
1191 Teuchos::RCP<const MV> ivec = problem_->getInitVec();
1192 TEUCHOS_TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::logic_error,
1193 "Anasazi::RTRBase::initialize(): Eigenproblem did not specify initial vectors to clone from.");
1194
1195 int initSize = MVT::GetNumberVecs(*ivec);
1196 if (initSize > blockSize_) {
1197 // we need only the first blockSize_ vectors from ivec; get a view of them
1198 initSize = blockSize_;
1199 std::vector<int> ind(blockSize_);
1200 for (int i=0; i<blockSize_; i++) ind[i] = i;
1201 ivec = MVT::CloneView(*ivec,ind);
1202 }
1203
1204 // assign ivec to first part of X
1205 if (initSize > 0) {
1206 std::vector<int> ind(initSize);
1207 for (int i=0; i<initSize; i++) ind[i] = i;
1208 MVT::SetBlock(*ivec,ind,*X);
1209 }
1210 // fill the rest of X with random
1211 if (blockSize_ > initSize) {
1212 std::vector<int> ind(blockSize_ - initSize);
1213 for (int i=0; i<blockSize_ - initSize; i++) ind[i] = initSize + i;
1214 Teuchos::RCP<MV> rX = MVT::CloneViewNonConst(*X,ind);
1215 MVT::MvRandom(*rX);
1216 rX = Teuchos::null;
1217 }
1218
1219 // put data in BX
1220 if (hasBOp_) {
1221#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1222 Teuchos::TimeMonitor lcltimer( *timerBOp_ );
1223#endif
1224 OPT::Apply(*BOp_,*X,*BX);
1225 counterBOp_ += blockSize_;
1226 }
1227 else {
1228 // the assignment BX==X would be redundant; take advantage of this opportunity to debug a little
1229 TEUCHOS_TEST_FOR_EXCEPTION(BX != X, std::logic_error, "Anasazi::RTRBase::initialize(): solver invariant not satisfied (BX==X).");
1230 }
1231
1232 // remove auxVecs from X and normalize it
1233 if (numAuxVecs_ > 0) {
1234#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1235 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1236#endif
1237 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
1238 int rank = orthman_->projectAndNormalizeMat(*X,auxVecs_,dummyC,Teuchos::null,BX);
1239 TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_, RTRInitFailure,
1240 "Anasazi::RTRBase::initialize(): Couldn't generate initial basis of full rank.");
1241 }
1242 else {
1243#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1244 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1245#endif
1246 int rank = orthman_->normalizeMat(*X,Teuchos::null,BX);
1247 TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_, RTRInitFailure,
1248 "Anasazi::RTRBase::initialize(): Couldn't generate initial basis of full rank.");
1249 }
1250
1251 // put data in AX
1252 if (isSkinny_) {
1253 AX_ = Z_;
1254 }
1255 {
1256#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1257 Teuchos::TimeMonitor lcltimer( *timerAOp_ );
1258#endif
1259 OPT::Apply(*AOp_,*X,*AX_);
1260 counterAOp_ += blockSize_;
1261 }
1262
1263 } // end if (newstate.X != Teuchos::null)
1264
1265
1266 // set up Ritz values
1267 if (newstate.T != Teuchos::null) {
1268 TEUCHOS_TEST_FOR_EXCEPTION( (signed int)(newstate.T->size()) < blockSize_,
1269 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.T must contain at least block size Ritz values.");
1270 for (int i=0; i<blockSize_; i++) {
1271 theta_[i] = (*newstate.T)[i];
1272 }
1273 }
1274 else {
1275 // get ritz vecs/vals
1276 Teuchos::SerialDenseMatrix<int,ScalarType> AA(blockSize_,blockSize_),
1277 BB(blockSize_,blockSize_),
1278 S(blockSize_,blockSize_);
1279 // project A
1280 {
1281#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1282 Teuchos::TimeMonitor lcltimer( *timerLocalProj_ );
1283#endif
1284 MVT::MvTransMv(ONE,*X,*AX_,AA);
1285 if (hasBOp_) {
1286 MVT::MvTransMv(ONE,*X,*BX,BB);
1287 }
1288 }
1289 nevLocal_ = blockSize_;
1290
1291 // solve the projected problem
1292 // project B
1293 int ret;
1294 if (hasBOp_) {
1295#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1296 Teuchos::TimeMonitor lcltimer( *timerDS_ );
1297#endif
1298 ret = Utils::directSolver(blockSize_,AA,Teuchos::rcpFromRef(BB),S,theta_,nevLocal_,1);
1299 }
1300 else {
1301#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1302 Teuchos::TimeMonitor lcltimer( *timerDS_ );
1303#endif
1304 ret = Utils::directSolver(blockSize_,AA,Teuchos::null,S,theta_,nevLocal_,10);
1305 }
1306 TEUCHOS_TEST_FOR_EXCEPTION(ret != 0,RTRInitFailure,
1307 "Anasazi::RTRBase::initialize(): failure solving projected eigenproblem after retraction. LAPACK returns " << ret);
1308 TEUCHOS_TEST_FOR_EXCEPTION(nevLocal_ != blockSize_,RTRInitFailure,"Anasazi::RTRBase::initialize(): retracted iterate failed in Ritz analysis.");
1309
1310 // We only have blockSize_ ritz pairs, ergo we do not need to select.
1311 // However, we still require them to be ordered correctly
1312 {
1313#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1314 Teuchos::TimeMonitor lcltimer( *timerSort_ );
1315#endif
1316
1317 std::vector<int> order(blockSize_);
1318 //
1319 // sort the first blockSize_ values in theta_
1320 sm_->sort(theta_, Teuchos::rcpFromRef(order), blockSize_); // don't catch exception
1321 //
1322 // apply the same ordering to the primitive ritz vectors
1323 Utils::permuteVectors(order,S);
1324 }
1325
1326 // compute ritz residual norms
1327 {
1328 Teuchos::BLAS<int,ScalarType> blas;
1329 Teuchos::SerialDenseMatrix<int,ScalarType> RR(blockSize_,blockSize_);
1330 // RR = AA*S - BB*S*diag(theta)
1331 int info;
1332 if (hasBOp_) {
1333 info = RR.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,BB,S,ZERO);
1334 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error, "Anasazi::RTRBase::initialize(): Logic error calling SerialDenseMatrix::multiply.");
1335 }
1336 else {
1337 RR.assign(S);
1338 }
1339 for (int i=0; i<blockSize_; i++) {
1340 blas.SCAL(blockSize_,theta_[i],RR[i],1);
1341 }
1342 info = RR.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,AA,S,-ONE);
1343 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error, "Anasazi::RTRBase::initialize(): Logic error calling SerialDenseMatrix::multiply.");
1344 for (int i=0; i<blockSize_; i++) {
1345 ritz2norms_[i] = blas.NRM2(blockSize_,RR[i],1);
1346 }
1347 }
1348
1349 // update the solution
1350 // use R as local work space
1351 // Z may already be in use as work space (holding AX if isSkinny)
1352 {
1353#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1354 Teuchos::TimeMonitor lcltimer( *timerLocalUpdate_ );
1355#endif
1356 // X <- X*S
1357 MVT::MvAddMv( ONE, *X, ZERO, *X, *R_ );
1358 MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *X );
1359 // AX <- AX*S
1360 MVT::MvAddMv( ONE, *AX_, ZERO, *AX_, *R_ );
1361 MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *AX_ );
1362 if (hasBOp_) {
1363 // BX <- BX*S
1364 MVT::MvAddMv( ONE, *BX, ZERO, *BX, *R_ );
1365 MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *BX );
1366 }
1367 }
1368 }
1369
1370 // done modifying X,BX; clear temp views and setup const views
1371 X = Teuchos::null;
1372 BX = Teuchos::null;
1373 {
1374 std::vector<int> ind(blockSize_);
1375 for (int i=0; i<blockSize_; ++i) ind[i] = numAuxVecs_+i;
1376 this->X_ = MVT::CloneView(static_cast<const MV&>(*V_),ind);
1377 if (hasBOp_) {
1378 this->BX_ = MVT::CloneView(static_cast<const MV&>(*BV_),ind);
1379 }
1380 else {
1381 this->BX_ = this->X_;
1382 }
1383 }
1384
1385
1386 // get objective function value
1387 fx_ = std::accumulate(theta_.begin(),theta_.end(),ZERO);
1388
1389 // set up R
1390 if (newstate.R != Teuchos::null) {
1391 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) < blockSize_,
1392 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): newstate.R must have blockSize number of vectors." );
1393 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
1394 std::invalid_argument, "Anasazi::RTRBase::initialize(newstate): vector length of newstate.R not correct." );
1395 MVT::SetBlock(*newstate.R,bsind,*R_);
1396 }
1397 else {
1398#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1399 Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
1400#endif
1401 // form R <- AX - BX*T
1402 MVT::MvAddMv(ZERO,*AX_,ONE,*AX_,*R_);
1403 Teuchos::SerialDenseMatrix<int,ScalarType> T(blockSize_,blockSize_);
1404 T.putScalar(ZERO);
1405 for (int i=0; i<blockSize_; i++) T(i,i) = theta_[i];
1406 MVT::MvTimesMatAddMv(-ONE,*BX_,T,ONE,*R_);
1407 }
1408
1409 // R has been updated; mark the norms as out-of-date
1410 Rnorms_current_ = false;
1411 R2norms_current_ = false;
1412
1413 // if isSkinny, then AX currently points to Z for temp storage
1414 // set AX back to null
1415 if (isSkinny_) {
1416 AX_ = Teuchos::null;
1417 }
1418
1419 // finally, we are initialized
1420 initialized_ = true;
1421
1422 if (om_->isVerbosity( Debug ) ) {
1423 // Check almost everything here
1424 CheckList chk;
1425 chk.checkX = true;
1426 chk.checkAX = true;
1427 chk.checkBX = true;
1428 chk.checkR = true;
1429 chk.checkQ = true;
1430 om_->print( Debug, accuracyCheck(chk, "after initialize()") );
1431 }
1432 }
1433
1434 template <class ScalarType, class MV, class OP>
1436 {
1438 initialize(empty);
1439 }
1440
1441
1442
1443
1445 // compute/return residual B-norms
1446 template <class ScalarType, class MV, class OP>
1447 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
1449 if (Rnorms_current_ == false) {
1450 // Update the residual norms
1451 orthman_->norm(*R_,Rnorms_);
1452 Rnorms_current_ = true;
1453 }
1454 return Rnorms_;
1455 }
1456
1457
1459 // compute/return residual 2-norms
1460 template <class ScalarType, class MV, class OP>
1461 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
1463 if (R2norms_current_ == false) {
1464 // Update the residual 2-norms
1465 MVT::MvNorm(*R_,R2norms_);
1466 R2norms_current_ = true;
1467 }
1468 return R2norms_;
1469 }
1470
1471
1472
1473
1475 // Check accuracy, orthogonality, and other debugging stuff
1476 //
1477 // bools specify which tests we want to run (instead of running more than we actually care about)
1478 //
1479 // we don't bother checking the following because they are computed explicitly:
1480 // AH == A*H
1481 //
1482 //
1483 // checkX : X orthonormal
1484 // orthogonal to auxvecs
1485 // checkAX : check AX == A*X
1486 // checkBX : check BX == B*X
1487 // checkEta : check that Eta'*B*X == 0
1488 // orthogonal to auxvecs
1489 // checkAEta : check that AEta = A*Eta
1490 // checkBEta : check that BEta = B*Eta
1491 // checkR : check R orthogonal to X
1492 // checkBR : check R B-orthogonal to X, valid in and immediately after solveTRSubproblem
1493 // checkQ : check that auxiliary vectors are actually orthonormal
1494 //
1495 // TODO:
1496 // add checkTheta
1497 //
1498 template <class ScalarType, class MV, class OP>
1499 std::string RTRBase<ScalarType,MV,OP>::accuracyCheck( const CheckList &chk, const std::string &where ) const
1500 {
1501 using std::setprecision;
1502 using std::scientific;
1503 using std::setw;
1504 using std::endl;
1505 std::stringstream os;
1506 MagnitudeType tmp;
1507
1508 os << " Debugging checks: " << where << endl;
1509
1510 // X and friends
1511 if (chk.checkX && initialized_) {
1512 tmp = orthman_->orthonormError(*X_);
1513 os << " >> Error in X^H B X == I : " << scientific << setprecision(10) << tmp << endl;
1514 for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1515 tmp = orthman_->orthogError(*X_,*auxVecs_[i]);
1516 os << " >> Error in X^H B Q[" << i << "] == 0 : " << scientific << setprecision(10) << tmp << endl;
1517 }
1518 }
1519 if (chk.checkAX && !isSkinny_ && initialized_) {
1520 tmp = Utils::errorEquality(*X_, *AX_, AOp_);
1521 os << " >> Error in AX == A*X : " << scientific << setprecision(10) << tmp << endl;
1522 }
1523 if (chk.checkBX && hasBOp_ && initialized_) {
1524 Teuchos::RCP<MV> BX = MVT::Clone(*this->X_,this->blockSize_);
1525 OPT::Apply(*BOp_,*this->X_,*BX);
1526 tmp = Utils::errorEquality(*BX, *BX_);
1527 os << " >> Error in BX == B*X : " << scientific << setprecision(10) << tmp << endl;
1528 }
1529
1530 // Eta and friends
1531 if (chk.checkEta && initialized_) {
1532 tmp = orthman_->orthogError(*X_,*eta_);
1533 os << " >> Error in X^H B Eta == 0 : " << scientific << setprecision(10) << tmp << endl;
1534 for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1535 tmp = orthman_->orthogError(*eta_,*auxVecs_[i]);
1536 os << " >> Error in Eta^H B Q[" << i << "]==0 : " << scientific << setprecision(10) << tmp << endl;
1537 }
1538 }
1539 if (chk.checkAEta && !isSkinny_ && initialized_) {
1540 tmp = Utils::errorEquality(*eta_, *Aeta_, AOp_);
1541 os << " >> Error in AEta == A*Eta : " << scientific << setprecision(10) << tmp << endl;
1542 }
1543 if (chk.checkBEta && !isSkinny_ && hasBOp_ && initialized_) {
1544 tmp = Utils::errorEquality(*eta_, *Beta_, BOp_);
1545 os << " >> Error in BEta == B*Eta : " << scientific << setprecision(10) << tmp << endl;
1546 }
1547
1548 // R: this is not B-orthogonality, but standard euclidean orthogonality
1549 if (chk.checkR && initialized_) {
1550 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
1551 MVT::MvTransMv(ONE,*X_,*R_,xTx);
1552 tmp = xTx.normFrobenius();
1553 os << " >> Error in X^H R == 0 : " << scientific << setprecision(10) << tmp << endl;
1554 }
1555
1556 // BR: this is B-orthogonality: this is only valid inside and immediately after solveTRSubproblem
1557 // only check if B != I, otherwise, it is equivalent to the above test
1558 if (chk.checkBR && hasBOp_ && initialized_) {
1559 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
1560 MVT::MvTransMv(ONE,*BX_,*R_,xTx);
1561 tmp = xTx.normFrobenius();
1562 os << " >> Error in X^H B R == 0 : " << scientific << setprecision(10) << tmp << endl;
1563 }
1564
1565 // Z: Z is preconditioned R, should be on tangent plane
1566 if (chk.checkZ && initialized_) {
1567 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
1568 MVT::MvTransMv(ONE,*BX_,*Z_,xTx);
1569 tmp = xTx.normFrobenius();
1570 os << " >> Error in X^H B Z == 0 : " << scientific << setprecision(10) << tmp << endl;
1571 }
1572
1573 // Q
1574 if (chk.checkQ) {
1575 for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1576 tmp = orthman_->orthonormError(*auxVecs_[i]);
1577 os << " >> Error in Q[" << i << "]^H B Q[" << i << "]==I: " << scientific << setprecision(10) << tmp << endl;
1578 for (Array_size_type j=i+1; j<auxVecs_.size(); j++) {
1579 tmp = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
1580 os << " >> Error in Q[" << i << "]^H B Q[" << j << "]==0: " << scientific << setprecision(10) << tmp << endl;
1581 }
1582 }
1583 }
1584 os << endl;
1585 return os.str();
1586 }
1587
1588
1590 // Print the current status of the solver
1591 template <class ScalarType, class MV, class OP>
1592 void
1594 {
1595 using std::setprecision;
1596 using std::scientific;
1597 using std::setw;
1598 using std::endl;
1599
1600 os <<endl;
1601 os <<"================================================================================" << endl;
1602 os << endl;
1603 os <<" RTRBase Solver Status" << endl;
1604 os << endl;
1605 os <<"The solver is "<<(initialized_ ? "initialized." : "not initialized.") << endl;
1606 os <<"The number of iterations performed is " << iter_ << endl;
1607 os <<"The current block size is " << blockSize_ << endl;
1608 os <<"The number of auxiliary vectors is " << numAuxVecs_ << endl;
1609 os <<"The number of operations A*x is " << counterAOp_ << endl;
1610 os <<"The number of operations B*x is " << counterBOp_ << endl;
1611 os <<"The number of operations Prec*x is " << counterPrec_ << endl;
1612 os <<"The most recent rho was " << scientific << setprecision(10) << rho_ << endl;
1613 os <<"The current value of f(x) is " << scientific << setprecision(10) << fx_ << endl;
1614
1615 if (initialized_) {
1616 os << endl;
1617 os <<"CURRENT EIGENVALUE ESTIMATES "<<endl;
1618 os << setw(20) << "Eigenvalue"
1619 << setw(20) << "Residual(B)"
1620 << setw(20) << "Residual(2)"
1621 << endl;
1622 os <<"--------------------------------------------------------------------------------"<<endl;
1623 for (int i=0; i<blockSize_; i++) {
1624 os << scientific << setprecision(10) << setw(20) << theta_[i];
1625 if (Rnorms_current_) os << scientific << setprecision(10) << setw(20) << Rnorms_[i];
1626 else os << scientific << setprecision(10) << setw(20) << "not current";
1627 if (R2norms_current_) os << scientific << setprecision(10) << setw(20) << R2norms_[i];
1628 else os << scientific << setprecision(10) << setw(20) << "not current";
1629 os << endl;
1630 }
1631 }
1632 os <<"================================================================================" << endl;
1633 os << endl;
1634 }
1635
1636
1638 // Inner product 1
1639 template <class ScalarType, class MV, class OP>
1640 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
1641 RTRBase<ScalarType,MV,OP>::ginner(const MV &xi) const
1642 {
1643 std::vector<MagnitudeType> d(MVT::GetNumberVecs(xi));
1644 MVT::MvNorm(xi,d);
1645 MagnitudeType ret = 0;
1646 for (vecMTiter i=d.begin(); i != d.end(); ++i) {
1647 ret += (*i)*(*i);
1648 }
1649 return ret;
1650 }
1651
1652
1654 // Inner product 2
1655 template <class ScalarType, class MV, class OP>
1656 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
1657 RTRBase<ScalarType,MV,OP>::ginner(const MV &xi, const MV &zeta) const
1658 {
1659 std::vector<ScalarType> d(MVT::GetNumberVecs(xi));
1660 MVT::MvDot(xi,zeta,d);
1661 return SCT::real(std::accumulate(d.begin(),d.end(),SCT::zero()));
1662 }
1663
1664
1666 // Inner product 1 without trace accumulation
1667 template <class ScalarType, class MV, class OP>
1668 void RTRBase<ScalarType,MV,OP>::ginnersep(
1669 const MV &xi,
1670 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &d) const
1671 {
1672 MVT::MvNorm(xi,d);
1673 for (vecMTiter i=d.begin(); i != d.end(); ++i) {
1674 (*i) = (*i)*(*i);
1675 }
1676 }
1677
1678
1680 // Inner product 2 without trace accumulation
1681 template <class ScalarType, class MV, class OP>
1682 void RTRBase<ScalarType,MV,OP>::ginnersep(
1683 const MV &xi, const MV &zeta,
1684 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &d) const
1685 {
1686 std::vector<ScalarType> dC(MVT::GetNumberVecs(xi));
1687 MVT::MvDot(xi,zeta,dC);
1688 vecMTiter iR=d.begin();
1689 vecSTiter iS=dC.begin();
1690 for (; iR != d.end(); iR++, iS++) {
1691 (*iR) = SCT::real(*iS);
1692 }
1693 }
1694
1695 template <class ScalarType, class MV, class OP>
1696 Teuchos::Array<Teuchos::RCP<const MV> > RTRBase<ScalarType,MV,OP>::getAuxVecs() const {
1697 return auxVecs_;
1698 }
1699
1700 template <class ScalarType, class MV, class OP>
1702 return(blockSize_);
1703 }
1704
1705 template <class ScalarType, class MV, class OP>
1707 return(*problem_);
1708 }
1709
1710 template <class ScalarType, class MV, class OP>
1712 return blockSize_;
1713 }
1714
1715 template <class ScalarType, class MV, class OP>
1717 {
1718 if (!initialized_) return 0;
1719 return nevLocal_;
1720 }
1721
1722 template <class ScalarType, class MV, class OP>
1723 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
1725 {
1726 std::vector<MagnitudeType> ret = ritz2norms_;
1727 ret.resize(nevLocal_);
1728 return ret;
1729 }
1730
1731 template <class ScalarType, class MV, class OP>
1732 std::vector<Value<ScalarType> >
1734 {
1735 std::vector<Value<ScalarType> > ret(nevLocal_);
1736 for (int i=0; i<nevLocal_; i++) {
1737 ret[i].realpart = theta_[i];
1738 ret[i].imagpart = ZERO;
1739 }
1740 return ret;
1741 }
1742
1743 template <class ScalarType, class MV, class OP>
1744 Teuchos::RCP<const MV>
1746 {
1747 return X_;
1748 }
1749
1750 template <class ScalarType, class MV, class OP>
1752 {
1753 iter_=0;
1754 }
1755
1756 template <class ScalarType, class MV, class OP>
1758 {
1759 return iter_;
1760 }
1761
1762 template <class ScalarType, class MV, class OP>
1764 {
1766 state.X = X_;
1767 state.AX = AX_;
1768 if (hasBOp_) {
1769 state.BX = BX_;
1770 }
1771 else {
1772 state.BX = Teuchos::null;
1773 }
1774 state.rho = rho_;
1775 state.R = R_;
1776 state.T = Teuchos::rcp(new std::vector<MagnitudeType>(theta_));
1777 return state;
1778 }
1779
1780 template <class ScalarType, class MV, class OP>
1782 {
1783 return initialized_;
1784 }
1785
1786 template <class ScalarType, class MV, class OP>
1788 {
1789 std::vector<int> ret(nevLocal_,0);
1790 return ret;
1791 }
1792
1793
1794} // end Anasazi namespace
1795
1796#endif // ANASAZI_RTRBASE_HPP
Pure virtual base class which describes the basic interface to the iterative eigensolver.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Declaration of basic traits for the multivector type.
Virtual base class which defines basic traits for the operator type.
Class which provides internal utilities for the Anasazi solvers.
Types and exceptions used within Anasazi solvers and interfaces.
An exception class parent to all Anasazi exceptions.
This class defines the interface required by an eigensolver and status test class to compute solution...
The Eigensolver is a templated virtual base class that defines the basic interface that any eigensolv...
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...
This class is an abstract base class for Implicit Riemannian Trust-Region based eigensolvers....
virtual void iterate()=0
This method performs RTR iterations until the status test indicates the need to stop or an error occu...
void setAuxVecs(const Teuchos::Array< Teuchos::RCP< const MV > > &auxvecs)
Set the auxiliary vectors for the solver.
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace. For RTR, this always returns getBlockSiz...
Teuchos::RCP< const MV > getRitzVectors()
Get the Ritz vectors from the previous iteration.
Teuchos::Array< Teuchos::RCP< const MV > > getAuxVecs() const
Get the current auxiliary vectors.
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > getStatusTest() const
Get the current StatusTest used by the solver.
void resetNumIters()
Reset the iteration count.
void setBlockSize(int blockSize)
Set the blocksize to be used by the iterative solver in solving this eigenproblem.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the eigenvalue problem.
std::vector< Value< ScalarType > > getRitzValues()
Get the Ritz values from the previous iteration.
virtual void currentStatus(std::ostream &os)
This method requests that the solver print out its current status to screen.
int getNumIters() const
Get the current iteration count.
void initialize()
Initialize the solver with the initial vectors from the eigenproblem or random data.
virtual ~RTRBase()
RTRBase destructor
std::vector< int > getRitzIndex()
Get the index used for extracting Ritz vectors from getRitzVectors().
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRitzRes2Norms()
Get the 2-norms of the Ritz residuals.
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current eigenvectors and eigenvalues.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this eigenproblem.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRes2Norms()
Get the current residual 2-norms.
RTRBase(const Teuchos::RCP< Eigenproblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< SortManager< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > &sorter, const Teuchos::RCP< OutputManager< ScalarType > > &printer, const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &tester, const Teuchos::RCP< GenOrthoManager< ScalarType, MV, OP > > &ortho, Teuchos::ParameterList &params, const std::string &solverLabel, bool skinnySolver)
RTRBase constructor with eigenproblem, solver utilities, and parameter list of solver options.
bool isInitialized() const
Indicates whether the solver has been initialized or not.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getResNorms()
Get the current residual norms.
void setStatusTest(Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test)
Set a new StatusTest for the solver.
RTRState< ScalarType, MV > getState() const
Get the current state of the eigensolver.
RTRInitFailure is thrown when the RTR solver is unable to generate an initial iterate in the RTRBase:...
RTROrthoFailure is thrown when an orthogonalization attempt fails.
RTRRitzFailure is thrown when the RTR solver is unable to continue a call to RTRBase::iterate() due t...
Anasazi's templated, static class providing utilities for the solvers.
Anasazi's templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
Common interface of stopping criteria for Anasazi's solvers.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
Structure to contain pointers to RTR state variables.
Teuchos::RCP< const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > T
The current Ritz values.
Teuchos::RCP< const MV > BX
The image of the current eigenvectors under B, or Teuchos::null if B was not specified.
Teuchos::ScalarTraits< ScalarType >::magnitudeType rho
The current rho value. This is only valid if the debugging level of verbosity is enabled.
Teuchos::RCP< const MV > R
The current residual vectors.
Teuchos::RCP< const MV > X
The current eigenvectors.
Teuchos::RCP< const MV > AX
The image of the current eigenvectors under A, or Teuchos::null is we implement a skinny solver.