Tpetra_Export.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_EXPORT_HPP
#define TPETRA_EXPORT_HPP

#include <Kokkos_DefaultNode.hpp>
#include <Teuchos_Describable.hpp>
#include <Teuchos_as.hpp>
#include "Tpetra_Map.hpp"
#include "Tpetra_Util.hpp"
#include "Tpetra_ImportExportData.hpp"
#include <iterator>

namespace Tpetra {

  template <class LocalOrdinal, class GlobalOrdinal = LocalOrdinal, class Node = Kokkos::DefaultNode::DefaultNodeType>
  class Export: public Teuchos::Describable {

  public:
    typedef Map<LocalOrdinal,GlobalOrdinal,Node> map_type;



    Export (const Teuchos::RCP<const map_type>& source, 
            const Teuchos::RCP<const map_type>& target);

    Export (const Teuchos::RCP<const map_type>& source,
            const Teuchos::RCP<const map_type>& target,
            const Teuchos::RCP<Teuchos::ParameterList>& plist);

    Export (const Export<LocalOrdinal,GlobalOrdinal,Node>& rhs);

    ~Export();





    inline size_t getNumSameIDs() const;

    inline size_t getNumPermuteIDs() const;

    inline ArrayView<const LocalOrdinal> getPermuteFromLIDs() const;

    inline ArrayView<const LocalOrdinal> getPermuteToLIDs() const;

    inline size_t getNumRemoteIDs() const;

    inline ArrayView<const LocalOrdinal> getRemoteLIDs() const;

    inline size_t getNumExportIDs() const;

    inline ArrayView<const LocalOrdinal> getExportLIDs() const;

    inline ArrayView<const int> getExportImageIDs() const;

    inline const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getSourceMap() const;

    inline const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getTargetMap() const;

    inline Distributor & getDistributor() const;

    Export<LocalOrdinal,GlobalOrdinal,Node>& 
    operator= (const Export<LocalOrdinal,GlobalOrdinal,Node>& rhs);




    virtual void print (std::ostream& os) const;


  private:

    RCP<ImportExportData<LocalOrdinal,GlobalOrdinal,Node> > ExportData_;



    //==============================================================================
    // sets up numSameIDs_, numPermuteIDs_, and the export IDs
    // these variables are already initialized to 0 by the ImportExportData ctr.
    // also sets up permuteToLIDs_, permuteFromLIDs_, exportGIDs_, and exportLIDs_
    void setupSamePermuteExport();
    void setupRemote();

  };


  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  Export<LocalOrdinal,GlobalOrdinal,Node>::
  Export (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& source,   
          const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& target)
  {
    using Teuchos::rcp;
    typedef ImportExportData<LocalOrdinal,GlobalOrdinal,Node> data_type;

    ExportData_ = rcp (new data_type (source, target));
    setupSamePermuteExport();
    if (source->isDistributed()) {
      setupRemote();
    }
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  Export<LocalOrdinal,GlobalOrdinal,Node>::
  Export (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& source,   
          const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& target,
          const Teuchos::RCP<Teuchos::ParameterList>& plist)
  {
    using Teuchos::rcp;
    typedef ImportExportData<LocalOrdinal,GlobalOrdinal,Node> data_type;

    ExportData_ = rcp (new data_type (source, target, plist));
    setupSamePermuteExport();
    if (source->isDistributed()) {
      setupRemote();
    }
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  Export<LocalOrdinal,GlobalOrdinal,Node>::Export(const Export<LocalOrdinal,GlobalOrdinal,Node> & rhs)
  : ExportData_(rhs.ExportData_)
  {}

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  Export<LocalOrdinal,GlobalOrdinal,Node>::~Export() 
  {}

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  size_t Export<LocalOrdinal,GlobalOrdinal,Node>::getNumSameIDs() const {
    return ExportData_->numSameIDs_;
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  size_t Export<LocalOrdinal,GlobalOrdinal,Node>::getNumPermuteIDs() const {
    return ExportData_->permuteFromLIDs_.size();
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  ArrayView<const LocalOrdinal> 
  Export<LocalOrdinal,GlobalOrdinal,Node>::getPermuteFromLIDs() const {
    return ExportData_->permuteFromLIDs_();
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  ArrayView<const LocalOrdinal>
  Export<LocalOrdinal,GlobalOrdinal,Node>::getPermuteToLIDs() const {
    return ExportData_->permuteToLIDs_();
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  size_t Export<LocalOrdinal,GlobalOrdinal,Node>::getNumRemoteIDs() const {
    return ExportData_->remoteLIDs_.size();
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  ArrayView<const LocalOrdinal> 
  Export<LocalOrdinal,GlobalOrdinal,Node>::getRemoteLIDs() const {
    return ExportData_->remoteLIDs_();
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  size_t Export<LocalOrdinal,GlobalOrdinal,Node>::getNumExportIDs() const {
    return ExportData_->exportLIDs_.size();
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  ArrayView<const LocalOrdinal> 
  Export<LocalOrdinal,GlobalOrdinal,Node>::getExportLIDs() const {
    return ExportData_->exportLIDs_();
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  ArrayView<const int> 
  Export<LocalOrdinal,GlobalOrdinal,Node>::getExportImageIDs() const {
    return ExportData_->exportImageIDs_();
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & 
  Export<LocalOrdinal,GlobalOrdinal,Node>::getSourceMap() const {
    return ExportData_->source_;
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & 
  Export<LocalOrdinal,GlobalOrdinal,Node>::getTargetMap() const {
    return ExportData_->target_;
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  Distributor & 
  Export<LocalOrdinal,GlobalOrdinal,Node>::getDistributor() const { 
    return ExportData_->distributor_;
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  Export<LocalOrdinal,GlobalOrdinal,Node>& 
  Export<LocalOrdinal,GlobalOrdinal,Node>::operator=(const Export<LocalOrdinal,GlobalOrdinal,Node> & rhs) {
    if (&rhs != this) {
      // It's bad form to clobber your own data in a self-assignment.
      // This can result in dangling pointers if some member data are
      // raw pointers that the class deallocates in the constructor.
      // It doesn't matter in this case, because ExportData_ is an
      // RCP, which defines self-assignment sensibly.  Nevertheless,
      // we include the check for self-assignment, because it's good
      // form and not expensive (just a raw pointer comparison).
      ExportData_ = rhs.ExportData_;
    }
    return *this;
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  void Export<LocalOrdinal,GlobalOrdinal,Node>::print(std::ostream& os) const {
    using Teuchos::getFancyOStream;
    using Teuchos::rcpFromRef;
    using std::endl;

    ArrayView<const LocalOrdinal> av;
    ArrayView<const int> avi;
    const RCP<const Comm<int> > & comm = getSourceMap()->getComm();
    const int myImageID = comm->getRank();
    const int numImages = comm->getSize();
    for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
      if (myImageID == imageCtr) {
        os << endl;
        if (myImageID == 0) { // this is the root node (only output this info once)
          os << "Export Data Members:" << endl;
        }
        os << "Image ID       : " << myImageID << endl;
        os << "permuteFromLIDs: {"; av = getPermuteFromLIDs(); std::copy(av.begin(),av.end(),std::ostream_iterator<LocalOrdinal>(os," ")); os << " }" << endl;
        os << "permuteToLIDs  : {"; av = getPermuteToLIDs();   std::copy(av.begin(),av.end(),std::ostream_iterator<LocalOrdinal>(os," ")); os << " }" << endl;
        os << "remoteLIDs     : {"; av = getRemoteLIDs();      std::copy(av.begin(),av.end(),std::ostream_iterator<LocalOrdinal>(os," ")); os << " }" << endl;
        os << "exportLIDs     : {"; av = getExportLIDs();      std::copy(av.begin(),av.end(),std::ostream_iterator<LocalOrdinal>(os," ")); os << " }" << endl;
        os << "exportImageIDs : {"; avi = getExportImageIDs();  std::copy(avi.begin(),avi.end(),std::ostream_iterator<int>(os," ")); os << " }" << endl;
        os << "numSameIDs     : " << getNumSameIDs() << endl;
        os << "numPermuteIDs  : " << getNumPermuteIDs() << endl;
        os << "numRemoteIDs   : " << getNumRemoteIDs() << endl;
        os << "numExportIDs   : " << getNumExportIDs() << endl;
      }
      // Do a few global ops to give I/O a chance to complete
      comm->barrier();
      comm->barrier();
      comm->barrier();
    }
    if (myImageID == 0) {
      os << endl << endl << "Source Map:" << endl << std::flush; 
    }
    comm->barrier();
    os << *getSourceMap();
    comm->barrier();

    if (myImageID == 0) {
      os << endl << endl << "Target Map:" << endl << std::flush; 
    }
    comm->barrier();
    os << *getTargetMap();
    comm->barrier();

    // It's also helpful for debugging to print the Distributor
    // object.  Epetra_Import::Print() does this (or _should_ do this,
    // but doesn't, as of 05 Jan 2012), so we can do a side-by-side
    // comparison.
    if (myImageID == 0) {
      os << endl << endl << "Distributor:" << endl << std::flush;
    }
    comm->barrier();
    getDistributor().describe (*(getFancyOStream (rcpFromRef (os))),
                               Teuchos::VERB_EXTREME);
  }


  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  void 
  Export<LocalOrdinal,GlobalOrdinal,Node>::setupSamePermuteExport() 
  {
    const Map<LocalOrdinal,GlobalOrdinal,Node> & source = *getSourceMap();
    const Map<LocalOrdinal,GlobalOrdinal,Node> & target = *getTargetMap();
    ArrayView<const GlobalOrdinal> sourceGIDs = source.getNodeElementList();
    ArrayView<const GlobalOrdinal> targetGIDs = target.getNodeElementList();

    // Compute numSameIDs_:
    //
    // Iterate through the source and target GID lists.  If the i-th
    // GID of both is the same, increment numSameIDs_ and try the
    // next.  As soon as you come to a nonmatching pair, give up.
    typename ArrayView<const GlobalOrdinal>::iterator sourceIter = sourceGIDs.begin(),
                                                      targetIter = targetGIDs.begin();
    while (sourceIter != sourceGIDs.end() && targetIter != targetGIDs.end() && *sourceIter == *targetIter) {
      ++ExportData_->numSameIDs_;
      ++sourceIter;
      ++targetIter;
    }
    // sourceIter should now point either to the GID of the first
    // non-same entry in sourceGIDs, or to the end of sourceGIDs (if
    // all the entries were the same).

    // -- compute numPermuteIDs --
    // -- fill permuteToLIDs_, permuteFromLIDs_ --
    //
    // Iterate through the remaining entries in sourceGIDs.  (There
    // may not be any, if all pairs of entries in sourceGIDs and
    // targetGIDs were the same.)  If target owns that GID, add
    // entries to permuteToLIDs_ and permuteFromLIDs_.  Otherwise, add
    // entries to exportGIDs_.
    //
    for (; sourceIter != sourceGIDs.end(); ++sourceIter) {
      if (target.isNodeGlobalElement(*sourceIter)) {
        // The current process owns this GID, for both the source and
        // the target Maps.  Determine the LIDs for this GID on both
        // Maps and add them to the permutation lists.
        ExportData_->permuteToLIDs_.push_back(  target.getLocalElement(*sourceIter));
        ExportData_->permuteFromLIDs_.push_back(source.getLocalElement(*sourceIter));
      }
      else {
        // The current GID is owned by this process in the source Map,
        // but is not owned by this process in the target Map.  Store
        // such GIDs.
        ExportData_->exportGIDs_.push_back(*sourceIter);
      }
    }

    // Above, we filled exportGIDs_ with all the GIDs which we own in
    // the source Map, but not in the target Map.  Now allocate
    // exportLIDs_, and fill it with the LIDs (from the source Map)
    // corresponding to those GIDs.
    if (ExportData_->exportGIDs_.size()) {
      ExportData_->exportLIDs_ = arcp<LocalOrdinal>(ExportData_->exportGIDs_.size());
    }
    {
      typename ArrayRCP<LocalOrdinal>::iterator liditer = ExportData_->exportLIDs_.begin();
      typename Array<GlobalOrdinal>::iterator   giditer = ExportData_->exportGIDs_.begin();
      for (; giditer != ExportData_->exportGIDs_.end(); ++liditer, ++giditer) {
        *liditer = source.getLocalElement(*giditer);
      }
    }

    TPETRA_ABUSE_WARNING( (getNumExportIDs() > 0) && (!source.isDistributed()), std::runtime_error,
        "::setupSamePermuteExport(): Source has export LIDs but Source is not distributed globally."
        << std::endl << "Importing to a submap of the target map.");

    // -- compute exportImageIDs_ --
    //
    // For each GID in exportGIDs_, find its corresponding owning
    // image (a.k.a. "process," "node") ID in the target Map.  Store
    // these image IDs in exportImageIDs_.  These are the image IDs to
    // which the Export needs to send data.
    //
    // We only need to do this if the source Map is distributed;
    // otherwise, the Export doesn't have to perform any
    // communication.
    if (source.isDistributed()) {
      ExportData_->exportImageIDs_ = arcp<int>(ExportData_->exportGIDs_.size());
      const LookupStatus lookup = target.getRemoteIndexList(ExportData_->exportGIDs_(), ExportData_->exportImageIDs_());
      TPETRA_ABUSE_WARNING( lookup == IDNotPresent, std::runtime_error, 
        "::setupSamePermuteExport(): The source Map has GIDs not found in the target Map.");

      // Get rid of IDs not in the Target Map
      typedef typename ArrayRCP<int>::difference_type size_type;
      if (lookup == IDNotPresent) {
        const size_type numInvalidExports = std::count_if( ExportData_->exportImageIDs_().begin(), 
                                                          ExportData_->exportImageIDs_().end(), 
                                                          std::bind1st(std::equal_to<int>(),-1) 
                                                        );
        // count number of valid and total number of exports
        const size_type totalNumExports = ExportData_->exportImageIDs_.size();
        if (numInvalidExports == totalNumExports) {
          // all exports are invalid; we have no exports; we can delete all exports
          ExportData_->exportGIDs_.resize(0);
          ExportData_->exportLIDs_ = null;
          ExportData_->exportImageIDs_ = null;
        }
        else {
          // some exports are valid; we need to keep the valid exports
          // pack and resize
          size_type numValidExports = 0;
          for (size_type e = 0; e < totalNumExports; ++e) {
            if (ExportData_->exportImageIDs_[e] != -1) {
              ExportData_->exportGIDs_[numValidExports]     = ExportData_->exportGIDs_[e];
              ExportData_->exportLIDs_[numValidExports]     = ExportData_->exportLIDs_[e];
              ExportData_->exportImageIDs_[numValidExports] = ExportData_->exportImageIDs_[e];
              ++numValidExports;
            }
          }
          ExportData_->exportGIDs_.resize(numValidExports);
          ExportData_->exportLIDs_.resize(numValidExports);
          ExportData_->exportImageIDs_.resize(numValidExports);
        }
      }
    }
  } // setupSamePermuteExport()


  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  void 
  Export<LocalOrdinal,GlobalOrdinal,Node>::setupRemote() 
  {
    const Map<LocalOrdinal,GlobalOrdinal,Node>& target = *getTargetMap();

    // Make sure the export IDs are ordered by image.  Sort
    // exportImageIDs_ in ascending order, and apply the same
    // permutation to exportGIDs_ and exportLIDs_.
    sort3 (ExportData_->exportImageIDs_.begin(), 
           ExportData_->exportImageIDs_.end(), 
           ExportData_->exportGIDs_.begin(), 
           ExportData_->exportLIDs_.begin());

    // Construct the list of entries that calling image needs to send
    // as a result of everyone asking for what it needs to receive.
    //
    // mfh 05 Jan 2012: I understand the above comment as follows:
    // Construct the communication plan from the list of image IDs to
    // which we need to send.
    size_t numRemoteIDs;
    numRemoteIDs = ExportData_->distributor_.createFromSends (ExportData_->exportImageIDs_());

    // Use the communication plan with ExportGIDs to find out who is
    // sending to us and get the proper ordering of GIDs for incoming
    // remote entries (these will be converted to LIDs when done).
    Array<GlobalOrdinal> remoteGIDs(numRemoteIDs);
    ExportData_->distributor_.doPostsAndWaits (ExportData_->exportGIDs_().getConst(), 1, remoteGIDs());

    // Remote IDs come in as GIDs; convert to LIDs.  LIDs tell us
    // where to store the incoming remote data.
    ExportData_->remoteLIDs_.resize(numRemoteIDs);
    {
      typename Array<GlobalOrdinal>::const_iterator i = remoteGIDs.begin();
      typename Array<LocalOrdinal>::iterator        j = ExportData_->remoteLIDs_.begin();
      while (i != remoteGIDs.end()) {
        *j++ = target.getLocalElement(*i++);
      }
    }
  }

  template <class LO, class GO, class Node> 
  RCP< const Export<LO,GO,Node> >
  createExport( const RCP<const Map<LO,GO,Node> > & src, 
                const RCP<const Map<LO,GO,Node> > & tgt )
  {
    if (src == tgt) return null;
#ifdef HAVE_TPETRA_DEBUG
    TEUCHOS_TEST_FOR_EXCEPTION(src == null || tgt == null, std::runtime_error,
        "Tpetra::createExport(): neither source nor target map may be null:\nsource: " << src << "\ntarget: " << tgt << "\n");
#endif
    return rcp(new Export<LO,GO,Node>(src,tgt));
  }

} // namespace Tpetra

#endif // TPETRA_EXPORT_HPP