Anasazi Version of the Day
Loading...
Searching...
No Matches
AnasaziTraceMinBase.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Anasazi: Block Eigensolvers Package
5// Copyright 2004 Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40// @HEADER
41
42// TODO: Modify code so R is not computed unless needed
43
48#ifndef ANASAZI_TRACEMIN_BASE_HPP
49#define ANASAZI_TRACEMIN_BASE_HPP
50
51#include "AnasaziBasicSort.hpp"
52#include "AnasaziConfigDefs.hpp"
62#include "AnasaziTypes.hpp"
63
64#include "Teuchos_ParameterList.hpp"
65#include "Teuchos_ScalarTraits.hpp"
66#include "Teuchos_SerialDenseMatrix.hpp"
67#include "Teuchos_SerialDenseSolver.hpp"
68#include "Teuchos_TimeMonitor.hpp"
69
70using Teuchos::RCP;
71using Teuchos::rcp;
72
73namespace Anasazi {
84namespace Experimental {
85
87
88
93 template <class ScalarType, class MV>
96 int curDim;
101 RCP<const MV> V;
103 RCP<const MV> KV;
105 RCP<const MV> MopV;
107 RCP<const MV> X;
109 RCP<const MV> KX;
111 RCP<const MV> MX;
113 RCP<const MV> R;
115 RCP<const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > T;
121 RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > KK;
123 RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > RV;
127 int NEV;
131 RCP< const std::vector<ScalarType> > ritzShifts;
132 TraceMinBaseState() : curDim(0), V(Teuchos::null), KV(Teuchos::null), MopV(Teuchos::null),
133 X(Teuchos::null), KX(Teuchos::null), MX(Teuchos::null), R(Teuchos::null),
134 T(Teuchos::null), KK(Teuchos::null), RV(Teuchos::null), isOrtho(false),
135 NEV(0), largestSafeShift(Teuchos::ScalarTraits<ScalarType>::zero()),
136 ritzShifts(Teuchos::null) {}
137 };
138
140
142
143
156 class TraceMinBaseInitFailure : public AnasaziError {public:
157 TraceMinBaseInitFailure(const std::string& what_arg) : AnasaziError(what_arg)
158 {}};
159
164 TraceMinBaseOrthoFailure(const std::string& what_arg) : AnasaziError(what_arg)
165 {}};
166
168
181 template <class ScalarType, class MV, class OP>
182 class TraceMinBase : public Eigensolver<ScalarType,MV,OP> {
183 public:
185
186
214 TraceMinBase( const RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
215 const RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
216 const RCP<OutputManager<ScalarType> > &printer,
217 const RCP<StatusTest<ScalarType,MV,OP> > &tester,
218 const RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
219 Teuchos::ParameterList &params
220 );
221
223 virtual ~TraceMinBase();
225
226
228
229
253 void iterate();
254
255 void harmonicIterate();
256
279 void initialize(TraceMinBaseState<ScalarType,MV>& newstate);
280
281 void harmonicInitialize(TraceMinBaseState<ScalarType,MV> newstate);
282
286 void initialize();
287
303 bool isInitialized() const;
304
313 TraceMinBaseState<ScalarType,MV> getState() const;
314
316
317
319
320
322 int getNumIters() const;
323
325 void resetNumIters();
326
334 RCP<const MV> getRitzVectors();
335
341 std::vector<Value<ScalarType> > getRitzValues();
342
343
352 std::vector<int> getRitzIndex();
353
354
360 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getResNorms();
361
362
368 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRes2Norms();
369
370
378 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRitzRes2Norms();
379
385 int getCurSubspaceDim() const;
386
388 int getMaxSubspaceDim() const;
389
391
392
394
395
397 void setStatusTest(RCP<StatusTest<ScalarType,MV,OP> > test);
398
400 RCP<StatusTest<ScalarType,MV,OP> > getStatusTest() const;
401
403 const Eigenproblem<ScalarType,MV,OP>& getProblem() const;
404
414 void setBlockSize(int blockSize);
415
417 int getBlockSize() const;
418
436 void setAuxVecs(const Teuchos::Array<RCP<const MV> > &auxvecs);
437
439 Teuchos::Array<RCP<const MV> > getAuxVecs() const;
440
442
444
445
452 void setSize(int blockSize, int numBlocks);
453
455
456
458 void currentStatus(std::ostream &os);
459
461
462 protected:
463 //
464 // Convenience typedefs
465 //
466 typedef SolverUtils<ScalarType,MV,OP> Utils;
467 typedef MultiVecTraits<ScalarType,MV> MVT;
468 typedef OperatorTraits<ScalarType,MV,OP> OPT;
469 typedef Teuchos::ScalarTraits<ScalarType> SCT;
470 typedef typename SCT::magnitudeType MagnitudeType;
471 typedef TraceMinRitzOp<ScalarType,MV,OP> tracemin_ritz_op_type;
472 typedef SaddleContainer<ScalarType,MV> saddle_container_type;
473 typedef SaddleOperator<ScalarType,MV,tracemin_ritz_op_type> saddle_op_type;
474 const MagnitudeType ONE;
475 const MagnitudeType ZERO;
476 const MagnitudeType NANVAL;
477 //
478 // Classes inputed through constructor that define the eigenproblem to be solved.
479 //
480 const RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
481 const RCP<SortManager<MagnitudeType> > sm_;
482 const RCP<OutputManager<ScalarType> > om_;
483 RCP<StatusTest<ScalarType,MV,OP> > tester_;
484 const RCP<MatOrthoManager<ScalarType,MV,OP> > orthman_;
485 //
486 // Information obtained from the eigenproblem
487 //
488 RCP<const OP> Op_;
489 RCP<const OP> MOp_;
490 RCP<const OP> Prec_;
491 bool hasM_;
492 //
493 // Internal timers
494 // TODO: Fix the timers
495 //
496 RCP<Teuchos::Time> timerOp_, timerMOp_, timerSaddle_, timerSortEval_, timerDS_,
497 timerLocal_, timerCompRes_, timerOrtho_, timerInit_;
498 //
499 // Internal structs
500 // TODO: Fix the checks
501 //
502 struct CheckList {
503 bool checkV, checkX, checkMX,
504 checkKX, checkQ, checkKK;
505 CheckList() : checkV(false),checkX(false),
506 checkMX(false),checkKX(false),
507 checkQ(false),checkKK(false) {};
508 };
509 //
510 // Internal methods
511 //
512 std::string accuracyCheck(const CheckList &chk, const std::string &where) const;
513
514 //
515 // Counters
516 //
517 int count_ApplyOp_, count_ApplyM_;
518
519 //
520 // Algorithmic parameters.
521 //
522 // blockSize_ is the solver block size; it controls the number of vectors added to the basis on each iteration.
523 int blockSize_;
524 // numBlocks_ is the size of the allocated space for the basis, in blocks.
525 int numBlocks_;
526
527 //
528 // Current solver state
529 //
530 // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
531 // is capable of running; _initialize is controlled by the initialize() member method
532 // For the implications of the state of initialized_, please see documentation for initialize()
533 bool initialized_;
534 //
535 // curDim_ reflects how much of the current basis is valid
536 // NOTE: 0 <= curDim_ <= blockSize_*numBlocks_
537 // this also tells us how many of the values in theta_ are valid Ritz values
538 int curDim_;
539 //
540 // State Multivecs
541 // V is the current basis
542 // X is the current set of eigenvectors
543 // R is the residual
544 RCP<MV> X_, KX_, MX_, KV_, MV_, R_, V_;
545 //
546 // Projected matrices
547 //
548 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > KK_, ritzVecs_;
549 //
550 // auxiliary vectors
551 Teuchos::Array<RCP<const MV> > auxVecs_, MauxVecs_;
552 int numAuxVecs_;
553 //
554 // Number of iterations that have been performed.
555 int iter_;
556 //
557 // Current eigenvalues, residual norms
558 std::vector<MagnitudeType> theta_, Rnorms_, R2norms_;
559 //
560 // are the residual norms current with the residual?
561 bool Rnorms_current_, R2norms_current_;
562
563 // TraceMin specific variables
564 RCP<tracemin_ritz_op_type> ritzOp_; // Operator which incorporates Ritz shifts
565 enum SaddleSolType saddleSolType_; // Type of saddle point solver being used
566 bool previouslyLeveled_; // True if the trace already leveled off
567 MagnitudeType previousTrace_; // Trace of X'KX at the previous iteration
568 bool posSafeToShift_, negSafeToShift_; // Whether there are multiple clusters
569 MagnitudeType largestSafeShift_; // The largest shift considered to be safe - generally the biggest converged eigenvalue
570 int NEV_; // Number of eigenvalues we seek - used in computation of trace
571 std::vector<ScalarType> ritzShifts_; // Current Ritz shifts
572
573 // This is only used if we're using the Schur complement solver
574 RCP<MV> Z_;
575
576 // User provided TraceMin parameters
577 enum WhenToShiftType whenToShift_; // What triggers a Ritz shift
578 enum HowToShiftType howToShift_; // Strategy for choosing the Ritz shifts
579 bool useMultipleShifts_; // Whether to use one Ritz shift or many
580 bool considerClusters_; // Don't shift if there is only one cluster
581 bool projectAllVecs_; // Use all vectors in projected minres, or just 1
582 bool projectLockedVecs_; // Project locked vectors too in minres; does nothing if projectAllVecs = false
583 bool computeAllRes_; // Compute all residuals, or just blocksize ones - helps make Ritz shifts safer
584 bool useRHSR_; // Use -R as the RHS of projected minres rather than AX
585 bool useHarmonic_;
586 MagnitudeType traceThresh_;
587 MagnitudeType alpha_;
588
589 // Variables that are only used if we're shifting when the residual gets small
590 // TODO: These are currently unused
591 std::string shiftNorm_; // Measure 2-norm or M-norm of residual
592 MagnitudeType shiftThresh_; // How small must the residual be?
593 bool useRelShiftThresh_; // Are we scaling the shift threshold by the eigenvalues?
594
595 // TraceMin specific functions
596 // Returns the trace of KK = X'KX
597 ScalarType getTrace() const;
598 // Returns true if the change in trace is very small (or was very small at one point)
599 // TODO: Check whether I want to return true if the trace WAS small
600 bool traceLeveled();
601 // Get the residuals of each cluster of eigenvalues
602 // TODO: Figure out whether I want to use these for all eigenvalues or just the first
603 std::vector<ScalarType> getClusterResids();
604 // Computes the Ritz shifts, which is not a trivial task
605 // TODO: Make it safer for indefinite matrices maybe?
606 void computeRitzShifts(const std::vector<ScalarType>& clusterResids);
607 // Compute the tolerance for the inner solves
608 // TODO: Come back to this and make sure it does what I want it to
609 std::vector<ScalarType> computeTol();
610 // Solves a saddle point problem
611 void solveSaddlePointProblem(RCP<MV> Delta);
612 // Solves a saddle point problem by using unpreconditioned projected minres
613 void solveSaddleProj(RCP<MV> Delta) const;
614 // Solves a saddle point problem by using preconditioned projected...Krylov...something
615 // TODO: Fix this. Replace minres with gmres or something.
616 void solveSaddleProjPrec(RCP<MV> Delta) const;
617 // Solves a saddle point problem by explicitly forming the inexact Schur complement
618 void solveSaddleSchur (RCP<MV> Delta) const;
619 // Solves a saddle point problem with a block diagonal preconditioner
620 void solveSaddleBDPrec (RCP<MV> Delta) const;
621 // Solves a saddle point problem with a Hermitian/non-Hermitian splitting preconditioner
622 void solveSaddleHSSPrec (RCP<MV> Delta) const;
623 // Computes KK = X'KX
624 void computeKK();
625 // Computes the eigenpairs of KK
626 void computeRitzPairs();
627 // Computes X = V evecs
628 void computeX();
629 // Updates KX and MX without a matvec by MAGIC (or by multiplying KV and MV by evecs)
630 void updateKXMX();
631 // Updates the residual
632 void updateResidual();
633 // Adds Delta to the basis
634 // TraceMin and TraceMinDavidson each do this differently, which is why this is pure virtual
635 virtual void addToBasis(const RCP<const MV> Delta) =0;
636 virtual void harmonicAddToBasis(const RCP<const MV> Delta) =0;
637 };
638
641 //
642 // Implementations
643 //
646
647
649 // Constructor
650 // TODO: Allow the users to supply their own linear solver
651 // TODO: Add additional checking for logic errors (like trying to use gmres with multiple shifts)
652 template <class ScalarType, class MV, class OP>
654 const RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
655 const RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
656 const RCP<OutputManager<ScalarType> > &printer,
657 const RCP<StatusTest<ScalarType,MV,OP> > &tester,
658 const RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
659 Teuchos::ParameterList &params
660 ) :
661 ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
662 ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
663 NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan()),
664 // problem, tools
665 problem_(problem),
666 sm_(sorter),
667 om_(printer),
668 tester_(tester),
669 orthman_(ortho),
670 // timers, counters
671#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
672 timerOp_(Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinBase::Operation Op*x")),
673 timerMOp_(Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinBase::Operation M*x")),
674 timerSaddle_(Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinBase::Solving saddle point problem")),
675 timerSortEval_(Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinBase::Sorting eigenvalues")),
676 timerDS_(Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinBase::Direct solve")),
677 timerLocal_(Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinBase::Local update")),
678 timerCompRes_(Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinBase::Computing residuals")),
679 timerOrtho_(Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinBase::Orthogonalization")),
680 timerInit_(Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinBase::Initialization")),
681#endif
682 count_ApplyOp_(0),
683 count_ApplyM_(0),
684 // internal data
685 blockSize_(0),
686 numBlocks_(0),
687 initialized_(false),
688 curDim_(0),
689 auxVecs_( Teuchos::Array<RCP<const MV> >(0) ),
690 MauxVecs_( Teuchos::Array<RCP<const MV> >(0) ),
691 numAuxVecs_(0),
692 iter_(0),
693 Rnorms_current_(false),
694 R2norms_current_(false),
695 // TraceMin variables
696 previouslyLeveled_(false),
697 previousTrace_(ZERO),
698 posSafeToShift_(false),
699 negSafeToShift_(false),
700 largestSafeShift_(ZERO),
701 Z_(Teuchos::null)
702 {
703 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument,
704 "Anasazi::TraceMinBase::constructor: user passed null problem pointer.");
705 TEUCHOS_TEST_FOR_EXCEPTION(sm_ == Teuchos::null,std::invalid_argument,
706 "Anasazi::TraceMinBase::constructor: user passed null sort manager pointer.");
707 TEUCHOS_TEST_FOR_EXCEPTION(om_ == Teuchos::null,std::invalid_argument,
708 "Anasazi::TraceMinBase::constructor: user passed null output manager pointer.");
709 TEUCHOS_TEST_FOR_EXCEPTION(tester_ == Teuchos::null,std::invalid_argument,
710 "Anasazi::TraceMinBase::constructor: user passed null status test pointer.");
711 TEUCHOS_TEST_FOR_EXCEPTION(orthman_ == Teuchos::null,std::invalid_argument,
712 "Anasazi::TraceMinBase::constructor: user passed null orthogonalization manager pointer.");
713 TEUCHOS_TEST_FOR_EXCEPTION(problem_->isHermitian() == false, std::invalid_argument,
714 "Anasazi::TraceMinBase::constructor: problem is not hermitian.");
715
716 // get the problem operators
717 Op_ = problem_->getOperator();
718 MOp_ = problem_->getM();
719 Prec_ = problem_->getPrec();
720 hasM_ = (MOp_ != Teuchos::null);
721
722 // Set the saddle point solver parameters
723 saddleSolType_ = params.get("Saddle Solver Type", PROJECTED_KRYLOV_SOLVER);
724 TEUCHOS_TEST_FOR_EXCEPTION(saddleSolType_ != PROJECTED_KRYLOV_SOLVER && saddleSolType_ != SCHUR_COMPLEMENT_SOLVER && saddleSolType_ != BD_PREC_MINRES && saddleSolType_ != HSS_PREC_GMRES, std::invalid_argument,
725 "Anasazi::TraceMin::constructor: Invalid value for \"Saddle Solver Type\"; valid options are PROJECTED_KRYLOV_SOLVER, SCHUR_COMPLEMENT_SOLVER, and BD_PREC_MINRES.");
726
727 // Set the Ritz shift parameters
728 whenToShift_ = params.get("When To Shift", ALWAYS_SHIFT);
729 TEUCHOS_TEST_FOR_EXCEPTION(whenToShift_ != NEVER_SHIFT && whenToShift_ != SHIFT_WHEN_TRACE_LEVELS && whenToShift_ != SHIFT_WHEN_RESID_SMALL && whenToShift_ != ALWAYS_SHIFT, std::invalid_argument,
730 "Anasazi::TraceMin::constructor: Invalid value for \"When To Shift\"; valid options are \"NEVER_SHIFT\", \"SHIFT_WHEN_TRACE_LEVELS\", \"SHIFT_WHEN_RESID_SMALL\", and \"ALWAYS_SHIFT\".");
731
732 traceThresh_ = params.get("Trace Threshold", 2e-2);
733 TEUCHOS_TEST_FOR_EXCEPTION(traceThresh_ < 0, std::invalid_argument,
734 "Anasazi::TraceMin::constructor: Invalid value for \"Trace Threshold\"; Must be positive.");
735
736 howToShift_ = params.get("How To Choose Shift", ADJUSTED_RITZ_SHIFT);
737 TEUCHOS_TEST_FOR_EXCEPTION(howToShift_ != LARGEST_CONVERGED_SHIFT && howToShift_ != ADJUSTED_RITZ_SHIFT && howToShift_ != RITZ_VALUES_SHIFT && howToShift_ != EXPERIMENTAL_SHIFT, std::invalid_argument,
738 "Anasazi::TraceMin::constructor: Invalid value for \"How To Choose Shift\"; valid options are \"LARGEST_CONVERGED_SHIFT\", \"ADJUSTED_RITZ_SHIFT\", \"RITZ_VALUES_SHIFT\".");
739
740 useMultipleShifts_ = params.get("Use Multiple Shifts", true);
741
742 shiftThresh_ = params.get("Shift Threshold", 1e-2);
743 useRelShiftThresh_ = params.get("Relative Shift Threshold", true);
744 shiftNorm_ = params.get("Shift Norm", "2");
745 TEUCHOS_TEST_FOR_EXCEPTION(shiftNorm_ != "2" && shiftNorm_ != "M", std::invalid_argument,
746 "Anasazi::TraceMin::constructor: Invalid value for \"Shift Norm\"; valid options are \"2\", \"M\".");
747
748 considerClusters_ = params.get("Consider Clusters", true);
749
750 projectAllVecs_ = params.get("Project All Vectors", true);
751 projectLockedVecs_ = params.get("Project Locked Vectors", true);
752 useRHSR_ = params.get("Use Residual as RHS", true);
753 useHarmonic_ = params.get("Use Harmonic Ritz Values", false);
754 computeAllRes_ = params.get("Compute All Residuals", true);
755
756 // set the block size and allocate data
757 int bs = params.get("Block Size", problem_->getNEV());
758 int nb = params.get("Num Blocks", 1);
759 setSize(bs,nb);
760
761 NEV_ = problem_->getNEV();
762
763 // Create the Ritz shift operator
764 ritzOp_ = rcp (new tracemin_ritz_op_type (Op_, MOp_, Prec_));
765
766 // Set the maximum number of inner iterations
767 const int innerMaxIts = params.get ("Maximum Krylov Iterations", 200);
768 ritzOp_->setMaxIts (innerMaxIts);
769
770 alpha_ = params.get ("HSS: alpha", ONE);
771 }
772
773
775 // Destructor
776 template <class ScalarType, class MV, class OP>
778
779
781 // Set the block size
782 // This simply calls setSize(), modifying the block size while retaining the number of blocks.
783 template <class ScalarType, class MV, class OP>
785 {
786 TEUCHOS_TEST_FOR_EXCEPTION(blockSize < 1, std::invalid_argument, "Anasazi::TraceMinBase::setSize(blocksize,numblocks): blocksize must be strictly positive.");
787 setSize(blockSize,numBlocks_);
788 }
789
790
792 // Return the current auxiliary vectors
793 template <class ScalarType, class MV, class OP>
794 Teuchos::Array<RCP<const MV> > TraceMinBase<ScalarType,MV,OP>::getAuxVecs() const {
795 return auxVecs_;
796 }
797
798
800 // return the current block size
801 template <class ScalarType, class MV, class OP>
803 return(blockSize_);
804 }
805
806
808 // return eigenproblem
809 template <class ScalarType, class MV, class OP>
810 const Eigenproblem<ScalarType,MV,OP>& TraceMinBase<ScalarType,MV,OP>::getProblem() const {
811 return(*problem_);
812 }
813
814
816 // return max subspace dim
817 template <class ScalarType, class MV, class OP>
819 return blockSize_*numBlocks_;
820 }
821
823 // return current subspace dim
824 template <class ScalarType, class MV, class OP>
826 if (!initialized_) return 0;
827 return curDim_;
828 }
829
830
832 // return ritz residual 2-norms
833 template <class ScalarType, class MV, class OP>
834 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
836 return getRes2Norms();
837 }
838
839
841 // return ritz index
842 template <class ScalarType, class MV, class OP>
844 std::vector<int> ret(curDim_,0);
845 return ret;
846 }
847
848
850 // return ritz values
851 template <class ScalarType, class MV, class OP>
852 std::vector<Value<ScalarType> > TraceMinBase<ScalarType,MV,OP>::getRitzValues() {
853 std::vector<Value<ScalarType> > ret(curDim_);
854 for (int i=0; i<curDim_; ++i) {
855 ret[i].realpart = theta_[i];
856 ret[i].imagpart = ZERO;
857 }
858 return ret;
859 }
860
861
863 // return pointer to ritz vectors
864 template <class ScalarType, class MV, class OP>
866 return X_;
867 }
868
869
871 // reset number of iterations
872 template <class ScalarType, class MV, class OP>
876
877
879 // return number of iterations
880 template <class ScalarType, class MV, class OP>
882 return(iter_);
883 }
884
885
887 // return state pointers
888 template <class ScalarType, class MV, class OP>
889 TraceMinBaseState<ScalarType,MV> TraceMinBase<ScalarType,MV,OP>::getState() const {
890 TraceMinBaseState<ScalarType,MV> state;
891 state.curDim = curDim_;
892 state.V = V_;
893 state.X = X_;
894 if (Op_ != Teuchos::null) {
895 state.KV = KV_;
896 state.KX = KX_;
897 }
898 else {
899 state.KV = Teuchos::null;
900 state.KX = Teuchos::null;
901 }
902 if (hasM_) {
903 state.MopV = MV_;
904 state.MX = MX_;
905 }
906 else {
907 state.MopV = Teuchos::null;
908 state.MX = Teuchos::null;
909 }
910 state.R = R_;
911 state.KK = KK_;
912 state.RV = ritzVecs_;
913 if (curDim_ > 0) {
914 state.T = rcp(new std::vector<MagnitudeType>(theta_.begin(),theta_.begin()+curDim_) );
915 } else {
916 state.T = rcp(new std::vector<MagnitudeType>(0));
917 }
918 state.ritzShifts = rcp(new std::vector<MagnitudeType>(ritzShifts_.begin(),ritzShifts_.begin()+blockSize_) );
919 state.isOrtho = true;
920 state.largestSafeShift = largestSafeShift_;
921
922 return state;
923 }
924
925
927 // Perform TraceMinBase iterations until the StatusTest tells us to stop.
928 template <class ScalarType, class MV, class OP>
930 {
931 if(useHarmonic_)
932 {
933 harmonicIterate();
934 return;
935 }
936
937 //
938 // Initialize solver state
939 if (initialized_ == false) {
940 initialize();
941 }
942
943 // as a data member, this would be redundant and require synchronization with
944 // blockSize_ and numBlocks_; we'll use a constant here.
945 const int searchDim = blockSize_*numBlocks_;
946
947 // Update obtained from solving saddle point problem
948 RCP<MV> Delta = MVT::Clone(*X_,blockSize_);
949
951 // iterate until the status test tells us to stop.
952 // also break if our basis is full
953 while (tester_->checkStatus(this) != Passed && (numBlocks_ == 1 || curDim_ < searchDim)) {
954
955 // Print information on current iteration
956 if (om_->isVerbosity(Debug)) {
957 currentStatus( om_->stream(Debug) );
958 }
959 else if (om_->isVerbosity(IterationDetails)) {
960 currentStatus( om_->stream(IterationDetails) );
961 }
962
963 ++iter_;
964
965 // Solve the saddle point problem
966 solveSaddlePointProblem(Delta);
967
968 // Insert Delta at the end of V
969 addToBasis(Delta);
970
971 if (om_->isVerbosity( Debug ) ) {
972 // Check almost everything here
973 CheckList chk;
974 chk.checkV = true;
975 om_->print( Debug, accuracyCheck(chk, ": after addToBasis(Delta)") );
976 }
977
978 // Compute KK := V'KV
979 computeKK();
980
981 if (om_->isVerbosity( Debug ) ) {
982 // Check almost everything here
983 CheckList chk;
984 chk.checkKK = true;
985 om_->print( Debug, accuracyCheck(chk, ": after computeKK()") );
986 }
987
988 // Compute the Ritz pairs
989 computeRitzPairs();
990
991 // Compute X := V RitzVecs
992 computeX();
993
994 if (om_->isVerbosity( Debug ) ) {
995 // Check almost everything here
996 CheckList chk;
997 chk.checkX = true;
998 om_->print( Debug, accuracyCheck(chk, ": after computeX()") );
999 }
1000
1001 // Compute KX := KV RitzVecs and MX := MV RitzVecs (if necessary)
1002 updateKXMX();
1003
1004 if (om_->isVerbosity( Debug ) ) {
1005 // Check almost everything here
1006 CheckList chk;
1007 chk.checkKX = true;
1008 chk.checkMX = true;
1009 om_->print( Debug, accuracyCheck(chk, ": after updateKXMX()") );
1010 }
1011
1012 // Update the residual vectors
1013 updateResidual();
1014 } // end while (statusTest == false)
1015
1016 } // end of iterate()
1017
1018
1019
1021 // Perform TraceMinBase iterations until the StatusTest tells us to stop.
1022 template <class ScalarType, class MV, class OP>
1024 {
1025 //
1026 // Initialize solver state
1027 if (initialized_ == false) {
1028 initialize();
1029 }
1030
1031 // as a data member, this would be redundant and require synchronization with
1032 // blockSize_ and numBlocks_; we'll use a constant here.
1033 const int searchDim = blockSize_*numBlocks_;
1034
1035 // Update obtained from solving saddle point problem
1036 RCP<MV> Delta = MVT::Clone(*X_,blockSize_);
1037
1039 // iterate until the status test tells us to stop.
1040 // also break if our basis is full
1041 while (tester_->checkStatus(this) != Passed && (numBlocks_ == 1 || curDim_ < searchDim)) {
1042
1043 // Print information on current iteration
1044 if (om_->isVerbosity(Debug)) {
1045 currentStatus( om_->stream(Debug) );
1046 }
1047 else if (om_->isVerbosity(IterationDetails)) {
1048 currentStatus( om_->stream(IterationDetails) );
1049 }
1050
1051 ++iter_;
1052
1053 // Solve the saddle point problem
1054 solveSaddlePointProblem(Delta);
1055
1056 // Insert Delta at the end of V
1057 harmonicAddToBasis(Delta);
1058
1059 if (om_->isVerbosity( Debug ) ) {
1060 // Check almost everything here
1061 CheckList chk;
1062 chk.checkV = true;
1063 om_->print( Debug, accuracyCheck(chk, ": after addToBasis(Delta)") );
1064 }
1065
1066 // Compute KK := V'KV
1067 computeKK();
1068
1069 if (om_->isVerbosity( Debug ) ) {
1070 // Check almost everything here
1071 CheckList chk;
1072 chk.checkKK = true;
1073 om_->print( Debug, accuracyCheck(chk, ": after computeKK()") );
1074 }
1075
1076 // Compute the Ritz pairs
1077 computeRitzPairs();
1078
1079 // Compute X := V RitzVecs
1080 computeX();
1081
1082 // Get norm of each vector in X
1083 int nvecs;
1084 if(computeAllRes_)
1085 nvecs = curDim_;
1086 else
1087 nvecs = blockSize_;
1088 std::vector<int> dimind(nvecs);
1089 for(int i=0; i<nvecs; i++)
1090 dimind[i] = i;
1091 RCP<MV> lclX = MVT::CloneViewNonConst(*X_,dimind);
1092 std::vector<ScalarType> normvec(nvecs);
1093 orthman_->normMat(*lclX,normvec);
1094
1095 // Scale X
1096 for(int i=0; i<nvecs; i++)
1097 normvec[i] = ONE/normvec[i];
1098 MVT::MvScale(*lclX,normvec);
1099
1100 // Scale eigenvalues
1101 for(int i=0; i<nvecs; i++)
1102 {
1103 theta_[i] = theta_[i] * normvec[i] * normvec[i];
1104 }
1105
1106 if (om_->isVerbosity( Debug ) ) {
1107 // Check almost everything here
1108 CheckList chk;
1109 chk.checkX = true;
1110 om_->print( Debug, accuracyCheck(chk, ": after computeX()") );
1111 }
1112
1113 // Compute KX := KV RitzVecs and MX := MV RitzVecs (if necessary)
1114 updateKXMX();
1115
1116 // Scale KX and MX
1117 if(Op_ != Teuchos::null)
1118 {
1119 RCP<MV> lclKX = MVT::CloneViewNonConst(*KX_,dimind);
1120 MVT::MvScale(*lclKX,normvec);
1121 }
1122 if(hasM_)
1123 {
1124 RCP<MV> lclMX = MVT::CloneViewNonConst(*MX_,dimind);
1125 MVT::MvScale(*lclMX,normvec);
1126 }
1127
1128 if (om_->isVerbosity( Debug ) ) {
1129 // Check almost everything here
1130 CheckList chk;
1131 chk.checkKX = true;
1132 chk.checkMX = true;
1133 om_->print( Debug, accuracyCheck(chk, ": after updateKXMX()") );
1134 }
1135
1136 // Update the residual vectors
1137 updateResidual();
1138 } // end while (statusTest == false)
1139
1140 } // end of harmonicIterate()
1141
1142
1144 // initialize the solver with default state
1145 template <class ScalarType, class MV, class OP>
1147 {
1148 TraceMinBaseState<ScalarType,MV> empty;
1149 initialize(empty);
1150 }
1151
1152
1154 /* Initialize the state of the solver
1155 *
1156 * POST-CONDITIONS:
1157 *
1158 * V_ is orthonormal, orthogonal to auxVecs_, for first curDim_ vectors
1159 * theta_ contains Ritz w.r.t. V_(1:curDim_)
1160 * X is Ritz vectors w.r.t. V_(1:curDim_)
1161 * KX = Op*X
1162 * MX = M*X if hasM_
1163 * R = KX - MX*diag(theta_)
1164 *
1165 */
1166 template <class ScalarType, class MV, class OP>
1167 void TraceMinBase<ScalarType,MV,OP>::initialize(TraceMinBaseState<ScalarType,MV>& newstate)
1168 {
1169 // TODO: Check for bad input
1170 // NOTE: memory has been allocated by setBlockSize(). Use setBlock below; do not Clone
1171 // NOTE: Overall time spent in this routine is counted to timerInit_; portions will also be counted towards other primitives
1172
1173#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1174 Teuchos::TimeMonitor inittimer( *timerInit_ );
1175#endif
1176
1177 previouslyLeveled_ = false;
1178
1179 if(useHarmonic_)
1180 {
1181 harmonicInitialize(newstate);
1182 return;
1183 }
1184
1185 std::vector<int> bsind(blockSize_);
1186 for (int i=0; i<blockSize_; ++i) bsind[i] = i;
1187
1188 // in TraceMinBase, V is primary
1189 // the order of dependence follows like so.
1190 // --init-> V,KK
1191 // --ritz analysis-> theta,X
1192 // --op apply-> KX,MX
1193 // --compute-> R
1194 //
1195 // if the user specifies all data for a level, we will accept it.
1196 // otherwise, we will generate the whole level, and all subsequent levels.
1197 //
1198 // the data members are ordered based on dependence, and the levels are
1199 // partitioned according to the amount of work required to produce the
1200 // items in a level.
1201 //
1202 // inconsistent multivectors widths and lengths will not be tolerated, and
1203 // will be treated with exceptions.
1204 //
1205 // for multivector pointers in newstate which point directly (as opposed to indirectly, via a view) to
1206 // multivectors in the solver, the copy will not be affected.
1207
1208 // set up V and KK: get them from newstate if user specified them
1209 // otherwise, set them manually
1210 RCP<MV> lclV, lclKV, lclMV;
1211 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclKK, lclRV;
1212
1213 // If V was supplied by the user...
1214 if (newstate.V != Teuchos::null) {
1215 om_->stream(Debug) << "Copying V from the user\n";
1216
1217 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.V) != MVT::GetGlobalLength(*V_), std::invalid_argument,
1218 "Anasazi::TraceMinBase::initialize(newstate): Vector length of V not correct." );
1219 TEUCHOS_TEST_FOR_EXCEPTION( newstate.curDim < blockSize_, std::invalid_argument,
1220 "Anasazi::TraceMinBase::initialize(newstate): Rank of new state must be at least blockSize().");
1221 TEUCHOS_TEST_FOR_EXCEPTION( newstate.curDim > blockSize_*numBlocks_, std::invalid_argument,
1222 "Anasazi::TraceMinBase::initialize(newstate): Rank of new state must be less than getMaxSubspaceDim().");
1223
1224 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < newstate.curDim, std::invalid_argument,
1225 "Anasazi::TraceMinBase::initialize(newstate): Multivector for basis in new state must be as large as specified state rank.");
1226
1227 curDim_ = newstate.curDim;
1228 // pick an integral amount
1229 curDim_ = (int)(curDim_ / blockSize_)*blockSize_;
1230
1231 TEUCHOS_TEST_FOR_EXCEPTION( curDim_ != newstate.curDim, std::invalid_argument,
1232 "Anasazi::TraceMinBase::initialize(newstate): Rank of new state must be a multiple of getBlockSize().");
1233
1234 // put data in V
1235 std::vector<int> nevind(curDim_);
1236 for (int i=0; i<curDim_; ++i) nevind[i] = i;
1237 if (newstate.V != V_) {
1238 if (curDim_ < MVT::GetNumberVecs(*newstate.V)) {
1239 newstate.V = MVT::CloneView(*newstate.V,nevind);
1240 }
1241 MVT::SetBlock(*newstate.V,nevind,*V_);
1242 }
1243 lclV = MVT::CloneViewNonConst(*V_,nevind);
1244 } // end if user specified V
1245 else
1246 {
1247 // get vectors from problem or generate something, projectAndNormalize
1248 RCP<const MV> ivec = problem_->getInitVec();
1249 TEUCHOS_TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::invalid_argument,
1250 "Anasazi::TraceMinBase::initialize(newstate): Eigenproblem did not specify initial vectors to clone from.");
1251 // clear newstate so we won't use any data from it below
1252 newstate.X = Teuchos::null;
1253 newstate.MX = Teuchos::null;
1254 newstate.KX = Teuchos::null;
1255 newstate.R = Teuchos::null;
1256 newstate.T = Teuchos::null;
1257 newstate.RV = Teuchos::null;
1258 newstate.KK = Teuchos::null;
1259 newstate.KV = Teuchos::null;
1260 newstate.MopV = Teuchos::null;
1261 newstate.isOrtho = false;
1262
1263 // If the user did not specify a current dimension, use that of the initial vectors they provided
1264 if(newstate.curDim > 0)
1265 curDim_ = newstate.curDim;
1266 else
1267 curDim_ = MVT::GetNumberVecs(*ivec);
1268
1269 // pick the largest multiple of blockSize_
1270 curDim_ = (int)(curDim_ / blockSize_)*blockSize_;
1271 if (curDim_ > blockSize_*numBlocks_) {
1272 // user specified too many vectors... truncate
1273 // this produces a full subspace, but that is okay
1274 curDim_ = blockSize_*numBlocks_;
1275 }
1276 bool userand = false;
1277 if (curDim_ == 0) {
1278 // we need at least blockSize_ vectors
1279 // use a random multivec: ignore everything from InitVec
1280 userand = true;
1281 curDim_ = blockSize_;
1282 }
1283
1284 std::vector<int> nevind(curDim_);
1285 for (int i=0; i<curDim_; ++i) nevind[i] = i;
1286
1287 // get a pointer into V
1288 // lclV has curDim vectors
1289 //
1290 // get pointer to first curDim vectors in V_
1291 lclV = MVT::CloneViewNonConst(*V_,nevind);
1292 if (userand)
1293 {
1294 // generate random vector data
1295 MVT::MvRandom(*lclV);
1296 }
1297 else
1298 {
1299 if(newstate.curDim > 0)
1300 {
1301 if(MVT::GetNumberVecs(*newstate.V) > curDim_) {
1302 RCP<const MV> helperMV = MVT::CloneView(*newstate.V,nevind);
1303 MVT::SetBlock(*helperMV,nevind,*lclV);
1304 }
1305 // assign ivec to first part of lclV
1306 MVT::SetBlock(*newstate.V,nevind,*lclV);
1307 }
1308 else
1309 {
1310 if(MVT::GetNumberVecs(*ivec) > curDim_) {
1311 ivec = MVT::CloneView(*ivec,nevind);
1312 }
1313 // assign ivec to first part of lclV
1314 MVT::SetBlock(*ivec,nevind,*lclV);
1315 }
1316 }
1317 } // end if user did not specify V
1318
1319 // A vector of relevant indices
1320 std::vector<int> dimind(curDim_);
1321 for (int i=0; i<curDim_; ++i) dimind[i] = i;
1322
1323 // Compute MV if necessary
1324 if(hasM_ && newstate.MopV == Teuchos::null)
1325 {
1326 om_->stream(Debug) << "Computing MV\n";
1327
1328#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1329 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1330#endif
1331 count_ApplyM_+= curDim_;
1332
1333 newstate.isOrtho = false;
1334 // Get a pointer to the relevant parts of MV
1335 lclMV = MVT::CloneViewNonConst(*MV_,dimind);
1336 OPT::Apply(*MOp_,*lclV,*lclMV);
1337 }
1338 // Copy MV if necessary
1339 else if(hasM_)
1340 {
1341 om_->stream(Debug) << "Copying MV\n";
1342
1343 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.MopV) != MVT::GetGlobalLength(*MV_), std::invalid_argument,
1344 "Anasazi::TraceMinBase::initialize(newstate): Vector length of MV not correct." );
1345
1346 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.MopV) < curDim_, std::invalid_argument,
1347 "Anasazi::TraceMinBase::initialize(newstate): Number of vectors in MV not correct.");
1348
1349 if(newstate.MopV != MV_) {
1350 if(curDim_ < MVT::GetNumberVecs(*newstate.MopV)) {
1351 newstate.MopV = MVT::CloneView(*newstate.MopV,dimind);
1352 }
1353 MVT::SetBlock(*newstate.MopV,dimind,*MV_);
1354 }
1355 lclMV = MVT::CloneViewNonConst(*MV_,dimind);
1356 }
1357 // There is no M, so there is no MV
1358 else
1359 {
1360 om_->stream(Debug) << "There is no MV\n";
1361
1362 lclMV = lclV;
1363 MV_ = V_;
1364 }
1365
1366 // Project and normalize so that V'V = I and Q'V=0
1367 if(newstate.isOrtho == false)
1368 {
1369 om_->stream(Debug) << "Project and normalize\n";
1370
1371#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1372 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1373#endif
1374
1375 // These things are now invalid
1376 newstate.X = Teuchos::null;
1377 newstate.MX = Teuchos::null;
1378 newstate.KX = Teuchos::null;
1379 newstate.R = Teuchos::null;
1380 newstate.T = Teuchos::null;
1381 newstate.RV = Teuchos::null;
1382 newstate.KK = Teuchos::null;
1383 newstate.KV = Teuchos::null;
1384
1385 int rank;
1386 if(auxVecs_.size() > 0)
1387 {
1388 rank = orthman_->projectAndNormalizeMat(*lclV, auxVecs_,
1389 Teuchos::tuple(RCP< Teuchos::SerialDenseMatrix< int, ScalarType > >(Teuchos::null)),
1390 Teuchos::null, lclMV, MauxVecs_);
1391 }
1392 else
1393 {
1394 rank = orthman_->normalizeMat(*lclV,Teuchos::null,lclMV);
1395 }
1396
1397 TEUCHOS_TEST_FOR_EXCEPTION(rank != curDim_,TraceMinBaseInitFailure,
1398 "Anasazi::TraceMinBase::initialize(): Couldn't generate initial basis of full rank.");
1399 }
1400
1401 // Compute KV
1402 if(Op_ != Teuchos::null && newstate.KV == Teuchos::null)
1403 {
1404 om_->stream(Debug) << "Computing KV\n";
1405
1406#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1407 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1408#endif
1409 count_ApplyOp_+= curDim_;
1410
1411 // These things are now invalid
1412 newstate.X = Teuchos::null;
1413 newstate.MX = Teuchos::null;
1414 newstate.KX = Teuchos::null;
1415 newstate.R = Teuchos::null;
1416 newstate.T = Teuchos::null;
1417 newstate.RV = Teuchos::null;
1418 newstate.KK = Teuchos::null;
1419 newstate.KV = Teuchos::null;
1420
1421 lclKV = MVT::CloneViewNonConst(*KV_,dimind);
1422 OPT::Apply(*Op_,*lclV,*lclKV);
1423 }
1424 // Copy KV
1425 else if(Op_ != Teuchos::null)
1426 {
1427 om_->stream(Debug) << "Copying MV\n";
1428
1429 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.KV) != MVT::GetGlobalLength(*KV_), std::invalid_argument,
1430 "Anasazi::TraceMinBase::initialize(newstate): Vector length of MV not correct." );
1431
1432 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.KV) < curDim_, std::invalid_argument,
1433 "Anasazi::TraceMinBase::initialize(newstate): Number of vectors in KV not correct.");
1434
1435 if (newstate.KV != KV_) {
1436 if (curDim_ < MVT::GetNumberVecs(*newstate.KV)) {
1437 newstate.KV = MVT::CloneView(*newstate.KV,dimind);
1438 }
1439 MVT::SetBlock(*newstate.KV,dimind,*KV_);
1440 }
1441 lclKV = MVT::CloneViewNonConst(*KV_,dimind);
1442 }
1443 else
1444 {
1445 om_->stream(Debug) << "There is no KV\n";
1446
1447 lclKV = lclV;
1448 KV_ = V_;
1449 }
1450
1451 // Compute KK
1452 if(newstate.KK == Teuchos::null)
1453 {
1454 om_->stream(Debug) << "Computing KK\n";
1455
1456 // These things are now invalid
1457 newstate.X = Teuchos::null;
1458 newstate.MX = Teuchos::null;
1459 newstate.KX = Teuchos::null;
1460 newstate.R = Teuchos::null;
1461 newstate.T = Teuchos::null;
1462 newstate.RV = Teuchos::null;
1463
1464 // Note: there used to be a bug here.
1465 // We can't just call computeKK because it only computes the new parts of KK
1466
1467 // Get a pointer to the part of KK we're interested in
1468 lclKK = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
1469
1470 // KK := V'KV
1471 MVT::MvTransMv(ONE,*lclV,*lclKV,*lclKK);
1472 }
1473 // Copy KK
1474 else
1475 {
1476 om_->stream(Debug) << "Copying KK\n";
1477
1478 // check size of KK
1479 TEUCHOS_TEST_FOR_EXCEPTION( newstate.KK->numRows() < curDim_ || newstate.KK->numCols() < curDim_, std::invalid_argument,
1480 "Anasazi::TraceMinBase::initialize(newstate): Projected matrix in new state must be as large as specified state rank.");
1481
1482 // put data into KK_
1483 lclKK = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
1484 if (newstate.KK != KK_) {
1485 if (newstate.KK->numRows() > curDim_ || newstate.KK->numCols() > curDim_) {
1486 newstate.KK = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*newstate.KK,curDim_,curDim_) );
1487 }
1488 lclKK->assign(*newstate.KK);
1489 }
1490 }
1491
1492 // Compute Ritz pairs
1493 if(newstate.T == Teuchos::null || newstate.RV == Teuchos::null)
1494 {
1495 om_->stream(Debug) << "Computing Ritz pairs\n";
1496
1497 // These things are now invalid
1498 newstate.X = Teuchos::null;
1499 newstate.MX = Teuchos::null;
1500 newstate.KX = Teuchos::null;
1501 newstate.R = Teuchos::null;
1502 newstate.T = Teuchos::null;
1503 newstate.RV = Teuchos::null;
1504
1505 computeRitzPairs();
1506 }
1507 // Copy Ritz pairs
1508 else
1509 {
1510 om_->stream(Debug) << "Copying Ritz pairs\n";
1511
1512 TEUCHOS_TEST_FOR_EXCEPTION((signed int)(newstate.T->size()) != curDim_,
1513 std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of T must be consistent with dimension of V.");
1514
1515 TEUCHOS_TEST_FOR_EXCEPTION( newstate.RV->numRows() < curDim_ || newstate.RV->numCols() < curDim_, std::invalid_argument,
1516 "Anasazi::TraceMinBase::initialize(newstate): Ritz vectors in new state must be as large as specified state rank.");
1517
1518 std::copy(newstate.T->begin(),newstate.T->end(),theta_.begin());
1519
1520 lclRV = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*ritzVecs_,curDim_,curDim_) );
1521 if (newstate.RV != ritzVecs_) {
1522 if (newstate.RV->numRows() > curDim_ || newstate.RV->numCols() > curDim_) {
1523 newstate.RV = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*newstate.RV,curDim_,curDim_) );
1524 }
1525 lclRV->assign(*newstate.RV);
1526 }
1527 }
1528
1529 // Compute X
1530 if(newstate.X == Teuchos::null)
1531 {
1532 om_->stream(Debug) << "Computing X\n";
1533
1534 // These things are now invalid
1535 newstate.MX = Teuchos::null;
1536 newstate.KX = Teuchos::null;
1537 newstate.R = Teuchos::null;
1538
1539 computeX();
1540 }
1541 // Copy X
1542 else
1543 {
1544 om_->stream(Debug) << "Copying X\n";
1545
1546 if(computeAllRes_ == false)
1547 {
1548 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.X) != blockSize_ || MVT::GetGlobalLength(*newstate.X) != MVT::GetGlobalLength(*X_),
1549 std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of X must be consistent with block size and length of V.");
1550
1551 if(newstate.X != X_) {
1552 MVT::SetBlock(*newstate.X,bsind,*X_);
1553 }
1554 }
1555 else
1556 {
1557 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.X) != curDim_ || MVT::GetGlobalLength(*newstate.X) != MVT::GetGlobalLength(*X_),
1558 std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of X must be consistent with current dimension and length of V.");
1559
1560 if(newstate.X != X_) {
1561 MVT::SetBlock(*newstate.X,dimind,*X_);
1562 }
1563 }
1564 }
1565
1566 // Compute KX and MX if necessary
1567 // TODO: These technically should be separate; it won't matter much in terms of running time though
1568 if((Op_ != Teuchos::null && newstate.KX == Teuchos::null) || (hasM_ && newstate.MX == Teuchos::null))
1569 {
1570 om_->stream(Debug) << "Computing KX and MX\n";
1571
1572 // These things are now invalid
1573 newstate.R = Teuchos::null;
1574
1575 updateKXMX();
1576 }
1577 // Copy KX and MX if necessary
1578 else
1579 {
1580 om_->stream(Debug) << "Copying KX and MX\n";
1581 if(Op_ != Teuchos::null)
1582 {
1583 if(computeAllRes_ == false)
1584 {
1585 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.KX) != blockSize_ || MVT::GetGlobalLength(*newstate.KX) != MVT::GetGlobalLength(*KX_),
1586 std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of KX must be consistent with block size and length of V.");
1587
1588 if(newstate.KX != KX_) {
1589 MVT::SetBlock(*newstate.KX,bsind,*KX_);
1590 }
1591 }
1592 else
1593 {
1594 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.KX) != curDim_ || MVT::GetGlobalLength(*newstate.KX) != MVT::GetGlobalLength(*KX_),
1595 std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of KX must be consistent with current dimension and length of V.");
1596
1597 if (newstate.KX != KX_) {
1598 MVT::SetBlock(*newstate.KX,dimind,*KX_);
1599 }
1600 }
1601 }
1602
1603 if(hasM_)
1604 {
1605 if(computeAllRes_ == false)
1606 {
1607 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.MX) != blockSize_ || MVT::GetGlobalLength(*newstate.MX) != MVT::GetGlobalLength(*MX_),
1608 std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of MX must be consistent with block size and length of V.");
1609
1610 if (newstate.MX != MX_) {
1611 MVT::SetBlock(*newstate.MX,bsind,*MX_);
1612 }
1613 }
1614 else
1615 {
1616 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.MX) != curDim_ || MVT::GetGlobalLength(*newstate.MX) != MVT::GetGlobalLength(*MX_),
1617 std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of MX must be consistent with current dimension and length of V.");
1618
1619 if (newstate.MX != MX_) {
1620 MVT::SetBlock(*newstate.MX,dimind,*MX_);
1621 }
1622 }
1623 }
1624 }
1625
1626 // Compute R
1627 if(newstate.R == Teuchos::null)
1628 {
1629 om_->stream(Debug) << "Computing R\n";
1630
1631 updateResidual();
1632 }
1633 // Copy R
1634 else
1635 {
1636 om_->stream(Debug) << "Copying R\n";
1637
1638 if(computeAllRes_ == false)
1639 {
1640 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*X_),
1641 std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): vector length of newstate.R not correct." );
1642 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.R) != blockSize_,
1643 std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): newstate.R must have at least block size vectors." );
1644 if (newstate.R != R_) {
1645 MVT::SetBlock(*newstate.R,bsind,*R_);
1646 }
1647 }
1648 else
1649 {
1650 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*X_),
1651 std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): vector length of newstate.R not correct." );
1652 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.R) != curDim_,
1653 std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): newstate.R must have at least curDim vectors." );
1654 if (newstate.R != R_) {
1655 MVT::SetBlock(*newstate.R,dimind,*R_);
1656 }
1657 }
1658 }
1659
1660 // R has been updated; mark the norms as out-of-date
1661 Rnorms_current_ = false;
1662 R2norms_current_ = false;
1663
1664 // Set the largest safe shift
1665 largestSafeShift_ = newstate.largestSafeShift;
1666
1667 // Copy over the Ritz shifts
1668 if(newstate.ritzShifts != Teuchos::null)
1669 {
1670 om_->stream(Debug) << "Copying Ritz shifts\n";
1671 std::copy(newstate.ritzShifts->begin(),newstate.ritzShifts->end(),ritzShifts_.begin());
1672 }
1673 else
1674 {
1675 om_->stream(Debug) << "Setting Ritz shifts to 0\n";
1676 for(size_t i=0; i<ritzShifts_.size(); i++)
1677 ritzShifts_[i] = ZERO;
1678 }
1679
1680 for(size_t i=0; i<ritzShifts_.size(); i++)
1681 om_->stream(Debug) << "Ritz shifts[" << i << "] = " << ritzShifts_[i] << std::endl;
1682
1683 // finally, we are initialized
1684 initialized_ = true;
1685
1686 if (om_->isVerbosity( Debug ) ) {
1687 // Check almost everything here
1688 CheckList chk;
1689 chk.checkV = true;
1690 chk.checkX = true;
1691 chk.checkKX = true;
1692 chk.checkMX = true;
1693 chk.checkQ = true;
1694 chk.checkKK = true;
1695 om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
1696 }
1697
1698 // Print information on current status
1699 if (om_->isVerbosity(Debug)) {
1700 currentStatus( om_->stream(Debug) );
1701 }
1702 else if (om_->isVerbosity(IterationDetails)) {
1703 currentStatus( om_->stream(IterationDetails) );
1704 }
1705 }
1706
1707
1708
1710 /* Initialize the state of the solver
1711 *
1712 * POST-CONDITIONS:
1713 *
1714 * V_ is orthonormal, orthogonal to auxVecs_, for first curDim_ vectors
1715 * theta_ contains Ritz w.r.t. V_(1:curDim_)
1716 * X is Ritz vectors w.r.t. V_(1:curDim_)
1717 * KX = Op*X
1718 * MX = M*X if hasM_
1719 * R = KX - MX*diag(theta_)
1720 *
1721 */
1722 template <class ScalarType, class MV, class OP>
1723 void TraceMinBase<ScalarType,MV,OP>::harmonicInitialize(TraceMinBaseState<ScalarType,MV> newstate)
1724 {
1725 // TODO: Check for bad input
1726 // NOTE: memory has been allocated by setBlockSize(). Use setBlock below; do not Clone
1727 // NOTE: Overall time spent in this routine is counted to timerInit_; portions will also be counted towards other primitives
1728
1729 std::vector<int> bsind(blockSize_);
1730 for (int i=0; i<blockSize_; ++i) bsind[i] = i;
1731
1732 // in TraceMinBase, V is primary
1733 // the order of dependence follows like so.
1734 // --init-> V,KK
1735 // --ritz analysis-> theta,X
1736 // --op apply-> KX,MX
1737 // --compute-> R
1738 //
1739 // if the user specifies all data for a level, we will accept it.
1740 // otherwise, we will generate the whole level, and all subsequent levels.
1741 //
1742 // the data members are ordered based on dependence, and the levels are
1743 // partitioned according to the amount of work required to produce the
1744 // items in a level.
1745 //
1746 // inconsistent multivectors widths and lengths will not be tolerated, and
1747 // will be treated with exceptions.
1748 //
1749 // for multivector pointers in newstate which point directly (as opposed to indirectly, via a view) to
1750 // multivectors in the solver, the copy will not be affected.
1751
1752 // set up V and KK: get them from newstate if user specified them
1753 // otherwise, set them manually
1754 RCP<MV> lclV, lclKV, lclMV;
1755 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclKK, lclRV;
1756
1757 // If V was supplied by the user...
1758 if (newstate.V != Teuchos::null) {
1759 om_->stream(Debug) << "Copying V from the user\n";
1760
1761 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.V) != MVT::GetGlobalLength(*V_), std::invalid_argument,
1762 "Anasazi::TraceMinBase::initialize(newstate): Vector length of V not correct." );
1763 TEUCHOS_TEST_FOR_EXCEPTION( newstate.curDim < blockSize_, std::invalid_argument,
1764 "Anasazi::TraceMinBase::initialize(newstate): Rank of new state must be at least blockSize().");
1765 TEUCHOS_TEST_FOR_EXCEPTION( newstate.curDim > blockSize_*numBlocks_, std::invalid_argument,
1766 "Anasazi::TraceMinBase::initialize(newstate): Rank of new state must be less than getMaxSubspaceDim().");
1767
1768 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < newstate.curDim, std::invalid_argument,
1769 "Anasazi::TraceMinBase::initialize(newstate): Multivector for basis in new state must be as large as specified state rank.");
1770
1771 curDim_ = newstate.curDim;
1772 // pick an integral amount
1773 curDim_ = (int)(curDim_ / blockSize_)*blockSize_;
1774
1775 TEUCHOS_TEST_FOR_EXCEPTION( curDim_ != newstate.curDim, std::invalid_argument,
1776 "Anasazi::TraceMinBase::initialize(newstate): Rank of new state must be a multiple of getBlockSize().");
1777
1778 // put data in V
1779 std::vector<int> nevind(curDim_);
1780 for (int i=0; i<curDim_; ++i) nevind[i] = i;
1781 if (newstate.V != V_) {
1782 if (curDim_ < MVT::GetNumberVecs(*newstate.V)) {
1783 newstate.V = MVT::CloneView(*newstate.V,nevind);
1784 }
1785 MVT::SetBlock(*newstate.V,nevind,*V_);
1786 }
1787 lclV = MVT::CloneViewNonConst(*V_,nevind);
1788 } // end if user specified V
1789 else
1790 {
1791 // get vectors from problem or generate something, projectAndNormalize
1792 RCP<const MV> ivec = problem_->getInitVec();
1793 TEUCHOS_TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::invalid_argument,
1794 "Anasazi::TraceMinBase::initialize(newstate): Eigenproblem did not specify initial vectors to clone from.");
1795 // clear newstate so we won't use any data from it below
1796 newstate.X = Teuchos::null;
1797 newstate.MX = Teuchos::null;
1798 newstate.KX = Teuchos::null;
1799 newstate.R = Teuchos::null;
1800 newstate.T = Teuchos::null;
1801 newstate.RV = Teuchos::null;
1802 newstate.KK = Teuchos::null;
1803 newstate.KV = Teuchos::null;
1804 newstate.MopV = Teuchos::null;
1805 newstate.isOrtho = false;
1806
1807 // If the user did not specify a current dimension, use that of the initial vectors they provided
1808 if(newstate.curDim > 0)
1809 curDim_ = newstate.curDim;
1810 else
1811 curDim_ = MVT::GetNumberVecs(*ivec);
1812
1813 // pick the largest multiple of blockSize_
1814 curDim_ = (int)(curDim_ / blockSize_)*blockSize_;
1815 if (curDim_ > blockSize_*numBlocks_) {
1816 // user specified too many vectors... truncate
1817 // this produces a full subspace, but that is okay
1818 curDim_ = blockSize_*numBlocks_;
1819 }
1820 bool userand = false;
1821 if (curDim_ == 0) {
1822 // we need at least blockSize_ vectors
1823 // use a random multivec: ignore everything from InitVec
1824 userand = true;
1825 curDim_ = blockSize_;
1826 }
1827
1828 std::vector<int> nevind(curDim_);
1829 for (int i=0; i<curDim_; ++i) nevind[i] = i;
1830
1831 // get a pointer into V
1832 // lclV has curDim vectors
1833 //
1834 // get pointer to first curDim vectors in V_
1835 lclV = MVT::CloneViewNonConst(*V_,nevind);
1836 if (userand)
1837 {
1838 // generate random vector data
1839 MVT::MvRandom(*lclV);
1840 }
1841 else
1842 {
1843 if(newstate.curDim > 0)
1844 {
1845 if(MVT::GetNumberVecs(*newstate.V) > curDim_) {
1846 RCP<const MV> helperMV = MVT::CloneView(*newstate.V,nevind);
1847 MVT::SetBlock(*helperMV,nevind,*lclV);
1848 }
1849 // assign ivec to first part of lclV
1850 MVT::SetBlock(*newstate.V,nevind,*lclV);
1851 }
1852 else
1853 {
1854 if(MVT::GetNumberVecs(*ivec) > curDim_) {
1855 ivec = MVT::CloneView(*ivec,nevind);
1856 }
1857 // assign ivec to first part of lclV
1858 MVT::SetBlock(*ivec,nevind,*lclV);
1859 }
1860 }
1861 } // end if user did not specify V
1862
1863 // Nuke everything from orbit
1864 // This is a temporary measure due to a bug in the code that I have not found yet
1865 // It adds a minimal amount of work
1866 newstate.X = Teuchos::null;
1867 newstate.MX = Teuchos::null;
1868 newstate.KX = Teuchos::null;
1869 newstate.R = Teuchos::null;
1870 newstate.T = Teuchos::null;
1871 newstate.RV = Teuchos::null;
1872 newstate.KK = Teuchos::null;
1873 newstate.KV = Teuchos::null;
1874 newstate.MopV = Teuchos::null;
1875 newstate.isOrtho = false;
1876
1877 // A vector of relevant indices
1878 std::vector<int> dimind(curDim_);
1879 for (int i=0; i<curDim_; ++i) dimind[i] = i;
1880
1881 // Project the auxVecs out of V
1882 if(auxVecs_.size() > 0)
1883 orthman_->projectMat(*lclV,auxVecs_);
1884
1885 // Compute KV
1886 if(Op_ != Teuchos::null && newstate.KV == Teuchos::null)
1887 {
1888 om_->stream(Debug) << "Computing KV\n";
1889
1890#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1891 Teuchos::TimeMonitor lcltimer( *timerOp_ );
1892#endif
1893 count_ApplyOp_+= curDim_;
1894
1895 // These things are now invalid
1896 newstate.X = Teuchos::null;
1897 newstate.MX = Teuchos::null;
1898 newstate.KX = Teuchos::null;
1899 newstate.R = Teuchos::null;
1900 newstate.T = Teuchos::null;
1901 newstate.RV = Teuchos::null;
1902 newstate.KK = Teuchos::null;
1903
1904 lclKV = MVT::CloneViewNonConst(*KV_,dimind);
1905 OPT::Apply(*Op_,*lclV,*lclKV);
1906 }
1907 // Copy KV
1908 else if(Op_ != Teuchos::null)
1909 {
1910 om_->stream(Debug) << "Copying KV\n";
1911
1912 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.KV) != MVT::GetGlobalLength(*KV_), std::invalid_argument,
1913 "Anasazi::TraceMinBase::initialize(newstate): Vector length of KV not correct." );
1914
1915 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.KV) < curDim_, std::invalid_argument,
1916 "Anasazi::TraceMinBase::initialize(newstate): Number of vectors in KV not correct.");
1917
1918 if (newstate.KV != KV_) {
1919 if (curDim_ < MVT::GetNumberVecs(*newstate.KV)) {
1920 newstate.KV = MVT::CloneView(*newstate.KV,dimind);
1921 }
1922 MVT::SetBlock(*newstate.KV,dimind,*KV_);
1923 }
1924 lclKV = MVT::CloneViewNonConst(*KV_,dimind);
1925 }
1926 else
1927 {
1928 om_->stream(Debug) << "There is no KV\n";
1929
1930 lclKV = lclV;
1931 KV_ = V_;
1932 }
1933
1934
1935
1936 // Project and normalize so that V'V = I and Q'V=0
1937 if(newstate.isOrtho == false)
1938 {
1939 om_->stream(Debug) << "Project and normalize\n";
1940
1941#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1942 Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
1943#endif
1944
1945 // These things are now invalid
1946 newstate.MopV = Teuchos::null;
1947 newstate.X = Teuchos::null;
1948 newstate.MX = Teuchos::null;
1949 newstate.KX = Teuchos::null;
1950 newstate.R = Teuchos::null;
1951 newstate.T = Teuchos::null;
1952 newstate.RV = Teuchos::null;
1953 newstate.KK = Teuchos::null;
1954
1955 // Normalize lclKV
1956 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > gamma = rcp(new Teuchos::SerialDenseMatrix<int,ScalarType>(curDim_,curDim_));
1957 int rank = orthman_->normalizeMat(*lclKV,gamma);
1958
1959 // lclV = lclV/gamma
1960 Teuchos::SerialDenseSolver<int,ScalarType> SDsolver;
1961 SDsolver.setMatrix(gamma);
1962 SDsolver.invert();
1963 RCP<MV> tempMV = MVT::CloneCopy(*lclV);
1964 MVT::MvTimesMatAddMv(ONE,*tempMV,*gamma,ZERO,*lclV);
1965
1966 TEUCHOS_TEST_FOR_EXCEPTION(rank != curDim_,TraceMinBaseInitFailure,
1967 "Anasazi::TraceMinBase::initialize(): Couldn't generate initial basis of full rank.");
1968 }
1969
1970 // Compute MV if necessary
1971 if(hasM_ && newstate.MopV == Teuchos::null)
1972 {
1973 om_->stream(Debug) << "Computing MV\n";
1974
1975#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1976 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
1977#endif
1978 count_ApplyM_+= curDim_;
1979
1980 // Get a pointer to the relevant parts of MV
1981 lclMV = MVT::CloneViewNonConst(*MV_,dimind);
1982 OPT::Apply(*MOp_,*lclV,*lclMV);
1983 }
1984 // Copy MV if necessary
1985 else if(hasM_)
1986 {
1987 om_->stream(Debug) << "Copying MV\n";
1988
1989 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.MopV) != MVT::GetGlobalLength(*MV_), std::invalid_argument,
1990 "Anasazi::TraceMinBase::initialize(newstate): Vector length of MV not correct." );
1991
1992 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.MopV) < curDim_, std::invalid_argument,
1993 "Anasazi::TraceMinBase::initialize(newstate): Number of vectors in MV not correct.");
1994
1995 if(newstate.MopV != MV_) {
1996 if(curDim_ < MVT::GetNumberVecs(*newstate.MopV)) {
1997 newstate.MopV = MVT::CloneView(*newstate.MopV,dimind);
1998 }
1999 MVT::SetBlock(*newstate.MopV,dimind,*MV_);
2000 }
2001 lclMV = MVT::CloneViewNonConst(*MV_,dimind);
2002 }
2003 // There is no M, so there is no MV
2004 else
2005 {
2006 om_->stream(Debug) << "There is no MV\n";
2007
2008 lclMV = lclV;
2009 MV_ = V_;
2010 }
2011
2012 // Compute KK
2013 if(newstate.KK == Teuchos::null)
2014 {
2015 om_->stream(Debug) << "Computing KK\n";
2016
2017 // These things are now invalid
2018 newstate.X = Teuchos::null;
2019 newstate.MX = Teuchos::null;
2020 newstate.KX = Teuchos::null;
2021 newstate.R = Teuchos::null;
2022 newstate.T = Teuchos::null;
2023 newstate.RV = Teuchos::null;
2024
2025 // Note: there used to be a bug here.
2026 // We can't just call computeKK because it only computes the new parts of KK
2027
2028 // Get a pointer to the part of KK we're interested in
2029 lclKK = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
2030
2031 // KK := V'KV
2032 MVT::MvTransMv(ONE,*lclV,*lclKV,*lclKK);
2033 }
2034 // Copy KK
2035 else
2036 {
2037 om_->stream(Debug) << "Copying KK\n";
2038
2039 // check size of KK
2040 TEUCHOS_TEST_FOR_EXCEPTION( newstate.KK->numRows() < curDim_ || newstate.KK->numCols() < curDim_, std::invalid_argument,
2041 "Anasazi::TraceMinBase::initialize(newstate): Projected matrix in new state must be as large as specified state rank.");
2042
2043 // put data into KK_
2044 lclKK = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
2045 if (newstate.KK != KK_) {
2046 if (newstate.KK->numRows() > curDim_ || newstate.KK->numCols() > curDim_) {
2047 newstate.KK = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*newstate.KK,curDim_,curDim_) );
2048 }
2049 lclKK->assign(*newstate.KK);
2050 }
2051 }
2052
2053 // Compute Ritz pairs
2054 if(newstate.T == Teuchos::null || newstate.RV == Teuchos::null)
2055 {
2056 om_->stream(Debug) << "Computing Ritz pairs\n";
2057
2058 // These things are now invalid
2059 newstate.X = Teuchos::null;
2060 newstate.MX = Teuchos::null;
2061 newstate.KX = Teuchos::null;
2062 newstate.R = Teuchos::null;
2063 newstate.T = Teuchos::null;
2064 newstate.RV = Teuchos::null;
2065
2066 computeRitzPairs();
2067 }
2068 // Copy Ritz pairs
2069 else
2070 {
2071 om_->stream(Debug) << "Copying Ritz pairs\n";
2072
2073 TEUCHOS_TEST_FOR_EXCEPTION((signed int)(newstate.T->size()) != curDim_,
2074 std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of T must be consistent with dimension of V.");
2075
2076 TEUCHOS_TEST_FOR_EXCEPTION( newstate.RV->numRows() < curDim_ || newstate.RV->numCols() < curDim_, std::invalid_argument,
2077 "Anasazi::TraceMinBase::initialize(newstate): Ritz vectors in new state must be as large as specified state rank.");
2078
2079 std::copy(newstate.T->begin(),newstate.T->end(),theta_.begin());
2080
2081 lclRV = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*ritzVecs_,curDim_,curDim_) );
2082 if (newstate.RV != ritzVecs_) {
2083 if (newstate.RV->numRows() > curDim_ || newstate.RV->numCols() > curDim_) {
2084 newstate.RV = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*newstate.RV,curDim_,curDim_) );
2085 }
2086 lclRV->assign(*newstate.RV);
2087 }
2088 }
2089
2090 // Compute X
2091 if(newstate.X == Teuchos::null)
2092 {
2093 om_->stream(Debug) << "Computing X\n";
2094
2095 // These things are now invalid
2096 newstate.MX = Teuchos::null;
2097 newstate.KX = Teuchos::null;
2098 newstate.R = Teuchos::null;
2099
2100 computeX();
2101 }
2102 // Copy X
2103 else
2104 {
2105 om_->stream(Debug) << "Copying X\n";
2106
2107 if(computeAllRes_ == false)
2108 {
2109 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.X) != blockSize_ || MVT::GetGlobalLength(*newstate.X) != MVT::GetGlobalLength(*X_),
2110 std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of X must be consistent with block size and length of V.");
2111
2112 if(newstate.X != X_) {
2113 MVT::SetBlock(*newstate.X,bsind,*X_);
2114 }
2115 }
2116 else
2117 {
2118 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.X) != curDim_ || MVT::GetGlobalLength(*newstate.X) != MVT::GetGlobalLength(*X_),
2119 std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of X must be consistent with current dimension and length of V.");
2120
2121 if(newstate.X != X_) {
2122 MVT::SetBlock(*newstate.X,dimind,*X_);
2123 }
2124 }
2125 }
2126
2127 // Compute KX and MX if necessary
2128 // TODO: These technically should be separate; it won't matter much in terms of running time though
2129 if((Op_ != Teuchos::null && newstate.KX == Teuchos::null) || (hasM_ && newstate.MX == Teuchos::null))
2130 {
2131 om_->stream(Debug) << "Computing KX and MX\n";
2132
2133 // These things are now invalid
2134 newstate.R = Teuchos::null;
2135
2136 updateKXMX();
2137 }
2138 // Copy KX and MX if necessary
2139 else
2140 {
2141 om_->stream(Debug) << "Copying KX and MX\n";
2142 if(Op_ != Teuchos::null)
2143 {
2144 if(computeAllRes_ == false)
2145 {
2146 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.KX) != blockSize_ || MVT::GetGlobalLength(*newstate.KX) != MVT::GetGlobalLength(*KX_),
2147 std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of KX must be consistent with block size and length of V.");
2148
2149 if(newstate.KX != KX_) {
2150 MVT::SetBlock(*newstate.KX,bsind,*KX_);
2151 }
2152 }
2153 else
2154 {
2155 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.KX) != curDim_ || MVT::GetGlobalLength(*newstate.KX) != MVT::GetGlobalLength(*KX_),
2156 std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of KX must be consistent with current dimension and length of V.");
2157
2158 if (newstate.KX != KX_) {
2159 MVT::SetBlock(*newstate.KX,dimind,*KX_);
2160 }
2161 }
2162 }
2163
2164 if(hasM_)
2165 {
2166 if(computeAllRes_ == false)
2167 {
2168 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.MX) != blockSize_ || MVT::GetGlobalLength(*newstate.MX) != MVT::GetGlobalLength(*MX_),
2169 std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of MX must be consistent with block size and length of V.");
2170
2171 if (newstate.MX != MX_) {
2172 MVT::SetBlock(*newstate.MX,bsind,*MX_);
2173 }
2174 }
2175 else
2176 {
2177 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.MX) != curDim_ || MVT::GetGlobalLength(*newstate.MX) != MVT::GetGlobalLength(*MX_),
2178 std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): Size of MX must be consistent with current dimension and length of V.");
2179
2180 if (newstate.MX != MX_) {
2181 MVT::SetBlock(*newstate.MX,dimind,*MX_);
2182 }
2183 }
2184 }
2185 }
2186
2187 // Scale X so each vector is of length 1
2188 {
2189 // Get norm of each vector in X
2190 const int nvecs = computeAllRes_ ? curDim_ : blockSize_;
2191 Teuchos::Range1D dimind2 (0, nvecs-1);
2192 RCP<MV> lclX = MVT::CloneViewNonConst(*X_, dimind2);
2193 std::vector<ScalarType> normvec(nvecs);
2194 orthman_->normMat(*lclX,normvec);
2195
2196 // Scale X, KX, and MX accordingly
2197 for (int i = 0; i < nvecs; ++i) {
2198 normvec[i] = ONE / normvec[i];
2199 }
2200 MVT::MvScale (*lclX, normvec);
2201 if (Op_ != Teuchos::null) {
2202 RCP<MV> lclKX = MVT::CloneViewNonConst (*KX_, dimind2);
2203 MVT::MvScale (*lclKX, normvec);
2204 }
2205 if (hasM_) {
2206 RCP<MV> lclMX = MVT::CloneViewNonConst (*MX_, dimind2);
2207 MVT::MvScale (*lclMX, normvec);
2208 }
2209
2210 // Scale eigenvalues
2211 for (int i = 0; i < nvecs; ++i) {
2212 theta_[i] = theta_[i] * normvec[i] * normvec[i];
2213 }
2214 }
2215
2216 // Compute R
2217 if(newstate.R == Teuchos::null)
2218 {
2219 om_->stream(Debug) << "Computing R\n";
2220
2221 updateResidual();
2222 }
2223 // Copy R
2224 else
2225 {
2226 om_->stream(Debug) << "Copying R\n";
2227
2228 if(computeAllRes_ == false)
2229 {
2230 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*X_),
2231 std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): vector length of newstate.R not correct." );
2232 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.R) != blockSize_,
2233 std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): newstate.R must have at least block size vectors." );
2234 if (newstate.R != R_) {
2235 MVT::SetBlock(*newstate.R,bsind,*R_);
2236 }
2237 }
2238 else
2239 {
2240 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*X_),
2241 std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): vector length of newstate.R not correct." );
2242 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*newstate.R) != curDim_,
2243 std::invalid_argument, "Anasazi::TraceMinBase::initialize(newstate): newstate.R must have at least curDim vectors." );
2244 if (newstate.R != R_) {
2245 MVT::SetBlock(*newstate.R,dimind,*R_);
2246 }
2247 }
2248 }
2249
2250 // R has been updated; mark the norms as out-of-date
2251 Rnorms_current_ = false;
2252 R2norms_current_ = false;
2253
2254 // Set the largest safe shift
2255 largestSafeShift_ = newstate.largestSafeShift;
2256
2257 // Copy over the Ritz shifts
2258 if(newstate.ritzShifts != Teuchos::null)
2259 {
2260 om_->stream(Debug) << "Copying Ritz shifts\n";
2261 std::copy(newstate.ritzShifts->begin(),newstate.ritzShifts->end(),ritzShifts_.begin());
2262 }
2263 else
2264 {
2265 om_->stream(Debug) << "Setting Ritz shifts to 0\n";
2266 for(size_t i=0; i<ritzShifts_.size(); i++)
2267 ritzShifts_[i] = ZERO;
2268 }
2269
2270 for(size_t i=0; i<ritzShifts_.size(); i++)
2271 om_->stream(Debug) << "Ritz shifts[" << i << "] = " << ritzShifts_[i] << std::endl;
2272
2273 // finally, we are initialized
2274 initialized_ = true;
2275
2276 if (om_->isVerbosity( Debug ) ) {
2277 // Check almost everything here
2278 CheckList chk;
2279 chk.checkV = true;
2280 chk.checkX = true;
2281 chk.checkKX = true;
2282 chk.checkMX = true;
2283 chk.checkQ = true;
2284 chk.checkKK = true;
2285 om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
2286 }
2287
2288 // Print information on current status
2289 if (om_->isVerbosity(Debug)) {
2290 currentStatus( om_->stream(Debug) );
2291 }
2292 else if (om_->isVerbosity(IterationDetails)) {
2293 currentStatus( om_->stream(IterationDetails) );
2294 }
2295 }
2296
2297
2299 // Return initialized state
2300 template <class ScalarType, class MV, class OP>
2301 bool TraceMinBase<ScalarType,MV,OP>::isInitialized() const { return initialized_; }
2302
2303
2305 // Set the block size and make necessary adjustments.
2306 template <class ScalarType, class MV, class OP>
2307 void TraceMinBase<ScalarType,MV,OP>::setSize (int blockSize, int numBlocks)
2308 {
2309 // This routine only allocates space; it doesn't not perform any computation
2310 // any change in size will invalidate the state of the solver.
2311
2312 TEUCHOS_TEST_FOR_EXCEPTION(blockSize < 1, std::invalid_argument, "Anasazi::TraceMinBase::setSize(blocksize,numblocks): blocksize must be strictly positive.");
2313
2314 if (blockSize == blockSize_ && numBlocks == numBlocks_) {
2315 // do nothing
2316 return;
2317 }
2318
2319 blockSize_ = blockSize;
2320 numBlocks_ = numBlocks;
2321
2322 RCP<const MV> tmp;
2323 // grab some Multivector to Clone
2324 // in practice, getInitVec() should always provide this, but it is possible to use a
2325 // Eigenproblem with nothing in getInitVec() by manually initializing with initialize();
2326 // in case of that strange scenario, we will try to Clone from X_ first, then resort to getInitVec()
2327 if (X_ != Teuchos::null) { // this is equivalent to blockSize_ > 0
2328 tmp = X_;
2329 }
2330 else {
2331 tmp = problem_->getInitVec();
2332 TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
2333 "Anasazi::TraceMinBase::setSize(): eigenproblem did not specify initial vectors to clone from.");
2334 }
2335
2336 TEUCHOS_TEST_FOR_EXCEPTION(numAuxVecs_+blockSize*static_cast<ptrdiff_t>(numBlocks) > MVT::GetGlobalLength(*tmp),std::invalid_argument,
2337 "Anasazi::TraceMinBase::setSize(): max subspace dimension and auxilliary subspace too large. Potentially impossible orthogonality constraints.");
2338
2339 // New subspace dimension
2340 int newsd = blockSize_*numBlocks_;
2341
2343 // blockSize dependent
2344 //
2345 ritzShifts_.resize(blockSize_,ZERO);
2346 if(computeAllRes_ == false)
2347 {
2348 // grow/allocate vectors
2349 Rnorms_.resize(blockSize_,NANVAL);
2350 R2norms_.resize(blockSize_,NANVAL);
2351 //
2352 // clone multivectors off of tmp
2353 //
2354 // free current allocation first, to make room for new allocation
2355 X_ = Teuchos::null;
2356 KX_ = Teuchos::null;
2357 MX_ = Teuchos::null;
2358 R_ = Teuchos::null;
2359 V_ = Teuchos::null;
2360 KV_ = Teuchos::null;
2361 MV_ = Teuchos::null;
2362
2363 om_->print(Debug," >> Allocating X_\n");
2364 X_ = MVT::Clone(*tmp,blockSize_);
2365 if(Op_ != Teuchos::null) {
2366 om_->print(Debug," >> Allocating KX_\n");
2367 KX_ = MVT::Clone(*tmp,blockSize_);
2368 }
2369 else {
2370 KX_ = X_;
2371 }
2372 if (hasM_) {
2373 om_->print(Debug," >> Allocating MX_\n");
2374 MX_ = MVT::Clone(*tmp,blockSize_);
2375 }
2376 else {
2377 MX_ = X_;
2378 }
2379 om_->print(Debug," >> Allocating R_\n");
2380 R_ = MVT::Clone(*tmp,blockSize_);
2381 }
2382 else
2383 {
2384 // grow/allocate vectors
2385 Rnorms_.resize(newsd,NANVAL);
2386 R2norms_.resize(newsd,NANVAL);
2387 //
2388 // clone multivectors off of tmp
2389 //
2390 // free current allocation first, to make room for new allocation
2391 X_ = Teuchos::null;
2392 KX_ = Teuchos::null;
2393 MX_ = Teuchos::null;
2394 R_ = Teuchos::null;
2395 V_ = Teuchos::null;
2396 KV_ = Teuchos::null;
2397 MV_ = Teuchos::null;
2398
2399 om_->print(Debug," >> Allocating X_\n");
2400 X_ = MVT::Clone(*tmp,newsd);
2401 if(Op_ != Teuchos::null) {
2402 om_->print(Debug," >> Allocating KX_\n");
2403 KX_ = MVT::Clone(*tmp,newsd);
2404 }
2405 else {
2406 KX_ = X_;
2407 }
2408 if (hasM_) {
2409 om_->print(Debug," >> Allocating MX_\n");
2410 MX_ = MVT::Clone(*tmp,newsd);
2411 }
2412 else {
2413 MX_ = X_;
2414 }
2415 om_->print(Debug," >> Allocating R_\n");
2416 R_ = MVT::Clone(*tmp,newsd);
2417 }
2418
2420 // blockSize*numBlocks dependent
2421 //
2422 theta_.resize(newsd,NANVAL);
2423 om_->print(Debug," >> Allocating V_\n");
2424 V_ = MVT::Clone(*tmp,newsd);
2425 KK_ = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(newsd,newsd) );
2426 ritzVecs_ = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(newsd,newsd) );
2427
2428 if(Op_ != Teuchos::null) {
2429 om_->print(Debug," >> Allocating KV_\n");
2430 KV_ = MVT::Clone(*tmp,newsd);
2431 }
2432 else {
2433 KV_ = V_;
2434 }
2435 if (hasM_) {
2436 om_->print(Debug," >> Allocating MV_\n");
2437 MV_ = MVT::Clone(*tmp,newsd);
2438 }
2439 else {
2440 MV_ = V_;
2441 }
2442
2443 om_->print(Debug," >> done allocating.\n");
2444
2445 initialized_ = false;
2446 curDim_ = 0;
2447 }
2448
2449
2451 // Set the auxiliary vectors
2452 template <class ScalarType, class MV, class OP>
2453 void TraceMinBase<ScalarType,MV,OP>::setAuxVecs(const Teuchos::Array<RCP<const MV> > &auxvecs) {
2454 typedef typename Teuchos::Array<RCP<const MV> >::iterator tarcpmv;
2455
2456 // set new auxiliary vectors
2457 auxVecs_ = auxvecs;
2458
2459 if(hasM_)
2460 MauxVecs_.resize(0);
2461 numAuxVecs_ = 0;
2462
2463 for (tarcpmv i=auxVecs_.begin(); i != auxVecs_.end(); ++i) {
2464 numAuxVecs_ += MVT::GetNumberVecs(**i);
2465
2466 if(hasM_)
2467 {
2468#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
2469 Teuchos::TimeMonitor lcltimer( *timerMOp_ );
2470#endif
2471 count_ApplyM_+= MVT::GetNumberVecs(**i);
2472
2473 RCP<MV> helperMV = MVT::Clone(**i,MVT::GetNumberVecs(**i));
2474 OPT::Apply(*MOp_,**i,*helperMV);
2475 MauxVecs_.push_back(helperMV);
2476 }
2477 }
2478
2479 // If the solver has been initialized, V is not necessarily orthogonal to new auxiliary vectors
2480 if (numAuxVecs_ > 0 && initialized_) {
2481 initialized_ = false;
2482 }
2483
2484 if (om_->isVerbosity( Debug ) ) {
2485 // Check almost everything here
2486 CheckList chk;
2487 chk.checkQ = true;
2488 om_->print( Debug, accuracyCheck(chk, ": after setAuxVecs()") );
2489 }
2490 }
2491
2492
2494 // compute/return residual M-norms
2495 template <class ScalarType, class MV, class OP>
2496 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
2498 if (Rnorms_current_ == false) {
2499 // Update the residual norms
2500 if(computeAllRes_)
2501 {
2502 std::vector<int> curind(curDim_);
2503 for(int i=0; i<curDim_; i++)
2504 curind[i] = i;
2505
2506 RCP<const MV> locR = MVT::CloneView(*R_,curind);
2507 std::vector<ScalarType> locNorms(curDim_);
2508 orthman_->norm(*locR,locNorms);
2509
2510 for(int i=0; i<curDim_; i++)
2511 Rnorms_[i] = locNorms[i];
2512 for(int i=curDim_+1; i<blockSize_*numBlocks_; i++)
2513 Rnorms_[i] = NANVAL;
2514
2515 Rnorms_current_ = true;
2516 locNorms.resize(blockSize_);
2517 return locNorms;
2518 }
2519 else
2520 orthman_->norm(*R_,Rnorms_);
2521 Rnorms_current_ = true;
2522 }
2523 else if(computeAllRes_)
2524 {
2525 std::vector<ScalarType> locNorms(blockSize_);
2526 for(int i=0; i<blockSize_; i++)
2527 locNorms[i] = Rnorms_[i];
2528 return locNorms;
2529 }
2530
2531 return Rnorms_;
2532 }
2533
2534
2536 // compute/return residual 2-norms
2537 template <class ScalarType, class MV, class OP>
2538 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>
2540 if (R2norms_current_ == false) {
2541 // Update the residual 2-norms
2542 if(computeAllRes_)
2543 {
2544 std::vector<int> curind(curDim_);
2545 for(int i=0; i<curDim_; i++)
2546 curind[i] = i;
2547
2548 RCP<const MV> locR = MVT::CloneView(*R_,curind);
2549 std::vector<ScalarType> locNorms(curDim_);
2550 MVT::MvNorm(*locR,locNorms);
2551
2552 for(int i=0; i<curDim_; i++)
2553 {
2554 R2norms_[i] = locNorms[i];
2555 }
2556 for(int i=curDim_+1; i<blockSize_*numBlocks_; i++)
2557 R2norms_[i] = NANVAL;
2558
2559 R2norms_current_ = true;
2560 locNorms.resize(blockSize_);
2561 return locNorms;
2562 }
2563 else
2564 MVT::MvNorm(*R_,R2norms_);
2565 R2norms_current_ = true;
2566 }
2567 else if(computeAllRes_)
2568 {
2569 std::vector<ScalarType> locNorms(blockSize_);
2570 for(int i=0; i<blockSize_; i++)
2571 locNorms[i] = R2norms_[i];
2572 return locNorms;
2573 }
2574
2575 return R2norms_;
2576 }
2577
2579 // Set a new StatusTest for the solver.
2580 template <class ScalarType, class MV, class OP>
2581 void TraceMinBase<ScalarType,MV,OP>::setStatusTest(RCP<StatusTest<ScalarType,MV,OP> > test) {
2582 TEUCHOS_TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument,
2583 "Anasazi::TraceMinBase::setStatusTest() was passed a null StatusTest.");
2584 tester_ = test;
2585 }
2586
2588 // Get the current StatusTest used by the solver.
2589 template <class ScalarType, class MV, class OP>
2590 RCP<StatusTest<ScalarType,MV,OP> > TraceMinBase<ScalarType,MV,OP>::getStatusTest() const {
2591 return tester_;
2592 }
2593
2595 // Print the current status of the solver
2596 template <class ScalarType, class MV, class OP>
2597 void
2599 {
2600 using std::endl;
2601
2602 os.setf(std::ios::scientific, std::ios::floatfield);
2603 os.precision(6);
2604 os <<endl;
2605 os <<"================================================================================" << endl;
2606 os << endl;
2607 os <<" TraceMinBase Solver Status" << endl;
2608 os << endl;
2609 os <<"The solver is "<<(initialized_ ? "initialized." : "not initialized.") << endl;
2610 os <<"The number of iterations performed is " <<iter_<<endl;
2611 os <<"The block size is " << blockSize_<<endl;
2612 os <<"The number of blocks is " << numBlocks_<<endl;
2613 os <<"The current basis size is " << curDim_<<endl;
2614 os <<"The number of auxiliary vectors is "<< numAuxVecs_ << endl;
2615 os <<"The number of operations Op*x is "<<count_ApplyOp_<<endl;
2616 os <<"The number of operations M*x is "<<count_ApplyM_<<endl;
2617
2618 os.setf(std::ios_base::right, std::ios_base::adjustfield);
2619
2620 if (initialized_) {
2621 os << endl;
2622 os <<"CURRENT EIGENVALUE ESTIMATES "<<endl;
2623 os << std::setw(20) << "Eigenvalue"
2624 << std::setw(20) << "Residual(M)"
2625 << std::setw(20) << "Residual(2)"
2626 << endl;
2627 os <<"--------------------------------------------------------------------------------"<<endl;
2628 for (int i=0; i<blockSize_; ++i) {
2629 os << std::setw(20) << theta_[i];
2630 if (Rnorms_current_) os << std::setw(20) << Rnorms_[i];
2631 else os << std::setw(20) << "not current";
2632 if (R2norms_current_) os << std::setw(20) << R2norms_[i];
2633 else os << std::setw(20) << "not current";
2634 os << endl;
2635 }
2636 }
2637 os <<"================================================================================" << endl;
2638 os << endl;
2639 }
2640
2642 template <class ScalarType, class MV, class OP>
2644 {
2645 ScalarType currentTrace = ZERO;
2646
2647 for(int i=0; i < blockSize_; i++)
2648 currentTrace += theta_[i];
2649
2650 return currentTrace;
2651 }
2652
2654 template <class ScalarType, class MV, class OP>
2655 bool TraceMinBase<ScalarType,MV,OP>::traceLeveled()
2656 {
2657 ScalarType ratioOfChange = traceThresh_;
2658
2659 if(previouslyLeveled_)
2660 {
2661 om_->stream(Debug) << "The trace already leveled, so we're not going to check it again\n";
2662 return true;
2663 }
2664
2665 ScalarType currentTrace = getTrace();
2666
2667 om_->stream(Debug) << "The current trace is " << currentTrace << std::endl;
2668
2669 // Compute the ratio of the change
2670 // We seek the point where the trace has leveled off
2671 // It should be reasonably safe to shift at this point
2672 if(previousTrace_ != ZERO)
2673 {
2674 om_->stream(Debug) << "The previous trace was " << previousTrace_ << std::endl;
2675
2676 ratioOfChange = std::abs(previousTrace_-currentTrace)/std::abs(previousTrace_);
2677 om_->stream(Debug) << "The ratio of change is " << ratioOfChange << std::endl;
2678 }
2679
2680 previousTrace_ = currentTrace;
2681
2682 if(ratioOfChange < traceThresh_)
2683 {
2684 previouslyLeveled_ = true;
2685 return true;
2686 }
2687 return false;
2688 }
2689
2691 // Compute the residual of each CLUSTER of eigenvalues
2692 // This is important for selecting the Ritz shifts
2693 template <class ScalarType, class MV, class OP>
2694 std::vector<ScalarType> TraceMinBase<ScalarType,MV,OP>::getClusterResids()
2695 {
2696 int nvecs;
2697
2698 if(computeAllRes_)
2699 nvecs = curDim_;
2700 else
2701 nvecs = blockSize_;
2702
2703 getRes2Norms();
2704
2705 std::vector<ScalarType> clusterResids(nvecs);
2706 std::vector<int> clusterIndices;
2707 if(considerClusters_)
2708 {
2709 for(int i=0; i < nvecs; i++)
2710 {
2711 // test for cluster
2712 if(clusterIndices.empty() || (theta_[i-1] + R2norms_[i-1] >= theta_[i] - R2norms_[i]))
2713 {
2714 // Add to cluster
2715 if(!clusterIndices.empty()) om_->stream(Debug) << theta_[i-1] << " is in a cluster with " << theta_[i] << " because " << theta_[i-1] + R2norms_[i-1] << " >= " << theta_[i] - R2norms_[i] << std::endl;
2716 clusterIndices.push_back(i);
2717 }
2718 // Cluster completed
2719 else
2720 {
2721 om_->stream(Debug) << theta_[i-1] << " is NOT in a cluster with " << theta_[i] << " because " << theta_[i-1] + R2norms_[i-1] << " < " << theta_[i] - R2norms_[i] << std::endl;
2722 ScalarType totalRes = ZERO;
2723 for(size_t j=0; j < clusterIndices.size(); j++)
2724 totalRes += (R2norms_[clusterIndices[j]]*R2norms_[clusterIndices[j]]);
2725
2726 // If the smallest magnitude value of this sign is in a cluster with the
2727 // largest magnitude cluster of this sign, it is not safe for the smallest
2728 // eigenvalue to use a shift
2729 if(theta_[clusterIndices[0]] < 0 && theta_[i] < 0)
2730 negSafeToShift_ = true;
2731 else if(theta_[clusterIndices[0]] > 0 && theta_[i] > 0)
2732 posSafeToShift_ = true;
2733
2734 for(size_t j=0; j < clusterIndices.size(); j++)
2735 clusterResids[clusterIndices[j]] = sqrt(totalRes);
2736
2737 clusterIndices.clear();
2738 clusterIndices.push_back(i);
2739 }
2740 }
2741
2742 // Handle last cluster
2743 ScalarType totalRes = ZERO;
2744 for(size_t j=0; j < clusterIndices.size(); j++)
2745 totalRes += R2norms_[clusterIndices[j]];
2746 for(size_t j=0; j < clusterIndices.size(); j++)
2747 clusterResids[clusterIndices[j]] = totalRes;
2748 }
2749 else
2750 {
2751 for(int j=0; j < nvecs; j++)
2752 clusterResids[j] = R2norms_[j];
2753 }
2754
2755 return clusterResids;
2756 }
2757
2759 // Compute the Ritz shifts based on the Ritz values and residuals
2760 // A note on shifting: if the matrix is indefinite, you NEED to use a large block size
2761 // TODO: resids[i] on its own is unsafe for the generalized EVP
2762 // See "A Parallel Implementation of the Trace Minimization Eigensolver"
2763 // by Eloy Romero and Jose E. Roman
2764 template <class ScalarType, class MV, class OP>
2765 void TraceMinBase<ScalarType,MV,OP>::computeRitzShifts(const std::vector<ScalarType>& clusterResids)
2766 {
2767 std::vector<ScalarType> thetaMag(theta_);
2768 bool traceHasLeveled = traceLeveled();
2769
2770 // Compute the magnitude of the eigenvalues
2771 for(int i=0; i<blockSize_; i++)
2772 thetaMag[i] = std::abs(thetaMag[i]);
2773
2774 // Test whether it is safe to shift
2775 // TODO: Add residual shift option
2776 // Note: If you choose single shift with an indefinite matrix, you're gonna have a bad time...
2777 // Note: This is not safe for indefinite matrices, and I don't even know that it CAN be made safe.
2778 // There just isn't any theoretical reason it should work.
2779 // TODO: It feels like this should be different for TraceMin-Base; we should be able to use all eigenvalues in the current subspace to determine whether we have a cluster
2780 if(whenToShift_ == ALWAYS_SHIFT || (whenToShift_ == SHIFT_WHEN_TRACE_LEVELS && traceHasLeveled))
2781 {
2782 // Set the shift to the largest safe shift
2783 if(howToShift_ == LARGEST_CONVERGED_SHIFT)
2784 {
2785 for(int i=0; i<blockSize_; i++)
2786 ritzShifts_[i] = largestSafeShift_;
2787 }
2788 // Set the shifts to the Ritz values
2789 else if(howToShift_ == RITZ_VALUES_SHIFT)
2790 {
2791 ritzShifts_[0] = theta_[0];
2792
2793 // If we're using mulitple shifts, set them to EACH Ritz value.
2794 // Otherwise, only use the smallest
2795 if(useMultipleShifts_)
2796 {
2797 for(int i=1; i<blockSize_; i++)
2798 ritzShifts_[i] = theta_[i];
2799 }
2800 else
2801 {
2802 for(int i=1; i<blockSize_; i++)
2803 ritzShifts_[i] = ritzShifts_[0];
2804 }
2805 }
2806 else if(howToShift_ == EXPERIMENTAL_SHIFT)
2807 {
2808 ritzShifts_[0] = std::max(largestSafeShift_,theta_[0]-clusterResids[0]);
2809 for(int i=1; i<blockSize_; i++)
2810 {
2811 ritzShifts_[i] = std::max(ritzShifts_[i-1],theta_[i]-clusterResids[i]);
2812 }
2813 }
2814 // Use Dr. Sameh's original shifting strategy
2815 else if(howToShift_ == ADJUSTED_RITZ_SHIFT)
2816 {
2817 om_->stream(Debug) << "\nSeeking a shift for theta[0]=" << thetaMag[0] << std::endl;
2818
2819 // This is my adjustment. If all eigenvalues are in a single cluster, it's probably a bad idea to shift the smallest one.
2820 // If all of your eigenvalues are in one cluster, it's either way to early to shift or your subspace is too small
2821 if((theta_[0] > 0 && posSafeToShift_) || (theta_[0] < 0 && negSafeToShift_) || considerClusters_ == false)
2822 {
2823 // Initialize with a conservative shift, either the biggest safe shift or the eigenvalue adjusted by its cluster's residual
2824 ritzShifts_[0] = std::max(largestSafeShift_,thetaMag[0]-clusterResids[0]);
2825
2826 om_->stream(Debug) << "Initializing with a conservative shift, either the most positive converged eigenvalue ("
2827 << largestSafeShift_ << ") or the eigenvalue adjusted by the residual (" << thetaMag[0] << "-"
2828 << clusterResids[0] << ").\n";
2829
2830 // If this eigenvalue is NOT in a cluster, do an aggressive shift
2831 if(R2norms_[0] == clusterResids[0])
2832 {
2833 ritzShifts_[0] = thetaMag[0];
2834 om_->stream(Debug) << "Since this eigenvalue is NOT in a cluster, we can use the eigenvalue itself as a shift: ritzShifts[0]=" << ritzShifts_[0] << std::endl;
2835 }
2836 else
2837 om_->stream(Debug) << "This eigenvalue is in a cluster, so it would not be safe to use the eigenvalue itself as a shift\n";
2838 }
2839 else
2840 {
2841 if(largestSafeShift_ > std::abs(ritzShifts_[0]))
2842 {
2843 om_->stream(Debug) << "Initializing with a conservative shift...the most positive converged eigenvalue: " << largestSafeShift_ << std::endl;
2844 ritzShifts_[0] = largestSafeShift_;
2845 }
2846 else
2847 om_->stream(Debug) << "Using the previous value of ritzShifts[0]=" << ritzShifts_[0];
2848
2849 }
2850
2851 om_->stream(Debug) << "ritzShifts[0]=" << ritzShifts_[0] << std::endl;
2852
2853 if(useMultipleShifts_)
2854 {
2856 // Compute shifts for other eigenvalues
2857 for(int i=1; i < blockSize_; i++)
2858 {
2859 om_->stream(Debug) << "\nSeeking a shift for theta[" << i << "]=" << thetaMag[i] << std::endl;
2860
2861 // If the previous shift was aggressive and we are not in a cluster, do an aggressive shift
2862 if(ritzShifts_[i-1] == thetaMag[i-1] && i < blockSize_-1 && thetaMag[i] < thetaMag[i+1] - clusterResids[i+1])
2863 {
2864 ritzShifts_[i] = thetaMag[i];
2865 om_->stream(Debug) << "Using an aggressive shift: ritzShifts_[" << i << "]=" << ritzShifts_[i] << std::endl;
2866 }
2867 else
2868 {
2869 if(ritzShifts_[0] > std::abs(ritzShifts_[i]))
2870 {
2871 om_->stream(Debug) << "It was unsafe to use the aggressive shift. Choose the shift used by theta[0]="
2872 << thetaMag[0] << ": ritzShifts[0]=" << ritzShifts_[0] << std::endl;
2873
2874 // Choose a conservative shift, that of the smallest positive eigenvalue
2875 ritzShifts_[i] = ritzShifts_[0];
2876 }
2877 else
2878 om_->stream(Debug) << "It was unsafe to use the aggressive shift. We will use the shift from the previous iteration: " << ritzShifts_[i] << std::endl;
2879
2880 om_->stream(Debug) << "Check whether any less conservative shifts would work (such as the biggest eigenvalue outside of the cluster, namely theta[ell] < "
2881 << thetaMag[i] << "-" << clusterResids[i] << " (" << thetaMag[i] - clusterResids[i] << ")\n";
2882
2883 // If possible, choose a less conservative shift, that of the biggest eigenvalue outside of the cluster
2884 for(int ell=0; ell < i; ell++)
2885 {
2886 if(thetaMag[ell] < thetaMag[i] - clusterResids[i])
2887 {
2888 ritzShifts_[i] = thetaMag[ell];
2889 om_->stream(Debug) << "ritzShifts_[" << i << "]=" << ritzShifts_[i] << " is valid\n";
2890 }
2891 else
2892 break;
2893 }
2894 } // end else
2895
2896 om_->stream(Debug) << "ritzShifts[" << i << "]=" << ritzShifts_[i] << std::endl;
2897 } // end for
2898 } // end if(useMultipleShifts_)
2899 else
2900 {
2901 for(int i=1; i<blockSize_; i++)
2902 ritzShifts_[i] = ritzShifts_[0];
2903 }
2904 } // end if(howToShift_ == "Adjusted Ritz Values")
2905 } // end if(whenToShift_ == "Always" || (whenToShift_ == "After Trace Levels" && traceHasLeveled))
2906
2907 // Set the correct sign
2908 for(int i=0; i<blockSize_; i++)
2909 {
2910 if(theta_[i] < 0)
2911 ritzShifts_[i] = -abs(ritzShifts_[i]);
2912 else
2913 ritzShifts_[i] = abs(ritzShifts_[i]);
2914 }
2915 }
2916
2918 template <class ScalarType, class MV, class OP>
2919 std::vector<ScalarType> TraceMinBase<ScalarType,MV,OP>::computeTol()
2920 {
2921 ScalarType temp1;
2922 std::vector<ScalarType> tolerances(blockSize_);
2923
2924 for(int i=0; i < blockSize_-1; i++)
2925 {
2926 if(std::abs(theta_[blockSize_-1]) != std::abs(ritzShifts_[i]))
2927 temp1 = std::abs(theta_[i]-ritzShifts_[i])/std::abs(std::abs(theta_[blockSize_-1])-std::abs(ritzShifts_[i]));
2928 else
2929 temp1 = ZERO;
2930
2931 // TODO: The min and max tolerances should not be hard coded
2932 // Neither should the maximum number of iterations
2933 tolerances[i] = std::min(temp1*temp1,0.5);
2934 }
2935
2936 if(blockSize_ > 1)
2937 tolerances[blockSize_-1] = tolerances[blockSize_-2];
2938
2939 return tolerances;
2940 }
2941
2943 template <class ScalarType, class MV, class OP>
2944 void TraceMinBase<ScalarType,MV,OP>::solveSaddlePointProblem(RCP<MV> Delta)
2945 {
2946#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
2947 Teuchos::TimeMonitor lcltimer( *timerSaddle_ );
2948#endif
2949
2950 // This case can arise when looking for the largest eigenpairs
2951 if(Op_ == Teuchos::null)
2952 {
2953 // dense solver
2954 Teuchos::SerialDenseSolver<int,ScalarType> My_Solver;
2955
2956 // Schur complement
2957 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclS = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(blockSize_,blockSize_) );
2958
2959 if(computeAllRes_)
2960 {
2961 // Get the valid indices of X
2962 std::vector<int> curind(blockSize_);
2963 for(int i=0; i<blockSize_; i++)
2964 curind[i] = i;
2965
2966 // Get a view of MX
2967 RCP<const MV> lclMX = MVT::CloneView(*MX_, curind);
2968
2969 // form S = X' M^2 X
2970 MVT::MvTransMv(ONE,*lclMX,*lclMX,*lclS);
2971
2972 // compute the inverse of S
2973 My_Solver.setMatrix(lclS);
2974 My_Solver.invert();
2975
2976 // Delta = X - MX/S
2977 RCP<const MV> lclX = MVT::CloneView(*X_, curind);
2978 MVT::Assign(*lclX,*Delta);
2979 MVT::MvTimesMatAddMv( -ONE, *lclMX, *lclS, ONE, *Delta);
2980 }
2981 else
2982 {
2983 // form S = X' M^2 X
2984 MVT::MvTransMv(ONE,*MX_,*MX_,*lclS);
2985
2986 // compute the inverse of S
2987 My_Solver.setMatrix(lclS);
2988 My_Solver.invert();
2989
2990 // Delta = X - MX/S
2991 MVT::Assign(*X_,*Delta);
2992 MVT::MvTimesMatAddMv( -ONE, *MX_, *lclS, ONE, *Delta);
2993 }
2994 }
2995 else
2996 {
2997 std::vector<int> order(curDim_);
2998 std::vector<ScalarType> tempvec(blockSize_);
2999// RCP<BasicSort<MagnitudeType> > sorter = rcp( new BasicSort<MagnitudeType>("SR") );
3000
3001 // Stores the residual of each CLUSTER of eigenvalues
3002 std::vector<ScalarType> clusterResids;
3003
3004/* // Sort the eigenvalues in ascending order for the Ritz shift selection
3005 sorter->sort(theta_, Teuchos::rcpFromRef(order), curDim_); // don't catch exception
3006
3007 // Apply the same ordering to the residual norms
3008 getRes2Norms();
3009 for (int i=0; i<blockSize_; i++)
3010 tempvec[i] = R2norms_[order[i]];
3011 R2norms_ = tempvec;*/
3012
3013 // Compute the residual of each CLUSTER of eigenvalues
3014 // This is important for selecting the Ritz shifts
3015 clusterResids = getClusterResids();
3016
3017/* // Sort the eigenvalues based on what the user wanted
3018 sm_->sort(theta_, Teuchos::rcpFromRef(order), blockSize_);
3019
3020 // Apply the same ordering to the residual norms and cluster residuals
3021 for (int i=0; i<blockSize_; i++)
3022 tempvec[i] = R2norms_[order[i]];
3023 R2norms_ = tempvec;
3024
3025 for (int i=0; i<blockSize_; i++)
3026 tempvec[i] = clusterResids[order[i]];
3027 clusterResids = tempvec;*/
3028
3029 // Compute the Ritz shifts
3030 computeRitzShifts(clusterResids);
3031
3032 // Compute the tolerances for the inner solves
3033 std::vector<ScalarType> tolerances = computeTol();
3034
3035 for(int i=0; i<blockSize_; i++)
3036 {
3037 om_->stream(IterationDetails) << "Choosing Ritz shifts...theta[" << i << "]="
3038 << theta_[i] << ", resids[" << i << "]=" << R2norms_[i] << ", clusterResids[" << i << "]=" << clusterResids[i]
3039 << ", ritzShifts[" << i << "]=" << ritzShifts_[i] << ", and tol[" << i << "]=" << tolerances[i] << std::endl;
3040 }
3041
3042 // Set the Ritz shifts for the solver
3043 ritzOp_->setRitzShifts(ritzShifts_);
3044
3045 // Set the inner stopping tolerance
3046 // This uses the Ritz values to determine when to stop
3047 ritzOp_->setInnerTol(tolerances);
3048
3049 // Solve the saddle point problem
3050 if(saddleSolType_ == PROJECTED_KRYLOV_SOLVER)
3051 {
3052 if(Prec_ != Teuchos::null)
3053 solveSaddleProjPrec(Delta);
3054 else
3055 solveSaddleProj(Delta);
3056 }
3057 else if(saddleSolType_ == SCHUR_COMPLEMENT_SOLVER)
3058 {
3059 if(Z_ == Teuchos::null || MVT::GetNumberVecs(*Z_) != blockSize_)
3060 {
3061 // We do NOT want Z to be 0, because that could result in stagnation
3062 // I know it's tempting to take out the MvRandom, but seriously, don't do it.
3063 Z_ = MVT::Clone(*X_,blockSize_);
3064 MVT::MvRandom(*Z_);
3065 }
3066 solveSaddleSchur(Delta);
3067 }
3068 else if(saddleSolType_ == BD_PREC_MINRES)
3069 {
3070 solveSaddleBDPrec(Delta);
3071// Delta->describe(*(Teuchos::VerboseObjectBase::getDefaultOStream()),Teuchos::VERB_EXTREME);
3072 }
3073 else if(saddleSolType_ == HSS_PREC_GMRES)
3074 {
3075 solveSaddleHSSPrec(Delta);
3076 }
3077 else
3078 std::cout << "Invalid saddle solver type\n";
3079 }
3080 }
3081
3083 // Solve the saddle point problem using projected minres
3084 // TODO: We should be able to choose KX or -R as RHS.
3085 template <class ScalarType, class MV, class OP>
3086 void TraceMinBase<ScalarType,MV,OP>::solveSaddleProj (RCP<MV> Delta) const
3087 {
3088 RCP<TraceMinProjRitzOp<ScalarType,MV,OP> > projOp;
3089
3090 if(computeAllRes_)
3091 {
3092 // Get the valid indices of X
3093 std::vector<int> curind(blockSize_);
3094 for(int i=0; i<blockSize_; i++)
3095 curind[i] = i;
3096
3097 RCP<const MV> locMX = MVT::CloneView(*MX_, curind);
3098
3099 // We should really project out the converged eigenvectors too
3100 if(projectAllVecs_)
3101 {
3102 if(projectLockedVecs_ && numAuxVecs_ > 0)
3103 projOp = rcp( new TraceMinProjRitzOp<ScalarType,MV,OP>(ritzOp_,locMX,orthman_,auxVecs_) );
3104 else
3105 projOp = rcp( new TraceMinProjRitzOp<ScalarType,MV,OP>(ritzOp_,locMX,orthman_) );
3106 }
3107 else
3108 projOp = rcp( new TraceMinProjRitzOp<ScalarType,MV,OP>(ritzOp_,locMX) );
3109
3110 // Remember, Delta0 must equal 0
3111 // This ensures B-orthogonality between Delta and X
3112 MVT::MvInit(*Delta);
3113
3114 if(useRHSR_)
3115 {
3116 RCP<const MV> locR = MVT::CloneView(*R_, curind);
3117 projOp->ApplyInverse(*locR, *Delta);
3118 }
3119 else
3120 {
3121 RCP<const MV> locKX = MVT::CloneView(*KX_, curind);
3122 projOp->ApplyInverse(*locKX, *Delta);
3123 }
3124 }
3125 else
3126 {
3127 // We should really project out the converged eigenvectors too
3128 if(projectAllVecs_)
3129 {
3130 if(projectLockedVecs_ && numAuxVecs_ > 0)
3131 projOp = rcp( new TraceMinProjRitzOp<ScalarType,MV,OP>(ritzOp_,MX_,orthman_,auxVecs_) );
3132 else
3133 projOp = rcp( new TraceMinProjRitzOp<ScalarType,MV,OP>(ritzOp_,MX_,orthman_) );
3134 }
3135 else
3136 projOp = rcp( new TraceMinProjRitzOp<ScalarType,MV,OP>(ritzOp_,MX_) );
3137
3138 // Remember, Delta0 must equal 0
3139 // This ensures B-orthogonality between Delta and X
3140 MVT::MvInit(*Delta);
3141
3142 if(useRHSR_) {
3143 projOp->ApplyInverse(*R_, *Delta);
3144 }
3145 else {
3146 projOp->ApplyInverse(*KX_, *Delta);
3147 }
3148 }
3149 }
3150
3152 // TODO: Fix preconditioning
3153 template <class ScalarType, class MV, class OP>
3154 void TraceMinBase<ScalarType,MV,OP>::solveSaddleProjPrec (RCP<MV> Delta) const
3155 {
3156 // If we don't have Belos installed, we can't use TraceMinProjRitzOpWithPrec
3157 // Of course, this problem will be detected in the constructor and an exception will be thrown
3158 // This is only here to make sure the code compiles properly
3159#ifdef HAVE_ANASAZI_BELOS
3160 RCP<TraceMinProjRitzOpWithPrec<ScalarType,MV,OP> > projOp;
3161
3162 if(computeAllRes_)
3163 {
3164 int dimension;
3165 if(projectAllVecs_)
3166 dimension = curDim_;
3167 else
3168 dimension = blockSize_;
3169
3170 // Get the valid indices of X
3171 std::vector<int> curind(dimension);
3172 for(int i=0; i<dimension; i++)
3173 curind[i] = i;
3174
3175 RCP<const MV> locMX = MVT::CloneView(*MX_, curind);
3176
3177 // We should really project out the converged eigenvectors too
3178 if(projectAllVecs_)
3179 {
3180 if(projectLockedVecs_ && numAuxVecs_ > 0)
3181 projOp = rcp( new TraceMinProjRitzOpWithPrec<ScalarType,MV,OP>(ritzOp_,locMX,orthman_,auxVecs_) );
3182 else
3183 projOp = rcp( new TraceMinProjRitzOpWithPrec<ScalarType,MV,OP>(ritzOp_,locMX,orthman_) );
3184 }
3185 else
3186 projOp = rcp( new TraceMinProjRitzOpWithPrec<ScalarType,MV,OP>(ritzOp_,locMX) );
3187
3188 // Remember, Delta0 must equal 0
3189 // This ensures B-orthogonality between Delta and X
3190 MVT::MvInit(*Delta);
3191
3192 std::vector<int> dimind(blockSize_);
3193 for(int i=0; i<blockSize_; i++)
3194 dimind[i] = i;
3195
3196 if(useRHSR_)
3197 {
3198 RCP<const MV> locR = MVT::CloneView(*R_, dimind);
3199 projOp->ApplyInverse(*locR, *Delta);
3200 MVT::MvScale(*Delta, -ONE);
3201 }
3202 else
3203 {
3204 RCP<const MV> locKX = MVT::CloneView(*KX_, dimind);
3205 projOp->ApplyInverse(*locKX, *Delta);
3206 }
3207 }
3208 else
3209 {
3210 // We should really project out the converged eigenvectors too
3211 if(projectAllVecs_)
3212 {
3213 if(projectLockedVecs_ && numAuxVecs_ > 0)
3214 projOp = rcp( new TraceMinProjRitzOpWithPrec<ScalarType,MV,OP>(ritzOp_,MX_,orthman_,auxVecs_) );
3215 else
3216 projOp = rcp( new TraceMinProjRitzOpWithPrec<ScalarType,MV,OP>(ritzOp_,MX_,orthman_) );
3217 }
3218 else
3219 projOp = rcp( new TraceMinProjRitzOpWithPrec<ScalarType,MV,OP>(ritzOp_,MX_) );
3220
3221 // Remember, Delta0 must equal 0
3222 // This ensures B-orthogonality between Delta and X
3223 MVT::MvInit(*Delta);
3224
3225 if(useRHSR_)
3226 {
3227 projOp->ApplyInverse(*R_, *Delta);
3228 MVT::MvScale(*Delta,-ONE);
3229 }
3230 else
3231 projOp->ApplyInverse(*KX_, *Delta);
3232 }
3233#endif
3234 }
3235
3237 // TODO: We can hold the Schur complement constant in later iterations
3238 // TODO: Make sure we're using the preconditioner correctly
3239 template <class ScalarType, class MV, class OP>
3240 void TraceMinBase<ScalarType,MV,OP>::solveSaddleSchur (RCP<MV> Delta) const
3241 {
3242 // dense solver
3243 Teuchos::SerialDenseSolver<int,ScalarType> My_Solver;
3244
3245 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclL;
3246 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclS = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(blockSize_,blockSize_) );
3247
3248 if(computeAllRes_)
3249 {
3250 // Get the valid indices of X
3251 std::vector<int> curind(blockSize_);
3252 for(int i=0; i<blockSize_; i++)
3253 curind[i] = i;
3254
3255 // Z = K \ MX
3256 // Why would I have wanted to set Z <- X? I want to leave Z's previous value alone...
3257 RCP<const MV> lclMX = MVT::CloneView(*MX_, curind);
3258
3259#ifdef USE_APPLY_INVERSE
3260 Op_->ApplyInverse(*lclMX,*Z_);
3261#else
3262 ritzOp_->ApplyInverse(*lclMX,*Z_);
3263#endif
3264
3265 // form S = X' M Z
3266 MVT::MvTransMv(ONE,*Z_,*lclMX,*lclS);
3267
3268 // solve S L = I
3269 My_Solver.setMatrix(lclS);
3270 My_Solver.invert();
3271 lclL = lclS;
3272
3273 // Delta = X - Z L
3274 RCP<const MV> lclX = MVT::CloneView(*X_, curind);
3275 MVT::Assign(*lclX,*Delta);
3276 MVT::MvTimesMatAddMv( -ONE, *Z_, *lclL, ONE, *Delta);
3277 }
3278 else
3279 {
3280 // Z = K \ MX
3281#ifdef USE_APPLY_INVERSE
3282 Op_->ApplyInverse(*MX_,*Z_);
3283#else
3284 ritzOp_->ApplyInverse(*MX_,*Z_);
3285#endif
3286
3287 // form S = X' M Z
3288 MVT::MvTransMv(ONE,*Z_,*MX_,*lclS);
3289
3290 // solve S L = I
3291 My_Solver.setMatrix(lclS);
3292 My_Solver.invert();
3293 lclL = lclS;
3294
3295 // Delta = X - Z L
3296 MVT::Assign(*X_,*Delta);
3297 MVT::MvTimesMatAddMv( -ONE, *Z_, *lclL, ONE, *Delta);
3298 }
3299 }
3300
3301
3303 // TODO: We can hold the Schur complement constant in later iterations
3304 template <class ScalarType, class MV, class OP>
3305 void TraceMinBase<ScalarType,MV,OP>::solveSaddleBDPrec (RCP<MV> Delta) const
3306 {
3307 RCP<MV> locKX, locMX;
3308 if(computeAllRes_)
3309 {
3310 std::vector<int> curind(blockSize_);
3311 for(int i=0; i<blockSize_; i++)
3312 curind[i] = i;
3313 locKX = MVT::CloneViewNonConst(*KX_, curind);
3314 locMX = MVT::CloneViewNonConst(*MX_, curind);
3315 }
3316 else
3317 {
3318 locKX = KX_;
3319 locMX = MX_;
3320 }
3321
3322 // Create the operator [A BX; X'B 0]
3323 RCP<saddle_op_type> sadOp = rcp(new saddle_op_type(ritzOp_,locMX));
3324
3325 // Create the RHS [AX; 0]
3326 RCP<saddle_container_type> sadRHS = rcp(new saddle_container_type(locKX));
3327
3328// locKX->describe(*(Teuchos::VerboseObjectBase::getDefaultOStream()),Teuchos::VERB_EXTREME);
3329
3330// locMX->describe(*(Teuchos::VerboseObjectBase::getDefaultOStream()),Teuchos::VERB_EXTREME);
3331
3332 // Create the solution vector [Delta; L]
3333 MVT::MvInit(*Delta);
3334 RCP<saddle_container_type> sadSol = rcp(new saddle_container_type(Delta));
3335
3336 // Create a minres solver
3337 RCP<PseudoBlockMinres<ScalarType,saddle_container_type,saddle_op_type > > sadSolver;
3338 if(Prec_ != Teuchos::null)
3339 {
3340 RCP<saddle_op_type> sadPrec = rcp(new saddle_op_type(ritzOp_->getPrec(),locMX,BD_PREC));
3341 sadSolver = rcp(new PseudoBlockMinres<ScalarType,saddle_container_type,saddle_op_type>(sadOp, sadPrec));
3342 }
3343 else {
3344 sadSolver = rcp(new PseudoBlockMinres<ScalarType,saddle_container_type,saddle_op_type>(sadOp));
3345 }
3346
3347 // Set the tolerance for the minres solver
3348 std::vector<ScalarType> tol;
3349 ritzOp_->getInnerTol(tol);
3350 sadSolver->setTol(tol);
3351
3352 // Set the maximum number of iterations
3353 sadSolver->setMaxIter(ritzOp_->getMaxIts());
3354
3355 // Set the solution vector
3356 sadSolver->setSol(sadSol);
3357
3358 // Set the RHS
3359 sadSolver->setRHS(sadRHS);
3360
3361 // Solve the saddle point problem
3362 sadSolver->solve();
3363 }
3364
3365
3366
3368 // TODO: We can hold the Schur complement constant in later iterations
3369 template <class ScalarType, class MV, class OP>
3370 void TraceMinBase<ScalarType,MV,OP>::solveSaddleHSSPrec (RCP<MV> Delta) const
3371 {
3372#ifdef HAVE_ANASAZI_BELOS
3373 typedef Belos::LinearProblem<ScalarType,saddle_container_type,saddle_op_type> LP;
3374 typedef Belos::PseudoBlockGmresSolMgr<ScalarType,saddle_container_type,saddle_op_type> GmresSolMgr;
3375
3376 RCP<MV> locKX, locMX;
3377 if(computeAllRes_)
3378 {
3379 std::vector<int> curind(blockSize_);
3380 for(int i=0; i<blockSize_; i++)
3381 curind[i] = i;
3382 locKX = MVT::CloneViewNonConst(*KX_, curind);
3383 locMX = MVT::CloneViewNonConst(*MX_, curind);
3384 }
3385 else
3386 {
3387 locKX = KX_;
3388 locMX = MX_;
3389 }
3390
3391 // Create the operator [A BX; X'B 0]
3392 RCP<saddle_op_type> sadOp = rcp(new saddle_op_type(ritzOp_,locMX,NONSYM));
3393
3394 // Create the RHS [AX; 0]
3395 RCP<saddle_container_type> sadRHS = rcp(new saddle_container_type(locKX));
3396
3397 // Create the solution vector [Delta; L]
3398 MVT::MvInit(*Delta);
3399 RCP<saddle_container_type> sadSol = rcp(new saddle_container_type(Delta));
3400
3401 // Create a parameter list for the gmres solver
3402 RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList());
3403
3404 // Set the tolerance for the gmres solver
3405 std::vector<ScalarType> tol;
3406 ritzOp_->getInnerTol(tol);
3407 pl->set("Convergence Tolerance",tol[0]);
3408
3409 // Set the maximum number of iterations
3410 // TODO: Come back to this
3411 pl->set("Maximum Iterations", ritzOp_->getMaxIts());
3412 pl->set("Num Blocks", ritzOp_->getMaxIts());
3413
3414 // Set the block size
3415 // TODO: Come back to this
3416 // TODO: This breaks the code right now, presumably because of a MVT cloneview issue.
3417 pl->set("Block Size", blockSize_);
3418
3419 // Set the verbosity of gmres
3420// pl->set("Verbosity", Belos::IterationDetails + Belos::StatusTestDetails + Belos::Debug);
3421// pl->set("Output Frequency", 1);
3422
3423 // Create the linear problem
3424 RCP<LP> problem = rcp(new LP(sadOp,sadSol,sadRHS));
3425
3426 // Set the preconditioner
3427 if(Prec_ != Teuchos::null)
3428 {
3429 RCP<saddle_op_type> sadPrec = rcp(new saddle_op_type(ritzOp_->getPrec(),locMX,HSS_PREC,alpha_));
3430 problem->setLeftPrec(sadPrec);
3431 }
3432
3433 // Set the problem
3434 problem->setProblem();
3435
3436 // Create a minres solver
3437 RCP<GmresSolMgr> sadSolver = rcp(new GmresSolMgr(problem,pl)) ;
3438
3439 // Solve the saddle point problem
3440 sadSolver->solve();
3441#else
3442 std::cout << "No Belos. This is bad\n";
3443#endif
3444 }
3445
3446
3447
3449 // Compute KK := V'KV
3450 // We only compute the NEW elements
3451 template <class ScalarType, class MV, class OP>
3452 void TraceMinBase<ScalarType,MV,OP>::computeKK()
3453 {
3454 // Get the valid indices of V
3455 std::vector<int> curind(curDim_);
3456 for(int i=0; i<curDim_; i++)
3457 curind[i] = i;
3458
3459 // Get a pointer to the valid parts of V
3460 RCP<const MV> lclV = MVT::CloneView(*V_,curind);
3461
3462 // Get the valid indices of KV
3463 curind.resize(blockSize_);
3464 for(int i=0; i<blockSize_; i++)
3465 curind[i] = curDim_-blockSize_+i;
3466 RCP<const MV> lclKV = MVT::CloneView(*KV_,curind);
3467
3468 // Get a pointer to the valid part of KK
3469 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > lclKK =
3470 rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,blockSize_,0,curDim_-blockSize_) );
3471
3472 // KK := V'KV
3473 MVT::MvTransMv(ONE,*lclV,*lclKV,*lclKK);
3474
3475 // We only constructed the upper triangular part of the matrix, but that's okay because KK is symmetric!
3476 for(int r=0; r<curDim_; r++)
3477 {
3478 for(int c=0; c<r; c++)
3479 {
3480 (*KK_)(r,c) = (*KK_)(c,r);
3481 }
3482 }
3483 }
3484
3486 // Compute the eigenpairs of KK, i.e. the Ritz pairs
3487 template <class ScalarType, class MV, class OP>
3488 void TraceMinBase<ScalarType,MV,OP>::computeRitzPairs()
3489 {
3490 // Get a pointer to the valid part of KK
3491 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > lclKK =
3492 rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*KK_,curDim_,curDim_) );
3493
3494 // Get a pointer to the valid part of ritzVecs
3495 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > lclRV =
3496 rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*ritzVecs_,curDim_,curDim_) );
3497
3498 // Compute Ritz pairs from KK
3499 {
3500#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
3501 Teuchos::TimeMonitor lcltimer( *timerDS_ );
3502#endif
3503 int rank = curDim_;
3504 Utils::directSolver(curDim_, *lclKK, Teuchos::null, *lclRV, theta_, rank, 10);
3505 // we want all ritz values back
3506 // TODO: This probably should not be an ortho failure
3507 TEUCHOS_TEST_FOR_EXCEPTION(rank != curDim_,TraceMinBaseOrthoFailure,
3508 "Anasazi::TraceMinBase::computeRitzPairs(): Failed to compute all eigenpairs of KK.");
3509 }
3510
3511 // Sort ritz pairs
3512 {
3513#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
3514 Teuchos::TimeMonitor lcltimer( *timerSortEval_ );
3515#endif
3516
3517 std::vector<int> order(curDim_);
3518 //
3519 // sort the first curDim_ values in theta_
3520 if(useHarmonic_)
3521 {
3523 sm.sort(theta_, Teuchos::rcpFromRef(order), curDim_);
3524 }
3525 else
3526 {
3527 sm_->sort(theta_, Teuchos::rcpFromRef(order), curDim_); // don't catch exception
3528 }
3529 //
3530 // apply the same ordering to the primitive ritz vectors
3531 Utils::permuteVectors(order,*lclRV);
3532 }
3533 }
3534
3536 // Compute X := V evecs
3537 template <class ScalarType, class MV, class OP>
3538 void TraceMinBase<ScalarType,MV,OP>::computeX()
3539 {
3540#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
3541 Teuchos::TimeMonitor lcltimer( *timerLocal_ );
3542#endif
3543
3544 // Get the valid indices of V
3545 std::vector<int> curind(curDim_);
3546 for(int i=0; i<curDim_; i++)
3547 curind[i] = i;
3548
3549 // Get a pointer to the valid parts of V
3550 RCP<const MV> lclV = MVT::CloneView(*V_,curind);
3551
3552 if(computeAllRes_)
3553 {
3554 // Capture the relevant eigenvectors
3555 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > relevantEvecs =
3556 rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*ritzVecs_,curDim_,curDim_) );
3557
3558 // X <- lclV*S
3559 RCP<MV> lclX = MVT::CloneViewNonConst(*X_,curind);
3560 MVT::MvTimesMatAddMv( ONE, *lclV, *relevantEvecs, ZERO, *lclX );
3561 }
3562 else
3563 {
3564 // Capture the relevant eigenvectors
3565 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > relevantEvecs =
3566 rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*ritzVecs_,curDim_,blockSize_) );
3567
3568 // X <- lclV*S
3569 MVT::MvTimesMatAddMv( ONE, *lclV, *relevantEvecs, ZERO, *X_ );
3570 }
3571 }
3572
3574 // Compute KX := KV evecs and (if necessary) MX := MV evecs
3575 template <class ScalarType, class MV, class OP>
3576 void TraceMinBase<ScalarType,MV,OP>::updateKXMX()
3577 {
3578#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
3579 Teuchos::TimeMonitor lcltimer( *timerLocal_ );
3580#endif
3581
3582 // Get the valid indices of V
3583 std::vector<int> curind(curDim_);
3584 for(int i=0; i<curDim_; i++)
3585 curind[i] = i;
3586
3587 // Get pointers to the valid parts of V, KV, and MV (if necessary)
3588 RCP<const MV> lclV = MVT::CloneView(*V_,curind);
3589 RCP<const MV> lclKV = MVT::CloneView(*KV_,curind);
3590
3591 if(computeAllRes_)
3592 {
3593 // Capture the relevant eigenvectors
3594 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > relevantEvecs =
3595 rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*ritzVecs_,curDim_,curDim_) );
3596
3597 // Update KX and MX
3598 RCP<MV> lclKX = MVT::CloneViewNonConst(*KX_,curind);
3599 MVT::MvTimesMatAddMv( ONE, *lclKV, *relevantEvecs, ZERO, *lclKX );
3600 if(hasM_)
3601 {
3602 RCP<const MV> lclMV = MVT::CloneView(*MV_,curind);
3603 RCP<MV> lclMX = MVT::CloneViewNonConst(*MX_,curind);
3604 MVT::MvTimesMatAddMv( ONE, *lclMV, *relevantEvecs, ZERO, *lclMX );
3605 }
3606 }
3607 else
3608 {
3609 // Capture the relevant eigenvectors
3610 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > relevantEvecs =
3611 rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*ritzVecs_,curDim_,blockSize_) );
3612
3613 // Update KX and MX
3614 MVT::MvTimesMatAddMv( ONE, *lclKV, *relevantEvecs, ZERO, *KX_ );
3615 if(hasM_)
3616 {
3617 RCP<const MV> lclMV = MVT::CloneView(*MV_,curind);
3618 MVT::MvTimesMatAddMv( ONE, *lclMV, *relevantEvecs, ZERO, *MX_ );
3619 }
3620 }
3621 }
3622
3624 // Update the residual R := KX - MX*T
3625 template <class ScalarType, class MV, class OP>
3626 void TraceMinBase<ScalarType,MV,OP>::updateResidual () {
3627#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
3628 Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
3629#endif
3630
3631 if(computeAllRes_)
3632 {
3633 // Get the valid indices of X
3634 std::vector<int> curind(curDim_);
3635 for(int i=0; i<curDim_; i++)
3636 curind[i] = i;
3637
3638 // Holds MX*T
3639 RCP<MV> MXT = MVT::CloneCopy(*MX_,curind);
3640
3641 // Holds the relevant part of theta
3642 std::vector<ScalarType> locTheta(curDim_);
3643 for(int i=0; i<curDim_; i++)
3644 locTheta[i] = theta_[i];
3645
3646 // Compute MX*T
3647 MVT::MvScale(*MXT,locTheta);
3648
3649 // form R <- KX - MX*T
3650 RCP<const MV> locKX = MVT::CloneView(*KX_,curind);
3651 RCP<MV> locR = MVT::CloneViewNonConst(*R_,curind);
3652 MVT::MvAddMv(ONE,*locKX,-ONE,*MXT,*locR);
3653 }
3654 else
3655 {
3656 // Holds MX*T
3657 RCP<MV> MXT = MVT::CloneCopy(*MX_);
3658
3659 // Holds the relevant part of theta
3660 std::vector<ScalarType> locTheta(blockSize_);
3661 for(int i=0; i<blockSize_; i++)
3662 locTheta[i] = theta_[i];
3663
3664 // Compute MX*T
3665 MVT::MvScale(*MXT,locTheta);
3666
3667 // form R <- KX - MX*T
3668 MVT::MvAddMv(ONE,*KX_,-ONE,*MXT,*R_);
3669 }
3670
3671 // R has been updated; mark the norms as out-of-date
3672 Rnorms_current_ = false;
3673 R2norms_current_ = false;
3674 }
3675
3676
3678 // Check accuracy, orthogonality, and other debugging stuff
3679 //
3680 // bools specify which tests we want to run (instead of running more than we actually care about)
3681 //
3682 // we don't bother checking the following because they are computed explicitly:
3683 // H == Prec*R
3684 // KH == K*H
3685 //
3686 //
3687 // checkV : V orthonormal
3688 // orthogonal to auxvecs
3689 // checkX : X orthonormal
3690 // orthogonal to auxvecs
3691 // checkMX: check MX == M*X
3692 // checkKX: check KX == K*X
3693 // checkH : H orthonormal
3694 // orthogonal to V and H and auxvecs
3695 // checkMH: check MH == M*H
3696 // checkR : check R orthogonal to X
3697 // checkQ : check that auxiliary vectors are actually orthonormal
3698 // checkKK: check that KK is symmetric in memory
3699 //
3700 // TODO:
3701 // add checkTheta
3702 //
3703 template <class ScalarType, class MV, class OP>
3704 std::string TraceMinBase<ScalarType,MV,OP>::accuracyCheck( const CheckList &chk, const std::string &where ) const
3705 {
3706 using std::endl;
3707
3708 std::stringstream os;
3709 os.precision(2);
3710 os.setf(std::ios::scientific, std::ios::floatfield);
3711
3712 os << " Debugging checks: iteration " << iter_ << where << endl;
3713
3714 // V and friends
3715 std::vector<int> lclind(curDim_);
3716 for (int i=0; i<curDim_; ++i) lclind[i] = i;
3717 RCP<const MV> lclV;
3718 if (initialized_) {
3719 lclV = MVT::CloneView(*V_,lclind);
3720 }
3721 if (chk.checkV && initialized_) {
3722 MagnitudeType err = orthman_->orthonormError(*lclV);
3723 os << " >> Error in V^H M V == I : " << err << endl;
3724 for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
3725 err = orthman_->orthogError(*lclV,*auxVecs_[i]);
3726 os << " >> Error in V^H M Q[" << i << "] == 0 : " << err << endl;
3727 }
3728 }
3729
3730 // X and friends
3731 RCP<const MV> lclX;
3732 if(initialized_)
3733 {
3734 if(computeAllRes_)
3735 lclX = MVT::CloneView(*X_,lclind);
3736 else
3737 lclX = X_;
3738 }
3739
3740 if (chk.checkX && initialized_) {
3741 MagnitudeType err = orthman_->orthonormError(*lclX);
3742 os << " >> Error in X^H M X == I : " << err << endl;
3743 for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
3744 err = orthman_->orthogError(*lclX,*auxVecs_[i]);
3745 os << " >> Error in X^H M Q[" << i << "] == 0 : " << err << endl;
3746 }
3747 }
3748 if (chk.checkMX && hasM_ && initialized_) {
3749 RCP<const MV> lclMX;
3750 if(computeAllRes_)
3751 lclMX = MVT::CloneView(*MX_,lclind);
3752 else
3753 lclMX = MX_;
3754
3755 MagnitudeType err = Utils::errorEquality(*lclX, *lclMX, MOp_);
3756 os << " >> Error in MX == M*X : " << err << endl;
3757 }
3758 if (Op_ != Teuchos::null && chk.checkKX && initialized_) {
3759 RCP<const MV> lclKX;
3760 if(computeAllRes_)
3761 lclKX = MVT::CloneView(*KX_,lclind);
3762 else
3763 lclKX = KX_;
3764
3765 MagnitudeType err = Utils::errorEquality(*lclX, *lclKX, Op_);
3766 os << " >> Error in KX == K*X : " << err << endl;
3767 }
3768
3769 // KK
3770 if (chk.checkKK && initialized_) {
3771 Teuchos::SerialDenseMatrix<int,ScalarType> curKK(curDim_,curDim_);
3772 if(Op_ != Teuchos::null) {
3773 RCP<MV> lclKV = MVT::Clone(*V_,curDim_);
3774 OPT::Apply(*Op_,*lclV,*lclKV);
3775 MVT::MvTransMv(ONE,*lclV,*lclKV,curKK);
3776 }
3777 else {
3778 MVT::MvTransMv(ONE,*lclV,*lclV,curKK);
3779 }
3780 Teuchos::SerialDenseMatrix<int,ScalarType> subKK(Teuchos::View,*KK_,curDim_,curDim_);
3781 curKK -= subKK;
3782 os << " >> Error in V^H K V == KK : " << curKK.normFrobenius() << endl;
3783
3784 Teuchos::SerialDenseMatrix<int,ScalarType> SDMerr(curDim_,curDim_);
3785 for (int j=0; j<curDim_; ++j) {
3786 for (int i=0; i<curDim_; ++i) {
3787 SDMerr(i,j) = subKK(i,j) - SCT::conjugate(subKK(j,i));
3788 }
3789 }
3790 os << " >> Error in KK - KK^H == 0 : " << SDMerr.normFrobenius() << endl;
3791 }
3792
3793 // Q
3794 if (chk.checkQ) {
3795 for (Array_size_type i=0; i<auxVecs_.size(); ++i) {
3796 MagnitudeType err = orthman_->orthonormError(*auxVecs_[i]);
3797 os << " >> Error in Q[" << i << "]^H M Q[" << i << "] == I : " << err << endl;
3798 for (Array_size_type j=i+1; j<auxVecs_.size(); ++j) {
3799 err = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
3800 os << " >> Error in Q[" << i << "]^H M Q[" << j << "] == 0 : " << err << endl;
3801 }
3802 }
3803 }
3804
3805 os << endl;
3806
3807 return os.str();
3808 }
3809
3810}} // End of namespace Anasazi
3811
3812#endif
3813
3814// End of file AnasaziTraceMinBase.hpp
Basic implementation of the Anasazi::SortManager class.
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
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.
Stores a set of vectors of the form [V; L] where V is a distributed multivector and L is a serialdens...
An operator of the form [A Y; Y' 0] where A is a sparse matrix and Y a multivector.
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.
An implementation of the Anasazi::SortManager that performs a collection of common sorting techniques...
void sort(std::vector< MagnitudeType > &evals, Teuchos::RCP< std::vector< int > > perm=Teuchos::null, int n=-1) const
Sort real eigenvalues, optionally returning the permutation vector.
The Eigensolver is a templated virtual base class that defines the basic interface that any eigensolv...
TraceMinBaseInitFailure is thrown when the TraceMinBase solver is unable to generate an initial itera...
TraceMinBaseOrthoFailure is thrown when the orthogonalization manager is unable to orthogonalize the ...
This is an abstract base class for the trace minimization eigensolvers.
TraceMinBase(const RCP< Eigenproblem< ScalarType, MV, OP > > &problem, const RCP< SortManager< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > &sorter, const RCP< OutputManager< ScalarType > > &printer, const RCP< StatusTest< ScalarType, MV, OP > > &tester, const RCP< MatOrthoManager< ScalarType, MV, OP > > &ortho, Teuchos::ParameterList &params)
TraceMinBase constructor with eigenproblem, solver utilities, and parameter list of solver options.
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...
std::vector< Value< ScalarType > > getRitzValues()
Get the Ritz values for the previous iteration.
void setStatusTest(RCP< StatusTest< ScalarType, MV, OP > > test)
Set a new StatusTest for the solver.
void resetNumIters()
Reset the iteration count.
RCP< StatusTest< ScalarType, MV, OP > > getStatusTest() const
Get the current StatusTest used by the solver.
RCP< const MV > getRitzVectors()
Get access to the current Ritz vectors.
std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > getRitzRes2Norms()
Get the 2-norms of the residuals.
TraceMinBaseState< 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 getBlockSize() const
Get the blocksize used by the iterative solver.
std::vector< int > getRitzIndex()
Get the index used for extracting individual Ritz vectors from getRitzVectors().
virtual ~TraceMinBase()
Anasazi::TraceMinBase destructor.
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace. For the trace minimization methods,...
void setAuxVecs(const Teuchos::Array< RCP< const MV > > &auxvecs)
Set the auxiliary vectors for the solver.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the eigenvalue problem.
void setBlockSize(int blockSize)
Set the blocksize.
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current eigenvectors and eigenvalues.
void setSize(int blockSize, int numBlocks)
Set the blocksize and number of blocks to be used by the iterative solver in solving this eigenproble...
bool isInitialized() const
Indicates whether the solver has been initialized or not.
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 initialize()
Initialize the solver with the initial vectors from the eigenproblem or random data.
Teuchos::Array< RCP< const MV > > getAuxVecs() const
Get the auxiliary vectors for the solver.
void iterate()
This method performs trace minimization iterations until the status test indicates the need to stop o...
int getNumIters() const
Get the current iteration count.
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 TraceMinBase state variables.
RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > KK
The current projected K matrix.
RCP< const MV > MX
The image of the current eigenvectors under M, or Teuchos::null if M was not specified.
RCP< const MV > KX
The image of the current eigenvectors under K.
ScalarType largestSafeShift
Largest safe shift.
int curDim
The current dimension of the solver.
RCP< const MV > MopV
The image of V under M, or Teuchos::null if M was not specified.
RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > RV
The current Ritz vectors.
bool isOrtho
Whether V has been projected and orthonormalized already.
RCP< const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > T
The current Ritz values. This vector is a copy of the internal data.
RCP< const MV > X
The current eigenvectors.
RCP< const MV > R
The current residual vectors.
int NEV
Number of unconverged eigenvalues.
RCP< const MV > KV
The image of V under K.
RCP< const std::vector< ScalarType > > ritzShifts
Current Ritz shifts.