BelosTsqrOrthoManager.hpp

Go to the documentation of this file.

//@HEADER
// ************************************************************************
//
//                 Belos: Block Linear Solvers Package
//                  Copyright 2004 Sandia Corporation
//
// Under the 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 __BelosTsqrOrthoManager_hpp
#define __BelosTsqrOrthoManager_hpp

#include "BelosTsqrOrthoManagerImpl.hpp"
// Belos doesn't have an SVQB orthomanager implemented yet, so we just
// use a default orthomanager for the case where the inner product
// matrix is nontrivial.
#include "BelosDGKSOrthoManager.hpp"


namespace Belos {

template<class Scalar, class MV>
class OutOfPlaceNormalizerMixin {
public:
  typedef Scalar scalar_type;
  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
  typedef MV multivector_type;

  typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
  typedef Teuchos::RCP<mat_type>                  mat_ptr;

  virtual int
  normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B) const = 0;

  virtual int
  projectAndNormalizeOutOfPlace (MV& X_in,
                                 MV& X_out,
                                 Teuchos::Array<mat_ptr> C,
                                 mat_ptr B,
                                 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const = 0;

  virtual ~OutOfPlaceNormalizerMixin () {}
};

template<class Scalar, class MV>
class TsqrOrthoManager :
    public OrthoManager<Scalar, MV>,
    public OutOfPlaceNormalizerMixin<Scalar, MV>
{
public:
  typedef Scalar scalar_type;
  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
  typedef MV multivector_type;

  typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
  typedef Teuchos::RCP<mat_type>                  mat_ptr;

  void setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& params) {
    impl_.setParameterList (params);
  }

  Teuchos::RCP<Teuchos::ParameterList> getNonconstParameterList () {
    return impl_.getNonconstParameterList ();
  }

  Teuchos::RCP<Teuchos::ParameterList> unsetParameterList () {
    return impl_.unsetParameterList ();
  }

  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const {
    return impl_.getValidParameters();
  }

  Teuchos::RCP<const Teuchos::ParameterList> getFastParameters() const {
    return impl_.getFastParameters();
  }

  TsqrOrthoManager (const Teuchos::RCP<Teuchos::ParameterList>& params,
                    const std::string& label = "Belos") :
    impl_ (params, label)
  {}

  TsqrOrthoManager (const std::string& label) :
    impl_ (label)
  {}

  virtual ~TsqrOrthoManager() {}

  void
  setReorthogonalizationCallback (const Teuchos::RCP<ReorthogonalizationCallback<Scalar> >& callback)
  {
    impl_.setReorthogonalizationCallback (callback);
  }

  void innerProd (const MV &X, const MV &Y, mat_type& Z) const {
    return impl_.innerProd (X, Y, Z);
  }

  void norm (const MV& X, std::vector<magnitude_type>& normVec) const {
    return impl_.norm (X, normVec);
  }

