Thyra_TsqrAdaptor.hpp

// @HEADER
// *****************************************************************************
//         Stratimikos: Thyra-based strategies for linear solvers
//
// Copyright 2006 NTESS and the Stratimikos contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
 
#ifndef __Thyra_TsqrAdaptor_hpp
#define __Thyra_TsqrAdaptor_hpp
 
#include "BelosConfigDefs.hpp"
 
// BelosThyraAdapter.hpp only includes this file if HAVE_BELOS_TSQR is
// defined.  Thus, it's OK to include TSQR header files here.
 
#include "Thyra_MultiVectorBase.hpp"
#include "Thyra_SpmdVectorSpaceBase.hpp"
 
#ifdef HAVE_MPI
#  include "Teuchos_DefaultMpiComm.hpp"
#endif // HAVE_MPI
#include "Teuchos_DefaultSerialComm.hpp"
#include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
 
#include <stdexcept>
 
 
namespace Thyra {
 
  template<class Scalar>
  class TsqrAdaptor : public Teuchos::ParameterListAcceptorDefaultBase {
  public:
    typedef Thyra::MultiVectorBase<Scalar> MV;
    typedef Scalar scalar_type;
    typedef int ordinal_type; // MultiVectorBase really does use int for this
    typedef Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> dense_matrix_type;
    typedef typename Teuchos::ScalarTraits<scalar_type>::magnitudeType magnitude_type;
 
    TsqrAdaptor (const Teuchos::RCP<Teuchos::ParameterList>& /* plist */)
    {
      throw std::logic_error ("Thyra adaptor for TSQR not implemented");
    }
 
    TsqrAdaptor ()
    {
      throw std::logic_error ("Thyra adaptor for TSQR not implemented");
    }
 
    Teuchos::RCP<const Teuchos::ParameterList>
    getValidParameters () const
    {
      throw std::logic_error ("Thyra adaptor for TSQR not implemented");
    }
 
    void
    setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& /* plist */)
    {
      throw std::logic_error ("Thyra adaptor for TSQR not implemented");
    }
 
    void
    factorExplicit (MV& /* A */,
                    MV& /* Q */,
                    dense_matrix_type& /* R */,
                    const bool /* forceNonnegativeDiagonal */ = false)
    {
      throw std::logic_error ("Thyra adaptor for TSQR not implemented");
    }
 
    int
    revealRank (MV& /* Q */,
                dense_matrix_type& /* R */,
                const magnitude_type& /* tol */)
    {
      throw std::logic_error ("Thyra adaptor for TSQR not implemented");
    }
 
  private:
    static Teuchos::RCP<const Teuchos::Comm<int> >
    getComm (const MV& X)
    {
      using Teuchos::RCP;
      using Teuchos::rcp;
      using Teuchos::rcp_dynamic_cast;
      using Teuchos::rcp_implicit_cast;
      typedef Thyra::VectorSpaceBase<Scalar> space_base_type;
      typedef Thyra::SpmdVectorSpaceBase<Scalar> space_type;
 
      // Thyra stores the communicator in the "vector space," but only
      // if that vector space is an SpmdVectorSpaceBase.
      RCP<const space_base_type> rangeBase = X.range ();
      TEUCHOS_TEST_FOR_EXCEPTION(rangeBase.is_null (), std::runtime_error, "X.range() is null.");
      RCP<const space_type> range = rcp_dynamic_cast<const space_type> (rangeBase);
      TEUCHOS_TEST_FOR_EXCEPTION(range.is_null (), std::runtime_error, "X.range() is not an SpmdVectorSpaceBase.");
 
      // Thyra annoyingly uses a (possibly) different template
      // parameter for its Teuchos::Comm than everybody else.  The
      // least hackish way to work around this is to convert the Comm
      // to one of two subclasses (MpiComm or SerialComm).  If it's an
      // MpiComm, we can extract the RCP<const OpaqueWrapper<MPI_Comm>
      // > and make a new MpiComm<int> from it.  If it's a SerialComm,
      // just create a new SerialComm<int>.  If it's neither of those,
      // then I have no idea what to do.  Note that MpiComm is only
      // defined if HAVE_MPI is defined.
      RCP<const Teuchos::Comm<Thyra::Ordinal> > thyraComm = range->getComm ();
#ifdef HAVE_MPI
      RCP<const Teuchos::MpiComm<Thyra::Ordinal> > thyraMpiComm =
        rcp_dynamic_cast<const Teuchos::MpiComm<Thyra::Ordinal> > (thyraComm);
      if (thyraMpiComm.is_null ()) {
        RCP<const Teuchos::SerialComm<Thyra::Ordinal> > thyraSerialComm =
          rcp_dynamic_cast<const Teuchos::SerialComm<Thyra::Ordinal> > (thyraComm);
        TEUCHOS_TEST_FOR_EXCEPTION(
          thyraSerialComm.is_null (), std::runtime_error,
          "Thyra's communicator is neither an MpiComm nor a SerialComm.  "
          "Sorry, I have no idea what to do with it in that case.");
        // It's a SerialComm.  Make a SerialComm of our own.
        // SerialComm instances are all the same, so there's no need
        // to keep the original one.
        return rcp_implicit_cast<const Teuchos::Comm<int> > (rcp (new Teuchos::SerialComm<int>));
      }
      else { // Yippie, we have an MpiComm.
        RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > rawMpiComm = thyraMpiComm->getRawMpiComm ();
        // NOTE (mfh 18 Jun 2013) Since the error handler is attached
        // to the MPI_Comm, not to the Teuchos widget, we don't have
        // to set the error handler again on the new MpiComm object.
        return rcp_implicit_cast<const Teuchos::Comm<int> > (rcp (new Teuchos::MpiComm<int> (rawMpiComm)));
      }
#else // NOT HAVE_MPI
      // Either it's a SerialComm or I don't know what to do with it.
      RCP<const Teuchos::SerialComm<Thyra::Ordinal> > thyraSerialComm =
        rcp_dynamic_cast<const Teuchos::SerialComm<Thyra::Ordinal> > (thyraComm);
      TEUCHOS_TEST_FOR_EXCEPTION(
        thyraSerialComm.is_null (), std::runtime_error,
        "Thyra's communicator is not a SerialComm, and MPI is not enabled, so "
        "it can't be an MpiComm either.  That means it must be some other "
        "subclass of Comm, about which I don't know.  "
        "Sorry, I have no idea what to do with it in that case.");
      // It's a SerialComm.  Make a SerialComm of our own.
      // SerialComm instances are all the same, so there's no need
      // to keep the original one.
      return rcp_implicit_cast<const Teuchos::Comm<int> > (rcp (new Teuchos::SerialComm<int>));
#endif // HAVE_MPI
    }
 
    void
    prepareDistTsqr (const MV& /* X */) {}
 
    void
    prepareTsqr (const MV& /* X */) {}
  };
 
} // namespace Tpetra
 
#endif // __Thyra_TsqrAdaptor_hpp