Anasazi Version of the Day
Loading...
Searching...
No Matches
AnasaziBlockKrylovSchur.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_KRYLOV_SCHUR_HPP
47#define ANASAZI_BLOCK_KRYLOV_SCHUR_HPP
48
49#include "AnasaziTypes.hpp"
50
54#include "Teuchos_ScalarTraits.hpp"
55#include "AnasaziHelperTraits.hpp"
56
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
79namespace Anasazi {
80
82
83
88 template <class ScalarType, class MulVec>
94 int curDim;
96 Teuchos::RCP<const MulVec> V;
102 Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > H;
104 Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > S;
106 Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > Q;
107 BlockKrylovSchurState() : curDim(0), V(Teuchos::null),
108 H(Teuchos::null), S(Teuchos::null),
109 Q(Teuchos::null) {}
110 };
111
113
115
116
130 BlockKrylovSchurInitFailure(const std::string& what_arg) : AnasaziError(what_arg)
131 {}};
132
140 BlockKrylovSchurOrthoFailure(const std::string& what_arg) : AnasaziError(what_arg)
141 {}};
142
144
145
146 template <class ScalarType, class MV, class OP>
147 class BlockKrylovSchur : public Eigensolver<ScalarType,MV,OP> {
148 public:
150
151
163 BlockKrylovSchur( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
164 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
165 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
166 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
167 const Teuchos::RCP<OrthoManager<ScalarType,MV> > &ortho,
168 Teuchos::ParameterList &params
169 );
170
172 virtual ~BlockKrylovSchur() {};
174
175
177
178
200 void iterate();
201
223 void initialize(BlockKrylovSchurState<ScalarType,MV>& state);
224
228 void initialize();
229
238 bool isInitialized() const { return initialized_; }
239
247 BlockKrylovSchurState<ScalarType,MV> getState() const {
248 BlockKrylovSchurState<ScalarType,MV> state;
249 state.curDim = curDim_;
250 state.V = V_;
251 state.H = H_;
252 state.Q = Q_;
253 state.S = schurH_;
254 return state;
255 }
256
258
259
261
262
264 int getNumIters() const { return(iter_); }
265
267 void resetNumIters() { iter_=0; }
268
276 Teuchos::RCP<const MV> getRitzVectors() { return ritzVectors_; }
277
285 std::vector<Value<ScalarType> > getRitzValues() {
286 std::vector<Value<ScalarType> > ret = ritzValues_;
287 ret.resize(ritzIndex_.size());
288 return ret;
289 }
290
296 std::vector<int> getRitzIndex() { return ritzIndex_; }
297
302 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getResNorms() {
303 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> ret(0);
304 return ret;
305 }
306
311 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRes2Norms() {
312 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> ret(0);
313 return ret;
314 }
315
320 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRitzRes2Norms() {
321 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> ret = ritzResiduals_;
322 ret.resize(ritzIndex_.size());
323 return ret;
324 }
325
327
329
330
332 void setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test);
333
335 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > getStatusTest() const;
336
338 const Eigenproblem<ScalarType,MV,OP>& getProblem() const { return(*problem_); };
339
346 void setSize(int blockSize, int numBlocks);
347
349 void setBlockSize(int blockSize);
350
352 void setStepSize(int stepSize);
353
355 void setNumRitzVectors(int numRitzVecs);
356
358 int getStepSize() const { return(stepSize_); }
359
361 int getBlockSize() const { return(blockSize_); }
362
364 int getNumRitzVectors() const { return(numRitzVecs_); }
365
371 int getCurSubspaceDim() const {
372 if (!initialized_) return 0;
373 return curDim_;
374 }
375
377 int getMaxSubspaceDim() const { return (problem_->isHermitian()?blockSize_*numBlocks_:blockSize_*numBlocks_+1); }
378
379
392 void setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs);
393
395 Teuchos::Array<Teuchos::RCP<const MV> > getAuxVecs() const {return auxVecs_;}
396
398
400
401
403 void currentStatus(std::ostream &os);
404
406
408
409
411 bool isRitzVecsCurrent() const { return ritzVecsCurrent_; }
412
414 bool isRitzValsCurrent() const { return ritzValsCurrent_; }
415
417 bool isSchurCurrent() const { return schurCurrent_; }
418
420
422
423
425 void computeRitzVectors();
426
428 void computeRitzValues();
429
431 void computeSchurForm( const bool sort = true );
432
434
435 private:
436 //
437 // Convenience typedefs
438 //
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 typedef typename std::vector<ScalarType>::iterator STiter;
444 typedef typename std::vector<MagnitudeType>::iterator MTiter;
445 const MagnitudeType MT_ONE;
446 const MagnitudeType MT_ZERO;
447 const MagnitudeType NANVAL;
448 const ScalarType ST_ONE;
449 const ScalarType ST_ZERO;
450 //
451 // Internal structs
452 //
453 struct CheckList {
454 bool checkV;
455 bool checkArn;
456 bool checkAux;
457 CheckList() : checkV(false), checkArn(false), checkAux(false) {};
458 };
459 //
460 // Internal methods
461 //
462 std::string accuracyCheck(const CheckList &chk, const std::string &where) const;
463 void sortSchurForm( Teuchos::SerialDenseMatrix<int,ScalarType>& H,
464 Teuchos::SerialDenseMatrix<int,ScalarType>& Q,
465 std::vector<int>& order );
466 //
467 // Classes inputed through constructor that define the eigenproblem to be solved.
468 //
469 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
470 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sm_;
471 const Teuchos::RCP<OutputManager<ScalarType> > om_;
472 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > tester_;
473 const Teuchos::RCP<OrthoManager<ScalarType,MV> > orthman_;
474 //
475 // Information obtained from the eigenproblem
476 //
477 Teuchos::RCP<const OP> Op_;
478 //
479 // Internal timers
480 //
481 Teuchos::RCP<Teuchos::Time> timerOp_, timerSortRitzVal_,
482 timerCompSF_, timerSortSF_,
483 timerCompRitzVec_, timerOrtho_;
484 //
485 // Counters
486 //
487 int count_ApplyOp_;
488
489 //
490 // Algorithmic parameters.
491 //
492 // blockSize_ is the solver block size; it controls the number of eigenvectors that
493 // we compute, the number of residual vectors that we compute, and therefore the number
494 // of vectors added to the basis on each iteration.
495 int blockSize_;
496 // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
497 int numBlocks_;
498 // stepSize_ dictates how many iterations are performed before eigenvectors and eigenvalues
499 // are computed again
500 int stepSize_;
501
502 //
503 // Current solver state
504 //
505 // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
506 // is capable of running; _initialize is controlled by the initialize() member method
507 // For the implications of the state of initialized_, please see documentation for initialize()
508 bool initialized_;
509 //
510 // curDim_ reflects how much of the current basis is valid
511 // NOTE: for Hermitian, 0 <= curDim_ <= blockSize_*numBlocks_
512 // for non-Hermitian, 0 <= curDim_ <= blockSize_*numBlocks_ + 1
513 // this also tells us how many of the values in _theta are valid Ritz values
514 int curDim_;
515 //
516 // State Multivecs
517 Teuchos::RCP<MV> ritzVectors_, V_;
518 int numRitzVecs_;
519 //
520 // Projected matrices
521 // H_ : Projected matrix from the Krylov-Schur factorization AV = VH + FB^T
522 //
523 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H_;
524 //
525 // Schur form of Projected matrices (these are only updated when the Ritz values/vectors are updated).
526 // schurH_: Schur form reduction of H
527 // Q_: Schur vectors of H
528 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > schurH_;
529 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Q_;
530 //
531 // Auxiliary vectors
532 Teuchos::Array<Teuchos::RCP<const MV> > auxVecs_;
533 int numAuxVecs_;
534 //
535 // Number of iterations that have been performed.
536 int iter_;
537 //
538 // State flags
539 bool ritzVecsCurrent_, ritzValsCurrent_, schurCurrent_;
540 //
541 // Current eigenvalues, residual norms
542 std::vector<Value<ScalarType> > ritzValues_;
543 std::vector<MagnitudeType> ritzResiduals_;
544 //
545 // Current index vector for Ritz values and vectors
546 std::vector<int> ritzIndex_; // computed by BKS
547 std::vector<int> ritzOrder_; // returned from sort manager
548 //
549 // Number of Ritz pairs to be printed upon output, if possible
550 int numRitzPrint_;
551 };
552
553
555 // Constructor
556 template <class ScalarType, class MV, class OP>
558 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
559 const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
560 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
561 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
562 const Teuchos::RCP<OrthoManager<ScalarType,MV> > &ortho,
563 Teuchos::ParameterList &params
564 ) :
565 MT_ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
566 MT_ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
567 NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan()),
568 ST_ONE(Teuchos::ScalarTraits<ScalarType>::one()),
569 ST_ZERO(Teuchos::ScalarTraits<ScalarType>::zero()),
570 // problem, tools
571 problem_(problem),
572 sm_(sorter),
573 om_(printer),
574 tester_(tester),
575 orthman_(ortho),
576 // timers, counters
577#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
578 timerOp_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchur::Operation Op*x")),
579 timerSortRitzVal_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchur::Sorting Ritz values")),
580 timerCompSF_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchur::Computing Schur form")),
581 timerSortSF_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchur::Sorting Schur form")),
582 timerCompRitzVec_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchur::Computing Ritz vectors")),
583 timerOrtho_(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockKrylovSchur::Orthogonalization")),
584#endif
585 count_ApplyOp_(0),
586 // internal data
587 blockSize_(0),
588 numBlocks_(0),
589 stepSize_(0),
590 initialized_(false),
591 curDim_(0),
592 numRitzVecs_(0),
593 auxVecs_( Teuchos::Array<Teuchos::RCP<const MV> >(0) ),
594 numAuxVecs_(0),
595 iter_(0),
596 ritzVecsCurrent_(false),
597 ritzValsCurrent_(false),
598 schurCurrent_(false),
599 numRitzPrint_(0)
600 {
601 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument,
602 "Anasazi::BlockKrylovSchur::constructor: user specified null problem pointer.");
603 TEUCHOS_TEST_FOR_EXCEPTION(sm_ == Teuchos::null,std::invalid_argument,
604 "Anasazi::BlockKrylovSchur::constructor: user passed null sort manager pointer.");
605 TEUCHOS_TEST_FOR_EXCEPTION(om_ == Teuchos::null,std::invalid_argument,
606 "Anasazi::BlockKrylovSchur::constructor: user passed null output manager pointer.");
607 TEUCHOS_TEST_FOR_EXCEPTION(tester_ == Teuchos::null,std::invalid_argument,
608 "Anasazi::BlockKrylovSchur::constructor: user passed null status test pointer.");
609 TEUCHOS_TEST_FOR_EXCEPTION(orthman_ == Teuchos::null,std::invalid_argument,
610 "Anasazi::BlockKrylovSchur::constructor: user passed null orthogonalization manager pointer.");
611 TEUCHOS_TEST_FOR_EXCEPTION(problem_->isProblemSet() == false, std::invalid_argument,
612 "Anasazi::BlockKrylovSchur::constructor: user specified problem is not set.");
613 TEUCHOS_TEST_FOR_EXCEPTION(sorter == Teuchos::null,std::invalid_argument,
614 "Anasazi::BlockKrylovSchur::constructor: user specified null sort manager pointer.");
615 TEUCHOS_TEST_FOR_EXCEPTION(printer == Teuchos::null,std::invalid_argument,
616 "Anasazi::BlockKrylovSchur::constructor: user specified null output manager pointer.");
617 TEUCHOS_TEST_FOR_EXCEPTION(tester == Teuchos::null,std::invalid_argument,
618 "Anasazi::BlockKrylovSchur::constructor: user specified null status test pointer.");
619 TEUCHOS_TEST_FOR_EXCEPTION(ortho == Teuchos::null,std::invalid_argument,
620 "Anasazi::BlockKrylovSchur::constructor: user specified null ortho manager pointer.");
621
622 // Get problem operator
623 Op_ = problem_->getOperator();
624
625 // get the step size
626 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Step Size"), std::invalid_argument,
627 "Anasazi::BlockKrylovSchur::constructor: mandatory parameter 'Step Size' is not specified.");
628 int ss = params.get("Step Size",numBlocks_);
629 setStepSize(ss);
630
631 // set the block size and allocate data
632 int bs = params.get("Block Size", 1);
633 int nb = params.get("Num Blocks", 3*problem_->getNEV());
634 setSize(bs,nb);
635
636 // get the number of Ritz vectors to compute and allocate data.
637 // --> if this parameter is not specified in the parameter list, then it's assumed that no Ritz vectors will be computed.
638 int numRitzVecs = params.get("Number of Ritz Vectors", 0);
639 setNumRitzVectors( numRitzVecs );
640
641 // get the number of Ritz values to print out when currentStatus is called.
642 numRitzPrint_ = params.get("Print Number of Ritz Values", bs);
643 }
644
645
647 // Set the block size
648 // This simply calls setSize(), modifying the block size while retaining the number of blocks.
649 template <class ScalarType, class MV, class OP>
651 {
652 setSize(blockSize,numBlocks_);
653 }
654
655
657 // Set the step size.
658 template <class ScalarType, class MV, class OP>
660 {
661 TEUCHOS_TEST_FOR_EXCEPTION(stepSize <= 0, std::invalid_argument, "Anasazi::BlockKrylovSchur::setStepSize(): new step size must be positive and non-zero.");
662 stepSize_ = stepSize;
663 }
664
666 // Set the number of Ritz vectors to compute.
667 template <class ScalarType, class MV, class OP>
669 {
670 // This routine only allocates space; it doesn't not perform any computation
671 // any change in size will invalidate the state of the solver.
672
673 TEUCHOS_TEST_FOR_EXCEPTION(numRitzVecs < 0, std::invalid_argument, "Anasazi::BlockKrylovSchur::setNumRitzVectors(): number of Ritz vectors to compute must be positive.");
674
675 // Check to see if the number of requested Ritz vectors has changed.
676 if (numRitzVecs != numRitzVecs_) {
677 if (numRitzVecs) {
678 ritzVectors_ = Teuchos::null;
679 ritzVectors_ = MVT::Clone(*V_, numRitzVecs);
680 } else {
681 ritzVectors_ = Teuchos::null;
682 }
683 numRitzVecs_ = numRitzVecs;
684 ritzVecsCurrent_ = false;
685 }
686 }
687
689 // Set the block size and make necessary adjustments.
690 template <class ScalarType, class MV, class OP>
691 void BlockKrylovSchur<ScalarType,MV,OP>::setSize (int blockSize, int numBlocks)
692 {
693 // This routine only allocates space; it doesn't not perform any computation
694 // any change in size will invalidate the state of the solver.
695
696 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0 || blockSize <= 0, std::invalid_argument, "Anasazi::BlockKrylovSchur::setSize was passed a non-positive argument.");
697 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks < 3, std::invalid_argument, "Anasazi::BlockKrylovSchur::setSize(): numBlocks must be at least three.");
698 if (blockSize == blockSize_ && numBlocks == numBlocks_) {
699 // do nothing
700 return;
701 }
702
703 blockSize_ = blockSize;
704 numBlocks_ = numBlocks;
705
706 Teuchos::RCP<const MV> tmp;
707 // grab some Multivector to Clone
708 // in practice, getInitVec() should always provide this, but it is possible to use a
709 // Eigenproblem with nothing in getInitVec() by manually initializing with initialize();
710 // in case of that strange scenario, we will try to Clone from V_; first resort to getInitVec(),
711 // because we would like to clear the storage associated with V_ so we have room for the new V_
712 if (problem_->getInitVec() != Teuchos::null) {
713 tmp = problem_->getInitVec();
714 }
715 else {
716 tmp = V_;
717 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
718 "Anasazi::BlockKrylovSchur::setSize(): eigenproblem did not specify initial vectors to clone from.");
719 }
720
721
723 // blockSize*numBlocks dependent
724 //
725 int newsd;
726 if (problem_->isHermitian()) {
727 newsd = blockSize_*numBlocks_;
728 } else {
729 newsd = blockSize_*numBlocks_+1;
730 }
731 // check that new size is valid
732 TEUCHOS_TEST_FOR_EXCEPTION(static_cast<ptrdiff_t>(newsd) > MVT::GetGlobalLength(*tmp),std::invalid_argument,
733 "Anasazi::BlockKrylovSchur::setSize(): maximum basis size is larger than problem dimension.");
734
735 ritzValues_.resize(newsd);
736 ritzResiduals_.resize(newsd,MT_ONE);
737 ritzOrder_.resize(newsd);
738 V_ = Teuchos::null;
739 V_ = MVT::Clone(*tmp,newsd+blockSize_);
740 H_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(newsd+blockSize_,newsd) );
741 Q_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(newsd,newsd) );
742
743 initialized_ = false;
744 curDim_ = 0;
745 }
746
747
749 // Set the auxiliary vectors
750 template <class ScalarType, class MV, class OP>
751 void BlockKrylovSchur<ScalarType,MV,OP>::setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs) {
752 typedef typename Teuchos::Array<Teuchos::RCP<const MV> >::iterator tarcpmv;
753
754 // set new auxiliary vectors
755 auxVecs_ = auxvecs;
756
757 if (om_->isVerbosity( Debug ) ) {
758 // Check almost everything here
759 CheckList chk;
760 chk.checkAux = true;
761 om_->print( Debug, accuracyCheck(chk, ": in setAuxVecs()") );
762 }
763
764 numAuxVecs_ = 0;
765 for (tarcpmv i=auxVecs_.begin(); i != auxVecs_.end(); i++) {
766 numAuxVecs_ += MVT::GetNumberVecs(**i);
767 }
768
769 // If the solver has been initialized, X and P are not necessarily orthogonal to new auxiliary vectors
770 if (numAuxVecs_ > 0 && initialized_) {
771 initialized_ = false;
772 }
773 }
774
776 /* Initialize the state of the solver
777 *
778 * POST-CONDITIONS:
779 *
780 * V_ is orthonormal, orthogonal to auxVecs_, for first curDim_ vectors
781 *
782 */
783
784 template <class ScalarType, class MV, class OP>
785 void BlockKrylovSchur<ScalarType,MV,OP>::initialize(BlockKrylovSchurState<ScalarType,MV>& newstate)
786 {
787 // NOTE: memory has been allocated by setBlockSize(). Use SetBlock below; do not Clone
788
789 std::vector<int> bsind(blockSize_);
790 for (int i=0; i<blockSize_; i++) bsind[i] = i;
791
792 // in BlockKrylovSchur, V and H are required.
793 // if either doesn't exist, then we will start with the initial vector.
794 //
795 // inconsistent multivectors widths and lengths will not be tolerated, and
796 // will be treated with exceptions.
797 //
798 std::string errstr("Anasazi::BlockKrylovSchur::initialize(): specified multivectors must have a consistent length and width.");
799
800 // set up V,H: if the user doesn't specify both of these these,
801 // we will start over with the initial vector.
802 if (newstate.V != Teuchos::null && newstate.H != Teuchos::null) {
803
804 // initialize V_,H_, and curDim_
805
806 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.V) != MVT::GetGlobalLength(*V_),
807 std::invalid_argument, errstr );
808 if (newstate.V != V_) {
809 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < blockSize_,
810 std::invalid_argument, errstr );
811 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) > getMaxSubspaceDim(),
812 std::invalid_argument, errstr );
813 }
814 TEUCHOS_TEST_FOR_EXCEPTION( newstate.curDim > getMaxSubspaceDim(),
815 std::invalid_argument, errstr );
816
817 curDim_ = newstate.curDim;
818 int lclDim = MVT::GetNumberVecs(*newstate.V);
819
820 // check size of H
821 TEUCHOS_TEST_FOR_EXCEPTION(newstate.H->numRows() < curDim_ || newstate.H->numCols() < curDim_, std::invalid_argument, errstr);
822
823 if (curDim_ == 0 && lclDim > blockSize_) {
824 om_->stream(Warnings) << "Anasazi::BlockKrylovSchur::initialize(): the solver was initialized with a kernel of " << lclDim << std::endl
825 << "The block size however is only " << blockSize_ << std::endl
826 << "The last " << lclDim - blockSize_ << " vectors of the kernel will be overwritten on the first call to iterate()." << std::endl;
827 }
828
829
830 // copy basis vectors from newstate into V
831 if (newstate.V != V_) {
832 std::vector<int> nevind(lclDim);
833 for (int i=0; i<lclDim; i++) nevind[i] = i;
834 MVT::SetBlock(*newstate.V,nevind,*V_);
835 }
836
837 // put data into H_, make sure old information is not still hanging around.
838 if (newstate.H != H_) {
839 H_->putScalar( ST_ZERO );
840 Teuchos::SerialDenseMatrix<int,ScalarType> newH(Teuchos::View,*newstate.H,curDim_+blockSize_,curDim_);
841 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclH;
842 lclH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*H_,curDim_+blockSize_,curDim_) );
843 lclH->assign(newH);
844
845 // done with local pointers
846 lclH = Teuchos::null;
847 }
848
849 }
850 else {
851 // user did not specify a basis V
852 // get vectors from problem or generate something, projectAndNormalize, call initialize() recursively
853 Teuchos::RCP<const MV> ivec = problem_->getInitVec();
854 TEUCHOS_TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::invalid_argument,
855 "Anasazi::BlockKrylovSchur::initialize(): eigenproblem did not specify initial vectors to clone from.");
856
857 int lclDim = MVT::GetNumberVecs(*ivec);
858 bool userand = false;
859 if (lclDim < blockSize_) {
860 // we need at least blockSize_ vectors
861 // use a random multivec
862 userand = true;
863 }
864
865 if (userand) {
866 // make an index
867 std::vector<int> dimind2(lclDim);
868 for (int i=0; i<lclDim; i++) { dimind2[i] = i; }
869
870 // alloc newV as a view of the first block of V
871 Teuchos::RCP<MV> newV1 = MVT::CloneViewNonConst(*V_,dimind2);
872
873 // copy the initial vectors into the first lclDim vectors of V
874 MVT::SetBlock(*ivec,dimind2,*newV1);
875
876 // resize / reinitialize the index vector
877 dimind2.resize(blockSize_-lclDim);
878 for (int i=0; i<blockSize_-lclDim; i++) { dimind2[i] = lclDim + i; }
879
880 // initialize the rest of the vectors with random vectors
881 Teuchos::RCP<MV> newV2 = MVT::CloneViewNonConst(*V_,dimind2);
882 MVT::MvRandom(*newV2);
883 }
884 else {
885 // alloc newV as a view of the first block of V
886 Teuchos::RCP<MV> newV1 = MVT::CloneViewNonConst(*V_,bsind);
887
888 // get a view of the first block of initial vectors
889 Teuchos::RCP<const MV> ivecV = MVT::CloneView(*ivec,bsind);
890
891 // assign ivec to first part of newV
892 MVT::SetBlock(*ivecV,bsind,*newV1);
893 }
894
895 // get pointer into first block of V
896 Teuchos::RCP<MV> newV = MVT::CloneViewNonConst(*V_,bsind);
897
898 // remove auxVecs from newV and normalize newV
899 if (auxVecs_.size() > 0) {
900#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
901 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
902#endif
903
904 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummy;
905 int rank = orthman_->projectAndNormalize(*newV,auxVecs_);
906 TEUCHOS_TEST_FOR_EXCEPTION( rank != blockSize_,BlockKrylovSchurInitFailure,
907 "Anasazi::BlockKrylovSchur::initialize(): couldn't generate initial basis of full rank." );
908 }
909 else {
910#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
911 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
912#endif
913
914 int rank = orthman_->normalize(*newV);
915 TEUCHOS_TEST_FOR_EXCEPTION( rank != blockSize_,BlockKrylovSchurInitFailure,
916 "Anasazi::BlockKrylovSchur::initialize(): couldn't generate initial basis of full rank." );
917 }
918
919 // set curDim
920 curDim_ = 0;
921
922 // clear pointer
923 newV = Teuchos::null;
924 }
925
926 // The Ritz vectors/values and Schur form are no longer current.
927 ritzVecsCurrent_ = false;
928 ritzValsCurrent_ = false;
929 schurCurrent_ = false;
930
931 // the solver is initialized
932 initialized_ = true;
933
934 if (om_->isVerbosity( Debug ) ) {
935 // Check almost everything here
936 CheckList chk;
937 chk.checkV = true;
938 chk.checkArn = true;
939 chk.checkAux = true;
940 om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
941 }
942
943 // Print information on current status
944 if (om_->isVerbosity(Debug)) {
945 currentStatus( om_->stream(Debug) );
946 }
947 else if (om_->isVerbosity(IterationDetails)) {
948 currentStatus( om_->stream(IterationDetails) );
949 }
950 }
951
952
954 // initialize the solver with default state
955 template <class ScalarType, class MV, class OP>
957 {
958 BlockKrylovSchurState<ScalarType,MV> empty;
959 initialize(empty);
960 }
961
962
964 // Perform BlockKrylovSchur iterations until the StatusTest tells us to stop.
965 template <class ScalarType, class MV, class OP>
967 {
968 //
969 // Allocate/initialize data structures
970 //
971 if (initialized_ == false) {
972 initialize();
973 }
974
975 // Compute the current search dimension.
976 // If the problem is non-Hermitian and the blocksize is one, let the solver use the extra vector.
977 int searchDim = blockSize_*numBlocks_;
978 if (problem_->isHermitian() == false) {
979 searchDim++;
980 }
981
983 // iterate until the status test tells us to stop.
984 //
985 // also break if our basis is full
986 //
987 while (tester_->checkStatus(this) != Passed && curDim_+blockSize_ <= searchDim) {
988
989 iter_++;
990
991 // F can be found at the curDim_ block, but the next block is at curDim_ + blockSize_.
992 int lclDim = curDim_ + blockSize_;
993
994 // Get the current part of the basis.
995 std::vector<int> curind(blockSize_);
996 for (int i=0; i<blockSize_; i++) { curind[i] = lclDim + i; }
997 Teuchos::RCP<MV> Vnext = MVT::CloneViewNonConst(*V_,curind);
998
999 // Get a view of the previous vectors
1000 // this is used for orthogonalization and for computing V^H K H
1001 for (int i=0; i<blockSize_; i++) { curind[i] = curDim_ + i; }
1002 Teuchos::RCP<const MV> Vprev = MVT::CloneView(*V_,curind);
1003 // om_->stream(Debug) << "Vprev: " << std::endl;
1004 // MVT::MvPrint(*Vprev,om_->stream(Debug));
1005
1006 // Compute the next vector in the Krylov basis: Vnext = Op*Vprev
1007 {
1008#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1009 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1010#endif
1011 OPT::Apply(*Op_,*Vprev,*Vnext);
1012 count_ApplyOp_ += blockSize_;
1013 }
1014 // om_->stream(Debug) << "Vnext: " << std::endl;
1015 // MVT::MvPrint(*Vnext,om_->stream(Debug));
1016 Vprev = Teuchos::null;
1017
1018 // Remove all previous Krylov-Schur basis vectors and auxVecs from Vnext
1019 {
1020#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1021 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1022#endif
1023
1024 // Get a view of all the previous vectors
1025 std::vector<int> prevind(lclDim);
1026 for (int i=0; i<lclDim; i++) { prevind[i] = i; }
1027 Vprev = MVT::CloneView(*V_,prevind);
1028 Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
1029
1030 // Get a view of the part of the Hessenberg matrix needed to hold the ortho coeffs.
1031 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
1032 subH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
1033 ( Teuchos::View,*H_,lclDim,blockSize_,0,curDim_ ) );
1034 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubH;
1035 AsubH.append( subH );
1036
1037 // Add the auxiliary vectors to the current basis vectors if any exist
1038 if (auxVecs_.size() > 0) {
1039 for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1040 AVprev.append( auxVecs_[i] );
1041 AsubH.append( Teuchos::null );
1042 }
1043 }
1044
1045 // Get a view of the part of the Hessenberg matrix needed to hold the norm coeffs.
1046 // om_->stream(Debug) << "V before ortho: " << std::endl;
1047 // MVT::MvPrint(*Vprev,om_->stream(Debug));
1048 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
1049 subR = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
1050 ( Teuchos::View,*H_,blockSize_,blockSize_,lclDim,curDim_ ) );
1051 int rank = orthman_->projectAndNormalize(*Vnext,AVprev,AsubH,subR);
1052 // om_->stream(Debug) << "Vnext after ortho: " << std::endl;
1053 // MVT::MvPrint(*Vnext,om_->stream(Debug));
1054 // om_->stream(Debug) << "subH: " << std::endl << *subH << std::endl;
1055 // om_->stream(Debug) << "subR: " << std::endl << *subR << std::endl;
1056 // om_->stream(Debug) << "H: " << std::endl << *H_ << std::endl;
1057 TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_,BlockKrylovSchurOrthoFailure,
1058 "Anasazi::BlockKrylovSchur::iterate(): couldn't generate basis of full rank.");
1059 }
1060 //
1061 // V has been extended, and H has been extended.
1062 //
1063 // Update basis dim and release all pointers.
1064 Vnext = Teuchos::null;
1065 curDim_ += blockSize_;
1066 // The Ritz vectors/values and Schur form are no longer current.
1067 ritzVecsCurrent_ = false;
1068 ritzValsCurrent_ = false;
1069 schurCurrent_ = false;
1070 //
1071 // Update Ritz values and residuals if needed
1072 if (!(iter_%stepSize_)) {
1073 computeRitzValues();
1074 }
1075
1076 // When required, monitor some orthogonalities
1077 if (om_->isVerbosity( Debug ) ) {
1078 // Check almost everything here
1079 CheckList chk;
1080 chk.checkV = true;
1081 chk.checkArn = true;
1082 om_->print( Debug, accuracyCheck(chk, ": after local update") );
1083 }
1084 else if (om_->isVerbosity( OrthoDetails ) ) {
1085 CheckList chk;
1086 chk.checkV = true;
1087 om_->print( OrthoDetails, accuracyCheck(chk, ": after local update") );
1088 }
1089
1090 // Print information on current iteration
1091 if (om_->isVerbosity(Debug)) {
1092 currentStatus( om_->stream(Debug) );
1093 }
1094 else if (om_->isVerbosity(IterationDetails)) {
1095 currentStatus( om_->stream(IterationDetails) );
1096 }
1097
1098 } // end while (statusTest == false)
1099
1100 } // end of iterate()
1101
1102
1104 // Check accuracy, orthogonality, and other debugging stuff
1105 //
1106 // bools specify which tests we want to run (instead of running more than we actually care about)
1107 //
1108 // checkV : V orthonormal
1109 // orthogonal to auxvecs
1110 // checkAux: check that auxiliary vectors are actually orthonormal
1111 //
1112 // checkArn: check the Arnoldi factorization
1113 //
1114 // NOTE: This method needs to check the current dimension of the subspace, since it is possible to
1115 // call this method when curDim_ = 0 (after initialization).
1116 template <class ScalarType, class MV, class OP>
1117 std::string BlockKrylovSchur<ScalarType,MV,OP>::accuracyCheck( const CheckList &chk, const std::string &where ) const
1118 {
1119 std::stringstream os;
1120 os.precision(2);
1121 os.setf(std::ios::scientific, std::ios::floatfield);
1122 MagnitudeType tmp;
1123
1124 os << " Debugging checks: iteration " << iter_ << where << std::endl;
1125
1126 // index vectors for V and F
1127 std::vector<int> lclind(curDim_);
1128 for (int i=0; i<curDim_; i++) lclind[i] = i;
1129 std::vector<int> bsind(blockSize_);
1130 for (int i=0; i<blockSize_; i++) { bsind[i] = curDim_ + i; }
1131
1132 Teuchos::RCP<const MV> lclV,lclF;
1133 Teuchos::RCP<MV> lclAV;
1134 if (curDim_)
1135 lclV = MVT::CloneView(*V_,lclind);
1136 lclF = MVT::CloneView(*V_,bsind);
1137
1138 if (chk.checkV) {
1139 if (curDim_) {
1140 tmp = orthman_->orthonormError(*lclV);
1141 os << " >> Error in V^H M V == I : " << tmp << std::endl;
1142 }
1143 tmp = orthman_->orthonormError(*lclF);
1144 os << " >> Error in F^H M F == I : " << tmp << std::endl;
1145 if (curDim_) {
1146 tmp = orthman_->orthogError(*lclV,*lclF);
1147 os << " >> Error in V^H M F == 0 : " << tmp << std::endl;
1148 }
1149 for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1150 if (curDim_) {
1151 tmp = orthman_->orthogError(*lclV,*auxVecs_[i]);
1152 os << " >> Error in V^H M Aux[" << i << "] == 0 : " << tmp << std::endl;
1153 }
1154 tmp = orthman_->orthogError(*lclF,*auxVecs_[i]);
1155 os << " >> Error in F^H M Aux[" << i << "] == 0 : " << tmp << std::endl;
1156 }
1157 }
1158
1159 if (chk.checkArn) {
1160
1161 if (curDim_) {
1162 // Compute AV
1163 lclAV = MVT::Clone(*V_,curDim_);
1164 {
1165#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1166 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1167#endif
1168 OPT::Apply(*Op_,*lclV,*lclAV);
1169 }
1170
1171 // Compute AV - VH
1172 Teuchos::SerialDenseMatrix<int,ScalarType> subH(Teuchos::View,*H_,curDim_,curDim_);
1173 MVT::MvTimesMatAddMv( -ST_ONE, *lclV, subH, ST_ONE, *lclAV );
1174
1175 // Compute FB_k^T - (AV-VH)
1176 Teuchos::SerialDenseMatrix<int,ScalarType> curB(Teuchos::View,*H_,
1177 blockSize_,curDim_, curDim_ );
1178 MVT::MvTimesMatAddMv( -ST_ONE, *lclF, curB, ST_ONE, *lclAV );
1179
1180 // Compute || FE_k^T - (AV-VH) ||
1181 std::vector<MagnitudeType> arnNorms( curDim_ );
1182 orthman_->norm( *lclAV, arnNorms );
1183
1184 for (int i=0; i<curDim_; i++) {
1185 os << " >> Error in Krylov-Schur factorization (R = AV-VS-FB^H), ||R[" << i << "]|| : " << arnNorms[i] << std::endl;
1186 }
1187 }
1188 }
1189
1190 if (chk.checkAux) {
1191 for (Array_size_type i=0; i<auxVecs_.size(); i++) {
1192 tmp = orthman_->orthonormError(*auxVecs_[i]);
1193 os << " >> Error in Aux[" << i << "]^H M Aux[" << i << "] == I : " << tmp << std::endl;
1194 for (Array_size_type j=i+1; j<auxVecs_.size(); j++) {
1195 tmp = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
1196 os << " >> Error in Aux[" << i << "]^H M Aux[" << j << "] == 0 : " << tmp << std::endl;
1197 }
1198 }
1199 }
1200
1201 os << std::endl;
1202
1203 return os.str();
1204 }
1205
1207 /* Get the current approximate eigenvalues, i.e. Ritz values.
1208 *
1209 * POST-CONDITIONS:
1210 *
1211 * ritzValues_ contains Ritz w.r.t. V, H
1212 * Q_ contains the Schur vectors w.r.t. H
1213 * schurH_ contains the Schur matrix w.r.t. H
1214 * ritzOrder_ contains the current ordering from sort manager
1215 */
1216
1217 template <class ScalarType, class MV, class OP>
1219 {
1220 // Can only call this if the solver is initialized
1221 if (initialized_) {
1222
1223 // This just updates the Ritz values and residuals.
1224 // --> ritzValsCurrent_ will be set to 'true' by this method.
1225 if (!ritzValsCurrent_) {
1226 // Compute the current Ritz values, through computing the Schur form
1227 // without updating the current projection matrix or sorting the Schur form.
1228 computeSchurForm( false );
1229 }
1230 }
1231 }
1232
1234 /* Get the current approximate eigenvectors, i.e. Ritz vectors.
1235 *
1236 * POST-CONDITIONS:
1237 *
1238 * ritzValues_ contains Ritz w.r.t. V, H
1239 * ritzVectors_ is first blockSize_ Ritz vectors w.r.t. V, H
1240 * Q_ contains the Schur vectors w.r.t. H
1241 * schurH_ contains the Schur matrix w.r.t. H
1242 * ritzOrder_ contains the current ordering from sort manager
1243 */
1244
1245 template <class ScalarType, class MV, class OP>
1247 {
1248#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1249 Teuchos::TimeMonitor LocalTimer(*timerCompRitzVec_);
1250#endif
1251
1252 TEUCHOS_TEST_FOR_EXCEPTION(numRitzVecs_==0, std::invalid_argument,
1253 "Anasazi::BlockKrylovSchur::computeRitzVectors(): no Ritz vectors were required from this solver.");
1254
1255 TEUCHOS_TEST_FOR_EXCEPTION(curDim_ < numRitzVecs_, std::invalid_argument,
1256 "Anasazi::BlockKrylovSchur::computeRitzVectors(): the current subspace is not large enough to compute the number of requested Ritz vectors.");
1257
1258
1259 // Check to see if the current subspace dimension is non-trivial and the solver is initialized
1260 if (curDim_ && initialized_) {
1261
1262 // Check to see if the Ritz vectors are current.
1263 if (!ritzVecsCurrent_) {
1264
1265 // Check to see if the Schur factorization of H (schurH_, Q) is current and sorted.
1266 if (!schurCurrent_) {
1267 // Compute the Schur factorization of the current H, which will not directly change H,
1268 // the factorization will be sorted and placed in (schurH_, Q)
1269 computeSchurForm( true );
1270 }
1271
1272 // After the Schur form is computed, then the Ritz values are current.
1273 // Thus, I can check the Ritz index vector to see if I have enough space for the Ritz vectors requested.
1274 TEUCHOS_TEST_FOR_EXCEPTION(ritzIndex_[numRitzVecs_-1]==1, std::logic_error,
1275 "Anasazi::BlockKrylovSchur::computeRitzVectors(): the number of required Ritz vectors splits a complex conjugate pair.");
1276
1277 Teuchos::LAPACK<int,ScalarType> lapack;
1278 Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
1279
1280 // Compute the Ritz vectors.
1281 // --> For a Hermitian problem this is simply the current basis times the first numRitzVecs_ Schur vectors
1282 //
1283 // --> For a non-Hermitian problem, this involves solving the projected eigenproblem, then
1284 // placing the product of the current basis times the first numRitzVecs_ Schur vectors times the
1285 // eigenvectors of interest into the Ritz vectors.
1286
1287 // Get a view of the current Krylov-Schur basis vectors and Schur vectors
1288 std::vector<int> curind( curDim_ );
1289 for (int i=0; i<curDim_; i++) { curind[i] = i; }
1290 Teuchos::RCP<const MV> Vtemp = MVT::CloneView( *V_, curind );
1291 if (problem_->isHermitian()) {
1292 // Get a view into the current Schur vectors
1293 Teuchos::SerialDenseMatrix<int,ScalarType> subQ( Teuchos::View, *Q_, curDim_, numRitzVecs_ );
1294
1295 // Compute the current Ritz vectors
1296 MVT::MvTimesMatAddMv( ST_ONE, *Vtemp, subQ, ST_ZERO, *ritzVectors_ );
1297
1298 // Double check that no complex Ritz values have snuck into the set of converged nev.
1299 bool complexRitzVals = false;
1300 for (int i=0; i<numRitzVecs_; i++) {
1301 if (ritzIndex_[i] != 0) {
1302 complexRitzVals = true;
1303 break;
1304 }
1305 }
1306 if (complexRitzVals)
1307 om_->stream(Warnings) << " Eigenproblem is Hermitian and complex eigenvalues have converged, corresponding eigenvectors will be incorrect!!!"
1308 << std::endl;
1309 } else {
1310
1311 // Get a view into the current Schur vectors.
1312 Teuchos::SerialDenseMatrix<int,ScalarType> subQ( Teuchos::View, *Q_, curDim_, curDim_ );
1313
1314 // Get a set of work vectors to hold the current Ritz vectors.
1315 Teuchos::RCP<MV> tmpritzVectors_ = MVT::Clone( *V_, curDim_ );
1316
1317 // Compute the current Krylov-Schur vectors.
1318 MVT::MvTimesMatAddMv( ST_ONE, *Vtemp, subQ, ST_ZERO, *tmpritzVectors_ );
1319
1320 // Now compute the eigenvectors of the Schur form
1321 // Reset the dense matrix and compute the eigenvalues of the Schur form.
1322 //
1323 // Allocate the work space. This space will be used below for calls to:
1324 // * TREVC (requires 3*N for real, 2*N for complex)
1325
1326 int lwork = 3*curDim_;
1327 std::vector<ScalarType> work( lwork );
1328 std::vector<MagnitudeType> rwork( curDim_ );
1329 char side = 'R';
1330 int mm, info = 0;
1331 const int ldvl = 1;
1332 ScalarType vl[ ldvl ];
1333 Teuchos::SerialDenseMatrix<int,ScalarType> copyQ( Teuchos::Copy, *Q_, curDim_, curDim_ );
1334 lapack.TREVC( side, curDim_, schurH_->values(), schurH_->stride(), vl, ldvl,
1335 copyQ.values(), copyQ.stride(), curDim_, &mm, &work[0], &rwork[0], &info );
1336 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
1337 "Anasazi::BlockKrylovSchur::computeRitzVectors(): TREVC(n==" << curDim_ << ") returned info " << info << " != 0.");
1338
1339 // Get a view into the eigenvectors of the Schur form
1340 Teuchos::SerialDenseMatrix<int,ScalarType> subCopyQ( Teuchos::View, copyQ, curDim_, numRitzVecs_ );
1341
1342 // Convert back to Ritz vectors of the operator.
1343 curind.resize( numRitzVecs_ ); // This is already initialized above
1344 Teuchos::RCP<MV> view_ritzVectors = MVT::CloneViewNonConst( *ritzVectors_, curind );
1345 MVT::MvTimesMatAddMv( ST_ONE, *tmpritzVectors_, subCopyQ, ST_ZERO, *view_ritzVectors );
1346
1347 // Compute the norm of the new Ritz vectors
1348 std::vector<MagnitudeType> ritzNrm( numRitzVecs_ );
1349 MVT::MvNorm( *view_ritzVectors, ritzNrm );
1350
1351 // Release memory used to compute Ritz vectors before scaling the current vectors.
1352 tmpritzVectors_ = Teuchos::null;
1353 view_ritzVectors = Teuchos::null;
1354
1355 // Scale the Ritz vectors to have Euclidean norm.
1356 ScalarType ritzScale = ST_ONE;
1357 for (int i=0; i<numRitzVecs_; i++) {
1358
1359 // If this is a conjugate pair then normalize by the real and imaginary parts.
1360 if (ritzIndex_[i] == 1 ) {
1361 ritzScale = ST_ONE/lapack_mag.LAPY2(ritzNrm[i],ritzNrm[i+1]);
1362 std::vector<int> newind(2);
1363 newind[0] = i; newind[1] = i+1;
1364 tmpritzVectors_ = MVT::CloneCopy( *ritzVectors_, newind );
1365 view_ritzVectors = MVT::CloneViewNonConst( *ritzVectors_, newind );
1366 MVT::MvAddMv( ritzScale, *tmpritzVectors_, ST_ZERO, *tmpritzVectors_, *view_ritzVectors );
1367
1368 // Increment counter for imaginary part
1369 i++;
1370 } else {
1371
1372 // This is a real Ritz value, normalize the vector
1373 std::vector<int> newind(1);
1374 newind[0] = i;
1375 tmpritzVectors_ = MVT::CloneCopy( *ritzVectors_, newind );
1376 view_ritzVectors = MVT::CloneViewNonConst( *ritzVectors_, newind );
1377 MVT::MvAddMv( ST_ONE/ritzNrm[i], *tmpritzVectors_, ST_ZERO, *tmpritzVectors_, *view_ritzVectors );
1378 }
1379 }
1380
1381 } // if (problem_->isHermitian())
1382
1383 // The current Ritz vectors have been computed.
1384 ritzVecsCurrent_ = true;
1385
1386 } // if (!ritzVecsCurrent_)
1387 } // if (curDim_)
1388 } // computeRitzVectors()
1389
1390
1392 // Set a new StatusTest for the solver.
1393 template <class ScalarType, class MV, class OP>
1394 void BlockKrylovSchur<ScalarType,MV,OP>::setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test) {
1395 TEUCHOS_TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument,
1396 "Anasazi::BlockKrylovSchur::setStatusTest() was passed a null StatusTest.");
1397 tester_ = test;
1398 }
1399
1401 // Get the current StatusTest used by the solver.
1402 template <class ScalarType, class MV, class OP>
1403 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > BlockKrylovSchur<ScalarType,MV,OP>::getStatusTest() const {
1404 return tester_;
1405 }
1406
1407
1409 /* Get the current approximate eigenvalues, i.e. Ritz values.
1410 *
1411 * POST-CONDITIONS:
1412 *
1413 * ritzValues_ contains Ritz w.r.t. V, H
1414 * Q_ contains the Schur vectors w.r.t. H
1415 * schurH_ contains the Schur matrix w.r.t. H
1416 * ritzOrder_ contains the current ordering from sort manager
1417 * schurCurrent_ = true if sort = true; i.e. the Schur form is sorted according to the index
1418 * vector returned by the sort manager.
1419 */
1420 template <class ScalarType, class MV, class OP>
1422 {
1423 // local timer
1424#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1425 Teuchos::TimeMonitor LocalTimer(*timerCompSF_);
1426#endif
1427
1428 // Check to see if the dimension of the factorization is greater than zero.
1429 if (curDim_) {
1430
1431 // Check to see if the Schur factorization is current.
1432 if (!schurCurrent_) {
1433
1434 // Check to see if the Ritz values are current
1435 // --> If they are then the Schur factorization is current but not sorted.
1436 if (!ritzValsCurrent_) {
1437 Teuchos::LAPACK<int,ScalarType> lapack;
1438 Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
1439 Teuchos::BLAS<int,ScalarType> blas;
1440 Teuchos::BLAS<int,MagnitudeType> blas_mag;
1441
1442 // Get a view into Q, the storage for H's Schur vectors.
1443 Teuchos::SerialDenseMatrix<int,ScalarType> subQ( Teuchos::View, *Q_, curDim_, curDim_ );
1444
1445 // Get a copy of H to compute/sort the Schur form.
1446 schurH_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *H_, curDim_, curDim_ ) );
1447 //
1448 //---------------------------------------------------
1449 // Compute the Schur factorization of subH
1450 // ---> Use driver GEES to first reduce to upper Hessenberg
1451 // form and then compute Schur form, outputting Ritz values
1452 //---------------------------------------------------
1453 //
1454 // Allocate the work space. This space will be used below for calls to:
1455 // * GEES (requires 3*N for real, 2*N for complex)
1456 // * TREVC (requires 3*N for real, 2*N for complex)
1457 // * TREXC (requires N for real, none for complex)
1458 // Furthermore, GEES requires a real array of length curDim_ (for complex datatypes)
1459 //
1460 int lwork = 3*curDim_;
1461 std::vector<ScalarType> work( lwork );
1462 std::vector<MagnitudeType> rwork( curDim_ );
1463 std::vector<MagnitudeType> tmp_rRitzValues( curDim_ );
1464 std::vector<MagnitudeType> tmp_iRitzValues( curDim_ );
1465 std::vector<int> bwork( curDim_ );
1466 int info = 0, sdim = 0;
1467 char jobvs = 'V';
1468 lapack.GEES( jobvs,curDim_, schurH_->values(), schurH_->stride(), &sdim, &tmp_rRitzValues[0],
1469 &tmp_iRitzValues[0], subQ.values(), subQ.stride(), &work[0], lwork,
1470 &rwork[0], &bwork[0], &info );
1471
1472 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
1473 "Anasazi::BlockKrylovSchur::computeSchurForm(): GEES(n==" << curDim_ << ") returned info " << info << " != 0.");
1474
1475 // Check if imaginary part is detected by the Schur factorization of subH for Hermitian eigenproblems
1476 // NOTE: Because of full orthogonalization, there will be small entries above the block tridiagonal in the block Hessenberg matrix.
1477 // The spectrum of this matrix may include imaginary eigenvalues with small imaginary part, which will mess up the Schur
1478 // form sorting below.
1479 bool hermImagDetected = false;
1480 if (problem_->isHermitian()) {
1481 for (int i=0; i<curDim_; i++)
1482 {
1483 if (tmp_iRitzValues[i] != MT_ZERO)
1484 {
1485 hermImagDetected = true;
1486 break;
1487 }
1488 }
1489 if (hermImagDetected) {
1490 // Warn the user that complex eigenvalues have been detected.
1491 om_->stream(Warnings) << " Eigenproblem is Hermitian and complex eigenvalues have been detected!!!" << std::endl;
1492 // Compute || H - H' || to indicate how bad the symmetry is in the projected eigenproblem
1493 Teuchos::SerialDenseMatrix<int,ScalarType> localH( Teuchos::View, *H_, curDim_, curDim_ );
1494 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > tLocalH;
1495 if (Teuchos::ScalarTraits<ScalarType>::isComplex)
1496 tLocalH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( localH, Teuchos::CONJ_TRANS ) );
1497 else
1498 tLocalH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( localH, Teuchos::TRANS ) );
1499 (*tLocalH) -= localH;
1500 MagnitudeType normF = tLocalH->normFrobenius();
1501 MagnitudeType norm1 = tLocalH->normOne();
1502 om_->stream(Warnings) << " Symmetry error in projected eigenproblem: || S - S' ||_F = " << normF
1503 << ", || S - S' ||_1 = " << norm1 << std::endl;
1504 }
1505 }
1506 //
1507 //---------------------------------------------------
1508 // Use the Krylov-Schur factorization to compute the current Ritz residuals
1509 // for ALL the eigenvalues estimates (Ritz values)
1510 // || Ax - x\theta || = || U_m+1*B_m+1^H*Q*s ||
1511 // = || B_m+1^H*Q*s ||
1512 //
1513 // where U_m+1 is the current Krylov-Schur basis, Q are the Schur vectors, and x = U_m+1*Q*s
1514 // NOTE: This means that s = e_i if the problem is hermitian, else the eigenvectors
1515 // of the Schur form need to be computed.
1516 //
1517 // First compute H_{m+1,m}*B_m^T, then determine what 's' is.
1518 //---------------------------------------------------
1519 //
1520 // Get current B_m+1
1521 Teuchos::SerialDenseMatrix<int,ScalarType> curB(Teuchos::View, *H_,
1522 blockSize_, curDim_, curDim_ );
1523 //
1524 // Compute B_m+1^H*Q
1525 Teuchos::SerialDenseMatrix<int,ScalarType> subB( blockSize_, curDim_ );
1526 blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, blockSize_, curDim_, curDim_, ST_ONE,
1527 curB.values(), curB.stride(), subQ.values(), subQ.stride(),
1528 ST_ZERO, subB.values(), subB.stride() );
1529 //
1530 // Determine what 's' is and compute Ritz residuals.
1531 //
1532 ScalarType* b_ptr = subB.values();
1533 if (problem_->isHermitian() && !hermImagDetected) {
1534 //
1535 // 's' is the i-th canonical basis vector.
1536 //
1537 for (int i=0; i<curDim_ ; i++) {
1538 ritzResiduals_[i] = blas.NRM2(blockSize_, b_ptr + i*blockSize_, 1);
1539 }
1540 } else {
1541 //
1542 // Compute S: the eigenvectors of the block upper triangular, Schur matrix.
1543 //
1544 char side = 'R';
1545 int mm;
1546 const int ldvl = 1;
1547 ScalarType vl[ ldvl ];
1548 Teuchos::SerialDenseMatrix<int,ScalarType> S( curDim_, curDim_ );
1549 lapack.TREVC( side, curDim_, schurH_->values(), schurH_->stride(), vl, ldvl,
1550 S.values(), S.stride(), curDim_, &mm, &work[0], &rwork[0], &info );
1551
1552 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
1553 "Anasazi::BlockKrylovSchur::computeSchurForm(): TREVC(n==" << curDim_ << ") returned info " << info << " != 0.");
1554 //
1555 // Scale the eigenvectors so that their Euclidean norms are all one.
1556 //
1557 HelperTraits<ScalarType>::scaleRitzVectors( tmp_iRitzValues, &S );
1558 //
1559 // Compute ritzRes = *B_m+1^H*Q*S where the i-th column of S is 's' for the i-th Ritz-value
1560 //
1561 Teuchos::SerialDenseMatrix<int,ScalarType> ritzRes( blockSize_, curDim_ );
1562 blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, blockSize_, curDim_, curDim_, ST_ONE,
1563 subB.values(), subB.stride(), S.values(), S.stride(),
1564 ST_ZERO, ritzRes.values(), ritzRes.stride() );
1565
1566 /* TO DO: There's be an incorrect assumption made in the computation of the Ritz residuals.
1567 This assumption is that the next vector in the Krylov subspace is Euclidean orthonormal.
1568 It may not be normalized using Euclidean norm.
1569 Teuchos::RCP<MV> ritzResVecs = MVT::Clone( *V_, curDim_ );
1570 std::vector<int> curind(blockSize_);
1571 for (int i=0; i<blockSize_; i++) { curind[i] = curDim_ + i; }
1572 Teuchos::RCP<MV> Vtemp = MVT::CloneView(*V_,curind);
1573
1574 MVT::MvTimesMatAddMv( ST_ONE, *Vtemp, ritzRes, ST_ZERO, *ritzResVecs );
1575 std::vector<MagnitudeType> ritzResNrms(curDim_);
1576 MVT::MvNorm( *ritzResVecs, ritzResNrms );
1577 i = 0;
1578 while( i < curDim_ ) {
1579 if ( tmp_ritzValues[curDim_+i] != MT_ZERO ) {
1580 ritzResiduals_[i] = lapack_mag.LAPY2( ritzResNrms[i], ritzResNrms[i+1] );
1581 ritzResiduals_[i+1] = ritzResiduals_[i];
1582 i = i+2;
1583 } else {
1584 ritzResiduals_[i] = ritzResNrms[i];
1585 i++;
1586 }
1587 }
1588 */
1589 //
1590 // Compute the Ritz residuals for each Ritz value.
1591 //
1592 HelperTraits<ScalarType>::computeRitzResiduals( tmp_iRitzValues, ritzRes, &ritzResiduals_ );
1593 }
1594 //
1595 // Sort the Ritz values.
1596 //
1597 {
1598#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1599 Teuchos::TimeMonitor LocalTimer2(*timerSortRitzVal_);
1600#endif
1601 int i=0;
1602 if (problem_->isHermitian() && !hermImagDetected) {
1603 //
1604 // Sort using just the real part of the Ritz values.
1605 sm_->sort(tmp_rRitzValues, Teuchos::rcpFromRef(ritzOrder_), curDim_); // don't catch exception
1606 ritzIndex_.clear();
1607 while ( i < curDim_ ) {
1608 // The Ritz value is not complex.
1609 ritzValues_[i].set(tmp_rRitzValues[i], MT_ZERO);
1610 ritzIndex_.push_back(0);
1611 i++;
1612 }
1613 }
1614 else {
1615 //
1616 // Sort using both the real and imaginary parts of the Ritz values.
1617 sm_->sort(tmp_rRitzValues, tmp_iRitzValues, Teuchos::rcpFromRef(ritzOrder_) , curDim_);
1618 HelperTraits<ScalarType>::sortRitzValues( tmp_rRitzValues, tmp_iRitzValues, &ritzValues_, &ritzOrder_, &ritzIndex_ );
1619 }
1620 //
1621 // Sort the ritzResiduals_ based on the ordering from the Sort Manager.
1622 std::vector<MagnitudeType> ritz2( curDim_ );
1623 for (i=0; i<curDim_; i++) { ritz2[i] = ritzResiduals_[ ritzOrder_[i] ]; }
1624 blas_mag.COPY( curDim_, &ritz2[0], 1, &ritzResiduals_[0], 1 );
1625
1626 // The Ritz values have now been updated.
1627 ritzValsCurrent_ = true;
1628 }
1629
1630 } // if (!ritzValsCurrent_) ...
1631 //
1632 //---------------------------------------------------
1633 // * The Ritz values and residuals have been updated at this point.
1634 //
1635 // * The Schur factorization of the projected matrix has been computed,
1636 // and is stored in (schurH_, Q_).
1637 //
1638 // Now the Schur factorization needs to be sorted.
1639 //---------------------------------------------------
1640 //
1641 // Sort the Schur form using the ordering from the Sort Manager.
1642 if (sort) {
1643 sortSchurForm( *schurH_, *Q_, ritzOrder_ );
1644 //
1645 // Indicate the Schur form in (schurH_, Q_) is current and sorted
1646 schurCurrent_ = true;
1647 }
1648 } // if (!schurCurrent_) ...
1649
1650 } // if (curDim_) ...
1651
1652 } // computeSchurForm( ... )
1653
1654
1656 // Sort the Schur form of H stored in (H,Q) using the ordering vector.
1657 template <class ScalarType, class MV, class OP>
1658 void BlockKrylovSchur<ScalarType,MV,OP>::sortSchurForm( Teuchos::SerialDenseMatrix<int,ScalarType>& H,
1659 Teuchos::SerialDenseMatrix<int,ScalarType>& Q,
1660 std::vector<int>& order )
1661 {
1662 // local timer
1663#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1664 Teuchos::TimeMonitor LocalTimer(*timerSortSF_);
1665#endif
1666 //
1667 //---------------------------------------------------
1668 // Reorder real Schur factorization, remember to add one to the indices for the
1669 // fortran call and determine offset. The offset is necessary since the TREXC
1670 // method reorders in a nonsymmetric fashion, thus we use the reordering in
1671 // a stack-like fashion. Also take into account conjugate pairs, which may mess
1672 // up the reordering, since the pair is moved if one of the pair is moved.
1673 //---------------------------------------------------
1674 //
1675 int i = 0, nevtemp = 0;
1676 char compq = 'V';
1677 std::vector<int> offset2( curDim_ );
1678 std::vector<int> order2( curDim_ );
1679
1680 // LAPACK objects.
1681 Teuchos::LAPACK<int,ScalarType> lapack;
1682 int lwork = 3*curDim_;
1683 std::vector<ScalarType> work( lwork );
1684
1685 while (i < curDim_) {
1686 if ( ritzIndex_[i] != 0 ) { // This is the first value of a complex conjugate pair
1687 offset2[nevtemp] = 0;
1688 for (int j=i; j<curDim_; j++) {
1689 if (order[j] > order[i]) { offset2[nevtemp]++; }
1690 }
1691 order2[nevtemp] = order[i];
1692 i = i+2;
1693 } else {
1694 offset2[nevtemp] = 0;
1695 for (int j=i; j<curDim_; j++) {
1696 if (order[j] > order[i]) { offset2[nevtemp]++; }
1697 }
1698 order2[nevtemp] = order[i];
1699 i++;
1700 }
1701 nevtemp++;
1702 }
1703 ScalarType *ptr_h = H.values();
1704 ScalarType *ptr_q = Q.values();
1705 int ldh = H.stride(), ldq = Q.stride();
1706 int info = 0;
1707 for (i=nevtemp-1; i>=0; i--) {
1708 int ifst = order2[i]+1+offset2[i];
1709 int ilst = 1;
1710 lapack.TREXC( compq, curDim_, ptr_h, ldh, ptr_q, ldq, &ifst,
1711 &ilst, &work[0], &info );
1712 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
1713 "Anasazi::BlockKrylovSchur::computeSchurForm(): TREXC(n==" << curDim_ << ") returned info " << info << " != 0.");
1714 }
1715 }
1716
1718 // Print the current status of the solver
1719 template <class ScalarType, class MV, class OP>
1721 {
1722 using std::endl;
1723
1724 os.setf(std::ios::scientific, std::ios::floatfield);
1725 os.precision(6);
1726 os <<"================================================================================" << endl;
1727 os << endl;
1728 os <<" BlockKrylovSchur Solver Status" << endl;
1729 os << endl;
1730 os <<"The solver is "<<(initialized_ ? "initialized." : "not initialized.") << endl;
1731 os <<"The number of iterations performed is " <<iter_<<endl;
1732 os <<"The block size is " << blockSize_<<endl;
1733 os <<"The number of blocks is " << numBlocks_<<endl;
1734 os <<"The current basis size is " << curDim_<<endl;
1735 os <<"The number of auxiliary vectors is " << numAuxVecs_ << endl;
1736 os <<"The number of operations Op*x is "<<count_ApplyOp_<<endl;
1737
1738 os.setf(std::ios_base::right, std::ios_base::adjustfield);
1739
1740 os << endl;
1741 if (initialized_) {
1742 os <<"CURRENT RITZ VALUES "<<endl;
1743 if (ritzIndex_.size() != 0) {
1744 int numPrint = (curDim_ < numRitzPrint_? curDim_: numRitzPrint_);
1745 if (problem_->isHermitian()) {
1746 os << std::setw(20) << "Ritz Value"
1747 << std::setw(20) << "Ritz Residual"
1748 << endl;
1749 os <<"--------------------------------------------------------------------------------"<<endl;
1750 for (int i=0; i<numPrint; i++) {
1751 os << std::setw(20) << ritzValues_[i].realpart
1752 << std::setw(20) << ritzResiduals_[i]
1753 << endl;
1754 }
1755 } else {
1756 os << std::setw(24) << "Ritz Value"
1757 << std::setw(30) << "Ritz Residual"
1758 << endl;
1759 os <<"--------------------------------------------------------------------------------"<<endl;
1760 for (int i=0; i<numPrint; i++) {
1761 // Print out the real eigenvalue.
1762 os << std::setw(15) << ritzValues_[i].realpart;
1763 if (ritzValues_[i].imagpart < MT_ZERO) {
1764 os << " - i" << std::setw(15) << Teuchos::ScalarTraits<MagnitudeType>::magnitude(ritzValues_[i].imagpart);
1765 } else {
1766 os << " + i" << std::setw(15) << ritzValues_[i].imagpart;
1767 }
1768 os << std::setw(20) << ritzResiduals_[i] << endl;
1769 }
1770 }
1771 } else {
1772 os << std::setw(20) << "[ NONE COMPUTED ]" << endl;
1773 }
1774 }
1775 os << endl;
1776 os <<"================================================================================" << endl;
1777 os << endl;
1778 }
1779
1780} // End of namespace Anasazi
1781
1782#endif
1783
1784// End of file AnasaziBlockKrylovSchur.hpp
Pure virtual base class which describes the basic interface to the iterative eigensolver.
Declaration of basic traits for the multivector type.
Virtual base class which defines basic traits for the operator type.
Templated virtual class for providing orthogonalization/orthonormalization methods.
Types and exceptions used within Anasazi solvers and interfaces.
An exception class parent to all Anasazi exceptions.
BlockKrylovSchurInitFailure is thrown when the BlockKrylovSchur solver is unable to generate an initi...
BlockKrylovSchurOrthoFailure is thrown when the orthogonalization manager is unable to generate ortho...
This class implements the block Krylov-Schur iteration, for solving linear eigenvalue problems.
void resetNumIters()
Reset the iteration count.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the eigenvalue problem.
void computeSchurForm(const bool sort=true)
Compute the Schur form of the projected eigenproblem from the current Krylov factorization.
void setNumRitzVectors(int numRitzVecs)
Set the number of Ritz vectors to compute.
BlockKrylovSchur(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< OrthoManager< ScalarType, MV > > &ortho, Teuchos::ParameterList &params)
BlockKrylovSchur constructor with eigenproblem, solver utilities, and parameter list of solver option...
BlockKrylovSchurState< ScalarType, MV > getState() const
Get the current state of the eigensolver.
void setAuxVecs(const Teuchos::Array< Teuchos::RCP< const MV > > &auxvecs)
Set the auxiliary vectors for the solver.
bool isInitialized() const
Indicates whether the solver has been initialized or not.
Teuchos::RCP< const MV > getRitzVectors()
Get the Ritz vectors.
void setStepSize(int stepSize)
Set the step size.
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
bool isSchurCurrent() const
Get the status of the Schur form currently stored in the eigensolver.
void iterate()
This method performs Block Krylov-Schur iterations until the status test indicates the need to stop o...
void computeRitzValues()
Compute the Ritz values using the current Krylov factorization.
Teuchos::Array< Teuchos::RCP< const MV > > getAuxVecs() const
Get the auxiliary vectors for the solver.
void setStatusTest(Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test)
Set a new StatusTest for the solver.
void setSize(int blockSize, int numBlocks)
Set the blocksize and number of blocks to be used by the iterative solver in solving this eigenproble...
void initialize()
Initialize the solver with the initial vectors from the eigenproblem or random data.
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.
bool isRitzVecsCurrent() const
Get the status of the Ritz vectors currently stored in the eigensolver.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRes2Norms()
Get the current residual 2-norms.
bool isRitzValsCurrent() const
Get the status of the Ritz values currently stored in the eigensolver.
int getNumIters() const
Get the current iteration count.
std::vector< int > getRitzIndex()
Get the Ritz index vector.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRitzRes2Norms()
Get the current Ritz residual 2-norms.
int getStepSize() const
Get the step size.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getResNorms()
Get the current residual norms.
std::vector< Value< ScalarType > > getRitzValues()
Get the Ritz values.
void computeRitzVectors()
Compute the Ritz vectors using the current Krylov factorization.
int getNumRitzVectors() const
Get the number of Ritz vectors to compute.
void currentStatus(std::ostream &os)
This method requests that the solver print out its current status to screen.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this eigenproblem.
virtual ~BlockKrylovSchur()
BlockKrylovSchur destructor.
The Eigensolver is a templated virtual base class that defines the basic interface that any eigensolv...
static void computeRitzResiduals(const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &iRV, const Teuchos::SerialDenseMatrix< int, ScalarType > &S, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > *RR)
Helper function for correctly computing the Ritz residuals of the projected eigenproblem.
static void scaleRitzVectors(const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &iRV, Teuchos::SerialDenseMatrix< int, ScalarType > *S)
Helper function for correctly scaling the eigenvectors of the projected eigenproblem.
static void sortRitzValues(const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &rRV, const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &iRV, std::vector< Value< ScalarType > > *RV, std::vector< int > *RO, std::vector< int > *RI)
Helper function for correctly storing the Ritz values when the eigenproblem is non-Hermitian.
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 BlockKrylovSchur state variables.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > S
The current Schur form reduction of the valid part of H.
int curDim
The current dimension of the reduction.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > Q
The current Schur vectors of the valid part of H.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
Teuchos::RCP< const MulVec > V
The current Krylov basis.