Tpetra_TsqrAdaptor_UQ_PCE.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//          Tpetra: Templated Linear Algebra Services Package
//
// Copyright 2008 NTESS and the Tpetra contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER

#ifndef TPETRA_TSQR_ADAPTOR_UQ_PCE_HPP
#define TPETRA_TSQR_ADAPTOR_UQ_PCE_HPP

#include <Tpetra_ConfigDefs.hpp> // HAVE_TPETRA_TSQR, etc.

#ifdef HAVE_TPETRA_TSQR

#include "Stokhos_Sacado_Kokkos_UQ_PCE.hpp"

#  include "Tsqr_NodeTsqrFactory.hpp" // create intranode TSQR object
#  include "Tsqr.hpp" // full (internode + intranode) TSQR
#  include "Tsqr_DistTsqr.hpp" // internode TSQR
// Subclass of TSQR::MessengerBase, implemented using Teuchos
// communicator template helper functions
#  include "Tsqr_TeuchosMessenger.hpp"
#  include "Tpetra_MultiVector.hpp"
#  include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
#  include <stdexcept>

// Base TsqrAdator template we will specialize
#  include "Tpetra_TsqrAdaptor.hpp"

namespace Tpetra {

  template <class Storage, class LO, class GO, class Node>
  class TsqrAdaptor< Tpetra::MultiVector< Sacado::UQ::PCE<Storage>,
                                          LO, GO, Node > > :
    public Teuchos::ParameterListAcceptorDefaultBase {
  public:
    typedef Tpetra::MultiVector< Sacado::UQ::PCE<Storage>, LO, GO, Node > MV;
    typedef typename MV::scalar_type mp_scalar_type;

    // For Sacado::UQ::PCE< Storage<Ordinal,Scalar,Device> > this is Scalar
    typedef typename mp_scalar_type::scalar_type scalar_type;
    typedef typename mp_scalar_type::ordinal_type mp_ordinal_type;
    typedef typename MV::local_ordinal_type ordinal_type;
    typedef Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> dense_matrix_type;
    typedef typename Teuchos::ScalarTraits<scalar_type>::magnitudeType magnitude_type;

  private:
    using node_tsqr_factory_type =
      TSQR::NodeTsqrFactory<scalar_type, ordinal_type,
                            typename MV::device_type>;
    using node_tsqr_type = TSQR::NodeTsqr<ordinal_type, scalar_type>;
    using dist_tsqr_type = TSQR::DistTsqr<ordinal_type, scalar_type>;
    using tsqr_type = TSQR::Tsqr<ordinal_type, scalar_type>;

  public:
    TsqrAdaptor (const Teuchos::RCP<Teuchos::ParameterList>& plist) :
      nodeTsqr_ (node_tsqr_factory_type::getNodeTsqr ()),
      distTsqr_ (new dist_tsqr_type),
      tsqr_ (new tsqr_type (nodeTsqr_, distTsqr_)),
      ready_ (false)
    {
      setParameterList (plist);
    }

    TsqrAdaptor () :
      nodeTsqr_ (new node_tsqr_factory_type::getNodeTsqr ()),
      distTsqr_ (new dist_tsqr_type),
      tsqr_ (new tsqr_type (nodeTsqr_, distTsqr_)),
      ready_ (false)
    {
      setParameterList (Teuchos::null);
    }

    Teuchos::RCP<const Teuchos::ParameterList>
    getValidParameters () const
    {
      using Teuchos::RCP;
      using Teuchos::rcp;
      using Teuchos::ParameterList;
      using Teuchos::parameterList;

      if (defaultParams_.is_null()) {
        RCP<ParameterList> params = parameterList ("TSQR implementation");
        params->set ("NodeTsqr", *(nodeTsqr_->getValidParameters ()));
        params->set ("DistTsqr", *(distTsqr_->getValidParameters ()));
        defaultParams_ = params;
      }
      return defaultParams_;
    }

    void
    setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& plist)
    {
      using Teuchos::ParameterList;
      using Teuchos::parameterList;
      using Teuchos::RCP;
      using Teuchos::sublist;

      RCP<ParameterList> params = plist.is_null() ?
        parameterList (*getValidParameters ()) : plist;
      nodeTsqr_->setParameterList (sublist (params, "NodeTsqr"));
      distTsqr_->setParameterList (sublist (params, "DistTsqr"));

      this->setMyParamList (params);
    }

