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