AnasaziStatusTestWithOrdering.hpp

Go to the documentation of this file.

// @HEADER
// ***********************************************************************
//
//                 Anasazi: Block Eigensolvers Package
//                 Copyright 2004 Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ***********************************************************************
// @HEADER
//

#ifndef ANASAZI_STATUS_TEST_ORDEREDRESNORM_HPP
#define ANASAZI_STATUS_TEST_ORDEREDRESNORM_HPP

#include "AnasaziStatusTest.hpp"
#include "Teuchos_ScalarTraits.hpp"
#include "Teuchos_LAPACK.hpp"

namespace Anasazi {


template <class ScalarType, class MV, class OP>
class StatusTestWithOrdering : public StatusTest<ScalarType,MV,OP> {

 private:
  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
  typedef Teuchos::ScalarTraits<MagnitudeType> MT;

 public:



  StatusTestWithOrdering(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test, Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sorter, int quorum = -1);

  virtual ~StatusTestWithOrdering() {};



  TestStatus checkStatus( Eigensolver<ScalarType,MV,OP>* solver );

  TestStatus getStatus() const { return state_; }


  std::vector<int> whichVecs() const {
    return ind_;
  }

  int howMany() const {
    return ind_.size();
  }




  void setQuorum(int quorum) {
    state_ = Undefined;
    quorum_ = quorum;
  }

  int getQuorum() const {
    return quorum_;
  }




  void reset() { 
    ind_.resize(0);
    state_ = Undefined;
    test_->reset();
  }


  void clearStatus() {
    ind_.resize(0);
    state_ = Undefined;
    test_->clearStatus();
  }

  void setAuxVals(const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &vals) {
    rvals_ = vals;
    ivals_.resize(rvals_.size(),MT::zero());
    state_ = Undefined;
  }

  void setAuxVals(const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &rvals, const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &ivals) {
    rvals_ = rvals;
    ivals_ = ivals;
    state_ = Undefined;
  }

  void getAuxVals(std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &rvals, std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &ivals) const {
    rvals = rvals_;
    ivals = ivals_;
  }



  
  std::ostream& print(std::ostream& os, int indent = 0) const;
 
  private:
    TestStatus state_;
    std::vector<int> ind_;
    int quorum_;
    std::vector<MagnitudeType> rvals_, ivals_;
    Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sorter_;
    Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test_;
};


template <class ScalarType, class MV, class OP>
StatusTestWithOrdering<ScalarType,MV,OP>::StatusTestWithOrdering(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test, Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sorter, int quorum)
  : state_(Undefined), ind_(0), quorum_(quorum), rvals_(0), ivals_(0), sorter_(sorter), test_(test)
{
  TEUCHOS_TEST_FOR_EXCEPTION(sorter_ == Teuchos::null, StatusTestError, "StatusTestWithOrdering::constructor() was passed null pointer for constituent SortManager.");
  TEUCHOS_TEST_FOR_EXCEPTION(test_ == Teuchos::null, StatusTestError, "StatusTestWithOrdering::constructor() was passed null pointer for constituent StatusTest.");
}

template <class ScalarType, class MV, class OP>
TestStatus StatusTestWithOrdering<ScalarType,MV,OP>::checkStatus( Eigensolver<ScalarType,MV,OP>* solver ) {

  // Algorithm
  // we PASS iff the "most signficant" values of "all values" PASS
  // "most significant" is measured by sorter
  // auxilliary values are automatically PASSED
  // constituent status test determines who else passes
  // "all values" mean {auxilliary values} UNION {solver's ritz values}
  //
  // test_->checkStatus()                   calls the constituent status test
  // cwhch = test_->whichVecs()             gets the passing indices from the constituent test
  // solval = solver->getRitzValues()       gets the solver's ritz values
  // allval = {solval auxval}               contains all values
  // sorter_->sort(allval,perm)             sort all values (we just need the perm vector)
  //
  // allpass = {cwhch numsolval+1 ... numsolval+numAux}
  // mostsig = {perm[1] ... perm[quorum]}
  // whichVecs = mostsig INTERSECT allpass
  // if size(whichVecs) >= quorum,
  //    PASS
  // else
  //    FAIL
  // 
  // finish: this needs to be better tested and revisited for efficiency.

  // call the constituent solver manager
  test_->checkStatus(solver);
  std::vector<int> cwhch( test_->whichVecs() );

  // get the ritzvalues from solver
  std::vector<Value<ScalarType> > solval = solver->getRitzValues();
  int numsolval = solval.size();
  int numauxval = rvals_.size();
  int numallval = numsolval + numauxval;

  if (numallval == 0) {
    ind_.resize(0);
    return Failed;
  }

  // make space for all values, real and imaginary parts
  std::vector<MagnitudeType> allvalr(numallval), allvali(numallval);
  // separate the real and imaginary parts of solver ritz values
  for (int i=0; i<numsolval; ++i) {
    allvalr[i] = solval[i].realpart;
    allvali[i] = solval[i].imagpart;
  }
  // put the auxiliary values at the end of the solver ritz values
  std::copy(rvals_.begin(),rvals_.end(),allvalr.begin()+numsolval);
  std::copy(ivals_.begin(),ivals_.end(),allvali.begin()+numsolval);

  // sort all values
  std::vector<int> perm(numallval);
  sorter_->sort(allvalr,allvali,Teuchos::rcpFromRef(perm),numallval);

  // make the set of passing values: allpass = {cwhch -1 ... -numauxval}
  std::vector<int> allpass(cwhch.size() + numauxval);
  std::copy(cwhch.begin(),cwhch.end(),allpass.begin());
  for (int i=0; i<numauxval; i++) {
    allpass[cwhch.size()+i] = -(i+1);
  }

  // make list, with length quorum, of most significant values, if there are that many
  int numsig = quorum_ < numallval ? quorum_ : numallval;
  // int numsig = cwhch.size() + numauxval;
  std::vector<int> mostsig(numsig);
  for (int i=0; i<numsig; ++i) {
    mostsig[i] = perm[i];
    // if perm[i] >= numsolval, then it corresponds to the perm[i]-numsolval aux val
    // the first aux val gets index -numauxval, the second -numauxval+1, and so forth
    if (mostsig[i] >= numsolval) {
      mostsig[i] = mostsig[i]-numsolval-numauxval;
    }
  }

  // who passed?
  // to use set_intersection, ind_ must have room for the result, which will have size() <= min( allpass.size(), mostsig.size() )
  // also, allpass and mostsig must be in ascending order; neither are
  // lastly, the return from set_intersection points to the last element in the intersection, which tells us how big the intersection is
  ind_.resize(numsig);
  std::vector<int>::iterator end;
  std::sort(mostsig.begin(),mostsig.end());
  std::sort(allpass.begin(),allpass.end());
  end = std::set_intersection(mostsig.begin(),mostsig.end(),allpass.begin(),allpass.end(),ind_.begin());
  ind_.resize(end - ind_.begin());

  // did we pass, overall
  if (ind_.size() >= (unsigned int)quorum_) {
    state_ = Passed;
  }
  else {
    state_ = Failed;
  }
  return state_;
}


template <class ScalarType, class MV, class OP>
std::ostream& StatusTestWithOrdering<ScalarType,MV,OP>::print(std::ostream& os, int indent) const {
  // build indent string
  std::string ind(indent,' ');
  // print header
  os << ind << "- StatusTestWithOrdering: ";
  switch (state_) {
  case Passed:
    os << "Passed" << std::endl;
    break;
  case Failed:
    os << "Failed" << std::endl;
    break;
  case Undefined:
    os << "Undefined" << std::endl;
    break;
  }
  // print parameters, namely, quorum
  os << ind << "  Quorum: " << quorum_ << std::endl;
  // print any auxilliary values
  os << ind << "  Auxiliary values: ";
  if (rvals_.size() > 0) {
    for (unsigned int i=0; i<rvals_.size(); i++) {
      os << "(" << rvals_[i] << ", " << ivals_[i] << ")  ";
    }
    os << std::endl;
  }
  else {
    os << "[empty]" << std::endl;
  }
  // print which vectors have passed
  if (state_ != Undefined) {
    os << ind << "  Which vectors: ";
    if (ind_.size() > 0) {
      for (unsigned int i=0; i<ind_.size(); i++) os << ind_[i] << " ";
      os << std::endl;
    }
    else {
      os << "[empty]" << std::endl;
    }
  }
  // print constituent test
  test_->print(os,indent+2);
  return os;
}


} // end of Anasazi namespace

#endif /* ANASAZI_STATUS_TEST_ORDEREDRESNORM_HPP */