  void
  project (MV &X,
           Teuchos::Array<mat_ptr> C,
           Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
  {
    return impl_.project (X, C, Q);
  }

  int
  normalize (MV &X, mat_ptr B) const
  {
    return impl_.normalize (X, B);
  }

protected:
  virtual int
  projectAndNormalizeImpl (MV &X,
                           Teuchos::Array<mat_ptr> C,
                           mat_ptr B,
                           Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
  {
    return impl_.projectAndNormalize (X, C, B, Q);
  }

public:
  int
  normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B) const
  {
    return impl_.normalizeOutOfPlace (X, Q, B);
  }

  int
  projectAndNormalizeOutOfPlace (MV& X_in,
                                 MV& X_out,
                                 Teuchos::Array<mat_ptr> C,
                                 mat_ptr B,
                                 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
  {
    return impl_.projectAndNormalizeOutOfPlace (X_in, X_out, C, B, Q);
  }

  magnitude_type orthonormError (const MV &X) const {
    return impl_.orthonormError (X);
  }

  magnitude_type orthogError (const MV &X1, const MV &X2) const {
    return impl_.orthogError (X1, X2);
  }

  void setLabel (const std::string& label) {
    impl_.setLabel (label);
  }

  const std::string& getLabel() const { return impl_.getLabel(); }

private:
  mutable TsqrOrthoManagerImpl<Scalar, MV> impl_;

  std::string label_;
};


template<class Scalar, class MV, class OP>
class TsqrMatOrthoManager :
    public MatOrthoManager<Scalar, MV, OP>,
    public OutOfPlaceNormalizerMixin<Scalar, MV>
{
public:
  typedef Scalar scalar_type;
  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
  typedef MV multivector_type;
  typedef OP operator_type;

  typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
  typedef Teuchos::RCP<mat_type>                  mat_ptr;

private:
  typedef MatOrthoManager<Scalar, MV, OP> base_type;

  typedef TsqrOrthoManagerImpl<Scalar, MV> tsqr_type;

  typedef DGKSOrthoManager<Scalar, MV, OP> dgks_type;

  typedef MultiVecTraits<Scalar, MV> MVT;

public:
  TsqrMatOrthoManager (const Teuchos::RCP<Teuchos::ParameterList>& params,
                       const std::string& label = "Belos",
                       Teuchos::RCP<const OP> Op = Teuchos::null) :
    MatOrthoManager<Scalar, MV, OP> (Op),
    tsqr_ (params, label),
    pDgks_ (Teuchos::null) // Initialized lazily
  {}

  TsqrMatOrthoManager (const std::string& label = "Belos",
                       Teuchos::RCP<const OP> Op = Teuchos::null) :
    MatOrthoManager<Scalar, MV, OP> (Op),
    tsqr_ (label),
    pDgks_ (Teuchos::null) // Initialized lazily
  {}

  virtual ~TsqrMatOrthoManager() {}

  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const {
    return tsqr_.getValidParameters ();
  }

  Teuchos::RCP<const Teuchos::ParameterList> getFastParameters() {
    return tsqr_.getFastParameters ();
  }

  void setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& params) {
    tsqr_.setParameterList (params);
  }

  const std::string& getLabel() const { return tsqr_.getLabel (); }

  void
  setOp (Teuchos::RCP<const OP> Op)
  {
    // We override the base class' setOp() so that the
    // DGKSOrthoManager instance gets the new op.
    //
    // The base_type notation helps C++ resolve the name for a
    // member function of a templated base class.
    base_type::setOp (Op); // base class gets a copy of the Op too

    if (! Op.is_null()) {
      ensureDgksInit (); // Make sure the DGKS object has been initialized
      pDgks_->setOp (Op);
    }
  }

  Teuchos::RCP<const OP> getOp () const {
    // The base_type notation helps C++ resolve the name for a
    // member function of a templated base class.
    return base_type::getOp();
  }

  void
  project (MV &X,
           Teuchos::RCP<MV> MX,
           Teuchos::Array<mat_ptr> C,
           Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
  {
    if (getOp().is_null()) {
      tsqr_.project (X, C, Q);
      if (! MX.is_null()) {
        // MX gets a copy of X; M is the identity operator.
        MVT::Assign (X, *MX);
      }
    } else {
      ensureDgksInit ();
      pDgks_->project (X, MX, C, Q);
    }
  }

  void
  project (MV &X,
           Teuchos::Array<mat_ptr> C,
           Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
  {
    project (X, Teuchos::null, C, Q);
  }

  int
  normalize (MV& X, Teuchos::RCP<MV> MX, mat_ptr B) const
  {
    if (getOp().is_null()) {
      const int rank = tsqr_.normalize (X, B);
      if (! MX.is_null()) {
        // MX gets a copy of X; M is the identity operator.
        MVT::Assign (X, *MX);
      }
      return rank;
    } else {
      ensureDgksInit ();
      return pDgks_->normalize (X, MX, B);
    }
  }

  int normalize (MV& X, mat_ptr B) const {
    return normalize (X, Teuchos::null, B);
  }

  // Attempted fix for a warning about hiding
  // OrthoManager::projectAndNormalize(). Unfortunately, this fix turns out
  // to produce a compilation error with cray++, see bug #6129,
  // <https://software.sandia.gov/bugzilla/show_bug.cgi?id=6129>.
  //using Belos::OrthoManager<Scalar, MV>::projectAndNormalize;

protected:
  virtual int
  projectAndNormalizeWithMxImpl (MV &X,
                                 Teuchos::RCP<MV> MX,
                                 Teuchos::Array<mat_ptr> C,
                                 mat_ptr B,
                                 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
  {
    if (getOp().is_null()) {
      const int rank = tsqr_.projectAndNormalize (X, C, B, Q);
      if (! MX.is_null()) {
        // MX gets a copy of X; M is the identity operator.
        MVT::Assign (X, *MX);
      }
      return rank;
    } else {
      ensureDgksInit ();
      return pDgks_->projectAndNormalize (X, MX, C, B, Q);
    }
  }

public:
  int
  normalizeOutOfPlace (MV& X, MV& Q, mat_ptr B) const
  {
    if (getOp().is_null()) {
      return tsqr_.normalizeOutOfPlace (X, Q, B);
    } else {
      // DGKS normalizes in place, so we have to copy.
      ensureDgksInit ();
      const int rank = pDgks_->normalize (X, B);
      MVT::Assign (X, Q);
      return rank;
    }
  }

  int
  projectAndNormalizeOutOfPlace (MV& X_in,
                                 MV& X_out,
                                 Teuchos::Array<mat_ptr> C,
                                 mat_ptr B,
                                 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
  {
    using Teuchos::null;

    if (getOp().is_null()) {
      return tsqr_.projectAndNormalizeOutOfPlace (X_in, X_out, C, B, Q);
    } else {
      // DGKS normalizes in place, so we have to copy.
      ensureDgksInit ();
      const int rank = pDgks_->projectAndNormalize (X_in, null, C, B, Q);
      MVT::Assign (X_in, X_out);
      return rank;
    }
  }

  magnitude_type
  orthonormError (const MV &X, Teuchos::RCP<const MV> MX) const
  {
    if (getOp().is_null()) {
      return tsqr_.orthonormError (X); // Ignore MX
    } else {
      ensureDgksInit ();
      return pDgks_->orthonormError (X, MX);
    }
  }

  magnitude_type orthonormError (const MV &X) const {
    return orthonormError (X, Teuchos::null);
  }

  magnitude_type orthogError (const MV &X1, const MV &X2) const {
    return orthogError (X1, Teuchos::null, X2);
  }

  magnitude_type
  orthogError (const MV &X1,
               Teuchos::RCP<const MV> MX1,
               const MV &X2) const
  {
    if (getOp ().is_null ()) {
      // Ignore MX1, since we don't need to write to it.
      return tsqr_.orthogError (X1, X2);
    } else {
      ensureDgksInit ();
      return pDgks_->orthogError (X1, MX1, X2);
    }
  }

  void
  setLabel (const std::string& label)
  {
    tsqr_.setLabel (label);

    // Make sure DGKS gets the new label, if it's initialized.
    // Otherwise, it will get the new label on initialization.
    if (! pDgks_.is_null ()) {
      pDgks_->setLabel (label);
    }
  }

private:
  void
  ensureDgksInit () const
  {
    // NOTE (mfh 11 Jan 2011) DGKS has a parameter that needs to be
    // set.  For now, we just use the default value of the parameter.
    if (pDgks_.is_null ()) {
      pDgks_ = Teuchos::rcp (new dgks_type (getLabel (), getOp ()));
    }
  }

  mutable tsqr_type tsqr_;

  mutable Teuchos::RCP<dgks_type> pDgks_;
};

} // namespace Belos

#endif // __BelosTsqrOrthoManager_hpp