Epetra_TsqrAdaptor.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 EPETRA_TSQRADAPTOR_HPP
#define EPETRA_TSQRADAPTOR_HPP
 
 
#include "Tpetra_ConfigDefs.hpp"
 
#if defined(TPETRA_ENABLE_DEPRECATED_CODE)
#if defined(TPETRA_DEPRECATED_DECLARATIONS)
#warning This file is deprecated due to Epetra removal and will be removed
#endif
#else
#error This file is deprecated due to Epetra removal and will be removed
#endif
 
#if defined(TPETRA_ENABLE_DEPRECATED_CODE) && defined(HAVE_TPETRA_EPETRA) && defined(HAVE_TPETRA_TSQR)
 
#include "Tsqr_NodeTsqrFactory.hpp"  // create intranode TSQR object
#include "Tsqr.hpp"                  // full (internode + intranode) TSQR
#include "Tsqr_DistTsqr.hpp"         // internode TSQR
#include "Epetra_Comm.h"
// Subclass of TSQR::MessengerBase, implemented using Teuchos
// communicator template helper functions
#include "Epetra_TsqrMessenger.hpp"
#include "Epetra_MultiVector.h"
#include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
#include <stdexcept>
 
namespace Epetra {
 
TPETRA_DEPRECATED_MSG("epetra removal")
class TsqrAdaptor : public Teuchos::ParameterListAcceptorDefaultBase {
 public:
  typedef Epetra_MultiVector MV;
 
  typedef double scalar_type;
 
  typedef int ordinal_type;
 
  using device_type =
      Kokkos::Device<Kokkos::DefaultHostExecutionSpace,
                     Kokkos::HostSpace>;
 
  using dense_matrix_type =
      Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>;
 
  using magnitude_type = double;
 
 private:
  using matview_type = TSQR::MatView<ordinal_type, scalar_type>;
  using node_tsqr_factory_type =
      TSQR::NodeTsqrFactory<scalar_type, ordinal_type, device_type>;
  // Don't need a "typename" here, because there are no template
  // parameters involved in the type definition.
  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_(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::ParameterList;
    using Teuchos::parameterList;
    using Teuchos::RCP;
    using Teuchos::rcp;
 
    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.
 
    scalar_type* const A_ptr   = A.Values();
    scalar_type* const Q_ptr   = Q.Values();
    scalar_type* const R_ptr   = R.values();
    const ordinal_type numRows = A.MyLength();
    const ordinal_type numCols = A.NumVectors();
    const ordinal_type lda     = A.Stride();
    const ordinal_type ldq     = Q.Stride();
    const ordinal_type ldr     = R.stride();
 
    const bool contiguousCacheBlocks = false;
    tsqr_->factorExplicitRaw(numRows, numCols, A_ptr, lda,
                             Q_ptr, ldq, R_ptr, ldr,
                             contiguousCacheBlocks,
                             forceNonnegativeDiagonal);
  }
 
  int revealRank(MV& Q,
                 dense_matrix_type& R,
                 const magnitude_type& tol) {
    TEUCHOS_TEST_FOR_EXCEPTION(!Q.ConstantStride(), std::invalid_argument,
                               "TsqrAdaptor::"
                               "revealRank: Input MultiVector Q must have constant stride.");
    prepareTsqr(Q);  // Finish initializing TSQR.
    // FIXME (mfh 25 Oct 2010) Check Epetra_Comm object in Q to make
    // sure it is the same communicator as the one we are using in
    // our dist_tsqr_type implementation.
    return tsqr_->revealRankRaw(Q.MyLength(), Q.NumVectors(),
                                Q.Values(), Q.Stride(),
                                R.values(), R.stride(), tol, false);
  }
 
 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;
    using TSQR::Epetra::makeTsqrMessenger;
    typedef TSQR::MessengerBase<scalar_type> base_mess_type;
 
    // If mv falls out of scope, its Epetra_Comm may become invalid.
    // Thus, we clone the input Epetra_Comm, so that the messenger
    // owns the object.
    RCP<const Epetra_Comm> comm  = rcp(mv.Comm().Clone());
    RCP<base_mess_type> messBase = makeTsqrMessenger<scalar_type>(comm);
    distTsqr_->init(messBase);
  }
};
 
}  // namespace Epetra
 
#endif  // defined(TPETRA_ENABLE_DEPRECATED_CODE) && defined(HAVE_TPETRA_EPETRA) && defined(HAVE_TPETRA_TSQR)
 
#endif  // EPETRA_TSQRADAPTOR_HPP