    void
    factorExplicit (MV& A,
                    MV& Q,
                    dense_matrix_type& R,
                    const bool forceNonnegativeDiagonal=false)
    {
      prepareTsqr (Q); // Finish initializing TSQR.

      ordinal_type numRows;
      ordinal_type numCols;
      ordinal_type LDA;
      ordinal_type LDQ;
      scalar_type* A_ptr;
      scalar_type* Q_ptr;

      getNonConstView (numRows, numCols, A_ptr, LDA, A);
      getNonConstView (numRows, numCols, Q_ptr, LDQ, Q);
      const bool contiguousCacheBlocks = false;
      tsqr_->factorExplicitRaw (numRows, numCols, A_ptr, LDA,
                                Q_ptr, LDQ, R.values (), R.stride (),
                                contiguousCacheBlocks,
                                forceNonnegativeDiagonal);
    }

    int
    revealRank (MV& Q,
                dense_matrix_type& R,
                const magnitude_type& tol)
    {
      prepareTsqr (Q); // Finish initializing TSQR.

      // FIXME (mfh 18 Oct 2010) Check Teuchos::Comm<int> object in Q
      // to make sure it is the same communicator as the one we are
      // using in our dist_tsqr_type implementation.

      ordinal_type numRows;
      ordinal_type numCols;
      scalar_type* Q_ptr;
      ordinal_type LDQ;
      getNonConstView (numRows, numCols, Q_ptr, LDQ, Q);
      const bool contiguousCacheBlocks = false;
      return tsqr_->revealRankRaw (numRows, numCols, Q_ptr, LDQ,
                                   R.values (), R.stride (), tol,
                                   contiguousCacheBlocks);
    }

  private:
    Teuchos::RCP<node_tsqr_type> nodeTsqr_;

    Teuchos::RCP<dist_tsqr_type> distTsqr_;

    Teuchos::RCP<tsqr_type> tsqr_;

    mutable Teuchos::RCP<const Teuchos::ParameterList> defaultParams_;

    bool ready_;

    void
    prepareTsqr (const MV& mv)
    {
      if (! ready_) {
        prepareDistTsqr (mv);
        ready_ = true;
      }
    }

    void
    prepareDistTsqr (const MV& mv)
    {
      using Teuchos::RCP;
      using Teuchos::rcp_implicit_cast;
      typedef TSQR::TeuchosMessenger<scalar_type> mess_type;
      typedef TSQR::MessengerBase<scalar_type> base_mess_type;

      RCP<const Teuchos::Comm<int> > comm = mv.getMap()->getComm();
      RCP<mess_type> mess (new mess_type (comm));
      RCP<base_mess_type> messBase = rcp_implicit_cast<base_mess_type> (mess);
      distTsqr_->init (messBase);
    }

    static void
    getNonConstView (ordinal_type& numRows,
                     ordinal_type& numCols,
                     scalar_type*& A_ptr,
                     ordinal_type& LDA,
                     const MV& A)
    {
      // FIXME (mfh 25 Oct 2010) We should be able to run TSQR even if
      // storage of A uses nonconstant stride internally.  We would
      // have to copy and pack into a matrix with constant stride, and
      // then unpack on exit.  For now we choose just to raise an
      // exception.
      TEUCHOS_TEST_FOR_EXCEPTION
        (! A.isConstantStride(), std::invalid_argument,
         "TSQR does not currently support Tpetra::MultiVector "
         "inputs that do not have constant stride.");

      // FIXME (mfh 16 Jan 2016) When I got here, I found strides[0]
      // instead of strides[1] for the stride.  I don't think this is
      // right.  However, I don't know about these Stokhos scalar
      // types so I'll just do what was here.
      //
      // STOKHOS' TYPES ARE NOT TESTED WITH TSQR REGULARLY SO IT IS
      // POSSIBLE THAT THE ORIGINAL CODE WAS WRONG.

      typedef typename MV::dual_view_type view_type;
      typedef typename view_type::t_dev::array_type flat_array_type;

      // Reinterpret the data as a longer array of the base scalar
      // type.  TSQR currently forbids MultiVector input with
      // nonconstant stride, so we need not worry about that here.

      flat_array_type flat_mv = A.getLocalViewDevice();

      numRows = static_cast<ordinal_type> (flat_mv.extent(0));
      numCols = static_cast<ordinal_type> (flat_mv.extent(1));
      A_ptr = flat_mv.data ();

      ordinal_type strides[2];
      flat_mv.stride (strides);
      LDA = strides[0];
    }
  };
} // namespace Tpetra

#endif // HAVE_TPETRA_TSQR

#endif // TPETRA_TSQR_ADAPTOR_UQ_PCE_HPP