Tpetra_MultiVector_decl.hpp

// @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 TPETRA_MULTIVECTOR_DECL_HPP
#define TPETRA_MULTIVECTOR_DECL_HPP

#include <Teuchos_LabeledObject.hpp>
#include <Teuchos_DataAccess.hpp>
#include <Teuchos_BLAS_types.hpp>
#include <Teuchos_Range1D.hpp>

#include <Kokkos_MultiVector.hpp>
#include <Kokkos_DefaultArithmetic.hpp>

#include "Tpetra_ConfigDefs.hpp"
#include "Tpetra_DistObject.hpp"
#include "Tpetra_Map.hpp"
#include "Tpetra_ViewAccepter.hpp"

// TODO: add principal use case instructions for memory management interfaces (view/copy extraction)
// TODO: expand user-visible documentation

namespace Tpetra {

#ifndef DOXYGEN_SHOULD_SKIP_THIS
  // forward declaration of Vector, needed to prevent circular inclusions
  template<class S, class LO, class GO, class N> class Vector;

  //template<class S, class LO, class GO, class N> class MultiVector;

  //template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  //RCP< MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
  //createMultiVectorFromView(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map,
  //                          const ArrayRCP<Scalar> &view, size_t LDA, size_t numVectors);
#endif


  template<class Scalar,
           class LocalOrdinal=int,
           class GlobalOrdinal=LocalOrdinal,
           class Node=Kokkos::DefaultNode::DefaultNodeType>
  class MultiVector :
    public DistObject<Scalar, LocalOrdinal, GlobalOrdinal, Node>
  {
  public:
    typedef Scalar        scalar_type;
    typedef LocalOrdinal  local_ordinal_type;
    typedef GlobalOrdinal global_ordinal_type;
    typedef Node          node_type;




    MultiVector (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& map,
                 size_t NumVectors,
                 bool zeroOut=true);

    MultiVector (const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &source);


    MultiVector (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& map,
                 const Teuchos::ArrayView<const Scalar>& A,
                 size_t LDA,
                 size_t NumVectors);


    MultiVector (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& map,
                 const Teuchos::ArrayView<const Teuchos::ArrayView<const Scalar> >&ArrayOfPtrs,
                 size_t NumVectors);

    virtual ~MultiVector();




    void
    replaceGlobalValue (GlobalOrdinal globalRow,
                        size_t vectorIndex,
                        const Scalar &value);

    void
    sumIntoGlobalValue (GlobalOrdinal globalRow,
                        size_t vectorIndex,
                        const Scalar &value);

    void
    replaceLocalValue (LocalOrdinal myRow,
                       size_t vectorIndex,
                       const Scalar &value);

    void
    sumIntoLocalValue (LocalOrdinal myRow,
                       size_t vectorIndex,
                       const Scalar &value);

    void putScalar (const Scalar &value);

    void randomize();

    void replaceMap(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map);

    void reduce();

    MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>&
    operator= (const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& source);




    Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
    subCopy (const Teuchos::Range1D &colRng) const;

    Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
    subCopy (const Teuchos::ArrayView<const size_t> &cols) const;

    Teuchos::RCP<const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
    subView (const Teuchos::Range1D &colRng) const;

    Teuchos::RCP<const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
    subView (const Teuchos::ArrayView<const size_t> &cols) const;

    Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
    subViewNonConst (const Teuchos::Range1D &colRng);

    Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
    subViewNonConst (const Teuchos::ArrayView<const size_t> &cols);


    Teuchos::RCP<const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
    offsetView (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& subMap,
                size_t offset) const;


    Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
    offsetViewNonConst (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &subMap,
                        size_t offset);

    Teuchos::RCP<const Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
    getVector (size_t j) const;

    Teuchos::RCP<Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
    getVectorNonConst (size_t j);

    Teuchos::ArrayRCP<const Scalar> getData(size_t j) const;

    Teuchos::ArrayRCP<Scalar> getDataNonConst(size_t j);

    void get1dCopy (Teuchos::ArrayView<Scalar> A, size_t LDA) const;

    void get2dCopy (Teuchos::ArrayView<const Teuchos::ArrayView<Scalar> > ArrayOfPtrs) const;

    Teuchos::ArrayRCP<const Scalar> get1dView() const;

    Teuchos::ArrayRCP<Teuchos::ArrayRCP<const Scalar> > get2dView() const;

    Teuchos::ArrayRCP<Scalar> get1dViewNonConst();

    Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar> > get2dViewNonConst();

    const Kokkos::MultiVector<Scalar,Node> & getLocalMV() const;

    Kokkos::MultiVector<Scalar,Node> & getLocalMVNonConst();




    void dot(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A, const Teuchos::ArrayView<Scalar> &dots) const;

    void abs(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A);

    void reciprocal(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A);

    void scale(const Scalar &alpha);

    void scale(Teuchos::ArrayView<const Scalar> alpha);

    void scale(const Scalar &alpha, const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A);

    void update(const Scalar &alpha, const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A, const Scalar &beta);

    void update(const Scalar &alpha, const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A, const Scalar &beta, const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &B, const Scalar &gamma);

    void norm1(const Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &norms) const;

    void norm2(const Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &norms) const;

    void normInf(const Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &norms) const;

    void normWeighted(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &weights, const Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &norms) const;

    void meanValue(const Teuchos::ArrayView<Scalar> &means) const;

    void
    multiply (Teuchos::ETransp transA,
              Teuchos::ETransp transB,
              const Scalar& alpha,
              const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
              const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& B,
              const Scalar& beta);


    void
    elementWiseMultiply (Scalar scalarAB,
                         const Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
                         const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& B,
                         Scalar scalarThis);



    size_t getNumVectors() const;

    size_t getLocalLength() const;

    global_size_t getGlobalLength() const;

    size_t getStride() const;

    bool isConstantStride() const;




    std::string description() const;

    void
    describe (Teuchos::FancyOStream& out,
              const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const;


  protected:

    typedef Kokkos::MultiVector<Scalar,Node> KMV;
    typedef Kokkos::DefaultArithmetic<KMV>   MVT;

    KMV lclMV_;

    Array<size_t> whichVectors_;



    template <class S,class LO,class GO,class N>
    friend RCP<MultiVector<S,LO,GO,N> >
    createMultiVectorFromView (const Teuchos::RCP<const Map<LO,GO,N> >&, const Teuchos::ArrayRCP<S>&, size_t, size_t);

    MultiVector (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& map,
                 const Teuchos::ArrayRCP<Scalar>& view,
                 size_t LDA,
                 size_t NumVectors,
                 EPrivateHostViewConstructor /* dummy */);

    inline bool vectorIndexOutOfRange(size_t VectorIndex) const {
      return (VectorIndex < 1 && VectorIndex != 0) || VectorIndex >= getNumVectors();
    }

    template <class T>
    ArrayRCP<T> getSubArrayRCP(ArrayRCP<T> arr, size_t j) const;

    MultiVector (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& map,
                 Teuchos::ArrayRCP<Scalar> data,
                 size_t LDA,
                 Teuchos::ArrayView<const size_t> whichVectors,
                 EPrivateComputeViewConstructor /* dummy */);

    MultiVector (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& map,
                 Teuchos::ArrayRCP<Scalar> data,
                 size_t LDA,
                 size_t NumVectors,
                 EPrivateComputeViewConstructor /* dummy */);




    bool
    checkSizes (const DistObject<Scalar,LocalOrdinal,GlobalOrdinal,Node>& sourceObj);

    void
    copyAndPermute (const DistObject<Scalar,LocalOrdinal,GlobalOrdinal,Node>& sourceObj,
                    size_t numSameIDs,
                    const ArrayView<const LocalOrdinal>& permuteToLIDs,
                    const ArrayView<const LocalOrdinal>& permuteFromLIDs);

    void
    packAndPrepare (const DistObject<Scalar,LocalOrdinal,GlobalOrdinal,Node>& sourceObj,
                    const ArrayView<const LocalOrdinal>& exportLIDs,
                    Array<Scalar>& exports,
                    const ArrayView<size_t>& numExportPacketsPerLID,
                    size_t& constantNumPackets,
                    Distributor& distor);

    void
    unpackAndCombine (const ArrayView<const LocalOrdinal>& importLIDs,
                      const ArrayView<const Scalar>& imports,
                      const ArrayView<size_t>& numPacketsPerLID,
                      size_t constantNumPackets,
                      Distributor& distor,
                      CombineMode CM);

    void createViews () const;
    void createViewsNonConst (Kokkos::ReadWriteOption rwo);
    void releaseViews () const;

    mutable ArrayRCP<Scalar> ncview_;

    mutable ArrayRCP<const Scalar> cview_;
  }; // class MultiVector

  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
  createMultiVector (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& map,
                     size_t numVectors)
  {
    using Teuchos::rcp;
    typedef MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;

    const bool initToZero = true;
    return rcp (new MV (map, numVectors, initToZero));
  }

  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
  createMultiVectorFromView (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& map,
                             const Teuchos::ArrayRCP<Scalar>& view,
                             size_t LDA,
                             size_t numVectors)
  {
    using Teuchos::rcp;
    typedef Tpetra::details::ViewAccepter<Node> VAN;
    typedef MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;

    // This uses a protected MultiVector constructor, but this
    // nonmember function was declared a friend of MultiVector.
    //
    // The ViewAccepter expression will fail to compile for
    // unsupported Kokkos Node types.
    return rcp (new MV (map, VAN::template acceptView<Scalar> (view),
                        LDA, numVectors, HOST_VIEW_CONSTRUCTOR));
  }

} // namespace Tpetra


#endif // TPETRA_MULTIVECTOR_DECL_HPP