Tpetra Matrix/Vector Services Version of the Day
Tpetra_Export.hpp
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_EXPORT_HPP
00043 #define TPETRA_EXPORT_HPP
00044 
00045 #include <Kokkos_DefaultNode.hpp>
00046 #include <Teuchos_Describable.hpp>
00047 #include <Teuchos_as.hpp>
00048 #include "Tpetra_Map.hpp"
00049 #include "Tpetra_Util.hpp"
00050 #include "Tpetra_ImportExportData.hpp"
00051 #include <iterator>
00052 
00053 namespace Tpetra {
00054 
00087   template <class LocalOrdinal, class GlobalOrdinal = LocalOrdinal, class Node = Kokkos::DefaultNode::DefaultNodeType>
00088   class Export: public Teuchos::Describable {
00089 
00090   public:
00092 
00093 
00101     Export (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& source, 
00102       const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& target);
00103 
00105     Export (const Export<LocalOrdinal,GlobalOrdinal,Node>& rhs);
00106 
00108     ~Export();
00109 
00111 
00113 
00114 
00116     inline size_t getNumSameIDs() const;
00117 
00119     inline size_t getNumPermuteIDs() const;
00120 
00122     inline ArrayView<const LocalOrdinal> getPermuteFromLIDs() const;
00123 
00125     inline ArrayView<const LocalOrdinal> getPermuteToLIDs() const;
00126 
00128     inline size_t getNumRemoteIDs() const;
00129 
00131     inline ArrayView<const LocalOrdinal> getRemoteLIDs() const;
00132 
00134     inline size_t getNumExportIDs() const;
00135 
00137     inline ArrayView<const LocalOrdinal> getExportLIDs() const;
00138 
00140     inline ArrayView<const int> getExportImageIDs() const;
00141 
00143     inline const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getSourceMap() const;
00144 
00146     inline const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getTargetMap() const;
00147 
00149     inline Distributor & getDistributor() const;
00150 
00152     Export<LocalOrdinal,GlobalOrdinal,Node>& 
00153     operator= (const Export<LocalOrdinal,GlobalOrdinal,Node>& rhs);
00154 
00156 
00158 
00159 
00175     virtual void print (std::ostream& os) const;
00176 
00178 
00179   private:
00180 
00181     RCP<ImportExportData<LocalOrdinal,GlobalOrdinal,Node> > ExportData_;
00182 
00184 
00185 
00186     //==============================================================================
00187     // sets up numSameIDs_, numPermuteIDs_, and the export IDs
00188     // these variables are already initialized to 0 by the ImportExportData ctr.
00189     // also sets up permuteToLIDs_, permuteFromLIDs_, exportGIDs_, and exportLIDs_
00190     void setupSamePermuteExport();
00191     void setupRemote();
00192 
00194   };
00195 
00196   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00197   Export<LocalOrdinal,GlobalOrdinal,Node>::Export(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & source,   
00198                                                   const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & target) {
00199     ExportData_ = rcp(new ImportExportData<LocalOrdinal,GlobalOrdinal,Node>(source, target));
00200     // call subfunctions
00201     setupSamePermuteExport();
00202     if (source->isDistributed()) {
00203       setupRemote();
00204     }
00205   }
00206 
00207   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00208   Export<LocalOrdinal,GlobalOrdinal,Node>::Export(const Export<LocalOrdinal,GlobalOrdinal,Node> & rhs)
00209   : ExportData_(rhs.ExportData_)
00210   {}
00211 
00212   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00213   Export<LocalOrdinal,GlobalOrdinal,Node>::~Export() 
00214   {}
00215 
00216   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00217   size_t Export<LocalOrdinal,GlobalOrdinal,Node>::getNumSameIDs() const {
00218     return ExportData_->numSameIDs_;
00219   }
00220 
00221   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00222   size_t Export<LocalOrdinal,GlobalOrdinal,Node>::getNumPermuteIDs() const {
00223     return ExportData_->permuteFromLIDs_.size();
00224   }
00225 
00226   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00227   ArrayView<const LocalOrdinal> 
00228   Export<LocalOrdinal,GlobalOrdinal,Node>::getPermuteFromLIDs() const {
00229     return ExportData_->permuteFromLIDs_();
00230   }
00231 
00232   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00233   ArrayView<const LocalOrdinal>
00234   Export<LocalOrdinal,GlobalOrdinal,Node>::getPermuteToLIDs() const {
00235     return ExportData_->permuteToLIDs_();
00236   }
00237 
00238   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00239   size_t Export<LocalOrdinal,GlobalOrdinal,Node>::getNumRemoteIDs() const {
00240     return ExportData_->remoteLIDs_.size();
00241   }
00242 
00243   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00244   ArrayView<const LocalOrdinal> 
00245   Export<LocalOrdinal,GlobalOrdinal,Node>::getRemoteLIDs() const {
00246     return ExportData_->remoteLIDs_();
00247   }
00248 
00249   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00250   size_t Export<LocalOrdinal,GlobalOrdinal,Node>::getNumExportIDs() const {
00251     return ExportData_->exportLIDs_.size();
00252   }
00253 
00254   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00255   ArrayView<const LocalOrdinal> 
00256   Export<LocalOrdinal,GlobalOrdinal,Node>::getExportLIDs() const {
00257     return ExportData_->exportLIDs_();
00258   }
00259 
00260   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00261   ArrayView<const int> 
00262   Export<LocalOrdinal,GlobalOrdinal,Node>::getExportImageIDs() const {
00263     return ExportData_->exportImageIDs_();
00264   }
00265 
00266   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00267   const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & 
00268   Export<LocalOrdinal,GlobalOrdinal,Node>::getSourceMap() const {
00269     return ExportData_->source_;
00270   }
00271 
00272   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00273   const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & 
00274   Export<LocalOrdinal,GlobalOrdinal,Node>::getTargetMap() const {
00275     return ExportData_->target_;
00276   }
00277 
00278   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00279   Distributor & 
00280   Export<LocalOrdinal,GlobalOrdinal,Node>::getDistributor() const { 
00281     return ExportData_->distributor_;
00282   }
00283 
00284   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00285   Export<LocalOrdinal,GlobalOrdinal,Node>& 
00286   Export<LocalOrdinal,GlobalOrdinal,Node>::operator=(const Export<LocalOrdinal,GlobalOrdinal,Node> & rhs) {
00287     if (&rhs != this) {
00288       // It's bad form to clobber your own data in a self-assignment.
00289       // This can result in dangling pointers if some member data are
00290       // raw pointers that the class deallocates in the constructor.
00291       // It doesn't matter in this case, because ExportData_ is an
00292       // RCP, which defines self-assignment sensibly.  Nevertheless,
00293       // we include the check for self-assignment, because it's good
00294       // form and not expensive (just a raw pointer comparison).
00295       ExportData_ = rhs.ExportData_;
00296     }
00297     return *this;
00298   }
00299 
00300   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00301   void Export<LocalOrdinal,GlobalOrdinal,Node>::print(std::ostream& os) const {
00302     using Teuchos::getFancyOStream;
00303     using Teuchos::rcpFromRef;
00304     using std::endl;
00305 
00306     ArrayView<const LocalOrdinal> av;
00307     ArrayView<const int> avi;
00308     const RCP<const Comm<int> > & comm = getSourceMap()->getComm();
00309     const int myImageID = comm->getRank();
00310     const int numImages = comm->getSize();
00311     for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
00312       if (myImageID == imageCtr) {
00313         os << endl;
00314         if (myImageID == 0) { // this is the root node (only output this info once)
00315           os << "Export Data Members:" << endl;
00316         }
00317         os << "Image ID       : " << myImageID << endl;
00318         os << "permuteFromLIDs: {"; av = getPermuteFromLIDs(); std::copy(av.begin(),av.end(),std::ostream_iterator<LocalOrdinal>(os," ")); os << " }" << endl;
00319         os << "permuteToLIDs  : {"; av = getPermuteToLIDs();   std::copy(av.begin(),av.end(),std::ostream_iterator<LocalOrdinal>(os," ")); os << " }" << endl;
00320         os << "remoteLIDs     : {"; av = getRemoteLIDs();      std::copy(av.begin(),av.end(),std::ostream_iterator<LocalOrdinal>(os," ")); os << " }" << endl;
00321         os << "exportLIDs     : {"; av = getExportLIDs();      std::copy(av.begin(),av.end(),std::ostream_iterator<LocalOrdinal>(os," ")); os << " }" << endl;
00322         os << "exportImageIDs : {"; avi = getExportImageIDs();  std::copy(avi.begin(),avi.end(),std::ostream_iterator<int>(os," ")); os << " }" << endl;
00323         os << "numSameIDs     : " << getNumSameIDs() << endl;
00324         os << "numPermuteIDs  : " << getNumPermuteIDs() << endl;
00325         os << "numRemoteIDs   : " << getNumRemoteIDs() << endl;
00326         os << "numExportIDs   : " << getNumExportIDs() << endl;
00327       }
00328       // Do a few global ops to give I/O a chance to complete
00329       comm->barrier();
00330       comm->barrier();
00331       comm->barrier();
00332     }
00333     if (myImageID == 0) {
00334       os << endl << endl << "Source Map:" << endl << std::flush; 
00335     }
00336     comm->barrier();
00337     os << *getSourceMap();
00338     comm->barrier();
00339 
00340     if (myImageID == 0) {
00341       os << endl << endl << "Target Map:" << endl << std::flush; 
00342     }
00343     comm->barrier();
00344     os << *getTargetMap();
00345     comm->barrier();
00346 
00347     // It's also helpful for debugging to print the Distributor
00348     // object.  Epetra_Import::Print() does this (or _should_ do this,
00349     // but doesn't, as of 05 Jan 2012), so we can do a side-by-side
00350     // comparison.
00351     if (myImageID == 0) {
00352       os << endl << endl << "Distributor:" << endl << std::flush;
00353     }
00354     comm->barrier();
00355     getDistributor().describe (*(getFancyOStream (rcpFromRef (os))),
00356              Teuchos::VERB_EXTREME);
00357   }
00358 
00359 
00360   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00361   void 
00362   Export<LocalOrdinal,GlobalOrdinal,Node>::setupSamePermuteExport() 
00363   {
00364     const Map<LocalOrdinal,GlobalOrdinal,Node> & source = *getSourceMap();
00365     const Map<LocalOrdinal,GlobalOrdinal,Node> & target = *getTargetMap();
00366     ArrayView<const GlobalOrdinal> sourceGIDs = source.getNodeElementList();
00367     ArrayView<const GlobalOrdinal> targetGIDs = target.getNodeElementList();
00368 
00369     // Compute numSameIDs_:
00370     //
00371     // Iterate through the source and target GID lists.  If the i-th
00372     // GID of both is the same, increment numSameIDs_ and try the
00373     // next.  As soon as you come to a nonmatching pair, give up.
00374     typename ArrayView<const GlobalOrdinal>::iterator sourceIter = sourceGIDs.begin(),
00375                                                       targetIter = targetGIDs.begin();
00376     while (sourceIter != sourceGIDs.end() && targetIter != targetGIDs.end() && *sourceIter == *targetIter) {
00377       ++ExportData_->numSameIDs_;
00378       ++sourceIter;
00379       ++targetIter;
00380     }
00381     // sourceIter should now point either to the GID of the first
00382     // non-same entry in sourceGIDs, or to the end of sourceGIDs (if
00383     // all the entries were the same).
00384 
00385     // -- compute numPermuteIDs --
00386     // -- fill permuteToLIDs_, permuteFromLIDs_ --
00387     //
00388     // Iterate through the remaining entries in sourceGIDs.  (There
00389     // may not be any, if all pairs of entries in sourceGIDs and
00390     // targetGIDs were the same.)  If target owns that GID, add
00391     // entries to permuteToLIDs_ and permuteFromLIDs_.  Otherwise, add
00392     // entries to exportGIDs_.
00393     //
00394     for (; sourceIter != sourceGIDs.end(); ++sourceIter) {
00395       if (target.isNodeGlobalElement(*sourceIter)) {
00396         // The current process owns this GID, for both the source and
00397   // the target Maps.  Determine the LIDs for this GID on both
00398   // Maps and add them to the permutation lists.
00399         ExportData_->permuteToLIDs_.push_back(  target.getLocalElement(*sourceIter));
00400         ExportData_->permuteFromLIDs_.push_back(source.getLocalElement(*sourceIter));
00401       }
00402       else {
00403         // The current GID is owned by this process in the source Map,
00404         // but is not owned by this process in the target Map.  Store
00405         // such GIDs.
00406         ExportData_->exportGIDs_.push_back(*sourceIter);
00407       }
00408     }
00409 
00410     // Above, we filled exportGIDs_ with all the GIDs which we own in
00411     // the source Map, but not in the target Map.  Now allocate
00412     // exportLIDs_, and fill it with the LIDs (from the source Map)
00413     // corresponding to those GIDs.
00414     if (ExportData_->exportGIDs_.size()) {
00415       ExportData_->exportLIDs_ = arcp<LocalOrdinal>(ExportData_->exportGIDs_.size());
00416     }
00417     {
00418       typename ArrayRCP<LocalOrdinal>::iterator liditer = ExportData_->exportLIDs_.begin();
00419       typename Array<GlobalOrdinal>::iterator   giditer = ExportData_->exportGIDs_.begin();
00420       for (; giditer != ExportData_->exportGIDs_.end(); ++liditer, ++giditer) {
00421         *liditer = source.getLocalElement(*giditer);
00422       }
00423     }
00424 
00425     TPETRA_ABUSE_WARNING( (getNumExportIDs() > 0) && (!source.isDistributed()), std::runtime_error,
00426         "::setupSamePermuteExport(): Source has export LIDs but Source is not distributed globally."
00427         << std::endl << "Importing to a submap of the target map.");
00428 
00429     // -- compute exportImageIDs_ --
00430     //
00431     // For each GID in exportGIDs_, find its corresponding owning
00432     // image (a.k.a. "process," "node") ID in the target Map.  Store
00433     // these image IDs in exportImageIDs_.  These are the image IDs to
00434     // which the Export needs to send data.
00435     //
00436     // We only need to do this if the source Map is distributed;
00437     // otherwise, the Export doesn't have to perform any
00438     // communication.
00439     if (source.isDistributed()) {
00440       ExportData_->exportImageIDs_ = arcp<int>(ExportData_->exportGIDs_.size());
00441       const LookupStatus lookup = target.getRemoteIndexList(ExportData_->exportGIDs_(), ExportData_->exportImageIDs_());
00442       TPETRA_ABUSE_WARNING( lookup == IDNotPresent, std::runtime_error, 
00443         "::setupSamePermuteExport(): The source Map has GIDs not found in the target Map.");
00444 
00445       // Get rid of IDs not in the Target Map
00446       typedef typename ArrayRCP<int>::difference_type size_type;
00447       if (lookup == IDNotPresent) {
00448         const size_type numInvalidExports = std::count_if( ExportData_->exportImageIDs_().begin(), 
00449                                                           ExportData_->exportImageIDs_().end(), 
00450                                                           std::bind1st(std::equal_to<int>(),-1) 
00451                                                         );
00452         // count number of valid and total number of exports
00453         const size_type totalNumExports = ExportData_->exportImageIDs_.size();
00454         if (numInvalidExports == totalNumExports) {
00455           // all exports are invalid; we have no exports; we can delete all exports
00456           ExportData_->exportGIDs_.resize(0);
00457           ExportData_->exportLIDs_ = null;
00458           ExportData_->exportImageIDs_ = null;
00459         }
00460         else {
00461           // some exports are valid; we need to keep the valid exports
00462           // pack and resize
00463           size_type numValidExports = 0;
00464           for (size_type e = 0; e < totalNumExports; ++e) {
00465             if (ExportData_->exportImageIDs_[e] != -1) {
00466               ExportData_->exportGIDs_[numValidExports]     = ExportData_->exportGIDs_[e];
00467               ExportData_->exportLIDs_[numValidExports]     = ExportData_->exportLIDs_[e];
00468               ExportData_->exportImageIDs_[numValidExports] = ExportData_->exportImageIDs_[e];
00469               ++numValidExports;
00470             }
00471           }
00472           ExportData_->exportGIDs_.resize(numValidExports);
00473           ExportData_->exportLIDs_.resize(numValidExports);
00474           ExportData_->exportImageIDs_.resize(numValidExports);
00475         }
00476       }
00477     }
00478   } // setupSamePermuteExport()
00479 
00480 
00481   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00482   void 
00483   Export<LocalOrdinal,GlobalOrdinal,Node>::setupRemote() 
00484   {
00485     const Map<LocalOrdinal,GlobalOrdinal,Node>& target = *getTargetMap();
00486 
00487     // Make sure the export IDs are ordered by image.  Sort
00488     // exportImageIDs_ in ascending order, and apply the same
00489     // permutation to exportGIDs_ and exportLIDs_.
00490     sort3 (ExportData_->exportImageIDs_.begin(), 
00491      ExportData_->exportImageIDs_.end(), 
00492      ExportData_->exportGIDs_.begin(), 
00493      ExportData_->exportLIDs_.begin());
00494 
00495     // Construct the list of entries that calling image needs to send
00496     // as a result of everyone asking for what it needs to receive.
00497     //
00498     // mfh 05 Jan 2012: I understand the above comment as follows:
00499     // Construct the communication plan from the list of image IDs to
00500     // which we need to send.
00501     size_t numRemoteIDs;
00502     numRemoteIDs = ExportData_->distributor_.createFromSends (ExportData_->exportImageIDs_());
00503 
00504     // Use the communication plan with ExportGIDs to find out who is
00505     // sending to us and get the proper ordering of GIDs for incoming
00506     // remote entries (these will be converted to LIDs when done).
00507     Array<GlobalOrdinal> remoteGIDs(numRemoteIDs);
00508     ExportData_->distributor_.doPostsAndWaits (ExportData_->exportGIDs_().getConst(), 1, remoteGIDs());
00509 
00510     // Remote IDs come in as GIDs; convert to LIDs.  LIDs tell us
00511     // where to store the incoming remote data.
00512     ExportData_->remoteLIDs_.resize(numRemoteIDs);
00513     {
00514       typename Array<GlobalOrdinal>::const_iterator i = remoteGIDs.begin();
00515       typename Array<LocalOrdinal>::iterator        j = ExportData_->remoteLIDs_.begin();
00516       while (i != remoteGIDs.end()) {
00517         *j++ = target.getLocalElement(*i++);
00518       }
00519     }
00520   }
00521 
00531   template <class LO, class GO, class Node> 
00532   RCP< const Export<LO,GO,Node> >
00533   createExport( const RCP<const Map<LO,GO,Node> > & src, 
00534                 const RCP<const Map<LO,GO,Node> > & tgt )
00535   {
00536     if (src == tgt) return null;
00537 #ifdef HAVE_TPETRA_DEBUG
00538     TEUCHOS_TEST_FOR_EXCEPTION(src == null || tgt == null, std::runtime_error,
00539         "Tpetra::createExport(): neither source nor target map may be null:\nsource: " << src << "\ntarget: " << tgt << "\n");
00540 #endif
00541     return rcp(new Export<LO,GO,Node>(src,tgt));
00542   }
00543 
00544 } // namespace Tpetra
00545 
00546 #endif // TPETRA_EXPORT_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines