Tpetra Matrix/Vector Services Version of the Day
Tpetra_Import.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_IMPORT_HPP
00030 #define TPETRA_IMPORT_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 "Tpetra_Distributor.hpp"
00039 #include <iterator>
00040 
00041 namespace Tpetra {
00042 
00044 
00055   template <class LocalOrdinal, class GlobalOrdinal = LocalOrdinal, class Node = Kokkos::DefaultNode::DefaultNodeType>
00056   class Import: public Teuchos::Describable {
00057 
00058   public:
00059 
00061 
00062 
00064     Import(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & source, 
00065            const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & target);
00066 
00068     Import(const Import<LocalOrdinal,GlobalOrdinal,Node> & import);
00069 
00071     ~Import();
00072 
00074 
00076 
00077 
00079     inline size_t getNumSameIDs() const;
00080 
00082     inline size_t getNumPermuteIDs() const;
00083 
00085     inline ArrayView<const LocalOrdinal> getPermuteFromLIDs() const;
00086 
00088     inline ArrayView<const LocalOrdinal> getPermuteToLIDs() const;
00089 
00091     inline size_t getNumRemoteIDs() const;
00092 
00094     inline ArrayView<const LocalOrdinal> getRemoteLIDs() const;
00095 
00097     inline size_t getNumExportIDs() const;
00098 
00100     inline ArrayView<const LocalOrdinal> getExportLIDs() const;
00101 
00103     inline ArrayView<const int> getExportImageIDs() const;
00104 
00106     inline const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getSourceMap() const;
00107 
00109     inline const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getTargetMap() const;
00110 
00111     inline Distributor & getDistributor() const;
00112 
00114     Import<LocalOrdinal,GlobalOrdinal,Node>& operator = (const Import<LocalOrdinal,GlobalOrdinal,Node> & Source);
00115 
00117 
00119 
00120 
00122     virtual void print(std::ostream& os) const;
00123 
00125 
00126   private:
00127 
00128     RCP<ImportExportData<LocalOrdinal,GlobalOrdinal,Node> > ImportData_;
00129     RCP<Array<GlobalOrdinal> > remoteGIDs_;
00130 
00131     // subfunctions used by constructor
00132     //==============================================================================
00133     // sets up numSameIDs_, numPermuteIDs_, and numRemoteIDs_
00134     // these variables are already initialized to 0 by the ImportExportData ctr.
00135     // also sets up permuteToLIDs_, permuteFromLIDs_, and remoteLIDs_
00136     void setupSamePermuteRemote();
00137     void setupExport();
00138   };
00139 
00140   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00141   Import<LocalOrdinal,GlobalOrdinal,Node>::Import(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & source, 
00142                                                   const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & target) {
00143     ImportData_ = rcp(new ImportExportData<LocalOrdinal,GlobalOrdinal,Node>(source, target));
00144     // call subfunctions
00145     setupSamePermuteRemote();
00146     if (source->isDistributed()) {
00147       setupExport();
00148     }
00149     // don't need remoteGIDs_ anymore
00150     remoteGIDs_ = null;
00151   }
00152 
00153   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00154   Import<LocalOrdinal,GlobalOrdinal,Node>::Import(const Import<LocalOrdinal,GlobalOrdinal,Node> & import)
00155   : ImportData_(import.ImportData_)
00156   {}
00157 
00158   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00159   Import<LocalOrdinal,GlobalOrdinal,Node>::~Import() 
00160   {}
00161 
00162   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00163   size_t Import<LocalOrdinal,GlobalOrdinal,Node>::getNumSameIDs() const {
00164     return ImportData_->numSameIDs_;
00165   }
00166 
00167   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00168   size_t Import<LocalOrdinal,GlobalOrdinal,Node>::getNumPermuteIDs() const {
00169     return ImportData_->permuteFromLIDs_.size();
00170   }
00171 
00172   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00173   ArrayView<const LocalOrdinal> 
00174   Import<LocalOrdinal,GlobalOrdinal,Node>::getPermuteFromLIDs() const {
00175     return ImportData_->permuteFromLIDs_();
00176   }
00177 
00178   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00179   ArrayView<const LocalOrdinal>
00180   Import<LocalOrdinal,GlobalOrdinal,Node>::getPermuteToLIDs() const {
00181     return ImportData_->permuteToLIDs_();
00182   }
00183 
00184   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00185   size_t Import<LocalOrdinal,GlobalOrdinal,Node>::getNumRemoteIDs() const {
00186     return ImportData_->remoteLIDs_.size();
00187   }
00188 
00189   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00190   ArrayView<const LocalOrdinal> 
00191   Import<LocalOrdinal,GlobalOrdinal,Node>::getRemoteLIDs() const {
00192     return ImportData_->remoteLIDs_();
00193   }
00194 
00195   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00196   size_t Import<LocalOrdinal,GlobalOrdinal,Node>::getNumExportIDs() const {
00197     return ImportData_->exportLIDs_.size();
00198   }
00199 
00200   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00201   ArrayView<const LocalOrdinal> 
00202   Import<LocalOrdinal,GlobalOrdinal,Node>::getExportLIDs() const {
00203     return ImportData_->exportLIDs_();
00204   }
00205 
00206   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00207   ArrayView<const int> 
00208   Import<LocalOrdinal,GlobalOrdinal,Node>::getExportImageIDs() const {
00209     return ImportData_->exportImageIDs_();
00210   }
00211 
00212   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00213   const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & 
00214   Import<LocalOrdinal,GlobalOrdinal,Node>::getSourceMap() const {
00215     return ImportData_->source_;
00216   }
00217 
00218   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00219   const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & 
00220   Import<LocalOrdinal,GlobalOrdinal,Node>::getTargetMap() const {
00221     return ImportData_->target_;
00222   }
00223 
00224   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00225   Distributor & 
00226   Import<LocalOrdinal,GlobalOrdinal,Node>::getDistributor() const {
00227     return ImportData_->distributor_;
00228   }
00229 
00230   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00231   Import<LocalOrdinal,GlobalOrdinal,Node>& 
00232   Import<LocalOrdinal,GlobalOrdinal,Node>::operator=(const Import<LocalOrdinal,GlobalOrdinal,Node> & source) {
00233     ImportData_ = source.ImportData_;
00234     return *this;
00235   }
00236 
00237   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00238   void Import<LocalOrdinal,GlobalOrdinal,Node>::print(std::ostream& os) const {
00239     using std::endl;
00240     ArrayView<const LocalOrdinal> av;
00241     ArrayView<const int> avi;
00242     const RCP<const Comm<int> > & comm = getSourceMap()->getComm();
00243     const int myImageID = comm->getRank();
00244     const int numImages = comm->getSize();
00245     for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
00246       if (myImageID == imageCtr) 
00247       {
00248         os << endl;
00249         if(myImageID == 0) { // this is the root node (only output this info once)
00250           os << "Import Data Members:" << endl;
00251         }
00252         os << "Image ID       : " << myImageID << endl;
00253         os << "permuteFromLIDs: {"; av = getPermuteFromLIDs(); std::copy(av.begin(),av.end(),std::ostream_iterator<LocalOrdinal>(os," ")); os << " }" << endl;
00254         os << "permuteToLIDs  : {"; av = getPermuteToLIDs();   std::copy(av.begin(),av.end(),std::ostream_iterator<LocalOrdinal>(os," ")); os << " }" << endl;
00255         os << "remoteLIDs     : {"; av = getRemoteLIDs();      std::copy(av.begin(),av.end(),std::ostream_iterator<LocalOrdinal>(os," ")); os << " }" << endl;
00256         os << "exportLIDs     : {"; av = getExportLIDs();      std::copy(av.begin(),av.end(),std::ostream_iterator<LocalOrdinal>(os," ")); os << " }" << endl;
00257         os << "exportImageIDs : {"; avi = getExportImageIDs();  std::copy(avi.begin(),avi.end(),std::ostream_iterator<int>(os," ")); os << " }" << endl;
00258         os << "numSameIDs     : " << getNumSameIDs() << endl;
00259         os << "numPermuteIDs  : " << getNumPermuteIDs() << endl;
00260         os << "numRemoteIDs   : " << getNumRemoteIDs() << endl;
00261         os << "numExportIDs   : " << getNumExportIDs() << endl;
00262       }
00263       // Do a few global ops to give I/O a chance to complete
00264       comm->barrier();
00265       comm->barrier();
00266       comm->barrier();
00267     }
00268     if (myImageID == 0) {
00269       os << "\nSource Map: " << endl; 
00270     }
00271     os << *getSourceMap();
00272     if (myImageID == 0) {
00273       os << "\nTarget Map: " << endl; 
00274     }
00275     os << *getTargetMap();
00276   }
00277 
00278 
00279   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00280   void Import<LocalOrdinal,GlobalOrdinal,Node>::setupSamePermuteRemote() {
00281     const Map<LocalOrdinal,GlobalOrdinal,Node> & source = *getSourceMap();
00282     const Map<LocalOrdinal,GlobalOrdinal,Node> & target = *getTargetMap();
00283     ArrayView<const GlobalOrdinal> sourceGIDs = source.getNodeElementList();
00284     ArrayView<const GlobalOrdinal> targetGIDs = target.getNodeElementList();
00285 
00286     // -- compute numSameIDs_ ---
00287     // go through GID lists of source and target. if the ith GID on both is the same, 
00288     // increment numSameIDs_ and try the next. as soon as you come to a pair that don't
00289     // match, give up.
00290     typename ArrayView<const GlobalOrdinal>::iterator sourceIter = sourceGIDs.begin(),
00291                                                       targetIter = targetGIDs.begin();
00292     while( sourceIter != sourceGIDs.end() && targetIter != targetGIDs.end() && *sourceIter == *targetIter )
00293     {
00294       ++ImportData_->numSameIDs_;
00295       ++sourceIter;
00296       ++targetIter;
00297     }
00298     // targetIter should now point to the GID of the first non-same entry or the end of targetGIDs
00299 
00300     // -- compute numPermuteIDs and numRemoteIDs --
00301     // -- fill permuteToLIDs_, permuteFromLIDs_, remoteGIDs_, and remoteLIDs_ --
00302     // go through remaining entries in targetGIDs. if source owns that GID, 
00303     // increment numPermuteIDs_, and add entries to permuteToLIDs_ and permuteFromLIDs_.
00304     // otherwise increment numRemoteIDs_ and add entries to remoteLIDs_ and remoteGIDs_.
00305     remoteGIDs_ = rcp( new Array<GlobalOrdinal>() );
00306     for (; targetIter != targetGIDs.end(); ++targetIter) {
00307       if (source.isNodeGlobalElement(*targetIter)) {
00308         // both source and target list this GID (*targetIter)
00309         // determine the LIDs for this GID on both Maps and add them to the permutation lists
00310         ImportData_->permuteToLIDs_.push_back(target.getLocalElement(*targetIter));
00311         ImportData_->permuteFromLIDs_.push_back(source.getLocalElement(*targetIter));
00312       }
00313       else {
00314         // this GID is on another processor; store it, along with its destination LID on this processor
00315         remoteGIDs_->push_back(*targetIter);
00316         ImportData_->remoteLIDs_.push_back(target.getLocalElement(*targetIter));
00317       }
00318     }
00319 
00320     TPETRA_ABUSE_WARNING((getNumRemoteIDs() > 0) && !source.isDistributed(),std::runtime_error,
00321         "::setupSamePermuteRemote(): Target has remote LIDs but Source is not distributed globally."
00322         << std::endl << "Importing to a submap of the target map.");
00323   }
00324 
00325 
00326   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00327   void Import<LocalOrdinal,GlobalOrdinal,Node>::setupExport() {
00328     typedef typename Array<int>::difference_type size_type;
00329     const Map<LocalOrdinal,GlobalOrdinal,Node> & source = *getSourceMap();
00330 
00331     // create remoteImageID list: for each entry remoteGIDs[i],
00332     // remoteImageIDs[i] will contain the ImageID of the image that owns that GID.
00333     // check for GIDs that exist in target but not in source: we see this if getRemoteIndexList returns true
00334     ArrayView<GlobalOrdinal> remoteGIDs = (*remoteGIDs_)();
00335     Array<int> remoteImageIDs(remoteGIDs.size());
00336     const LookupStatus lookup = source.getRemoteIndexList(remoteGIDs, remoteImageIDs());
00337     TPETRA_ABUSE_WARNING( lookup == IDNotPresent, std::runtime_error, "::setupExport(): Target has GIDs not found in Source." );
00338 
00339     // get rid of ids that don't exist in the source map
00340     if ( lookup == IDNotPresent ) {
00341       const size_type numInvalidRemote = std::count_if( remoteImageIDs.begin(), remoteImageIDs.end(), std::bind1st(std::equal_to<int>(),-1) );
00342       // if all of them are invalid, we can delete the whole array
00343       const size_type totalNumRemote = getNumRemoteIDs();
00344       if (numInvalidRemote == totalNumRemote) {
00345         // all remotes are invalid; we have no remotes; we can delete the remotes
00346         remoteImageIDs.clear();
00347         (*remoteGIDs_).clear();
00348         ImportData_->remoteLIDs_.clear();
00349       }
00350       else {
00351         // some remotes are valid; we need to keep the valid ones
00352         // pack and resize
00353         size_type numValidRemote = 0;
00354         for (size_type r = 0; r < totalNumRemote; ++r) {
00355           if (remoteImageIDs[r] != -1) {
00356             remoteImageIDs[numValidRemote] = remoteImageIDs[r];
00357             (*remoteGIDs_)[numValidRemote] = (*remoteGIDs_)[r];
00358             ImportData_->remoteLIDs_[numValidRemote] = ImportData_->remoteLIDs_[r];
00359             ++numValidRemote;
00360           }
00361         }
00362         TEST_FOR_EXCEPTION( numValidRemote != totalNumRemote - numInvalidRemote, std::logic_error,
00363             typeName(*this) << "::setupExport(): internal logic error. Please contact Tpetra team.")
00364         remoteImageIDs.resize(numValidRemote);
00365         (*remoteGIDs_).resize(numValidRemote);
00366         ImportData_->remoteLIDs_.resize(numValidRemote);
00367       }
00368       remoteGIDs = (*remoteGIDs_)();
00369     }
00370 
00371     // sort remoteImageIDs in ascending order
00372     // apply same permutation to remoteGIDs_
00373     sort2(remoteImageIDs.begin(), remoteImageIDs.end(), remoteGIDs.begin());
00374 
00375     // call Distributor.createFromRecvs()
00376     // takes in remoteGIDs and remoteImageIDs_
00377     // returns exportLIDs_, exportImageIDs_ 
00378     ArrayRCP<GlobalOrdinal> exportGIDs;
00379     ImportData_->distributor_.createFromRecvs(remoteGIDs().getConst(), remoteImageIDs, exportGIDs, ImportData_->exportImageIDs_);
00380     // -- exportGIDs and exportImageIDs_ allocated by createFromRecvs (the former contains GIDs, we will convert to LIDs below) --
00381 
00382     // convert exportGIDs from GIDs to LIDs
00383     if (exportGIDs != null) {
00384       ImportData_->exportLIDs_ = arcp<LocalOrdinal>(exportGIDs.size());
00385     }
00386     typename ArrayRCP<LocalOrdinal>::iterator dst = ImportData_->exportLIDs_.begin();
00387     typename ArrayRCP<GlobalOrdinal>::const_iterator src = exportGIDs.begin();
00388     while (src != exportGIDs.end())
00389     {
00390       (*dst++) = source.getLocalElement(*src++);
00391     }
00392   }
00393 
00403   template <class LO, class GO, class Node> 
00404   RCP< const Import<LO,GO,Node> >
00405   createImport( const RCP<const Map<LO,GO,Node> > & src, 
00406                 const RCP<const Map<LO,GO,Node> > & tgt )
00407   {
00408     if (src == tgt) return null;
00409 #ifdef HAVE_TPETRA_DEBUG
00410     TEST_FOR_EXCEPTION(src == null || tgt == null, std::runtime_error,
00411         "Tpetra::createImport(): neither source nor target map may be null:\nsource: " << src << "\ntarget: " << tgt << "\n");
00412 #endif
00413     return rcp(new Import<LO,GO,Node>(src,tgt));
00414   }
00415 
00418   template <class LO, class GO, class Node> 
00419   RCP< const Import<LO,GO,Node> >
00420   TPETRA_DEPRECATED makeImport( const RCP<const Map<LO,GO,Node> > & src, 
00421                                 const RCP<const Map<LO,GO,Node> > & tgt )
00422   {
00423     return createImport<LO,GO,Node>(src,tgt);
00424   }
00425 
00426 } // namespace Tpetra
00427 
00428 #endif // TPETRA_IMPORT_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines