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 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 #ifndef ANASAZI_STATUS_TEST_SPECTRANS_HPP
43 #define ANASAZI_STATUS_TEST_SPECTRANS_HPP
44 
45 #include "AnasaziTypes.hpp"
47 #include "AnasaziTraceMinBase.hpp"
48 
49 using Teuchos::RCP;
50 
51 namespace Anasazi {
52 namespace Experimental {
53 
54  template<class ScalarType, class MV, class OP>
55  class StatusTestSpecTrans : public StatusTest<ScalarType,MV,OP> {
56 
57  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
58  typedef MultiVecTraits<ScalarType,MV> MVT;
59  typedef OperatorTraits<ScalarType,MV,OP> OPT;
60 
61  public:
62 
63  // Constructor
64  StatusTestSpecTrans(MagnitudeType tol, int quorum = -1, ResType whichNorm = RES_2NORM, bool scaled = true, bool throwExceptionOnNan = true, const RCP<const OP> Mop = Teuchos::null);
65 
66  // Destructor
67  virtual ~StatusTestSpecTrans() {};
68 
69  // Check whether the test passed or failed
70  TestStatus checkStatus(Eigensolver<ScalarType,MV,OP> *solver);
71 
72  // Return the result of the most recent checkStatus call
73  TestStatus getStatus() const { return state_; }
74 
75  // Get the indices for the vectors that passed the test
76  std::vector<int> whichVecs() const { return ind_; }
77 
78  // Get the number of vectors that passed the test
79  int howMany() const { return ind_.size(); }
80 
81  void setQuorum (int quorum) {
82  state_ = Undefined;
83  quorum_ = quorum;
84  }
85 
86  int getQuorum() const { return quorum_; }
87 
88  void setTolerance(MagnitudeType tol)
89  {
90  state_ = Undefined;
91  tol_ = tol;
92  }
93 
94  MagnitudeType getTolerance() const { return tol_; }
95 
96  void setWhichNorm(ResType whichNorm)
97  {
98  state_ = Undefined;
99  whichNorm_ = whichNorm;
100  }
101 
102  ResType getWhichNorm() const { return whichNorm_; }
103 
104  void setScale(bool relscale)
105  {
106  state_ = Undefined;
107  scaled_ = relscale;
108  }
109 
110  bool getScale() const { return scaled_; }
111 
112  // Informs the status test that it should reset its internal configuration to the uninitialized state
113  void reset()
114  {
115  ind_.resize(0);
116  state_ = Undefined;
117  }
118 
119  // Clears the results of the last status test
120  void clearStatus() { reset(); };
121 
122  // Output formatted description of stopping test to output stream
123  std::ostream & print(std::ostream &os, int indent=0) const;
124 
125  private:
126  TestStatus state_;
127  MagnitudeType tol_;
128  std::vector<int> ind_;
129  int quorum_;
130  bool scaled_;
131  enum ResType whichNorm_;
132  bool throwExceptionOnNaN_;
133  RCP<const OP> M_;
134 
135  const MagnitudeType ONE;
136  };
137 
138 
139 
140  template <class ScalarType, class MV, class OP>
141  StatusTestSpecTrans<ScalarType,MV,OP>::StatusTestSpecTrans(MagnitudeType tol, int quorum, ResType whichNorm, bool scaled, bool throwExceptionOnNaN, const RCP<const OP> Mop)
142  : state_(Undefined),
143  tol_(tol),
144  quorum_(quorum),
145  scaled_(scaled),
146  whichNorm_(whichNorm),
147  throwExceptionOnNaN_(throwExceptionOnNaN),
148  M_(Mop),
149  ONE(Teuchos::ScalarTraits<MagnitudeType>::one())
150  {}
151 
152 
153 
154  template <class ScalarType, class MV, class OP>
155  TestStatus StatusTestSpecTrans<ScalarType,MV,OP>::checkStatus( Eigensolver<ScalarType,MV,OP>* solver )
156  {
157  typedef Teuchos::ScalarTraits<MagnitudeType> MT;
158  typedef TraceMinBase<ScalarType,MV,OP> TS;
159 
160  // Cast the eigensolver to a TraceMin solver
161  // Throw an exception if this fails
162  TS* tm_solver = dynamic_cast<TS*>(solver);
163  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!");
164 
165  // Get the residual Ax-\lambda Bx, which is not computed when there's a spectral transformation
166  // TraceMin computes Bx-1/\lambda Ax
167  TraceMinBaseState<ScalarType,MV> state = tm_solver->getState();
168 
169  size_t nvecs = state.ritzShifts->size();
170  std::vector<int> curind(nvecs);
171  for(size_t i=0; i<nvecs; i++)
172  curind[i] = i;
173 
174  RCP<const MV> locKX, locMX, locX;
175  RCP<MV> R;
176  locX = MVT::CloneView(*state.X,curind);
177  if(state.KX != Teuchos::null)
178  locKX = MVT::CloneView(*state.KX,curind);
179  else
180  locKX = locX;
181  if(state.MX != Teuchos::null)
182  locMX = MVT::CloneView(*state.MX,curind);
183  else
184  locMX = locX;
185  R = MVT::CloneCopy(*locKX,curind);
186 
187  std::vector<MagnitudeType> evals(nvecs);
188  for(size_t i=0; i<nvecs; i++)
189  evals[i] = ONE/(*state.T)[i];
190  MVT::MvScale(*R,evals);
191  MVT::MvAddMv(-ONE,*R,ONE,*locMX,*R);
192 
193  // Compute the norms
194  std::vector<MagnitudeType> res(nvecs);
195  switch (whichNorm_) {
196  case RES_2NORM:
197  {
198  MVT::MvNorm(*R,res);
199  break;
200  }
201  case RES_ORTH:
202  {
203  RCP<MV> MR = MVT::Clone(*R,nvecs);
204  OPT::Apply(*M_,*R,*MR);
205  MVT::MvDot(*R,*MR,res);
206  for(size_t i=0; i<nvecs; i++)
207  res[i] = MT::squareroot(res[i]);
208  break;
209  }
210  case RITZRES_2NORM:
211  {
212  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "The trace minimization solvers do not define a Ritz residual. Please choose a different residual type");
213  break;
214  }
215  }
216 
217  // if appropriate, scale the norms by the magnitude of the eigenvalue estimate
218  if(scaled_)
219  {
220  for(size_t i=0; i<nvecs; i++)
221  res[i] /= std::abs(evals[i]);
222  }
223 
224  // test the norms
225  ind_.resize(0);
226  for(size_t i=0; i<nvecs; i++)
227  {
228  TEUCHOS_TEST_FOR_EXCEPTION( MT::isnaninf(res[i]), ResNormNaNError,
229  "StatusTestSpecTrans::checkStatus(): residual norm is nan or inf" );
230  if(res[i] < tol_)
231  ind_.push_back(i);
232  }
233  int have = ind_.size();
234  int need = (quorum_ == -1) ? nvecs : quorum_;
235  state_ = (have >= need) ? Passed : Failed;
236  return state_;
237  }
238 
239 
240 
241  template <class ScalarType, class MV, class OP>
242  std::ostream& StatusTestSpecTrans<ScalarType,MV,OP>::print(std::ostream& os, int indent) const
243  {
244  std::string ind(indent,' ');
245  os << ind << "- StatusTestSpecTrans: ";
246  switch (state_) {
247  case Passed:
248  os << "Passed\n";
249  break;
250  case Failed:
251  os << "Failed\n";
252  break;
253  case Undefined:
254  os << "Undefined\n";
255  break;
256  }
257  os << ind << " (Tolerance, WhichNorm,Scaled,Quorum): "
258  << "(" << tol_;
259  switch (whichNorm_) {
260  case RES_ORTH:
261  os << ",RES_ORTH";
262  break;
263  case RES_2NORM:
264  os << ",RES_2NORM";
265  break;
266  case RITZRES_2NORM:
267  os << ",RITZRES_2NORM";
268  break;
269  }
270  os << "," << (scaled_ ? "true" : "false")
271  << "," << quorum_
272  << ")\n";
273 
274  if (state_ != Undefined) {
275  os << ind << " Which vectors: ";
276  if (ind_.size() > 0) {
277  for(size_t i=0; i<ind_.size(); i++) os << ind_[i] << " ";
278  os << std::endl;
279  }
280  else
281  os << "[empty]\n";
282  }
283  return os;
284  }
285 
286 }} // end of namespace
287 
288 #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.