Anasazi  Version of the Day
AnasaziStatusTestSpecTrans.hpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Anasazi: Block Eigensolvers Package
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25 //
26 // ***********************************************************************
27 // @HEADER
28 
29 #ifndef ANASAZI_STATUS_TEST_SPECTRANS_HPP
30 #define ANASAZI_STATUS_TEST_SPECTRANS_HPP
31 
32 #include "AnasaziTypes.hpp"
34 #include "AnasaziTraceMinBase.hpp"
35 
36 using Teuchos::RCP;
37 
38 namespace Anasazi {
39 namespace Experimental {
40 
41  template<class ScalarType, class MV, class OP>
42  class StatusTestSpecTrans : public StatusTest<ScalarType,MV,OP> {
43 
44  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
45  typedef MultiVecTraits<ScalarType,MV> MVT;
46  typedef OperatorTraits<ScalarType,MV,OP> OPT;
47 
48  public:
49 
50  // Constructor
51  StatusTestSpecTrans(MagnitudeType tol, int quorum = -1, ResType whichNorm = RES_2NORM, bool scaled = true, bool throwExceptionOnNan = true, const RCP<const OP> Mop = Teuchos::null);
52 
53  // Destructor
54  virtual ~StatusTestSpecTrans() {};
55 
56  // Check whether the test passed or failed
57  TestStatus checkStatus(Eigensolver<ScalarType,MV,OP> *solver);
58 
59  // Return the result of the most recent checkStatus call
60  TestStatus getStatus() const { return state_; }
61 
62  // Get the indices for the vectors that passed the test
63  std::vector<int> whichVecs() const { return ind_; }
64 
65  // Get the number of vectors that passed the test
66  int howMany() const { return ind_.size(); }
67 
68  void setQuorum (int quorum) {
69  state_ = Undefined;
70  quorum_ = quorum;
71  }
72 
73  int getQuorum() const { return quorum_; }
74 
75  void setTolerance(MagnitudeType tol)
76  {
77  state_ = Undefined;
78  tol_ = tol;
79  }
80 
81  MagnitudeType getTolerance() const { return tol_; }
82 
83  void setWhichNorm(ResType whichNorm)
84  {
85  state_ = Undefined;
86  whichNorm_ = whichNorm;
87  }
88 
89  ResType getWhichNorm() const { return whichNorm_; }
90 
91  void setScale(bool relscale)
92  {
93  state_ = Undefined;
94  scaled_ = relscale;
95  }
96 
97  bool getScale() const { return scaled_; }
98 
99  // Informs the status test that it should reset its internal configuration to the uninitialized state
100  void reset()
101  {
102  ind_.resize(0);
103  state_ = Undefined;
104  }
105 
106  // Clears the results of the last status test
107  void clearStatus() { reset(); };
108 
109  // Output formatted description of stopping test to output stream
110  std::ostream & print(std::ostream &os, int indent=0) const;
111 
112  private:
113  TestStatus state_;
114  MagnitudeType tol_;
115  std::vector<int> ind_;
116  int quorum_;
117  bool scaled_;
118  enum ResType whichNorm_;
119  bool throwExceptionOnNaN_;
120  RCP<const OP> M_;
121 
122  const MagnitudeType ONE;
123  };
124 
125 
126 
127  template <class ScalarType, class MV, class OP>
128  StatusTestSpecTrans<ScalarType,MV,OP>::StatusTestSpecTrans(MagnitudeType tol, int quorum, ResType whichNorm, bool scaled, bool throwExceptionOnNaN, const RCP<const OP> Mop)
129  : state_(Undefined),
130  tol_(tol),
131  quorum_(quorum),
132  scaled_(scaled),
133  whichNorm_(whichNorm),
134  throwExceptionOnNaN_(throwExceptionOnNaN),
135  M_(Mop),
136  ONE(Teuchos::ScalarTraits<MagnitudeType>::one())
137  {}
138 
139 
140 
141  template <class ScalarType, class MV, class OP>
142  TestStatus StatusTestSpecTrans<ScalarType,MV,OP>::checkStatus( Eigensolver<ScalarType,MV,OP>* solver )
143  {
144  typedef Teuchos::ScalarTraits<MagnitudeType> MT;
145  typedef TraceMinBase<ScalarType,MV,OP> TS;
146 
147  // Cast the eigensolver to a TraceMin solver
148  // Throw an exception if this fails
149  TS* tm_solver = dynamic_cast<TS*>(solver);
150  TEUCHOS_TEST_FOR_EXCEPTION(tm_solver == 0, std::invalid_argument, "The status test for spectral transformations currently only works for the trace minimization eigensolvers. Sorry!");
151 
152  // Get the residual Ax-\lambda Bx, which is not computed when there's a spectral transformation
153  // TraceMin computes Bx-1/\lambda Ax
154  TraceMinBaseState<ScalarType,MV> state = tm_solver->getState();
155 
156  size_t nvecs = state.ritzShifts->size();
157  std::vector<int> curind(nvecs);
158  for(size_t i=0; i<nvecs; i++)
159  curind[i] = i;
160 
161  RCP<const MV> locKX, locMX, locX;
162  RCP<MV> R;
163  locX = MVT::CloneView(*state.X,curind);
164  if(state.KX != Teuchos::null)
165  locKX = MVT::CloneView(*state.KX,curind);
166  else
167  locKX = locX;
168  if(state.MX != Teuchos::null)
169  locMX = MVT::CloneView(*state.MX,curind);
170  else
171  locMX = locX;
172  R = MVT::CloneCopy(*locKX,curind);
173 
174  std::vector<MagnitudeType> evals(nvecs);
175  for(size_t i=0; i<nvecs; i++)
176  evals[i] = ONE/(*state.T)[i];
177  MVT::MvScale(*R,evals);
178  MVT::MvAddMv(-ONE,*R,ONE,*locMX,*R);
179 
180  // Compute the norms
181  std::vector<MagnitudeType> res(nvecs);
182  switch (whichNorm_) {
183  case RES_2NORM:
184  {
185  MVT::MvNorm(*R,res);
186  break;
187  }
188  case RES_ORTH:
189  {
190  RCP<MV> MR = MVT::Clone(*R,nvecs);
191  OPT::Apply(*M_,*R,*MR);
192  MVT::MvDot(*R,*MR,res);
193  for(size_t i=0; i<nvecs; i++)
194  res[i] = MT::squareroot(res[i]);
195  break;
196  }
197  case RITZRES_2NORM:
198  {
199  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "The trace minimization solvers do not define a Ritz residual. Please choose a different residual type");
200  break;
201  }
202  }
203 
204  // if appropriate, scale the norms by the magnitude of the eigenvalue estimate
205  if(scaled_)
206  {
207  for(size_t i=0; i<nvecs; i++)
208  res[i] /= std::abs(evals[i]);
209  }
210 
211  // test the norms
212  ind_.resize(0);
213  for(size_t i=0; i<nvecs; i++)
214  {
215  TEUCHOS_TEST_FOR_EXCEPTION( MT::isnaninf(res[i]), ResNormNaNError,
216  "StatusTestSpecTrans::checkStatus(): residual norm is nan or inf" );
217  if(res[i] < tol_)
218  ind_.push_back(i);
219  }
220  int have = ind_.size();
221  int need = (quorum_ == -1) ? nvecs : quorum_;
222  state_ = (have >= need) ? Passed : Failed;
223  return state_;
224  }
225 
226 
227 
228  template <class ScalarType, class MV, class OP>
229  std::ostream& StatusTestSpecTrans<ScalarType,MV,OP>::print(std::ostream& os, int indent) const
230  {
231  std::string ind(indent,' ');
232  os << ind << "- StatusTestSpecTrans: ";
233  switch (state_) {
234  case Passed:
235  os << "Passed\n";
236  break;
237  case Failed:
238  os << "Failed\n";
239  break;
240  case Undefined:
241  os << "Undefined\n";
242  break;
243  }
244  os << ind << " (Tolerance, WhichNorm,Scaled,Quorum): "
245  << "(" << tol_;
246  switch (whichNorm_) {
247  case RES_ORTH:
248  os << ",RES_ORTH";
249  break;
250  case RES_2NORM:
251  os << ",RES_2NORM";
252  break;
253  case RITZRES_2NORM:
254  os << ",RITZRES_2NORM";
255  break;
256  }
257  os << "," << (scaled_ ? "true" : "false")
258  << "," << quorum_
259  << ")\n";
260 
261  if (state_ != Undefined) {
262  os << ind << " Which vectors: ";
263  if (ind_.size() > 0) {
264  for(size_t i=0; i<ind_.size(); i++) os << ind_[i] << " ";
265  os << std::endl;
266  }
267  else
268  os << "[empty]\n";
269  }
270  return os;
271  }
272 
273 }} // end of namespace
274 
275 #endif
Abstract base class for trace minimization eigensolvers.
ResType
Enumerated type used to specify which residual norm used by residual norm status tests.
TestStatus
Enumerated type used to pass back information from a StatusTest.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
A status test for testing the norm of the eigenvectors residuals.
Types and exceptions used within Anasazi solvers and interfaces.