AnasaziStatusTestResNorm.hpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Anasazi: Block Eigensolvers Package
00005 //                 Copyright (2004) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
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     // get the eigenvector/ritz residuals norms (using the appropriate norm)
00236     // get the eigenvalues/ritzvalues and ritz index as well
00237     std::vector<Value<ScalarType> > vals = solver->getRitzValues();
00238     switch (whichNorm_) {
00239       case RES_2NORM:
00240         res = solver->getRes2Norms();
00241         // we want only the ritz values corresponding to our eigenvector residuals
00242         vals.resize(res.size());
00243         break;
00244       case RES_ORTH:
00245         res = solver->getResNorms();
00246         // we want only the ritz values corresponding to our eigenvector residuals
00247         vals.resize(res.size());
00248         break;
00249       case RITZRES_2NORM:
00250         res = solver->getRitzRes2Norms();
00251         break;
00252     }
00253 
00254     // if appropriate, scale the norms by the magnitude of the eigenvalue estimate
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         // scale by the newly computed magnitude of the ritz values
00261         if ( tmp != MT::zero() ) {
00262           res[i] /= tmp;
00263         }
00264       }
00265     }
00266 
00267     // test the norms
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 } // end of Anasazi namespace
00333 
00334 #endif /* ANASAZI_STATUS_TEST_RESNORM_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends
Generated on Wed Apr 13 09:56:59 2011 for Anasazi by  doxygen 1.6.3