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 
00059 namespace Anasazi {
00060 
00061 
00062 template <class ScalarType, class MV, class OP>
00063 class StatusTestResNorm : public StatusTest<ScalarType,MV,OP> {
00064 
00065   typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00066 
00067  public:
00068 
00070 
00071 
00075   enum ResType {
00076     RES_ORTH,
00077     RES_2NORM,
00078     RITZRES_2NORM
00079   };
00080 
00082 
00084 
00085 
00087   StatusTestResNorm(MagnitudeType tol, int quorum = -1, ResType whichNorm = RES_ORTH, bool scaled = true);
00088 
00090   virtual ~StatusTestResNorm() {};
00092 
00094 
00095 
00098   TestStatus checkStatus( Eigensolver<ScalarType,MV,OP>* solver );
00099 
00101   TestStatus getStatus() const { return state_; }
00103 
00105 
00106 
00112   void setQuorum(int quorum) {
00113     state_ = Undefined;
00114     quorum_ = quorum;
00115   }
00116 
00119   int getQuorum() {
00120     return quorum_;
00121   }
00122 
00126   void setTolerance(MagnitudeType tol) {
00127     state_ = Undefined;
00128     tol_ = tol;
00129   }
00130 
00132   MagnitudeType getTolerance() {return tol_;}
00133 
00138   void setWhichNorm(ResType whichNorm) {
00139     state_ = Undefined;
00140     whichNorm_ = whichNorm;
00141   }
00142 
00144   ResType getWhichNorm() {return whichNorm_;}
00145 
00149   void setScale(bool relscale) {
00150     state_ = Undefined;
00151     scaled_ = relscale;
00152   }
00153 
00155   bool getScale() {return scaled_;}
00156 
00158   std::vector<int> whichVecs() {
00159     return ind_;
00160   }
00161 
00163   int howMany() {
00164     return ind_.size();
00165   }
00166 
00168 
00170 
00171 
00172 
00177   void reset() { 
00178     state_ = Undefined;
00179   }
00180 
00182 
00187   void clearStatus() {
00188     state_ = Undefined;
00189   }
00190 
00192 
00194 
00195   
00197   ostream& print(ostream& os, int indent = 0) const;
00198  
00200   private:
00201     TestStatus state_;
00202     MagnitudeType tol_;
00203     std::vector<int> ind_;
00204     int quorum_;
00205     bool scaled_;
00206     ResType whichNorm_;
00207 };
00208 
00209 
00210 template <class ScalarType, class MV, class OP>
00211 StatusTestResNorm<ScalarType,MV,OP>::StatusTestResNorm(MagnitudeType tol, int quorum, ResType whichNorm, bool scaled)
00212   : state_(Undefined), tol_(tol), quorum_(quorum), scaled_(scaled), whichNorm_(whichNorm) {}
00213 
00214 template <class ScalarType, class MV, class OP>
00215 TestStatus StatusTestResNorm<ScalarType,MV,OP>::checkStatus( Eigensolver<ScalarType,MV,OP>* solver ) {
00216   typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00217   
00218   std::vector<MagnitudeType> res;
00219 
00220   // get the eigenvector/ritz residuals norms (using the appropriate norm)
00221   // get the eigenvalues/ritzvalues and ritz index as well
00222   std::vector<Value<ScalarType> > vals = solver->getRitzValues();
00223   switch (whichNorm_) {
00224     case RES_2NORM:
00225       res = solver->getRes2Norms();
00226       // we want only the ritz values corresponding to our eigenvector residuals
00227       vals.resize(res.size());
00228       break;
00229     case RES_ORTH:
00230       res = solver->getResNorms();
00231       // we want only the ritz values corresponding to our eigenvector residuals
00232       vals.resize(res.size());
00233       break;
00234     case RITZRES_2NORM:
00235       res = solver->getRitzRes2Norms();
00236       break;
00237   }
00238 
00239   // if appropriate, scale the norms by the magnitude of the eigenvalue estimate
00240   if (scaled_) {
00241     Teuchos::LAPACK<int,MagnitudeType> lapack;
00242 
00243     for (unsigned int i=0; i<res.size(); i++) {
00244       MagnitudeType tmp = lapack.LAPY2(vals[i].realpart,vals[i].imagpart);
00245       // scale by the newly computed magnitude of the ritz values
00246       if ( tmp != MT::zero() ) {
00247         res[i] /= tmp;
00248       }
00249     }
00250   }
00251 
00252   // test the norms
00253   int have = 0;
00254   ind_.resize(res.size());
00255   for (unsigned int i=0; i<res.size(); i++) {
00256     TEST_FOR_EXCEPTION( MT::isnaninf(res[i]), StatusTestError, "StatusTestResNorm::checkStatus(): residual norm is nan or inf" );
00257     if (res[i] < tol_) {
00258       ind_[have] = i;
00259       have++;
00260     }
00261   }
00262   ind_.resize(have);
00263   int need = (quorum_ == -1) ? res.size() : quorum_;
00264   state_ = (have >= need) ? Passed : Failed;
00265   return state_;
00266 }
00267 
00268 
00269 template <class ScalarType, class MV, class OP>
00270 ostream& StatusTestResNorm<ScalarType,MV,OP>::print(ostream& os, int indent) const {
00271   string ind(indent,' ');
00272   os << ind << "- StatusTestResNorm: ";
00273   switch (state_) {
00274   case Passed:
00275     os << "Passed" << endl;
00276     break;
00277   case Failed:
00278     os << "Failed" << endl;
00279     break;
00280   case Undefined:
00281     os << "Undefined" << endl;
00282     break;
00283   }
00284   os << ind << "  (Tolerance,WhichNorm,Scaled,Quorum): " 
00285             << "(" << tol_;
00286   switch (whichNorm_) {
00287   case RES_ORTH:
00288     os << ",RES_ORTH";
00289     break;
00290   case RES_2NORM:
00291     os << ",RES_2NORM";
00292     break;
00293   case RITZRES_2NORM:
00294     os << ",RITZRES_2NORM";
00295     break;
00296   }
00297   os        << "," << (scaled_   ? "true" : "false")
00298             << "," << quorum_ 
00299             << ")" << endl;
00300 
00301   if (state_ != Undefined) {
00302     os << ind << "  Which vectors: ";
00303     if (ind_.size() > 0) {
00304       for (unsigned int i=0; i<ind_.size(); i++) os << ind_[i] << " ";
00305       os << endl;
00306     }
00307     else {
00308       os << "[empty]" << endl;
00309     }
00310   }
00311   return os;
00312 }
00313 
00314 
00315 } // end of Anasazi namespace
00316 
00317 #endif /* ANASAZI_STATUS_TEST_RESNORM_HPP */

Generated on Thu Sep 18 12:31:38 2008 for Anasazi by doxygen 1.3.9.1