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_RESNORM_HPP
00031 #define ANASAZI_STATUS_TEST_RESNORM_HPP
00032
00039 #include "AnasaziStatusTest.hpp"
00040 #include "Teuchos_ScalarTraits.hpp"
00041 namespace Anasazi {
00042
00044
00045
00054 class ResNormNaNError : public AnasaziError {public:
00055 ResNormNaNError(const std::string& what_arg) : AnasaziError(what_arg)
00056 {}};
00057
00059
00076 template <class ScalarType, class MV, class OP>
00077 class StatusTestResNorm : public StatusTest<ScalarType,MV,OP> {
00078
00079 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00080
00081 public:
00082
00086 enum ResType {
00087 RES_ORTH,
00088 RES_2NORM,
00089 RITZRES_2NORM
00090 };
00091
00093
00094
00096 StatusTestResNorm(typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol, int quorum = -1, ResType whichNorm = RES_ORTH, bool scaled = true, bool throwExceptionOnNaN = true);
00097
00099 virtual ~StatusTestResNorm() {};
00101
00103
00104
00108 TestStatus checkStatus( Eigensolver<ScalarType,MV,OP>* solver );
00109
00111 TestStatus getStatus() const { return state_; }
00112
00114 std::vector<int> whichVecs() const {
00115 return ind_;
00116 }
00117
00119 int howMany() const {
00120 return ind_.size();
00121 }
00122
00124
00126
00127
00133 void setQuorum(int quorum) {
00134 state_ = Undefined;
00135 quorum_ = quorum;
00136 }
00137
00140 int getQuorum() const {
00141 return quorum_;
00142 }
00143
00147 void setTolerance(typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol) {
00148 state_ = Undefined;
00149 tol_ = tol;
00150 }
00151
00153 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType getTolerance() {return tol_;}
00154
00159 void setWhichNorm(ResType whichNorm) {
00160 state_ = Undefined;
00161 whichNorm_ = whichNorm;
00162 }
00163
00165 ResType getWhichNorm() {return whichNorm_;}
00166
00170 void setScale(bool relscale) {
00171 state_ = Undefined;
00172 scaled_ = relscale;
00173 }
00174
00176 bool getScale() {return scaled_;}
00178
00180
00181
00182
00187 void reset() {
00188 ind_.resize(0);
00189 state_ = Undefined;
00190 }
00191
00193
00198 void clearStatus() {
00199 ind_.resize(0);
00200 state_ = Undefined;
00201 }
00202
00204
00206
00207
00209 std::ostream& print(std::ostream& os, int indent = 0) const;
00210
00212 private:
00213 TestStatus state_;
00214 MagnitudeType tol_;
00215 std::vector<int> ind_;
00216 int quorum_;
00217 bool scaled_;
00218 ResType whichNorm_;
00219 bool throwExceptionOnNaN_;
00220 };
00221
00222
00223 template <class ScalarType, class MV, class OP>
00224 StatusTestResNorm<ScalarType,MV,OP>::StatusTestResNorm(typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol, int quorum, ResType whichNorm, bool scaled, bool throwExceptionOnNaN)
00225 : state_(Undefined), tol_(tol), quorum_(quorum), scaled_(scaled), whichNorm_(whichNorm), throwExceptionOnNaN_(throwExceptionOnNaN)
00226 {}
00227
00228 template <class ScalarType, class MV, class OP>
00229 TestStatus StatusTestResNorm<ScalarType,MV,OP>::checkStatus( Eigensolver<ScalarType,MV,OP>* solver )
00230 {
00231 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00232
00233 std::vector<MagnitudeType> res;
00234
00235
00236
00237 std::vector<Value<ScalarType> > vals = solver->getRitzValues();
00238 switch (whichNorm_) {
00239 case RES_2NORM:
00240 res = solver->getRes2Norms();
00241
00242 vals.resize(res.size());
00243 break;
00244 case RES_ORTH:
00245 res = solver->getResNorms();
00246
00247 vals.resize(res.size());
00248 break;
00249 case RITZRES_2NORM:
00250 res = solver->getRitzRes2Norms();
00251 break;
00252 }
00253
00254
00255 if (scaled_) {
00256 Teuchos::LAPACK<int,MagnitudeType> lapack;
00257
00258 for (unsigned int i=0; i<res.size(); i++) {
00259 MagnitudeType tmp = lapack.LAPY2(vals[i].realpart,vals[i].imagpart);
00260
00261 if ( tmp != MT::zero() ) {
00262 res[i] /= tmp;
00263 }
00264 }
00265 }
00266
00267
00268 int have = 0;
00269 ind_.resize(res.size());
00270 for (unsigned int i=0; i<res.size(); i++) {
00271 TEST_FOR_EXCEPTION( MT::isnaninf(res[i]), ResNormNaNError,
00272 "StatusTestResNorm::checkStatus(): residual norm is nan or inf" );
00273 if (res[i] < tol_) {
00274 ind_[have] = i;
00275 have++;
00276 }
00277 }
00278 ind_.resize(have);
00279 int need = (quorum_ == -1) ? res.size() : quorum_;
00280 state_ = (have >= need) ? Passed : Failed;
00281 return state_;
00282 }
00283
00284
00285 template <class ScalarType, class MV, class OP>
00286 std::ostream& StatusTestResNorm<ScalarType,MV,OP>::print(std::ostream& os, int indent) const
00287 {
00288 std::string ind(indent,' ');
00289 os << ind << "- StatusTestResNorm: ";
00290 switch (state_) {
00291 case Passed:
00292 os << "Passed" << std::endl;
00293 break;
00294 case Failed:
00295 os << "Failed" << std::endl;
00296 break;
00297 case Undefined:
00298 os << "Undefined" << std::endl;
00299 break;
00300 }
00301 os << ind << " (Tolerance,WhichNorm,Scaled,Quorum): "
00302 << "(" << tol_;
00303 switch (whichNorm_) {
00304 case RES_ORTH:
00305 os << ",RES_ORTH";
00306 break;
00307 case RES_2NORM:
00308 os << ",RES_2NORM";
00309 break;
00310 case RITZRES_2NORM:
00311 os << ",RITZRES_2NORM";
00312 break;
00313 }
00314 os << "," << (scaled_ ? "true" : "false")
00315 << "," << quorum_
00316 << ")" << std::endl;
00317
00318 if (state_ != Undefined) {
00319 os << ind << " Which vectors: ";
00320 if (ind_.size() > 0) {
00321 for (unsigned int i=0; i<ind_.size(); i++) os << ind_[i] << " ";
00322 os << std::endl;
00323 }
00324 else {
00325 os << "[empty]" << std::endl;
00326 }
00327 }
00328 return os;
00329 }
00330
00331
00332 }
00333
00334 #endif