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 terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef TPETRA_EXPORT_HPP
00030 #define TPETRA_EXPORT_HPP
00031 
00032 #include <Kokkos_DefaultNode.hpp>
00033 #include <Teuchos_Describable.hpp>
00034 #include <Teuchos_as.hpp>
00035 #include "Tpetra_Map.hpp"
00036 #include "Tpetra_Util.hpp"
00037 #include "Tpetra_ImportExportData.hpp"
00038 #include <iterator>
00039 
00040 namespace Tpetra {
00041 
00043 
00054   template <class LocalOrdinal, class GlobalOrdinal = LocalOrdinal, class Node = Kokkos::DefaultNode::DefaultNodeType>
00055   class Export: public Teuchos::Describable {
00056 
00057   public:
00058 
00060 
00061 
00063     Export(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & source, 
00064            const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & target);
00065 
00067     Export(const Export<LocalOrdinal,GlobalOrdinal,Node> & import);
00068 
00070     ~Export();
00071 
00073 
00075 
00076 
00078     inline size_t getNumSameIDs() const;
00079 
00081     inline size_t getNumPermuteIDs() const;
00082 
00084     inline ArrayView<const LocalOrdinal> getPermuteFromLIDs() const;
00085 
00087     inline ArrayView<const LocalOrdinal> getPermuteToLIDs() const;
00088 
00090     inline size_t getNumRemoteIDs() const;
00091 
00093     inline ArrayView<const LocalOrdinal> getRemoteLIDs() const;
00094 
00096     inline size_t getNumExportIDs() const;
00097 
00099     inline ArrayView<const LocalOrdinal> getExportLIDs() const;
00100 
00102     inline ArrayView<const int> getExportImageIDs() const;
00103 
00105     inline const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getSourceMap() const;
00106 
00108     inline const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getTargetMap() const;
00109 
00110     inline Distributor & getDistributor() const;
00111 
00113     Export<LocalOrdinal,GlobalOrdinal,Node>& operator = (const Export<LocalOrdinal,GlobalOrdinal,Node> & Source);
00114 
00116 
00118 
00119 
00121     virtual void print(std::ostream& os) const;
00122 
00124 
00125   private:
00126 
00127     RCP<ImportExportData<LocalOrdinal,GlobalOrdinal,Node> > ExportData_;
00128 
00129     // subfunctions used by constructor
00130     //==============================================================================
00131     // sets up numSameIDs_, numPermuteIDs_, and the export IDs
00132     // these variables are already initialized to 0 by the ImportExportData ctr.
00133     // also sets up permuteToLIDs_, permuteFromLIDs_, exportGIDs_, and exportLIDs_
00134     void setupSamePermuteExport();
00135     void setupRemote();
00136   };
00137 
00138   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00139   Export<LocalOrdinal,GlobalOrdinal,Node>::Export(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & source,   
00140                                                   const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & target) {
00141     ExportData_ = rcp(new ImportExportData<LocalOrdinal,GlobalOrdinal,Node>(source, target));
00142     // call subfunctions
00143     setupSamePermuteExport();
00144     if (source->isDistributed()) {
00145       setupRemote();
00146     }
00147   }
00148 
00149   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00150   Export<LocalOrdinal,GlobalOrdinal,Node>::Export(const Export<LocalOrdinal,GlobalOrdinal,Node> & rhs)
00151   : ExportData_(rhs.ExportData_)
00152   {}
00153 
00154   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00155   Export<LocalOrdinal,GlobalOrdinal,Node>::~Export() 
00156   {}
00157 
00158   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00159   size_t Export<LocalOrdinal,GlobalOrdinal,Node>::getNumSameIDs() const {
00160     return ExportData_->numSameIDs_;
00161   }
00162 
00163   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00164   size_t Export<LocalOrdinal,GlobalOrdinal,Node>::getNumPermuteIDs() const {
00165     return ExportData_->permuteFromLIDs_.size();
00166   }
00167 
00168   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00169   ArrayView<const LocalOrdinal> 
00170   Export<LocalOrdinal,GlobalOrdinal,Node>::getPermuteFromLIDs() const {
00171     return ExportData_->permuteFromLIDs_();
00172   }
00173 
00174   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00175   ArrayView<const LocalOrdinal>
00176   Export<LocalOrdinal,GlobalOrdinal,Node>::getPermuteToLIDs() const {
00177     return ExportData_->permuteToLIDs_();
00178   }
00179 
00180   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00181   size_t Export<LocalOrdinal,GlobalOrdinal,Node>::getNumRemoteIDs() const {
00182     return ExportData_->remoteLIDs_.size();
00183   }
00184 
00185   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00186   ArrayView<const LocalOrdinal> 
00187   Export<LocalOrdinal,GlobalOrdinal,Node>::getRemoteLIDs() const {
00188     return ExportData_->remoteLIDs_();
00189   }
00190 
00191   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00192   size_t Export<LocalOrdinal,GlobalOrdinal,Node>::getNumExportIDs() const {
00193     return ExportData_->exportLIDs_.size();
00194   }
00195 
00196   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00197   ArrayView<const LocalOrdinal> 
00198   Export<LocalOrdinal,GlobalOrdinal,Node>::getExportLIDs() const {
00199     return ExportData_->exportLIDs_();
00200   }
00201 
00202   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00203   ArrayView<const int> 
00204   Export<LocalOrdinal,GlobalOrdinal,Node>::getExportImageIDs() const {
00205     return ExportData_->exportImageIDs_();
00206   }
00207 
00208   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00209   const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & 
00210   Export<LocalOrdinal,GlobalOrdinal,Node>::getSourceMap() const {
00211     return ExportData_->source_;
00212   }
00213 
00214   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00215   const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & 
00216   Export<LocalOrdinal,GlobalOrdinal,Node>::getTargetMap() const {
00217     return ExportData_->target_;
00218   }
00219 
00220   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00221   Distributor & 
00222   Export<LocalOrdinal,GlobalOrdinal,Node>::getDistributor() const { 
00223     return ExportData_->distributor_;
00224   }
00225 
00226   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00227   Export<LocalOrdinal,GlobalOrdinal,Node>& 
00228   Export<LocalOrdinal,GlobalOrdinal,Node>::operator=(const Export<LocalOrdinal,GlobalOrdinal,Node> & source) {
00229     ExportData_ = source.ExportData_;
00230     return *this;
00231   }
00232 
00233   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00234   void Export<LocalOrdinal,GlobalOrdinal,Node>::print(std::ostream& os) const {
00235     using std::endl;
00236     ArrayView<const LocalOrdinal> av;
00237     ArrayView<const int> avi;
00238     const RCP<const Comm<int> > & comm = getSourceMap()->getComm();
00239     const int myImageID = comm->getRank();
00240     const int numImages = comm->getSize();
00241     for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
00242       if (myImageID == imageCtr) {
00243         os << endl;
00244         if(myImageID == 0) { // this is the root node (only output this info once)
00245           os << "Export Data Members:" << endl;
00246         }
00247         os << "Image ID       : " << myImageID << endl;
00248         os << "permuteFromLIDs: {"; av = getPermuteFromLIDs(); std::copy(av.begin(),av.end(),std::ostream_iterator<LocalOrdinal>(os," ")); os << " }" << endl;
00249         os << "permuteToLIDs  : {"; av = getPermuteToLIDs();   std::copy(av.begin(),av.end(),std::ostream_iterator<LocalOrdinal>(os," ")); os << " }" << endl;
00250         os << "remoteLIDs     : {"; av = getRemoteLIDs();      std::copy(av.begin(),av.end(),std::ostream_iterator<LocalOrdinal>(os," ")); os << " }" << endl;
00251         os << "exportLIDs     : {"; av = getExportLIDs();      std::copy(av.begin(),av.end(),std::ostream_iterator<LocalOrdinal>(os," ")); os << " }" << endl;
00252         os << "exportImageIDs : {"; avi = getExportImageIDs();  std::copy(avi.begin(),avi.end(),std::ostream_iterator<int>(os," ")); os << " }" << endl;
00253         os << "numSameIDs     : " << getNumSameIDs() << endl;
00254         os << "numPermuteIDs  : " << getNumPermuteIDs() << endl;
00255         os << "numRemoteIDs   : " << getNumRemoteIDs() << endl;
00256         os << "numExportIDs   : " << getNumExportIDs() << endl;
00257       }
00258       // Do a few global ops to give I/O a chance to complete
00259       comm->barrier();
00260       comm->barrier();
00261       comm->barrier();
00262     }
00263     if (myImageID == 0) {
00264       os << "\nSource Map: " << endl; 
00265     }
00266     os << *getSourceMap();
00267     if (myImageID == 0) {
00268       os << "\nTarget Map: " << endl; 
00269     }
00270     os << *getTargetMap();
00271   }
00272 
00273 
00274   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00275   void Export<LocalOrdinal,GlobalOrdinal,Node>::setupSamePermuteExport() {
00276     const Map<LocalOrdinal,GlobalOrdinal,Node> & source = *getSourceMap();
00277     const Map<LocalOrdinal,GlobalOrdinal,Node> & target = *getTargetMap();
00278     ArrayView<const GlobalOrdinal> sourceGIDs = source.getNodeElementList();
00279     ArrayView<const GlobalOrdinal> targetGIDs = target.getNodeElementList();
00280 
00281     // -- compute numSameIDs_ ---
00282     // go through GID lists of source and target. if the ith GID on both is the same, 
00283     // increment numSameIDs_ and try the next. as soon as you come to a pair that don't
00284     // match, give up.
00285     typename ArrayView<const GlobalOrdinal>::iterator sourceIter = sourceGIDs.begin(),
00286                                                       targetIter = targetGIDs.begin();
00287     while( sourceIter != sourceGIDs.end() && targetIter != targetGIDs.end() && *sourceIter == *targetIter )
00288     {
00289       ++ExportData_->numSameIDs_;
00290       ++sourceIter;
00291       ++targetIter;
00292     }
00293     // sourceIter should now point to the GID of the first non-same entry or the end of targetGIDs
00294 
00295     // -- compute numPermuteIDs --
00296     // -- fill permuteToLIDs_, permuteFromLIDs_ --
00297     // go through remaining entries in sourceGIDs. if target owns that GID, add entries to permuteToLIDs_ and permuteFromLIDs_
00298     // otherwise add entries to exportGIDs_
00299     //
00300     for (; sourceIter != sourceGIDs.end(); ++sourceIter) {
00301       if (target.isNodeGlobalElement(*sourceIter)) {
00302         // both source and target list this GID (*targetIter)
00303         // determine the LIDs for this GID on both Maps and add them to the permutation lists
00304         ExportData_->permuteToLIDs_.push_back(  target.getLocalElement(*sourceIter));
00305         ExportData_->permuteFromLIDs_.push_back(source.getLocalElement(*sourceIter));
00306       }
00307       else {
00308         // this GID is on another processor; store it
00309         ExportData_->exportGIDs_.push_back(*sourceIter);
00310       }
00311     }
00312 
00313     // allocate and assign exportLIDs_
00314     if (ExportData_->exportGIDs_.size()) {
00315       ExportData_->exportLIDs_     = arcp<LocalOrdinal>(ExportData_->exportGIDs_.size());
00316     }
00317     {
00318       typename ArrayRCP<LocalOrdinal>::iterator liditer = ExportData_->exportLIDs_.begin();
00319       typename Array<GlobalOrdinal>::iterator   giditer = ExportData_->exportGIDs_.begin();
00320       for (; giditer != ExportData_->exportGIDs_.end(); ++liditer, ++giditer) {
00321         *liditer = source.getLocalElement(*giditer);
00322       }
00323     }
00324 
00325     TPETRA_ABUSE_WARNING( (getNumExportIDs() > 0) && (!source.isDistributed()), std::runtime_error,
00326         "::setupSamePermuteExport(): Source has export LIDs but Source is not distributed globally."
00327         << std::endl << "Importing to a submap of the target map.");
00328 
00329     // -- compute exportImageIDs_ --
00330     // get list of images that own the GIDs in exportGIDs_ (in the target Map)
00331     if (source.isDistributed()) {
00332       ExportData_->exportImageIDs_ = arcp<int>(ExportData_->exportGIDs_.size());
00333       const LookupStatus lookup = target.getRemoteIndexList(ExportData_->exportGIDs_(), ExportData_->exportImageIDs_());
00334       TPETRA_ABUSE_WARNING( lookup == IDNotPresent, std::runtime_error, 
00335           "::setupSamePermuteExport(): Source has GIDs not found in Target.");
00336 
00337       // Get rid of IDs not in the Target Map
00338       typedef typename ArrayRCP<int>::difference_type size_type;
00339       if (lookup == IDNotPresent) {
00340         const size_type numInvalidExports = std::count_if( ExportData_->exportImageIDs_().begin(), 
00341                                                           ExportData_->exportImageIDs_().end(), 
00342                                                           std::bind1st(std::equal_to<int>(),-1) 
00343                                                         );
00344         // count number of valid and total number of exports
00345         const size_type totalNumExports = ExportData_->exportImageIDs_.size();
00346         if (numInvalidExports == totalNumExports) {
00347           // all exports are invalid; we have no exports; we can delete all exports
00348           ExportData_->exportGIDs_.resize(0);
00349           ExportData_->exportLIDs_ = null;
00350           ExportData_->exportImageIDs_ = null;
00351         }
00352         else {
00353           // some exports are valid; we need to keep the valid exports
00354           // pack and resize
00355           size_type numValidExports = 0;
00356           for (size_type e = 0; e < totalNumExports; ++e) {
00357             if (ExportData_->exportImageIDs_[e] != -1) {
00358               ExportData_->exportGIDs_[numValidExports]     = ExportData_->exportGIDs_[e];
00359               ExportData_->exportLIDs_[numValidExports]     = ExportData_->exportLIDs_[e];
00360               ExportData_->exportImageIDs_[numValidExports] = ExportData_->exportImageIDs_[e];
00361               ++numValidExports;
00362             }
00363           }
00364           ExportData_->exportGIDs_.resize(numValidExports);
00365           ExportData_->exportLIDs_.resize(numValidExports);
00366           ExportData_->exportImageIDs_.resize(numValidExports);
00367         }
00368       }
00369     }
00370   } // setupSamePermuteExport()
00371 
00372 
00373   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00374   void Export<LocalOrdinal,GlobalOrdinal,Node>::setupRemote() {
00375     const Map<LocalOrdinal,GlobalOrdinal,Node> & target = *getTargetMap();
00376 
00377     // make sure export IDs are ordered by image
00378     // sort exportImageIDs_ in ascending order,
00379     // and apply the same permutation to exportGIDs_ and exportLIDs_.
00380     sort3(ExportData_->exportImageIDs_.begin(), ExportData_->exportImageIDs_.end(), ExportData_->exportGIDs_.begin(), ExportData_->exportLIDs_.begin());
00381 
00382     // Construct list of entries that calling image needs to send as a result
00383     // of everyone asking for what it needs to receive.
00384     size_t numRemoteIDs;
00385     numRemoteIDs = ExportData_->distributor_.createFromSends(ExportData_->exportImageIDs_());
00386 
00387     // Use comm plan with ExportGIDs to find out who is sending to us and
00388     // get proper ordering of GIDs for remote entries 
00389     // (these will be converted to LIDs when done).
00390     Array<GlobalOrdinal> remoteGIDs(numRemoteIDs);
00391     ExportData_->distributor_.doPostsAndWaits(ExportData_->exportGIDs_().getConst(),1,remoteGIDs());
00392 
00393     // Remote IDs come in as GIDs, convert to LIDs
00394     ExportData_->remoteLIDs_.resize(numRemoteIDs);
00395     {
00396       typename Array<GlobalOrdinal>::const_iterator i = remoteGIDs.begin();
00397       typename Array<LocalOrdinal>::iterator       j = ExportData_->remoteLIDs_.begin();
00398       while (i != remoteGIDs.end()) 
00399       {
00400         *j++ = target.getLocalElement(*i++);
00401       }
00402     }
00403   }
00404 
00414   template <class LO, class GO, class Node> 
00415   RCP< const Export<LO,GO,Node> >
00416   createExport( const RCP<const Map<LO,GO,Node> > & src, 
00417                 const RCP<const Map<LO,GO,Node> > & tgt )
00418   {
00419     if (src == tgt) return null;
00420 #ifdef HAVE_TPETRA_DEBUG
00421     TEST_FOR_EXCEPTION(src == null || tgt == null, std::runtime_error,
00422         "Tpetra::createExport(): neither source nor target map may be null:\nsource: " << src << "\ntarget: " << tgt << "\n");
00423 #endif
00424     return rcp(new Export<LO,GO,Node>(src,tgt));
00425   }
00426 
00427 } // namespace Tpetra
00428 
00429 #endif // TPETRA_EXPORT_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines