Anasazi Version of the Day
Loading...
Searching...
No Matches
AnasaziBlockDavidson.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_BLOCK_DAVIDSON_HPP
47#define ANASAZI_BLOCK_DAVIDSON_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
80namespace Anasazi {
81
83
84
89 template <class ScalarType, class MV>
95 int curDim;
100 Teuchos::RCP<const MV> V;
102 Teuchos::RCP<const MV> X;
104 Teuchos::RCP<const MV> KX;
106 Teuchos::RCP<const MV> MX;
108 Teuchos::RCP<const MV> R;
113 Teuchos::RCP<const MV> H;
115 Teuchos::RCP<const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > T;
121 Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > KK;
122 BlockDavidsonState() : curDim(0), V(Teuchos::null),
123 X(Teuchos::null), KX(Teuchos::null), MX(Teuchos::null),
124 R(Teuchos::null), H(Teuchos::null),
125 T(Teuchos::null), KK(Teuchos::null) {}
126 };
127
129
131
132
146 BlockDavidsonInitFailure(const std::string& what_arg) : AnasaziError(what_arg)
147 {}};
148
157 BlockDavidsonOrthoFailure(const std::string& what_arg) : AnasaziError(what_arg)
158 {}};
159
161
162
163 template <class ScalarType, class MV, class OP>
164 class BlockDavidson : public Eigensolver<ScalarType,MV,OP> {
165 public:
167
168
176 BlockDavidson( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
177 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
178 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
179 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
180 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
181 Teuchos::ParameterList &params
182 );
183
185 virtual ~BlockDavidson();
187
188
190
191
215 void iterate();
216
250 void initialize(BlockDavidsonState<ScalarType,MV>& newstate);
251
255 void initialize();
256
272 bool isInitialized() const;
273
286 BlockDavidsonState<ScalarType,MV> getState() const;
287
289
290
292
293
295 int getNumIters() const;
296
298 void resetNumIters();
299
307 Teuchos::RCP<const MV> getRitzVectors();
308
314 std::vector<Value<ScalarType> > getRitzValues();
315
316
324 std::vector<int> getRitzIndex();
325
326
332 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getResNorms();
333
334
340 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRes2Norms();
341
342
350 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRitzRes2Norms();
351
357 int getCurSubspaceDim() const;
358
360 int getMaxSubspaceDim() const;
361
363
364
366
367
369 void setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test);
370
372 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > getStatusTest() const;
373
375 const Eigenproblem<ScalarType,MV,OP>& getProblem() const;
376
386 void setBlockSize(int blockSize);
387
389 int getBlockSize() const;
390
403 void setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs);
404
406 Teuchos::Array<Teuchos::RCP<const MV> > getAuxVecs() const;
407
409
411
412
422 void setSize(int blockSize, int numBlocks);
423
425
427
428
430 void currentStatus(std::ostream &os);
431
433
434 private:
435 //
436 // Convenience typedefs
437 //
438 typedef SolverUtils<ScalarType,MV,OP> Utils;
439 typedef MultiVecTraits<ScalarType,MV> MVT;
440 typedef OperatorTraits<ScalarType,MV,OP> OPT;
441 typedef Teuchos::ScalarTraits<ScalarType> SCT;
442 typedef typename SCT::magnitudeType MagnitudeType;
443 const MagnitudeType ONE;
444 const MagnitudeType ZERO;
445 const MagnitudeType NANVAL;
446 //
447 // Internal structs
448 //
449 struct CheckList {
450 bool checkV;
451 bool checkX, checkMX, checkKX;
452 bool checkH, checkMH, checkKH;
453 bool checkR, checkQ;
454 bool checkKK;
455 CheckList() : checkV(false),
456 checkX(false),checkMX(false),checkKX(false),
457 checkH(false),checkMH(false),checkKH(false),
458 checkR(false),checkQ(false),checkKK(false) {};
459 };
460 //
461 // Internal methods
462 //
463 std::string accuracyCheck(const CheckList &chk, const std::string &where) const;
464 //
465 // Classes inputed through constructor that define the eigenproblem to be solved.
466 //
467 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
468 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sm_;
469 const Teuchos::RCP<OutputManager<ScalarType> > om_;
470 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > tester_;
471 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > orthman_;
472 //
473 // Information obtained from the eigenproblem
474 //
475 Teuchos::RCP<const OP> Op_;
476 Teuchos::RCP<const OP> MOp_;
477 Teuchos::RCP<const OP> Prec_;
478 bool hasM_;
479 //
480 // Internal timers
481 //
482 Teuchos::RCP<Teuchos::Time> timerOp_, timerMOp_, timerPrec_,
483 timerSortEval_, timerDS_,
484 timerLocal_, timerCompRes_,
485 timerOrtho_, timerInit_;
486 //
487 // Counters
488 //
489 int count_ApplyOp_, count_ApplyM_, count_ApplyPrec_;
490
491 //
492 // Algorithmic parameters.
493 //
494 // blockSize_ is the solver block size; it controls the number of eigenvectors that
495 // we compute, the number of residual vectors that we compute, and therefore the number
496 // of vectors added to the basis on each iteration.
497 int blockSize_;
498 // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
499 int numBlocks_;
500
501 //
502 // Current solver state
503 //
504 // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
505 // is capable of running; _initialize is controlled by the initialize() member method
506 // For the implications of the state of initialized_, please see documentation for initialize()
507 bool initialized_;
508 //
509 // curDim_ reflects how much of the current basis is valid
510 // NOTE: 0 <= curDim_ <= blockSize_*numBlocks_
511 // this also tells us how many of the values in theta_ are valid Ritz values
512 int curDim_;
513 //
514 // State Multivecs
515 // H_,KH_,MH_ will not own any storage
516 // H_ will occasionally point at the current block of vectors in the basis V_
517 // MH_,KH_ will occasionally point at MX_,KX_ when they are used as temporary storage
518 Teuchos::RCP<MV> X_, KX_, MX_, R_,
519 H_, KH_, MH_,
520 V_;
521 //
522 // Projected matrices
523 //
524 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > KK_;
525 //
526 // auxiliary vectors
527 Teuchos::Array<Teuchos::RCP<const MV> > auxVecs_;
528 int numAuxVecs_;
529 //
530 // Number of iterations that have been performed.
531 int iter_;
532 //
533 // Current eigenvalues, residual norms
534 std::vector<MagnitudeType> theta_, Rnorms_, R2norms_;
535 //
536 // are the residual norms current with the residual?
537 bool Rnorms_current_, R2norms_current_;
538
539 };
540
543 //
544 // Implementations
545 //
548
549
551 // Constructor
552 template <class ScalarType, class MV, class OP>
554 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
555 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
556 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
557 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
558 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
559 Teuchos::ParameterList &params
560 ) :
561 ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
562 ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
563 NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan()),
564 // problem, tools
565 problem_(problem),
566 sm_(sorter),
567 om_(printer),
568 tester_(tester),
569 orthman_(ortho),
570 // timers, counters
571#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
572 timerOp_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Operation Op*x")),
573 timerMOp_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Operation M*x")),
574 timerPrec_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Operation Prec*x")),
575 timerSortEval_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Sorting eigenvalues")),
576 timerDS_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Direct solve")),
577 timerLocal_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Local update")),
578 timerCompRes_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Computing residuals")),
579 timerOrtho_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Orthogonalization")),
580 timerInit_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidson::Initialization")),
581#endif
582 count_ApplyOp_(0),
583 count_ApplyM_(0),
584 count_ApplyPrec_(0),
585 // internal data
586 blockSize_(0),
587 numBlocks_(0),
588 initialized_(false),
589 curDim_(0),
590 auxVecs_( Teuchos::Array<Teuchos::RCP<const MV> >(0) ),
591 numAuxVecs_(0),
592 iter_(0),
593 Rnorms_current_(false),
594 R2norms_current_(false)
595 {
596 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument,
597 "Anasazi::BlockDavidson::constructor: user passed null problem pointer.");
598 TEUCHOS_TEST_FOR_EXCEPTION(sm_ == Teuchos::null,std::invalid_argument,
599 "Anasazi::BlockDavidson::constructor: user passed null sort manager pointer.");
600 TEUCHOS_TEST_FOR_EXCEPTION(om_ == Teuchos::null,std::invalid_argument,
601 "Anasazi::BlockDavidson::constructor: user passed null output manager pointer.");
602 TEUCHOS_TEST_FOR_EXCEPTION(tester_ == Teuchos::null,std::invalid_argument,
603 "Anasazi::BlockDavidson::constructor: user passed null status test pointer.");
604 TEUCHOS_TEST_FOR_EXCEPTION(orthman_ == Teuchos::null,std::invalid_argument,
605 "Anasazi::BlockDavidson::constructor: user passed null orthogonalization manager pointer.");
606 TEUCHOS_TEST_FOR_EXCEPTION(problem_->isProblemSet() == false, std::invalid_argument,
607 "Anasazi::BlockDavidson::constructor: problem is not set.");
608 TEUCHOS_TEST_FOR_EXCEPTION(problem_->isHermitian() == false, std::invalid_argument,
609 "Anasazi::BlockDavidson::constructor: problem is not hermitian.");
610
611 // get the problem operators
612 Op_ = problem_->getOperator();
613 TEUCHOS_TEST_FOR_EXCEPTION(Op_ == Teuchos::null, std::invalid_argument,
614 "Anasazi::BlockDavidson::constructor: problem provides no operator.");
615 MOp_ = problem_->getM();
616 Prec_ = problem_->getPrec();
617 hasM_ = (MOp_ != Teuchos::null);
618
619 // set the block size and allocate data
620 int bs = params.get("Block Size", problem_->getNEV());
621 int nb = params.get("Num Blocks", 2);
622 setSize(bs,nb);
623 }
624
625
627 // Destructor
628 template <class ScalarType, class MV, class OP>
630
631
633 // Set the block size
634 // This simply calls setSize(), modifying the block size while retaining the number of blocks.
635 template <class ScalarType, class MV, class OP>
637 {
638 setSize(blockSize,numBlocks_);
639 }
640
641
643 // Return the current auxiliary vectors
644 template <class ScalarType, class MV, class OP>
645 Teuchos::Array<Teuchos::RCP<const MV> > BlockDavidson<ScalarType,MV,OP>::getAuxVecs() const {
646 return auxVecs_;
647 }
648
649
650
652 // return the current block size
653 template <class ScalarType, class MV, class OP>
655 return(blockSize_);
656 }
657
658
660 // return eigenproblem
661 template <class ScalarType, class MV, class OP>
662 const Eigenproblem<ScalarType,MV,OP>& BlockDavidson<ScalarType,MV,OP>::getProblem() const {
663 return(*problem_);
664 }
665
666
668 // return max subspace dim
669 template <class ScalarType, class MV, class OP>
671 return blockSize_*numBlocks_;
672 }
673
675 // return current subspace dim
676 template <class ScalarType, class MV, class OP>
678 if (!initialized_) return 0;
679 return curDim_;
680 }
681
682
684 // return ritz residual 2-norms
685 template <class ScalarType, class MV, class OP>
686 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
688 return this->getRes2Norms();
689 }
690
691
693 // return ritz index
694 template <class ScalarType, class MV, class OP>
696 std::vector<int> ret(curDim_,0);
697 return ret;
698 }
699
700
702 // return ritz values
703 template <class ScalarType, class MV, class OP>
704 std::vector<Value<ScalarType> > BlockDavidson<ScalarType,MV,OP>::getRitzValues() {
705 std::vector<Value<ScalarType> > ret(curDim_);
706 for (int i=0; i<curDim_; ++i) {
707 ret[i].realpart = theta_[i];
708 ret[i].imagpart = ZERO;
709 }
710 return ret;
711 }
712
713
715 // return pointer to ritz vectors
716 template <class ScalarType, class MV, class OP>
718 return X_;
719 }
720
721
723 // reset number of iterations
724 template <class ScalarType, class MV, class OP>
728
729
731 // return number of iterations
732 template <class ScalarType, class MV, class OP>
734 return(iter_);
735 }
736
737
739 // return state pointers
740 template <class ScalarType, class MV, class OP>
741 BlockDavidsonState<ScalarType,MV> BlockDavidson<ScalarType,MV,OP>::getState() const {
742 BlockDavidsonState<ScalarType,MV> state;
743 state.curDim = curDim_;
744 state.V = V_;
745 state.X = X_;
746 state.KX = KX_;
747 if (hasM_) {
748 state.MX = MX_;
749 }
750 else {
751 state.MX = Teuchos::null;
752 }
753 state.R = R_;
754 state.H = H_;
755 state.KK = KK_;
756 if (curDim_ > 0) {
757 state.T = Teuchos::rcp(new std::vector<MagnitudeType>(theta_.begin(),theta_.begin()+curDim_) );
758 } else {
759 state.T = Teuchos::rcp(new std::vector<MagnitudeType>(0));
760 }
761 return state;
762 }
763
764
766 // Return initialized state
767 template <class ScalarType, class MV, class OP>
768 bool BlockDavidson<ScalarType,MV,OP>::isInitialized() const { return initialized_; }
769
770
772 // Set the block size and make necessary adjustments.
773 template <class ScalarType, class MV, class OP>
774 void BlockDavidson<ScalarType,MV,OP>::setSize (int blockSize, int numBlocks)
775 {
776 // time spent here counts towards timerInit_
777#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
778 Teuchos::TimeMonitor initimer( *timerInit_ );
779#endif
780
781 // This routine only allocates space; it doesn't not perform any computation
782 // any change in size will invalidate the state of the solver.
783
784 TEUCHOS_TEST_FOR_EXCEPTION(blockSize < 1, std::invalid_argument, "Anasazi::BlockDavidson::setSize(blocksize,numblocks): blocksize must be strictly positive.");
785 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks < 2, std::invalid_argument, "Anasazi::BlockDavidson::setSize(blocksize,numblocks): numblocks must be greater than one.");
786 if (blockSize == blockSize_ && numBlocks == numBlocks_) {
787 // do nothing
788 return;
789 }
790
791 blockSize_ = blockSize;
792 numBlocks_ = numBlocks;
793
794 Teuchos::RCP<const MV> tmp;
795 // grab some Multivector to Clone
796 // in practice, getInitVec() should always provide this, but it is possible to use a
797 // Eigenproblem with nothing in getInitVec() by manually initializing with initialize();
798 // in case of that strange scenario, we will try to Clone from X_ first, then resort to getInitVec()
799 if (X_ != Teuchos::null) { // this is equivalent to blockSize_ > 0
800 tmp = X_;
801 }
802 else {
803 tmp = problem_->getInitVec();
804 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
805 "Anasazi::BlockDavidson::setSize(): eigenproblem did not specify initial vectors to clone from.");
806 }
807
808 TEUCHOS_TEST_FOR_EXCEPTION(numAuxVecs_+blockSize*static_cast<ptrdiff_t>(numBlocks) > MVT::GetGlobalLength(*tmp),std::invalid_argument,
809 "Anasazi::BlockDavidson::setSize(): max subspace dimension and auxilliary subspace too large.");
810
811
813 // blockSize dependent
814 //
815 // grow/allocate vectors
816 Rnorms_.resize(blockSize_,NANVAL);
817 R2norms_.resize(blockSize_,NANVAL);
818 //
819 // clone multivectors off of tmp
820 //
821 // free current allocation first, to make room for new allocation
822 X_ = Teuchos::null;
823 KX_ = Teuchos::null;
824 MX_ = Teuchos::null;
825 R_ = Teuchos::null;
826 V_ = Teuchos::null;
827
828 om_->print(Debug," >> Allocating X_\n");
829 X_ = MVT::Clone(*tmp,blockSize_);
830 om_->print(Debug," >> Allocating KX_\n");
831 KX_ = MVT::Clone(*tmp,blockSize_);
832 if (hasM_) {
833 om_->print(Debug," >> Allocating MX_\n");
834 MX_ = MVT::Clone(*tmp,blockSize_);
835 }
836 else {
837 MX_ = X_;
838 }
839 om_->print(Debug," >> Allocating R_\n");
840 R_ = MVT::Clone(*tmp,blockSize_);
841
843 // blockSize*numBlocks dependent
844 //
845 int newsd = blockSize_*numBlocks_;
846 theta_.resize(blockSize_*numBlocks_,NANVAL);
847 om_->print(Debug," >> Allocating V_\n");
848 V_ = MVT::Clone(*tmp,newsd);
849 KK_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(newsd,newsd) );
850
851 om_->print(Debug," >> done allocating.\n");
852
853 initialized_ = false;
854 curDim_ = 0;
855 }
856
857
859 // Set the auxiliary vectors
860 template <class ScalarType, class MV, class OP>
861 void BlockDavidson<ScalarType,MV,OP>::setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs) {
862 typedef typename Teuchos::Array<Teuchos::RCP<const MV> >::iterator tarcpmv;
863
864 // set new auxiliary vectors
865 auxVecs_ = auxvecs;
866 numAuxVecs_ = 0;
867 for (tarcpmv i=auxVecs_.begin(); i != auxVecs_.end(); ++i) {
868 numAuxVecs_ += MVT::GetNumberVecs(**i);
869 }
870
871 // If the solver has been initialized, V is not necessarily orthogonal to new auxiliary vectors
872 if (numAuxVecs_ > 0 && initialized_) {
873 initialized_ = false;
874 }
875
876 if (om_->isVerbosity( Debug ) ) {
877 CheckList chk;
878 chk.checkQ = true;
879 om_->print( Debug, accuracyCheck(chk, ": in setAuxVecs()") );
880 }
881 }
882
883
885 /* Initialize the state of the solver
886 *
887 * POST-CONDITIONS:
888 *
889 * V_ is orthonormal, orthogonal to auxVecs_, for first curDim_ vectors
890 * theta_ contains Ritz w.r.t. V_(1:curDim_)
891 * X is Ritz vectors w.r.t. V_(1:curDim_)
892 * KX = Op*X
893 * MX = M*X if hasM_
894 * R = KX - MX*diag(theta_)
895 *
896 */
897 template <class ScalarType, class MV, class OP>
898 void BlockDavidson<ScalarType,MV,OP>::initialize(BlockDavidsonState<ScalarType,MV>& newstate)
899 {
900 // NOTE: memory has been allocated by setBlockSize(). Use setBlock below; do not Clone
901 // NOTE: Overall time spent in this routine is counted to timerInit_; portions will also be counted towards other primitives
902
903#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
904 Teuchos::TimeMonitor inittimer( *timerInit_ );
905#endif
906
907 std::vector<int> bsind(blockSize_);
908 for (int i=0; i<blockSize_; ++i) bsind[i] = i;
909
910 Teuchos::BLAS<int,ScalarType> blas;
911
912 // in BlockDavidson, V is primary
913 // the order of dependence follows like so.
914 // --init-> V,KK
915 // --ritz analysis-> theta,X
916 // --op apply-> KX,MX
917 // --compute-> R
918 //
919 // if the user specifies all data for a level, we will accept it.
920 // otherwise, we will generate the whole level, and all subsequent levels.
921 //
922 // the data members are ordered based on dependence, and the levels are
923 // partitioned according to the amount of work required to produce the
924 // items in a level.
925 //
926 // inconsistent multivectors widths and lengths will not be tolerated, and
927 // will be treated with exceptions.
928 //
929 // for multivector pointers in newstate which point directly (as opposed to indirectly, via a view) to
930 // multivectors in the solver, the copy will not be affected.
931
932 // set up V and KK: get them from newstate if user specified them
933 // otherwise, set them manually
934 Teuchos::RCP<MV> lclV;
935 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclKK;
936
937 if (newstate.V != Teuchos::null && newstate.KK != Teuchos::null) {
938 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.V) != MVT::GetGlobalLength(*V_), std::invalid_argument,
939 "Anasazi::BlockDavidson::initialize(newstate): Vector length of V not correct." );
940 TEUCHOS_TEST_FOR_EXCEPTION( newstate.curDim < blockSize_, std::invalid_argument,
941 "Anasazi::BlockDavidson::initialize(newstate): Rank of new state must be at least blockSize().");
942 TEUCHOS_TEST_FOR_EXCEPTION( newstate.curDim > blockSize_*numBlocks_, std::invalid_argument,
943 "Anasazi::BlockDavidson::initialize(newstate): Rank of new state must be less than getMaxSubspaceDim().");
944 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < newstate.curDim, std::invalid_argument,
945 "Anasazi::BlockDavidson::initialize(newstate): Multivector for basis in new state must be as large as specified state rank.");
946
947 curDim_ = newstate.curDim;
948 // pick an integral amount
949 curDim_ = (int)(curDim_ / blockSize_)*blockSize_;
950
951 TEUCHOS_TEST_FOR_EXCEPTION( curDim_ != newstate.curDim, std::invalid_argument,
952 "Anasazi::BlockDavidson::initialize(newstate): Rank of new state must be a multiple of getBlockSize().");
953
954 // check size of KK
955 TEUCHOS_TEST_FOR_EXCEPTION( newstate.KK->numRows() < curDim_ || newstate.KK->numCols() < curDim_, std::invalid_argument,
956 "Anasazi::BlockDavidson::initialize(newstate): Projected matrix in new state must be as large as specified state rank.");
957
958 // put data in V
959 std::vector<int> nevind(curDim_);
960 for (int i=0; i<curDim_; ++i) nevind[i] = i;
961 if (newstate.V != V_) {
962 if (curDim_ < MVT::GetNumberVecs(*newstate.V)) {
963 newstate.V = MVT::CloneView(*newstate.V,nevind);
964 }
965 MVT::SetBlock(*newstate.V,nevind,*V_);
966 }
967 lclV = MVT::CloneViewNonConst(*V_,nevind);
968
969 // put data into KK_
970 lclKK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
971 if (newstate.KK != KK_) {
972 if (newstate.KK->numRows() > curDim_ || newstate.KK->numCols() > curDim_) {
973 newstate.KK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*newstate.KK,curDim_,curDim_) );
974 }
975 lclKK->assign(*newstate.KK);
976 }
977 //
978 // make lclKK Hermitian in memory (copy the upper half to the lower half)
979 for (int j=0; j<curDim_-1; ++j) {
980 for (int i=j+1; i<curDim_; ++i) {
981 (*lclKK)(i,j) = SCT::conjugate((*lclKK)(j,i));
982 }
983 }
984 }
985 else {
986 // user did not specify one of V or KK
987 // get vectors from problem or generate something, projectAndNormalize
988 Teuchos::RCP<const MV> ivec = problem_->getInitVec();
989 TEUCHOS_TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::invalid_argument,
990 "Anasazi::BlockDavdison::initialize(newstate): Eigenproblem did not specify initial vectors to clone from.");
991 // clear newstate so we won't use any data from it below
992 newstate.X = Teuchos::null;
993 newstate.MX = Teuchos::null;
994 newstate.KX = Teuchos::null;
995 newstate.R = Teuchos::null;
996 newstate.H = Teuchos::null;
997 newstate.T = Teuchos::null;
998 newstate.KK = Teuchos::null;
999 newstate.V = Teuchos::null;
1000 newstate.curDim = 0;
1001
1002 curDim_ = MVT::GetNumberVecs(*ivec);
1003 // pick the largest multiple of blockSize_
1004 curDim_ = (int)(curDim_ / blockSize_)*blockSize_;
1005 if (curDim_ > blockSize_*numBlocks_) {
1006 // user specified too many vectors... truncate
1007 // this produces a full subspace, but that is okay
1008 curDim_ = blockSize_*numBlocks_;
1009 }
1010 bool userand = false;
1011 if (curDim_ == 0) {
1012 // we need at least blockSize_ vectors
1013 // use a random multivec: ignore everything from InitVec
1014 userand = true;
1015 curDim_ = blockSize_;
1016 }
1017
1018 // get pointers into V,KV,MV
1019 // tmpVecs will be used below for M*V and K*V (not simultaneously)
1020 // lclV has curDim vectors
1021 // if there is space for lclV and tmpVecs in V_, point tmpVecs into V_
1022 // otherwise, we must allocate space for these products
1023 //
1024 // get pointer to first curDim vector in V_
1025 std::vector<int> dimind(curDim_);
1026 for (int i=0; i<curDim_; ++i) dimind[i] = i;
1027 lclV = MVT::CloneViewNonConst(*V_,dimind);
1028 if (userand) {
1029 // generate random vector data
1030 MVT::MvRandom(*lclV);
1031 }
1032 else {
1033 if (MVT::GetNumberVecs(*ivec) > curDim_) {
1034 ivec = MVT::CloneView(*ivec,dimind);
1035 }
1036 // assign ivec to first part of lclV
1037 MVT::SetBlock(*ivec,dimind,*lclV);
1038 }
1039 Teuchos::RCP<MV> tmpVecs;
1040 if (curDim_*2 <= blockSize_*numBlocks_) {
1041 // partition V_ = [lclV tmpVecs _leftover_]
1042 std::vector<int> block2(curDim_);
1043 for (int i=0; i<curDim_; ++i) block2[i] = curDim_+i;
1044 tmpVecs = MVT::CloneViewNonConst(*V_,block2);
1045 }
1046 else {
1047 // allocate space for tmpVecs
1048 tmpVecs = MVT::Clone(*V_,curDim_);
1049 }
1050
1051 // compute M*lclV if hasM_
1052 if (hasM_) {
1053#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1054 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1055#endif
1056 OPT::Apply(*MOp_,*lclV,*tmpVecs);
1057 count_ApplyM_ += curDim_;
1058 }
1059
1060 // remove auxVecs from lclV and normalize it
1061 if (auxVecs_.size() > 0) {
1062#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1063 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1064#endif
1065
1066 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
1067 int rank = orthman_->projectAndNormalizeMat(*lclV,auxVecs_,dummyC,Teuchos::null,tmpVecs);
1068 TEUCHOS_TEST_FOR_EXCEPTION(rank != curDim_,BlockDavidsonInitFailure,
1069 "Anasazi::BlockDavidson::initialize(): Couldn't generate initial basis of full rank.");
1070 }
1071 else {
1072#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1073 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1074#endif
1075
1076 int rank = orthman_->normalizeMat(*lclV,Teuchos::null,tmpVecs);
1077 TEUCHOS_TEST_FOR_EXCEPTION(rank != curDim_,BlockDavidsonInitFailure,
1078 "Anasazi::BlockDavidson::initialize(): Couldn't generate initial basis of full rank.");
1079 }
1080
1081 // compute K*lclV: we are re-using tmpVecs to store the result
1082 {
1083#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1084 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1085#endif
1086 OPT::Apply(*Op_,*lclV,*tmpVecs);
1087 count_ApplyOp_ += curDim_;
1088 }
1089
1090 // generate KK
1091 lclKK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
1092 MVT::MvTransMv(ONE,*lclV,*tmpVecs,*lclKK);
1093
1094 // clear tmpVecs
1095 tmpVecs = Teuchos::null;
1096 }
1097
1098 // X,theta require Ritz analysis; if we have to generate one of these, we might as well generate both
1099 if (newstate.X != Teuchos::null && newstate.T != Teuchos::null) {
1100 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.X) != blockSize_ || MVT::GetGlobalLength(*newstate.X) != MVT::GetGlobalLength(*X_),
1101 std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): Size of X must be consistent with block size and length of V.");
1102 TEUCHOS_TEST_FOR_EXCEPTION((signed int)(newstate.T->size()) != curDim_,
1103 std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): Size of T must be consistent with dimension of V.");
1104
1105 if (newstate.X != X_) {
1106 MVT::SetBlock(*newstate.X,bsind,*X_);
1107 }
1108
1109 std::copy(newstate.T->begin(),newstate.T->end(),theta_.begin());
1110 }
1111 else {
1112 // compute ritz vecs/vals
1113 Teuchos::SerialDenseMatrix<int,ScalarType> S(curDim_,curDim_);
1114 {
1115#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1116 Teuchos::TimeMonitor lcltimer( *timerDS_ );
1117#endif
1118 int rank = curDim_;
1119 Utils::directSolver(curDim_, *lclKK, Teuchos::null, S, theta_, rank, 10);
1120 // we want all ritz values back
1121 TEUCHOS_TEST_FOR_EXCEPTION(rank != curDim_,BlockDavidsonInitFailure,
1122 "Anasazi::BlockDavidson::initialize(newstate): Not enough Ritz vectors to initialize algorithm.");
1123 }
1124 // sort ritz pairs
1125 {
1126#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1127 Teuchos::TimeMonitor lcltimer( *timerSortEval_ );
1128#endif
1129
1130 std::vector<int> order(curDim_);
1131 //
1132 // sort the first curDim_ values in theta_
1133 sm_->sort(theta_, Teuchos::rcpFromRef(order), curDim_); // don't catch exception
1134 //
1135 // apply the same ordering to the primitive ritz vectors
1136 Utils::permuteVectors(order,S);
1137 }
1138
1139 // compute eigenvectors
1140 Teuchos::SerialDenseMatrix<int,ScalarType> S1(Teuchos::View,S,curDim_,blockSize_);
1141 {
1142#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1143 Teuchos::TimeMonitor lcltimer( *timerLocal_ );
1144#endif
1145
1146 // X <- lclV*S
1147 MVT::MvTimesMatAddMv( ONE, *lclV, S1, ZERO, *X_ );
1148 }
1149 // we generated theta,X so we don't want to use the user's KX,MX
1150 newstate.KX = Teuchos::null;
1151 newstate.MX = Teuchos::null;
1152 }
1153
1154 // done with local pointers
1155 lclV = Teuchos::null;
1156 lclKK = Teuchos::null;
1157
1158 // set up KX
1159 if ( newstate.KX != Teuchos::null ) {
1160 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.KX) != blockSize_,
1161 std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): vector length of newstate.KX not correct." );
1162 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetGlobalLength(*newstate.KX) != MVT::GetGlobalLength(*X_),
1163 std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): newstate.KX must have at least block size vectors." );
1164 if (newstate.KX != KX_) {
1165 MVT::SetBlock(*newstate.KX,bsind,*KX_);
1166 }
1167 }
1168 else {
1169 // generate KX
1170 {
1171#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1172 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1173#endif
1174 OPT::Apply(*Op_,*X_,*KX_);
1175 count_ApplyOp_ += blockSize_;
1176 }
1177 // we generated KX; we will generate R as well
1178 newstate.R = Teuchos::null;
1179 }
1180
1181 // set up MX
1182 if (hasM_) {
1183 if ( newstate.MX != Teuchos::null ) {
1184 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.MX) != blockSize_,
1185 std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): vector length of newstate.MX not correct." );
1186 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetGlobalLength(*newstate.MX) != MVT::GetGlobalLength(*X_),
1187 std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): newstate.MX must have at least block size vectors." );
1188 if (newstate.MX != MX_) {
1189 MVT::SetBlock(*newstate.MX,bsind,*MX_);
1190 }
1191 }
1192 else {
1193 // generate MX
1194 {
1195#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1196 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1197#endif
1198 OPT::Apply(*MOp_,*X_,*MX_);
1199 count_ApplyOp_ += blockSize_;
1200 }
1201 // we generated MX; we will generate R as well
1202 newstate.R = Teuchos::null;
1203 }
1204 }
1205 else {
1206 // the assignment MX_==X_ would be redundant; take advantage of this opportunity to debug a little
1207 TEUCHOS_TEST_FOR_EXCEPTION(MX_ != X_, std::logic_error, "Anasazi::BlockDavidson::initialize(): solver invariant not satisfied (MX==X).");
1208 }
1209
1210 // set up R
1211 if (newstate.R != Teuchos::null) {
1212 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.R) != blockSize_,
1213 std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): vector length of newstate.R not correct." );
1214 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*X_),
1215 std::invalid_argument, "Anasazi::BlockDavidson::initialize(newstate): newstate.R must have at least block size vectors." );
1216 if (newstate.R != R_) {
1217 MVT::SetBlock(*newstate.R,bsind,*R_);
1218 }
1219 }
1220 else {
1221#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1222 Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
1223#endif
1224
1225 // form R <- KX - MX*T
1226 MVT::MvAddMv(ZERO,*KX_,ONE,*KX_,*R_);
1227 Teuchos::SerialDenseMatrix<int,ScalarType> T(blockSize_,blockSize_);
1228 T.putScalar(ZERO);
1229 for (int i=0; i<blockSize_; ++i) T(i,i) = theta_[i];
1230 MVT::MvTimesMatAddMv(-ONE,*MX_,T,ONE,*R_);
1231
1232 }
1233
1234 // R has been updated; mark the norms as out-of-date
1235 Rnorms_current_ = false;
1236 R2norms_current_ = false;
1237
1238 // finally, we are initialized
1239 initialized_ = true;
1240
1241 if (om_->isVerbosity( Debug ) ) {
1242 // Check almost everything here
1243 CheckList chk;
1244 chk.checkV = true;
1245 chk.checkX = true;
1246 chk.checkKX = true;
1247 chk.checkMX = true;
1248 chk.checkR = true;
1249 chk.checkQ = true;
1250 chk.checkKK = true;
1251 om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
1252 }
1253
1254 // Print information on current status
1255 if (om_->isVerbosity(Debug)) {
1256 currentStatus( om_->stream(Debug) );
1257 }
1258 else if (om_->isVerbosity(IterationDetails)) {
1259 currentStatus( om_->stream(IterationDetails) );
1260 }
1261 }
1262
1263
1265 // initialize the solver with default state
1266 template <class ScalarType, class MV, class OP>
1268 {
1269 BlockDavidsonState<ScalarType,MV> empty;
1270 initialize(empty);
1271 }
1272
1273
1275 // Perform BlockDavidson iterations until the StatusTest tells us to stop.
1276 template <class ScalarType, class MV, class OP>
1278 {
1279 //
1280 // Initialize solver state
1281 if (initialized_ == false) {
1282 initialize();
1283 }
1284
1285 // as a data member, this would be redundant and require synchronization with
1286 // blockSize_ and numBlocks_; we'll use a constant here.
1287 const int searchDim = blockSize_*numBlocks_;
1288
1289 Teuchos::BLAS<int,ScalarType> blas;
1290
1291 //
1292 // The projected matrices are part of the state, but the eigenvectors are defined locally.
1293 // S = Local eigenvectors (size: searchDim * searchDim
1294 Teuchos::SerialDenseMatrix<int,ScalarType> S( searchDim, searchDim );
1295
1296
1298 // iterate until the status test tells us to stop.
1299 // also break if our basis is full
1300 while (tester_->checkStatus(this) != Passed && curDim_ < searchDim) {
1301
1302 // Print information on current iteration
1303 if (om_->isVerbosity(Debug)) {
1304 currentStatus( om_->stream(Debug) );
1305 }
1306 else if (om_->isVerbosity(IterationDetails)) {
1307 currentStatus( om_->stream(IterationDetails) );
1308 }
1309
1310 ++iter_;
1311
1312 // get the current part of the basis
1313 std::vector<int> curind(blockSize_);
1314 for (int i=0; i<blockSize_; ++i) curind[i] = curDim_ + i;
1315 H_ = MVT::CloneViewNonConst(*V_,curind);
1316
1317 // Apply the preconditioner on the residuals: H <- Prec*R
1318 // H = Prec*R
1319 if (Prec_ != Teuchos::null) {
1320#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1321 Teuchos::TimeMonitor lcltimer( *timerPrec_ );
1322#endif
1323 OPT::Apply( *Prec_, *R_, *H_ ); // don't catch the exception
1324 count_ApplyPrec_ += blockSize_;
1325 }
1326 else {
1327 std::vector<int> bsind(blockSize_);
1328 for (int i=0; i<blockSize_; ++i) { bsind[i] = i; }
1329 MVT::SetBlock(*R_,bsind,*H_);
1330 }
1331
1332 // Apply the mass matrix on H
1333 if (hasM_) {
1334 // use memory at MX_ for temporary storage
1335 MH_ = MX_;
1336#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1337 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1338#endif
1339 OPT::Apply( *MOp_, *H_, *MH_); // don't catch the exception
1340 count_ApplyM_ += blockSize_;
1341 }
1342 else {
1343 MH_ = H_;
1344 }
1345
1346 // Get a view of the previous vectors
1347 // this is used for orthogonalization and for computing V^H K H
1348 std::vector<int> prevind(curDim_);
1349 for (int i=0; i<curDim_; ++i) prevind[i] = i;
1350 Teuchos::RCP<const MV> Vprev = MVT::CloneView(*V_,prevind);
1351
1352 // Orthogonalize H against the previous vectors and the auxiliary vectors, and normalize
1353 {
1354#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1355 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1356#endif
1357
1358 Teuchos::Array<Teuchos::RCP<const MV> > against = auxVecs_;
1359 against.push_back(Vprev);
1360 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
1361 int rank = orthman_->projectAndNormalizeMat(*H_,against,
1362 dummyC,
1363 Teuchos::null,MH_);
1364 TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_,BlockDavidsonOrthoFailure,
1365 "Anasazi::BlockDavidson::iterate(): unable to compute orthonormal basis for H.");
1366 }
1367
1368 // Apply the stiffness matrix to H
1369 {
1370 // use memory at KX_ for temporary storage
1371 KH_ = KX_;
1372#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1373 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1374#endif
1375 OPT::Apply( *Op_, *H_, *KH_); // don't catch the exception
1376 count_ApplyOp_ += blockSize_;
1377 }
1378
1379 if (om_->isVerbosity( Debug ) ) {
1380 CheckList chk;
1381 chk.checkH = true;
1382 chk.checkMH = true;
1383 chk.checkKH = true;
1384 om_->print( Debug, accuracyCheck(chk, ": after ortho H") );
1385 }
1386 else if (om_->isVerbosity( OrthoDetails ) ) {
1387 CheckList chk;
1388 chk.checkH = true;
1389 chk.checkMH = true;
1390 chk.checkKH = true;
1391 om_->print( OrthoDetails, accuracyCheck(chk,": after ortho H") );
1392 }
1393
1394 // compute next part of the projected matrices
1395 // this this in two parts
1396 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > nextKK;
1397 // Vprev*K*H
1398 nextKK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,blockSize_,0,curDim_) );
1399 MVT::MvTransMv(ONE,*Vprev,*KH_,*nextKK);
1400 // H*K*H
1401 nextKK = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,blockSize_,blockSize_,curDim_,curDim_) );
1402 MVT::MvTransMv(ONE,*H_,*KH_,*nextKK);
1403 //
1404 // make sure that KK_ is Hermitian in memory
1405 nextKK = Teuchos::null;
1406 for (int i=curDim_; i<curDim_+blockSize_; ++i) {
1407 for (int j=0; j<i; ++j) {
1408 (*KK_)(i,j) = SCT::conjugate((*KK_)(j,i));
1409 }
1410 }
1411
1412 // V has been extended, and KK has been extended. Update basis dim and release all pointers.
1413 curDim_ += blockSize_;
1414 H_ = KH_ = MH_ = Teuchos::null;
1415 Vprev = Teuchos::null;
1416
1417 if (om_->isVerbosity( Debug ) ) {
1418 CheckList chk;
1419 chk.checkKK = true;
1420 om_->print( Debug, accuracyCheck(chk, ": after expanding KK") );
1421 }
1422
1423 // Get pointer to complete basis
1424 curind.resize(curDim_);
1425 for (int i=0; i<curDim_; ++i) curind[i] = i;
1426 Teuchos::RCP<const MV> curV = MVT::CloneView(*V_,curind);
1427
1428 // Perform spectral decomposition
1429 {
1430#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1431 Teuchos::TimeMonitor lcltimer(*timerDS_);
1432#endif
1433 int nevlocal = curDim_;
1434 int info = Utils::directSolver(curDim_,*KK_,Teuchos::null,S,theta_,nevlocal,10);
1435 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,std::logic_error,"Anasazi::BlockDavidson::iterate(): direct solve returned error code.");
1436 // we did not ask directSolver to perform deflation, so nevLocal better be curDim_
1437 TEUCHOS_TEST_FOR_EXCEPTION(nevlocal != curDim_,std::logic_error,"Anasazi::BlockDavidson::iterate(): direct solve did not compute all eigenvectors."); // this should never happen
1438 }
1439
1440 // Sort ritz pairs
1441 {
1442#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1443 Teuchos::TimeMonitor lcltimer( *timerSortEval_ );
1444#endif
1445
1446 std::vector<int> order(curDim_);
1447 //
1448 // sort the first curDim_ values in theta_
1449 sm_->sort(theta_, Teuchos::rcp(&order,false), curDim_); // don't catch exception
1450 //
1451 // apply the same ordering to the primitive ritz vectors
1452 Teuchos::SerialDenseMatrix<int,ScalarType> curS(Teuchos::View,S,curDim_,curDim_);
1453 Utils::permuteVectors(order,curS);
1454 }
1455
1456 // Create a view matrix of the first blockSize_ vectors
1457 Teuchos::SerialDenseMatrix<int,ScalarType> S1( Teuchos::View, S, curDim_, blockSize_ );
1458
1459 // Compute the new Ritz vectors
1460 {
1461#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1462 Teuchos::TimeMonitor lcltimer( *timerLocal_ );
1463#endif
1464 MVT::MvTimesMatAddMv(ONE,*curV,S1,ZERO,*X_);
1465 }
1466
1467 // Apply the stiffness matrix for the Ritz vectors
1468 {
1469#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1470 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1471#endif
1472 OPT::Apply( *Op_, *X_, *KX_); // don't catch the exception
1473 count_ApplyOp_ += blockSize_;
1474 }
1475 // Apply the mass matrix for the Ritz vectors
1476 if (hasM_) {
1477#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1478 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1479#endif
1480 OPT::Apply(*MOp_,*X_,*MX_);
1481 count_ApplyM_ += blockSize_;
1482 }
1483 else {
1484 MX_ = X_;
1485 }
1486
1487 // Compute the residual
1488 // R = KX - MX*diag(theta)
1489 {
1490#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1491 Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
1492#endif
1493
1494 MVT::MvAddMv( ONE, *KX_, ZERO, *KX_, *R_ );
1495 Teuchos::SerialDenseMatrix<int,ScalarType> T( blockSize_, blockSize_ );
1496 for (int i = 0; i < blockSize_; ++i) {
1497 T(i,i) = theta_[i];
1498 }
1499 MVT::MvTimesMatAddMv( -ONE, *MX_, T, ONE, *R_ );
1500 }
1501
1502 // R has been updated; mark the norms as out-of-date
1503 Rnorms_current_ = false;
1504 R2norms_current_ = false;
1505
1506
1507 // When required, monitor some orthogonalities
1508 if (om_->isVerbosity( Debug ) ) {
1509 // Check almost everything here
1510 CheckList chk;
1511 chk.checkV = true;
1512 chk.checkX = true;
1513 chk.checkKX = true;
1514 chk.checkMX = true;
1515 chk.checkR = true;
1516 om_->print( Debug, accuracyCheck(chk, ": after local update") );
1517 }
1518 else if (om_->isVerbosity( OrthoDetails )) {
1519 CheckList chk;
1520 chk.checkX = true;
1521 chk.checkKX = true;
1522 chk.checkMX = true;
1523 chk.checkR = true;
1524 om_->print( OrthoDetails, accuracyCheck(chk, ": after local update") );
1525 }
1526 } // end while (statusTest == false)
1527
1528 } // end of iterate()
1529
1530
1531
1533 // compute/return residual M-norms
1534 template <class ScalarType, class MV, class OP>
1535 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
1537 if (Rnorms_current_ == false) {
1538 // Update the residual norms
1539 orthman_->norm(*R_,Rnorms_);
1540 Rnorms_current_ = true;
1541 }
1542 return Rnorms_;
1543 }
1544
1545
1547 // compute/return residual 2-norms
1548 template <class ScalarType, class MV, class OP>
1549 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
1551 if (R2norms_current_ == false) {
1552 // Update the residual 2-norms
1553 MVT::MvNorm(*R_,R2norms_);
1554 R2norms_current_ = true;
1555 }
1556 return R2norms_;
1557 }
1558
1560 // Set a new StatusTest for the solver.
1561 template <class ScalarType, class MV, class OP>
1562 void BlockDavidson<ScalarType,MV,OP>::setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test) {
1563 TEUCHOS_TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument,
1564 "Anasazi::BlockDavidson::setStatusTest() was passed a null StatusTest.");
1565 tester_ = test;
1566 }
1567
1569 // Get the current StatusTest used by the solver.
1570 template <class ScalarType, class MV, class OP>
1571 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > BlockDavidson<ScalarType,MV,OP>::getStatusTest() const {
1572 return tester_;
1573 }
1574
1575
1577 // Check accuracy, orthogonality, and other debugging stuff
1578 //
1579 // bools specify which tests we want to run (instead of running more than we actually care about)
1580 //
1581 // we don't bother checking the following because they are computed explicitly:
1582 // H == Prec*R
1583 // KH == K*H
1584 //
1585 //
1586 // checkV : V orthonormal
1587 // orthogonal to auxvecs
1588 // checkX : X orthonormal
1589 // orthogonal to auxvecs
1590 // checkMX: check MX == M*X
1591 // checkKX: check KX == K*X
1592 // checkH : H orthonormal
1593 // orthogonal to V and H and auxvecs
1594 // checkMH: check MH == M*H
1595 // checkR : check R orthogonal to X
1596 // checkQ : check that auxiliary vectors are actually orthonormal
1597 // checkKK: check that KK is symmetric in memory
1598 //
1599 // TODO:
1600 // add checkTheta
1601 //
1602 template <class ScalarType, class MV, class OP>
1603 std::string BlockDavidson<ScalarType,MV,OP>::accuracyCheck( const CheckList &chk, const std::string &where ) const
1604 {
1605 using std::endl;
1606
1607 std::stringstream os;
1608 os.precision(2);
1609 os.setf(std::ios::scientific, std::ios::floatfield);
1610
1611 os << " Debugging checks: iteration " << iter_ << where << endl;
1612
1613 // V and friends
1614 std::vector<int> lclind(curDim_);
1615 for (int i=0; i<curDim_; ++i) lclind[i] = i;
1616 Teuchos::RCP<const MV> lclV;
1617 if (initialized_) {
1618 lclV = MVT::CloneView(*V_,lclind);
1619 }
1620 if (chk.checkV && initialized_) {
1621 MagnitudeType err = orthman_->orthonormError(*lclV);
1622 os << " >> Error in V^H M V == I : " << err << endl;
1623 for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
1624 err = orthman_->orthogError(*lclV,*auxVecs_[i]);
1625 os << " >> Error in V^H M Q[" << i << "] == 0 : " << err << endl;
1626 }
1627 Teuchos::SerialDenseMatrix<int,ScalarType> curKK(curDim_,curDim_);
1628 Teuchos::RCP<MV> lclKV = MVT::Clone(*V_,curDim_);
1629 OPT::Apply(*Op_,*lclV,*lclKV);
1630 MVT::MvTransMv(ONE,*lclV,*lclKV,curKK);
1631 Teuchos::SerialDenseMatrix<int,ScalarType> subKK(Teuchos::View,*KK_,curDim_,curDim_);
1632 curKK -= subKK;
1633 // dup the lower tri part
1634 for (int j=0; j<curDim_; ++j) {
1635 for (int i=j+1; i<curDim_; ++i) {
1636 curKK(i,j) = curKK(j,i);
1637 }
1638 }
1639 os << " >> Error in V^H K V == KK : " << curKK.normFrobenius() << endl;
1640 }
1641
1642 // X and friends
1643 if (chk.checkX && initialized_) {
1644 MagnitudeType err = orthman_->orthonormError(*X_);
1645 os << " >> Error in X^H M X == I : " << err << endl;
1646 for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
1647 err = orthman_->orthogError(*X_,*auxVecs_[i]);
1648 os << " >> Error in X^H M Q[" << i << "] == 0 : " << err << endl;
1649 }
1650 }
1651 if (chk.checkMX && hasM_ && initialized_) {
1652 MagnitudeType err = Utils::errorEquality(*X_, *MX_, MOp_);
1653 os << " >> Error in MX == M*X : " << err << endl;
1654 }
1655 if (chk.checkKX && initialized_) {
1656 MagnitudeType err = Utils::errorEquality(*X_, *KX_, Op_);
1657 os << " >> Error in KX == K*X : " << err << endl;
1658 }
1659
1660 // H and friends
1661 if (chk.checkH && initialized_) {
1662 MagnitudeType err = orthman_->orthonormError(*H_);
1663 os << " >> Error in H^H M H == I : " << err << endl;
1664 err = orthman_->orthogError(*H_,*lclV);
1665 os << " >> Error in H^H M V == 0 : " << err << endl;
1666 err = orthman_->orthogError(*H_,*X_);
1667 os << " >> Error in H^H M X == 0 : " << err << endl;
1668 for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
1669 err = orthman_->orthogError(*H_,*auxVecs_[i]);
1670 os << " >> Error in H^H M Q[" << i << "] == 0 : " << err << endl;
1671 }
1672 }
1673 if (chk.checkKH && initialized_) {
1674 MagnitudeType err = Utils::errorEquality(*H_, *KH_, Op_);
1675 os << " >> Error in KH == K*H : " << err << endl;
1676 }
1677 if (chk.checkMH && hasM_ && initialized_) {
1678 MagnitudeType err = Utils::errorEquality(*H_, *MH_, MOp_);
1679 os << " >> Error in MH == M*H : " << err << endl;
1680 }
1681
1682 // R: this is not M-orthogonality, but standard Euclidean orthogonality
1683 if (chk.checkR && initialized_) {
1684 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
1685 MVT::MvTransMv(ONE,*X_,*R_,xTx);
1686 MagnitudeType err = xTx.normFrobenius();
1687 os << " >> Error in X^H R == 0 : " << err << endl;
1688 }
1689
1690 // KK
1691 if (chk.checkKK && initialized_) {
1692 Teuchos::SerialDenseMatrix<int,ScalarType> SDMerr(curDim_,curDim_), lclKK(Teuchos::View,*KK_,curDim_,curDim_);
1693 for (int j=0; j<curDim_; ++j) {
1694 for (int i=0; i<curDim_; ++i) {
1695 SDMerr(i,j) = lclKK(i,j) - SCT::conjugate(lclKK(j,i));
1696 }
1697 }
1698 os << " >> Error in KK - KK^H == 0 : " << SDMerr.normFrobenius() << endl;
1699 }
1700
1701 // Q
1702 if (chk.checkQ) {
1703 for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
1704 MagnitudeType err = orthman_->orthonormError(*auxVecs_[i]);
1705 os << " >> Error in Q[" << i << "]^H M Q[" << i << "] == I : " << err << endl;
1706 for (Array_size_type j=i+1; j<auxVecs_.size(); ++j) {
1707 err = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
1708 os << " >> Error in Q[" << i << "]^H M Q[" << j << "] == 0 : " << err << endl;
1709 }
1710 }
1711 }
1712
1713 os << endl;
1714
1715 return os.str();
1716 }
1717
1718
1720 // Print the current status of the solver
1721 template <class ScalarType, class MV, class OP>
1722 void
1724 {
1725 using std::endl;
1726
1727 os.setf(std::ios::scientific, std::ios::floatfield);
1728 os.precision(6);
1729 os <<endl;
1730 os <<"================================================================================" << endl;
1731 os << endl;
1732 os <<" BlockDavidson Solver Status" << endl;
1733 os << endl;
1734 os <<"The solver is "<<(initialized_ ? "initialized." : "not initialized.") << endl;
1735 os <<"The number of iterations performed is " <<iter_<<endl;
1736 os <<"The block size is " << blockSize_<<endl;
1737 os <<"The number of blocks is " << numBlocks_<<endl;
1738 os <<"The current basis size is " << curDim_<<endl;
1739 os <<"The number of auxiliary vectors is "<< numAuxVecs_ << endl;
1740 os <<"The number of operations Op*x is "<<count_ApplyOp_<<endl;
1741 os <<"The number of operations M*x is "<<count_ApplyM_<<endl;
1742 os <<"The number of operations Prec*x is "<<count_ApplyPrec_<<endl;
1743
1744 os.setf(std::ios_base::right, std::ios_base::adjustfield);
1745
1746 if (initialized_) {
1747 os << endl;
1748 os <<"CURRENT EIGENVALUE ESTIMATES "<<endl;
1749 os << std::setw(20) << "Eigenvalue"
1750 << std::setw(20) << "Residual(M)"
1751 << std::setw(20) << "Residual(2)"
1752 << endl;
1753 os <<"--------------------------------------------------------------------------------"<<endl;
1754 for (int i=0; i<blockSize_; ++i) {
1755 os << std::setw(20) << theta_[i];
1756 if (Rnorms_current_) os << std::setw(20) << Rnorms_[i];
1757 else os << std::setw(20) << "not current";
1758 if (R2norms_current_) os << std::setw(20) << R2norms_[i];
1759 else os << std::setw(20) << "not current";
1760 os << endl;
1761 }
1762 }
1763 os <<"================================================================================" << endl;
1764 os << endl;
1765 }
1766
1767} // End of namespace Anasazi
1768
1769#endif
1770
1771// End of file AnasaziBlockDavidson.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.
BlockDavidsonInitFailure is thrown when the BlockDavidson solver is unable to generate an initial ite...
BlockDavidsonOrthoFailure is thrown when the orthogonalization manager is unable to orthogonalize the...
This class implements a Block Davidson iteration, a preconditioned iteration for solving linear Hermi...
Teuchos::Array< Teuchos::RCP< const MV > > getAuxVecs() const
Get the auxiliary vectors for the solver.
void resetNumIters()
Reset the iteration count.
std::vector< int > getRitzIndex()
Get the index used for extracting individual Ritz vectors from getRitzVectors().
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current eigenvectors and eigenvalues.
void setBlockSize(int blockSize)
Set the blocksize.
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > getStatusTest() const
Get the current StatusTest used by the solver.
std::vector< Value< ScalarType > > getRitzValues()
Get the Ritz values for the previous iteration.
void setStatusTest(Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test)
Set a new StatusTest for the solver.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRitzRes2Norms()
Get the 2-norms of the residuals.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the eigenvalue problem.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRes2Norms()
Get the current residual 2-norms, computing the norms if they are not up-to-date with the current res...
void setAuxVecs(const Teuchos::Array< Teuchos::RCP< const MV > > &auxvecs)
Set the auxiliary vectors for the solver.
void initialize()
Initialize the solver with the initial vectors from the eigenproblem or random data.
BlockDavidsonState< ScalarType, MV > getState() const
Get access to the current state of the eigensolver.
void currentStatus(std::ostream &os)
This method requests that the solver print out its current status to the given output stream.
int getNumIters() const
Get the current iteration count.
BlockDavidson(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< MatOrthoManager< ScalarType, MV, OP > > &ortho, Teuchos::ParameterList &params)
BlockDavidson constructor with eigenproblem, solver utilities, and parameter list of solver options.
void iterate()
This method performs BlockDavidson iterations until the status test indicates the need to stop or an ...
void setSize(int blockSize, int numBlocks)
Set the blocksize and number of blocks to be used by the iterative solver in solving this eigenproble...
virtual ~BlockDavidson()
Anasazi::BlockDavidson destructor.
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, computing the norms if they are not up-to-date with the current resid...
Teuchos::RCP< const MV > getRitzVectors()
Get access to the current Ritz vectors.
int getBlockSize() const
Get the blocksize used by the iterative solver.
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace. For BlockDavidson, this always returns n...
The Eigensolver is a templated virtual base class that defines the basic interface that any eigensolv...
Anasazi's templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
Structure to contain pointers to BlockDavidson state variables.
int curDim
The current dimension of the solver.
Teuchos::RCP< const MV > V
The basis for the Krylov space.
Teuchos::RCP< const MV > H
The current preconditioned residual vectors.
Teuchos::RCP< const MV > X
The current eigenvectors.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > KK
The current projected K matrix.
Teuchos::RCP< const MV > R
The current residual vectors.
Teuchos::RCP< const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > T
The current Ritz values. This vector is a copy of the internal data.
Teuchos::RCP< const MV > KX
The image of the current eigenvectors under K.
Teuchos::RCP< const MV > MX
The image of the current eigenvectors under M, or Teuchos::null if M was not specified.