Anasazi  Version of the Day
AnasaziStatusTestWithOrdering.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, 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 
30 #ifndef ANASAZI_STATUS_TEST_ORDEREDRESNORM_HPP
31 #define ANASAZI_STATUS_TEST_ORDEREDRESNORM_HPP
32 
40 #include "AnasaziStatusTest.hpp"
41 #include "Teuchos_ScalarTraits.hpp"
42 #include "Teuchos_LAPACK.hpp"
43 
67 namespace Anasazi {
68 
69 
70 template <class ScalarType, class MV, class OP>
71 class StatusTestWithOrdering : public StatusTest<ScalarType,MV,OP> {
72 
73  private:
74  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
75  typedef Teuchos::ScalarTraits<MagnitudeType> MT;
76 
77  public:
78 
80 
81 
83  StatusTestWithOrdering(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test, Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sorter, int quorum = -1);
84 
88 
90 
91 
96 
98  TestStatus getStatus() const { return state_; }
99 
101 
106  std::vector<int> whichVecs() const {
107  return ind_;
108  }
109 
111  int howMany() const {
112  return ind_.size();
113  }
114 
116 
118 
119 
125  void setQuorum(int quorum) {
126  state_ = Undefined;
127  quorum_ = quorum;
128  }
129 
132  int getQuorum() const {
133  return quorum_;
134  }
135 
137 
139 
140 
146  void reset() {
147  ind_.resize(0);
148  state_ = Undefined;
149  test_->reset();
150  }
151 
153 
158  void clearStatus() {
159  ind_.resize(0);
160  state_ = Undefined;
161  test_->clearStatus();
162  }
163 
168  void setAuxVals(const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &vals) {
169  rvals_ = vals;
170  ivals_.resize(rvals_.size(),MT::zero());
171  state_ = Undefined;
172  }
173 
178  void setAuxVals(const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &rvals, const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &ivals) {
179  rvals_ = rvals;
180  ivals_ = ivals;
181  state_ = Undefined;
182  }
183 
188  void getAuxVals(std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &rvals, std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &ivals) const {
189  rvals = rvals_;
190  ivals = ivals_;
191  }
192 
194 
196 
197 
199  std::ostream& print(std::ostream& os, int indent = 0) const;
200 
202  private:
203  TestStatus state_;
204  std::vector<int> ind_;
205  int quorum_;
206  std::vector<MagnitudeType> rvals_, ivals_;
207  Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sorter_;
208  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test_;
209 };
210 
211 
212 template <class ScalarType, class MV, class OP>
213 StatusTestWithOrdering<ScalarType,MV,OP>::StatusTestWithOrdering(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test, Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sorter, int quorum)
214  : state_(Undefined), ind_(0), quorum_(quorum), rvals_(0), ivals_(0), sorter_(sorter), test_(test)
215 {
216  TEUCHOS_TEST_FOR_EXCEPTION(sorter_ == Teuchos::null, StatusTestError, "StatusTestWithOrdering::constructor() was passed null pointer for constituent SortManager.");
217  TEUCHOS_TEST_FOR_EXCEPTION(test_ == Teuchos::null, StatusTestError, "StatusTestWithOrdering::constructor() was passed null pointer for constituent StatusTest.");
218 }
219 
220 template <class ScalarType, class MV, class OP>
222 
223  // Algorithm
224  // we PASS iff the "most signficant" values of "all values" PASS
225  // "most significant" is measured by sorter
226  // auxilliary values are automatically PASSED
227  // constituent status test determines who else passes
228  // "all values" mean {auxilliary values} UNION {solver's ritz values}
229  //
230  // test_->checkStatus() calls the constituent status test
231  // cwhch = test_->whichVecs() gets the passing indices from the constituent test
232  // solval = solver->getRitzValues() gets the solver's ritz values
233  // allval = {solval auxval} contains all values
234  // sorter_->sort(allval,perm) sort all values (we just need the perm vector)
235  //
236  // allpass = {cwhch numsolval+1 ... numsolval+numAux}
237  // mostsig = {perm[1] ... perm[quorum]}
238  // whichVecs = mostsig INTERSECT allpass
239  // if size(whichVecs) >= quorum,
240  // PASS
241  // else
242  // FAIL
243  //
244  // finish: this needs to be better tested and revisited for efficiency.
245 
246  // call the constituent solver manager
247  test_->checkStatus(solver);
248  std::vector<int> cwhch( test_->whichVecs() );
249 
250  // get the ritzvalues from solver
251  std::vector<Value<ScalarType> > solval = solver->getRitzValues();
252  int numsolval = solval.size();
253  int numauxval = rvals_.size();
254  int numallval = numsolval + numauxval;
255 
256  if (numallval == 0) {
257  ind_.resize(0);
258  return Failed;
259  }
260 
261  // make space for all values, real and imaginary parts
262  std::vector<MagnitudeType> allvalr(numallval), allvali(numallval);
263  // separate the real and imaginary parts of solver ritz values
264  for (int i=0; i<numsolval; ++i) {
265  allvalr[i] = solval[i].realpart;
266  allvali[i] = solval[i].imagpart;
267  }
268  // put the auxiliary values at the end of the solver ritz values
269  std::copy(rvals_.begin(),rvals_.end(),allvalr.begin()+numsolval);
270  std::copy(ivals_.begin(),ivals_.end(),allvali.begin()+numsolval);
271 
272  // sort all values
273  std::vector<int> perm(numallval);
274  sorter_->sort(allvalr,allvali,Teuchos::rcpFromRef(perm),numallval);
275 
276  // make the set of passing values: allpass = {cwhch -1 ... -numauxval}
277  std::vector<int> allpass(cwhch.size() + numauxval);
278  std::copy(cwhch.begin(),cwhch.end(),allpass.begin());
279  for (int i=0; i<numauxval; i++) {
280  allpass[cwhch.size()+i] = -(i+1);
281  }
282 
283  // make list, with length quorum, of most significant values, if there are that many
284  int numsig = quorum_ < numallval ? quorum_ : numallval;
285  // int numsig = cwhch.size() + numauxval;
286  std::vector<int> mostsig(numsig);
287  for (int i=0; i<numsig; ++i) {
288  mostsig[i] = perm[i];
289  // if perm[i] >= numsolval, then it corresponds to the perm[i]-numsolval aux val
290  // the first aux val gets index -numauxval, the second -numauxval+1, and so forth
291  if (mostsig[i] >= numsolval) {
292  mostsig[i] = mostsig[i]-numsolval-numauxval;
293  }
294  }
295 
296  // who passed?
297  // to use set_intersection, ind_ must have room for the result, which will have size() <= min( allpass.size(), mostsig.size() )
298  // also, allpass and mostsig must be in ascending order; neither are
299  // lastly, the return from set_intersection points to the last element in the intersection, which tells us how big the intersection is
300  ind_.resize(numsig);
301  std::vector<int>::iterator end;
302  std::sort(mostsig.begin(),mostsig.end());
303  std::sort(allpass.begin(),allpass.end());
304  end = std::set_intersection(mostsig.begin(),mostsig.end(),allpass.begin(),allpass.end(),ind_.begin());
305  ind_.resize(end - ind_.begin());
306 
307  // did we pass, overall
308  if (ind_.size() >= (unsigned int)quorum_) {
309  state_ = Passed;
310  }
311  else {
312  state_ = Failed;
313  }
314  return state_;
315 }
316 
317 
318 template <class ScalarType, class MV, class OP>
319 std::ostream& StatusTestWithOrdering<ScalarType,MV,OP>::print(std::ostream& os, int indent) const {
320  // build indent string
321  std::string ind(indent,' ');
322  // print header
323  os << ind << "- StatusTestWithOrdering: ";
324  switch (state_) {
325  case Passed:
326  os << "Passed" << std::endl;
327  break;
328  case Failed:
329  os << "Failed" << std::endl;
330  break;
331  case Undefined:
332  os << "Undefined" << std::endl;
333  break;
334  }
335  // print parameters, namely, quorum
336  os << ind << " Quorum: " << quorum_ << std::endl;
337  // print any auxilliary values
338  os << ind << " Auxiliary values: ";
339  if (rvals_.size() > 0) {
340  for (unsigned int i=0; i<rvals_.size(); i++) {
341  os << "(" << rvals_[i] << ", " << ivals_[i] << ") ";
342  }
343  os << std::endl;
344  }
345  else {
346  os << "[empty]" << std::endl;
347  }
348  // print which vectors have passed
349  if (state_ != Undefined) {
350  os << ind << " Which vectors: ";
351  if (ind_.size() > 0) {
352  for (unsigned int i=0; i<ind_.size(); i++) os << ind_[i] << " ";
353  os << std::endl;
354  }
355  else {
356  os << "[empty]" << std::endl;
357  }
358  }
359  // print constituent test
360  test_->print(os,indent+2);
361  return os;
362 }
363 
364 
365 } // end of Anasazi namespace
366 
367 #endif /* ANASAZI_STATUS_TEST_ORDEREDRESNORM_HPP */
virtual std::vector< Value< ScalarType > > getRitzValues()=0
Get the Ritz values from the previous iteration.
TestStatus getStatus() const
Return the result of the most recent checkStatus call, or undefined if it has not been run...
std::ostream & print(std::ostream &os, int indent=0) const
Output formatted description of stopping test to output stream.
Exception thrown to signal error in a status test during Anasazi::StatusTest::checkStatus().
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...
TestStatus checkStatus(Eigensolver< ScalarType, MV, OP > *solver)
void getAuxVals(std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &rvals, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &ivals) const
Get the auxiliary eigenvalues.
std::vector< int > whichVecs() const
Get the indices for the vectors that passed the test.
StatusTestWithOrdering(Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test, Teuchos::RCP< SortManager< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > sorter, int quorum=-1)
Constructor.
void setAuxVals(const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &rvals, const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &ivals)
Set the auxiliary eigenvalues.
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
int howMany() const
Get the number of vectors that passed the test.
void reset()
Informs the status test that it should reset its internal configuration to the uninitialized state...
void setAuxVals(const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &vals)
Set the auxiliary eigenvalues.
Anasazi&#39;s templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
void clearStatus()
Clears the results of the last status test.
Common interface of stopping criteria for Anasazi&#39;s solvers.
The Eigensolver is a templated virtual base class that defines the basic interface that any eigensolv...
Declaration and definition of Anasazi::StatusTest.