Tpetra Matrix/Vector Services: Tpetra_MultiVector

00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //          Tpetra: Templated Linear Algebra Services Package
00005 //                 Copyright (2008) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 // @HEADER
00041 
00042 #ifndef TPETRA_MULTIVECTOR_DECL_HPP
00043 #define TPETRA_MULTIVECTOR_DECL_HPP
00044 
00045 #include <Teuchos_LabeledObject.hpp>
00046 #include <Teuchos_DataAccess.hpp>
00047 #include <Teuchos_BLAS_types.hpp>
00048 #include <Teuchos_Range1D.hpp>
00049 
00050 #include <Kokkos_MultiVector.hpp>
00051 #include <Kokkos_DefaultArithmetic.hpp>
00052 
00053 #include "Tpetra_ConfigDefs.hpp"
00054 #include "Tpetra_DistObject.hpp"
00055 #include "Tpetra_Map.hpp"
00056 #include "Tpetra_ViewAccepter.hpp"
00057 
00058 // TODO: add principal use case instructions for memory management interfaces (view/copy extraction)
00059 // TODO: expand user-visible documentation
00060 
00061 namespace Tpetra {
00062 
00063 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00064   // forward declaration of Vector, needed to prevent circular inclusions
00065   template<class S, class LO, class GO, class N> class Vector;
00066 
00067   //template<class S, class LO, class GO, class N> class MultiVector;
00068 
00069   //template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00070   //RCP< MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
00071   //createMultiVectorFromView(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map,
00072   //                          const ArrayRCP<Scalar> &view, size_t LDA, size_t numVectors);
00073 #endif
00074 
00323   template<class Scalar,
00324            class LocalOrdinal=int,
00325            class GlobalOrdinal=LocalOrdinal,
00326            class Node=Kokkos::DefaultNode::DefaultNodeType>
00327   class MultiVector :
00328     public DistObject<Scalar, LocalOrdinal, GlobalOrdinal, Node>
00329   {
00330   public:
00332 
00333 
00335     typedef Scalar        scalar_type;
00337     typedef LocalOrdinal  local_ordinal_type;
00339     typedef GlobalOrdinal global_ordinal_type;
00341     typedef Node          node_type;
00342 
00344 
00345 
00346 
00353     MultiVector (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& map,
00354                  size_t NumVectors,
00355                  bool zeroOut=true);
00356 
00358     MultiVector (const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &source);
00359 
00375     MultiVector (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& map,
00376                  const Teuchos::ArrayView<const Scalar>& A,
00377                  size_t LDA,
00378                  size_t NumVectors);
00379 
00393     MultiVector (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& map,
00394                  const Teuchos::ArrayView<const Teuchos::ArrayView<const Scalar> >&ArrayOfPtrs,
00395                  size_t NumVectors);
00396 
00398     virtual ~MultiVector();
00399 
00401 
00402 
00403 
00412     void
00413     replaceGlobalValue (GlobalOrdinal globalRow,
00414                         size_t vectorIndex,
00415                         const Scalar &value);
00416 
00425     void
00426     sumIntoGlobalValue (GlobalOrdinal globalRow,
00427                         size_t vectorIndex,
00428                         const Scalar &value);
00429 
00438     void
00439     replaceLocalValue (LocalOrdinal myRow,
00440                        size_t vectorIndex,
00441                        const Scalar &value);
00442 
00451     void
00452     sumIntoLocalValue (LocalOrdinal myRow,
00453                        size_t vectorIndex,
00454                        const Scalar &value);
00455 
00457     void putScalar (const Scalar &value);
00458 
00474     void randomize();
00475 
00501     void replaceMap(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map);
00502 
00509     void reduce();
00510 
00525     MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>&
00526     operator= (const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& source);
00527 
00529 
00531 
00539 
00541     Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
00542     subCopy (const Teuchos::Range1D &colRng) const;
00543 
00545     Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
00546     subCopy (const Teuchos::ArrayView<const size_t> &cols) const;
00547 
00549     Teuchos::RCP<const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
00550     subView (const Teuchos::Range1D &colRng) const;
00551 
00553     Teuchos::RCP<const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
00554     subView (const Teuchos::ArrayView<const size_t> &cols) const;
00555 
00557     Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
00558     subViewNonConst (const Teuchos::Range1D &colRng);
00559 
00561     Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
00562     subViewNonConst (const Teuchos::ArrayView<const size_t> &cols);
00563 
00565 
00576     Teuchos::RCP<const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
00577     offsetView (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& subMap,
00578                 size_t offset) const;
00579 
00581 
00592     Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
00593     offsetViewNonConst (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &subMap,
00594                         size_t offset);
00595 
00597     Teuchos::RCP<const Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
00598     getVector (size_t j) const;
00599 
00601     Teuchos::RCP<Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
00602     getVectorNonConst (size_t j);
00603 
00605     Teuchos::ArrayRCP<const Scalar> getData(size_t j) const;
00606 
00608     Teuchos::ArrayRCP<Scalar> getDataNonConst(size_t j);
00609 
00616     void get1dCopy (Teuchos::ArrayView<Scalar> A, size_t LDA) const;
00617 
00623     void get2dCopy (Teuchos::ArrayView<const Teuchos::ArrayView<Scalar> > ArrayOfPtrs) const;
00624 
00630     Teuchos::ArrayRCP<const Scalar> get1dView() const;
00631 
00633     Teuchos::ArrayRCP<Teuchos::ArrayRCP<const Scalar> > get2dView() const;
00634 
00640     Teuchos::ArrayRCP<Scalar> get1dViewNonConst();
00641 
00643     Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar> > get2dViewNonConst();
00644 
00646     const Kokkos::MultiVector<Scalar,Node> & getLocalMV() const;
00647 
00649     Kokkos::MultiVector<Scalar,Node> & getLocalMVNonConst();
00650 
00652 
00654 
00655 
00657     void dot(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A, const Teuchos::ArrayView<Scalar> &dots) const;
00658 
00660     void abs(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A);
00661 
00663     void reciprocal(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A);
00664 
00666     void scale(const Scalar &alpha);
00667 
00669     void scale(Teuchos::ArrayView<const Scalar> alpha);
00670 
00672     void scale(const Scalar &alpha, const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A);
00673 
00675     void update(const Scalar &alpha, const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &A, const Scalar &beta);
00676 
00678     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);
00679 
00681     void norm1(const Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &norms) const;
00682 
00684     void norm2(const Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &norms) const;
00685 
00687     void normInf(const Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &norms) const;
00688 
00690     void normWeighted(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &weights, const Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &norms) const;
00691 
00693     void meanValue(const Teuchos::ArrayView<Scalar> &means) const;
00694 
00696     void
00697     multiply (Teuchos::ETransp transA,
00698               Teuchos::ETransp transB,
00699               const Scalar& alpha,
00700               const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
00701               const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& B,
00702               const Scalar& beta);
00703 
00705 
00713     void
00714     elementWiseMultiply (Scalar scalarAB,
00715                          const Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
00716                          const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& B,
00717                          Scalar scalarThis);
00719 
00721 
00722 
00724     size_t getNumVectors() const;
00725 
00727     size_t getLocalLength() const;
00728 
00730     global_size_t getGlobalLength() const;
00731 
00737     size_t getStride() const;
00738 
00742     bool isConstantStride() const;
00743 
00745 
00747 
00748 
00750     std::string description() const;
00751 
00780     void
00781     describe (Teuchos::FancyOStream& out,
00782               const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const;
00783 
00785 
00786   protected:
00787 
00788     typedef Kokkos::MultiVector<Scalar,Node> KMV;
00789     typedef Kokkos::DefaultArithmetic<KMV>   MVT;
00790 
00792     KMV lclMV_;
00793 
00806     Array<size_t> whichVectors_;
00807 
00809 
00810 
00811     template <class S,class LO,class GO,class N>
00812     friend RCP<MultiVector<S,LO,GO,N> >
00813     createMultiVectorFromView (const Teuchos::RCP<const Map<LO,GO,N> >&, const Teuchos::ArrayRCP<S>&, size_t, size_t);
00814 
00824     MultiVector (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& map,
00825                  const Teuchos::ArrayRCP<Scalar>& view,
00826                  size_t LDA,
00827                  size_t NumVectors,
00828                  EPrivateHostViewConstructor /* dummy */);
00829 
00830     inline bool vectorIndexOutOfRange(size_t VectorIndex) const {
00831       return (VectorIndex < 1 && VectorIndex != 0) || VectorIndex >= getNumVectors();
00832     }
00833 
00839     template <class T>
00840     ArrayRCP<T> getSubArrayRCP(ArrayRCP<T> arr, size_t j) const;
00841 
00843     MultiVector (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& map,
00844                  Teuchos::ArrayRCP<Scalar> data,
00845                  size_t LDA,
00846                  Teuchos::ArrayView<const size_t> whichVectors,
00847                  EPrivateComputeViewConstructor /* dummy */);
00848 
00850     MultiVector (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& map,
00851                  Teuchos::ArrayRCP<Scalar> data,
00852                  size_t LDA,
00853                  size_t NumVectors,
00854                  EPrivateComputeViewConstructor /* dummy */);
00855 
00857 
00858 
00859 
00864     bool
00865     checkSizes (const DistObject<Scalar,LocalOrdinal,GlobalOrdinal,Node>& sourceObj);
00866 
00867     void
00868     copyAndPermute (const DistObject<Scalar,LocalOrdinal,GlobalOrdinal,Node>& sourceObj,
00869                     size_t numSameIDs,
00870                     const ArrayView<const LocalOrdinal>& permuteToLIDs,
00871                     const ArrayView<const LocalOrdinal>& permuteFromLIDs);
00872 
00873     void
00874     packAndPrepare (const DistObject<Scalar,LocalOrdinal,GlobalOrdinal,Node>& sourceObj,
00875                     const ArrayView<const LocalOrdinal>& exportLIDs,
00876                     Array<Scalar>& exports,
00877                     const ArrayView<size_t>& numExportPacketsPerLID,
00878                     size_t& constantNumPackets,
00879                     Distributor& distor);
00880 
00881     void
00882     unpackAndCombine (const ArrayView<const LocalOrdinal>& importLIDs,
00883                       const ArrayView<const Scalar>& imports,
00884                       const ArrayView<size_t>& numPacketsPerLID,
00885                       size_t constantNumPackets,
00886                       Distributor& distor,
00887                       CombineMode CM);
00888 
00889     void createViews () const;
00890     void createViewsNonConst (Kokkos::ReadWriteOption rwo);
00891     void releaseViews () const;
00892 
00894     mutable ArrayRCP<Scalar> ncview_;
00895 
00897     mutable ArrayRCP<const Scalar> cview_;
00898 
00900 
00901 
00902 
00904     mutable bool releaseViewsRaisedEfficiencyWarning_;
00906     mutable bool createViewsRaisedEfficiencyWarning_;
00908     mutable bool createViewsNonConstRaisedEfficiencyWarning_;
00909 
00911   }; // class MultiVector
00912 
00920   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00921   Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
00922   createMultiVector (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& map,
00923                      size_t numVectors)
00924   {
00925     using Teuchos::rcp;
00926     typedef MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
00927 
00928     const bool initToZero = true;
00929     return rcp (new MV (map, numVectors, initToZero));
00930   }
00931 
00957   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00958   Teuchos::RCP<MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
00959   createMultiVectorFromView (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& map,
00960                              const Teuchos::ArrayRCP<Scalar>& view,
00961                              size_t LDA,
00962                              size_t numVectors)
00963   {
00964     using Teuchos::rcp;
00965     typedef Tpetra::details::ViewAccepter<Node> VAN;
00966     typedef MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
00967 
00968     // This uses a protected MultiVector constructor, but this
00969     // nonmember function was declared a friend of MultiVector.
00970     //
00971     // The ViewAccepter expression will fail to compile for
00972     // unsupported Kokkos Node types.
00973     return rcp (new MV (map, VAN::template acceptView<Scalar> (view),
00974                         LDA, numVectors, HOST_VIEW_CONSTRUCTOR));
00975   }
00976 
00977 } // namespace Tpetra
00978 
00979 
00980 #endif // TPETRA_MULTIVECTOR_DECL_HPP