Anasazi  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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  {
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.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
TestStatus
Enumerated type used to pass back information from a StatusTest.
A status test for testing the norm of the eigenvectors residuals.
Types and exceptions used within Anasazi solvers and interfaces.