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