00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030 #ifndef ANASAZI_STATUS_TEST_ORDEREDRESNORM_HPP
00031 #define ANASAZI_STATUS_TEST_ORDEREDRESNORM_HPP
00032
00040 #include "AnasaziStatusTest.hpp"
00041 #include "Teuchos_ScalarTraits.hpp"
00042 #include "Teuchos_LAPACK.hpp"
00043
00067 namespace Anasazi {
00068
00069
00070 template <class ScalarType, class MV, class OP>
00071 class StatusTestWithOrdering : public StatusTest<ScalarType,MV,OP> {
00072
00073 private:
00074 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00075 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00076
00077 public:
00078
00080
00081
00083 StatusTestWithOrdering(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test, Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sorter, int quorum = -1);
00084
00086 virtual ~StatusTestWithOrdering() {};
00088
00090
00091
00095 TestStatus checkStatus( Eigensolver<ScalarType,MV,OP>* solver );
00096
00098 TestStatus getStatus() const { return state_; }
00099
00101
00106 std::vector<int> whichVecs() const {
00107 return ind_;
00108 }
00109
00111 int howMany() const {
00112 return ind_.size();
00113 }
00114
00116
00118
00119
00125 void setQuorum(int quorum) {
00126 state_ = Undefined;
00127 quorum_ = quorum;
00128 }
00129
00132 int getQuorum() const {
00133 return quorum_;
00134 }
00135
00137
00139
00140
00141
00146 void reset() {
00147 ind_.resize(0);
00148 state_ = Undefined;
00149 test_->reset();
00150 }
00151
00153
00158 void clearStatus() {
00159 ind_.resize(0);
00160 state_ = Undefined;
00161 test_->clearStatus();
00162 }
00163
00168 void setAuxVals(const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &vals) {
00169 rvals_ = vals;
00170 ivals_.resize(rvals_.size(),MT::zero());
00171 state_ = Undefined;
00172 }
00173
00178 void setAuxVals(const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &rvals, const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &ivals) {
00179 rvals_ = rvals;
00180 ivals_ = ivals;
00181 state_ = Undefined;
00182 }
00183
00188 void getAuxVals(std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &rvals, std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &ivals) const {
00189 rvals = rvals_;
00190 ivals = ivals_;
00191 }
00192
00194
00196
00197
00199 std::ostream& print(std::ostream& os, int indent = 0) const;
00200
00202 private:
00203 TestStatus state_;
00204 std::vector<int> ind_;
00205 int quorum_;
00206 std::vector<MagnitudeType> rvals_, ivals_;
00207 Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sorter_;
00208 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test_;
00209 };
00210
00211
00212 template <class ScalarType, class MV, class OP>
00213 StatusTestWithOrdering<ScalarType,MV,OP>::StatusTestWithOrdering(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test, Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sorter, int quorum)
00214 : state_(Undefined), ind_(0), quorum_(quorum), rvals_(0), ivals_(0), sorter_(sorter), test_(test)
00215 {
00216 TEST_FOR_EXCEPTION(sorter_ == Teuchos::null, StatusTestError, "StatusTestWithOrdering::constructor() was passed null pointer for constituent SortManager.");
00217 TEST_FOR_EXCEPTION(test_ == Teuchos::null, StatusTestError, "StatusTestWithOrdering::constructor() was passed null pointer for constituent StatusTest.");
00218 }
00219
00220 template <class ScalarType, class MV, class OP>
00221 TestStatus StatusTestWithOrdering<ScalarType,MV,OP>::checkStatus( Eigensolver<ScalarType,MV,OP>* solver ) {
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247 test_->checkStatus(solver);
00248 std::vector<int> cwhch( test_->whichVecs() );
00249
00250
00251 std::vector<Value<ScalarType> > solval = solver->getRitzValues();
00252 int numsolval = solval.size();
00253 int numauxval = rvals_.size();
00254 int numallval = numsolval + numauxval;
00255
00256 if (numallval == 0) {
00257 ind_.resize(0);
00258 return Failed;
00259 }
00260
00261
00262 std::vector<MagnitudeType> allvalr(numallval), allvali(numallval);
00263
00264 for (int i=0; i<numsolval; ++i) {
00265 allvalr[i] = solval[i].realpart;
00266 allvali[i] = solval[i].imagpart;
00267 }
00268
00269 std::copy(rvals_.begin(),rvals_.end(),allvalr.begin()+numsolval);
00270 std::copy(ivals_.begin(),ivals_.end(),allvali.begin()+numsolval);
00271
00272
00273 std::vector<int> perm(numallval);
00274 sorter_->sort(allvalr,allvali,Teuchos::rcpFromRef(perm),numallval);
00275
00276
00277 std::vector<int> allpass(cwhch.size() + numauxval);
00278 std::copy(cwhch.begin(),cwhch.end(),allpass.begin());
00279 for (int i=0; i<numauxval; i++) {
00280 allpass[cwhch.size()+i] = -(i+1);
00281 }
00282
00283
00284 int numsig = quorum_ < numallval ? quorum_ : numallval;
00285 std::vector<int> mostsig(numsig);
00286 for (int i=0; i<numsig; ++i) {
00287 mostsig[i] = perm[i];
00288
00289
00290 if (mostsig[i] >= numsolval) {
00291 mostsig[i] = mostsig[i]-numsolval-numauxval;
00292 }
00293 }
00294
00295
00296
00297
00298
00299 ind_.resize(numsig);
00300 std::vector<int>::iterator end;
00301 std::sort(mostsig.begin(),mostsig.end());
00302 std::sort(allpass.begin(),allpass.end());
00303 end = std::set_intersection(mostsig.begin(),mostsig.end(),allpass.begin(),allpass.end(),ind_.begin());
00304 ind_.resize(end - ind_.begin());
00305
00306
00307 if (ind_.size() >= (unsigned int)quorum_) {
00308 state_ = Passed;
00309 }
00310 else {
00311 state_ = Failed;
00312 }
00313 return state_;
00314 }
00315
00316
00317 template <class ScalarType, class MV, class OP>
00318 std::ostream& StatusTestWithOrdering<ScalarType,MV,OP>::print(std::ostream& os, int indent) const {
00319
00320 std::string ind(indent,' ');
00321
00322 os << ind << "- StatusTestWithOrdering: ";
00323 switch (state_) {
00324 case Passed:
00325 os << "Passed" << std::endl;
00326 break;
00327 case Failed:
00328 os << "Failed" << std::endl;
00329 break;
00330 case Undefined:
00331 os << "Undefined" << std::endl;
00332 break;
00333 }
00334
00335 os << ind << " Quorum: " << quorum_ << std::endl;
00336
00337 os << ind << " Auxiliary values: ";
00338 if (rvals_.size() > 0) {
00339 for (unsigned int i=0; i<rvals_.size(); i++) {
00340 os << "(" << rvals_[i] << ", " << ivals_[i] << ") ";
00341 }
00342 os << std::endl;
00343 }
00344 else {
00345 os << "[empty]" << std::endl;
00346 }
00347
00348 if (state_ != Undefined) {
00349 os << ind << " Which vectors: ";
00350 if (ind_.size() > 0) {
00351 for (unsigned int i=0; i<ind_.size(); i++) os << ind_[i] << " ";
00352 os << std::endl;
00353 }
00354 else {
00355 os << "[empty]" << std::endl;
00356 }
00357 }
00358
00359 test_->print(os,indent+2);
00360 return os;
00361 }
00362
00363
00364 }
00365
00366 #endif