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  {
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.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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.