Epetra_TsqrAdaptor.hpp

Go to the documentation of this file.

// @HEADER
// ***********************************************************************
// 
//          Tpetra: Templated Linear Algebra Services Package
//                 Copyright (2008) 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 __Epetra_TsqrAdaptor_hpp
#define __Epetra_TsqrAdaptor_hpp


#include <Tpetra_ConfigDefs.hpp>

#if defined(HAVE_TPETRA_EPETRA) && defined(HAVE_TPETRA_TSQR)

#include <Kokkos_DefaultNode.hpp> // Include minimal Kokkos Node types
#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 {

  class TsqrAdaptor : public Teuchos::ParameterListAcceptorDefaultBase {
  public:
    typedef Epetra_MultiVector MV;

    typedef double scalar_type;

    typedef int ordinal_type;

    typedef KokkosClassic::SerialNode node_type;

    typedef Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> dense_matrix_type;

    typedef double magnitude_type;

  private:
    typedef TSQR::MatView<ordinal_type, scalar_type> matview_type;
    typedef TSQR::NodeTsqrFactory<node_type, scalar_type, ordinal_type> node_tsqr_factory_type;
    // Don't need a "typename" here, because there are no template
    // parameters involved in the type definition.
    typedef node_tsqr_factory_type::node_tsqr_type node_tsqr_type;
    typedef TSQR::DistTsqr<ordinal_type, scalar_type> dist_tsqr_type;
    typedef TSQR::Tsqr<ordinal_type, scalar_type, node_tsqr_type> tsqr_type;

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

    TsqrAdaptor () : 
      nodeTsqr_ (new node_tsqr_type),
      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)
    {
      typedef KokkosClassic::MultiVector<scalar_type, node_type> KMV;

      prepareTsqr (Q); // Finish initializing TSQR.
      KMV A_view = getNonConstView (A);
      KMV Q_view = getNonConstView (Q);
      tsqr_->factorExplicit (A_view, Q_view, R, false,
                             forceNonnegativeDiagonal);
    }

    int
    revealRank (MV& Q,
                dense_matrix_type& R,
                const magnitude_type& tol)
    {
      typedef KokkosClassic::MultiVector<scalar_type, node_type> KMV;

      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.
      KMV Q_view = getNonConstView (Q);
      return tsqr_->revealRank (Q_view, R, 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);
        prepareNodeTsqr (mv);
        ready_ = true;
      }
    }

    void
    prepareNodeTsqr (const MV& mv)
    {
      (void) mv; // Epetra objects don't have a Kokkos Node.

      // KokkosClassic::SerialNode wants an empty ParameterList.
      Teuchos::ParameterList plist;
      Teuchos::RCP<node_type> node (new node_type (plist));
      node_tsqr_factory_type::prepareNodeTsqr (nodeTsqr_, node);
    }

    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);
    }

    static KokkosClassic::MultiVector<scalar_type, node_type>
    getNonConstView (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.ConstantStride(), std::invalid_argument,
                                 "TSQR does not currently support Epetra_MultiVector "
                                 "inputs that do not have constant stride.");
      const int numRows = A.MyLength();
      const int numCols = A.NumVectors();
      const int stride  = A.Stride();
      // A_ptr does _not_ own the data.  TSQR only operates within the
      // scope of the multivector objects on which it operates, so it
      // doesn't need ownership of the data.
      Teuchos::ArrayRCP<double> A_ptr (A.Values(), 0, numRows*stride, false);

      typedef KokkosClassic::MultiVector<scalar_type, node_type> KMV;
      // KMV objects want a Kokkos Node instance.  Epetra objects
      // don't have a Kokkos Node, so we make a temporary node just
      // for the KMV.
      //
      // KokkosClassic::SerialNode wants an empty ParameterList.
      Teuchos::ParameterList plist;
      Teuchos::RCP<node_type> node (new node_type (plist));
      KMV A_kmv (node);
      A_kmv.initializeValues (numRows, numCols, A_ptr, stride);
      return A_kmv;
    }
  };

} // namespace Epetra

#endif // defined(HAVE_TPETRA_EPETRA) && defined(HAVE_TPETRA_TSQR)

#endif // __Epetra_TsqrAdaptor_hpp