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 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_IMPORT_HPP
00043 #define TPETRA_IMPORT_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 "Tpetra_Distributor.hpp"
00052 #include <iterator>
00053 
00054 namespace Tpetra {
00055 
00087   template <class LocalOrdinal, class GlobalOrdinal = LocalOrdinal, class Node = Kokkos::DefaultNode::DefaultNodeType>
00088   class Import: public Teuchos::Describable {
00089 
00090   public:
00091 
00093 
00094 
00096     Import(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & source, 
00097            const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & target);
00098 
00100     Import(const Import<LocalOrdinal,GlobalOrdinal,Node> & import);
00101 
00103     ~Import();
00104 
00106 
00108 
00109 
00111     inline size_t getNumSameIDs() const;
00112 
00114     inline size_t getNumPermuteIDs() const;
00115 
00117     inline ArrayView<const LocalOrdinal> getPermuteFromLIDs() const;
00118 
00120     inline ArrayView<const LocalOrdinal> getPermuteToLIDs() const;
00121 
00123     inline size_t getNumRemoteIDs() const;
00124 
00126     inline ArrayView<const LocalOrdinal> getRemoteLIDs() const;
00127 
00129     inline size_t getNumExportIDs() const;
00130 
00132     inline ArrayView<const LocalOrdinal> getExportLIDs() const;
00133 
00135     inline ArrayView<const int> getExportImageIDs() const;
00136 
00138     inline const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getSourceMap() const;
00139 
00141     inline const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getTargetMap() const;
00142 
00143     inline Distributor & getDistributor() const;
00144 
00146     Import<LocalOrdinal,GlobalOrdinal,Node>& operator = (const Import<LocalOrdinal,GlobalOrdinal,Node> & Source);
00147 
00149 
00151 
00152 
00154     virtual void print(std::ostream& os) const;
00155 
00157 
00158   private:
00159 
00160     RCP<ImportExportData<LocalOrdinal,GlobalOrdinal,Node> > ImportData_;
00161     RCP<Array<GlobalOrdinal> > remoteGIDs_;
00162 
00164 
00165 
00166     //==============================================================================
00167     // sets up numSameIDs_, numPermuteIDs_, and numRemoteIDs_
00168     // these variables are already initialized to 0 by the ImportExportData ctr.
00169     // also sets up permuteToLIDs_, permuteFromLIDs_, and remoteLIDs_
00170 
00201     void setupSamePermuteRemote();
00202 
00228     void setupExport();
00229 
00231   };
00232 
00233   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00234   Import<LocalOrdinal,GlobalOrdinal,Node>::Import(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & source, 
00235                                                   const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & target) {
00236     ImportData_ = rcp(new ImportExportData<LocalOrdinal,GlobalOrdinal,Node>(source, target));
00237     // call subfunctions
00238     setupSamePermuteRemote();
00239     if (source->isDistributed()) {
00240       setupExport();
00241     }
00242     // don't need remoteGIDs_ anymore
00243     remoteGIDs_ = null;
00244   }
00245 
00246   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00247   Import<LocalOrdinal,GlobalOrdinal,Node>::Import(const Import<LocalOrdinal,GlobalOrdinal,Node> & import)
00248   : ImportData_(import.ImportData_)
00249   {}
00250 
00251   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00252   Import<LocalOrdinal,GlobalOrdinal,Node>::~Import() 
00253   {}
00254 
00255   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00256   size_t Import<LocalOrdinal,GlobalOrdinal,Node>::getNumSameIDs() const {
00257     return ImportData_->numSameIDs_;
00258   }
00259 
00260   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00261   size_t Import<LocalOrdinal,GlobalOrdinal,Node>::getNumPermuteIDs() const {
00262     return ImportData_->permuteFromLIDs_.size();
00263   }
00264 
00265   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00266   ArrayView<const LocalOrdinal> 
00267   Import<LocalOrdinal,GlobalOrdinal,Node>::getPermuteFromLIDs() const {
00268     return ImportData_->permuteFromLIDs_();
00269   }
00270 
00271   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00272   ArrayView<const LocalOrdinal>
00273   Import<LocalOrdinal,GlobalOrdinal,Node>::getPermuteToLIDs() const {
00274     return ImportData_->permuteToLIDs_();
00275   }
00276 
00277   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00278   size_t Import<LocalOrdinal,GlobalOrdinal,Node>::getNumRemoteIDs() const {
00279     return ImportData_->remoteLIDs_.size();
00280   }
00281 
00282   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00283   ArrayView<const LocalOrdinal> 
00284   Import<LocalOrdinal,GlobalOrdinal,Node>::getRemoteLIDs() const {
00285     return ImportData_->remoteLIDs_();
00286   }
00287 
00288   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00289   size_t Import<LocalOrdinal,GlobalOrdinal,Node>::getNumExportIDs() const {
00290     return ImportData_->exportLIDs_.size();
00291   }
00292 
00293   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00294   ArrayView<const LocalOrdinal> 
00295   Import<LocalOrdinal,GlobalOrdinal,Node>::getExportLIDs() const {
00296     return ImportData_->exportLIDs_();
00297   }
00298 
00299   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00300   ArrayView<const int> 
00301   Import<LocalOrdinal,GlobalOrdinal,Node>::getExportImageIDs() const {
00302     return ImportData_->exportImageIDs_();
00303   }
00304 
00305   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00306   const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & 
00307   Import<LocalOrdinal,GlobalOrdinal,Node>::getSourceMap() const {
00308     return ImportData_->source_;
00309   }
00310 
00311   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00312   const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & 
00313   Import<LocalOrdinal,GlobalOrdinal,Node>::getTargetMap() const {
00314     return ImportData_->target_;
00315   }
00316 
00317   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00318   Distributor & 
00319   Import<LocalOrdinal,GlobalOrdinal,Node>::getDistributor() const {
00320     return ImportData_->distributor_;
00321   }
00322 
00323   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00324   Import<LocalOrdinal,GlobalOrdinal,Node>& 
00325   Import<LocalOrdinal,GlobalOrdinal,Node>::operator=(const Import<LocalOrdinal,GlobalOrdinal,Node> & source) {
00326     ImportData_ = source.ImportData_;
00327     return *this;
00328   }
00329 
00330   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00331   void Import<LocalOrdinal,GlobalOrdinal,Node>::print(std::ostream& os) const {
00332     using Teuchos::getFancyOStream;
00333     using Teuchos::rcpFromRef;
00334     using std::endl;
00335 
00336     ArrayView<const LocalOrdinal> av;
00337     ArrayView<const int> avi;
00338     const RCP<const Comm<int> > & comm = getSourceMap()->getComm();
00339     const int myImageID = comm->getRank();
00340     const int numImages = comm->getSize();
00341     for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
00342       if (myImageID == imageCtr) {
00343         os << endl;
00344         if (myImageID == 0) { // this is the root node (only output this info once)
00345           os << "Import Data Members:" << endl;
00346         }
00347         os << "Image ID       : " << myImageID << endl;
00348         os << "permuteFromLIDs: {"; av = getPermuteFromLIDs(); std::copy(av.begin(),av.end(),std::ostream_iterator<LocalOrdinal>(os," ")); os << "}" << endl;
00349         os << "permuteToLIDs  : {"; av = getPermuteToLIDs();   std::copy(av.begin(),av.end(),std::ostream_iterator<LocalOrdinal>(os," ")); os << "}" << endl;
00350         os << "remoteLIDs     : {"; av = getRemoteLIDs();      std::copy(av.begin(),av.end(),std::ostream_iterator<LocalOrdinal>(os," ")); os << "}" << endl;
00351         os << "exportLIDs     : {"; av = getExportLIDs();      std::copy(av.begin(),av.end(),std::ostream_iterator<LocalOrdinal>(os," ")); os << "}" << endl;
00352         os << "exportImageIDs : {"; avi = getExportImageIDs();  std::copy(avi.begin(),avi.end(),std::ostream_iterator<int>(os," ")); os << "}" << endl;
00353         os << "numSameIDs     : " << getNumSameIDs() << endl;
00354         os << "numPermuteIDs  : " << getNumPermuteIDs() << endl;
00355         os << "numRemoteIDs   : " << getNumRemoteIDs() << endl;
00356         os << "numExportIDs   : " << getNumExportIDs() << endl;
00357       }
00358       // Do a few global ops to give I/O a chance to complete
00359       comm->barrier();
00360       comm->barrier();
00361       comm->barrier();
00362     }
00363     if (myImageID == 0) {
00364       os << endl << endl << "Source Map:" << endl << std::flush; 
00365     }
00366     comm->barrier();
00367     os << *getSourceMap();
00368     comm->barrier();
00369 
00370     if (myImageID == 0) {
00371       os << endl << endl << "Target Map:" << endl << std::flush; 
00372     }
00373     comm->barrier();
00374     os << *getTargetMap();
00375     comm->barrier();
00376 
00377     // It's also helpful for debugging to print the Distributor
00378     // object.  Epetra_Import::Print() does this, so we can do a
00379     // side-by-side comparison.
00380     if (myImageID == 0) {
00381       os << endl << endl << "Distributor:" << endl << std::flush;
00382     }
00383     comm->barrier();
00384     getDistributor().describe (*(getFancyOStream (rcpFromRef (os))),
00385              Teuchos::VERB_EXTREME);
00386   }
00387 
00388 
00389   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00390   void Import<LocalOrdinal,GlobalOrdinal,Node>::setupSamePermuteRemote() {
00391     const Map<LocalOrdinal,GlobalOrdinal,Node> & source = *getSourceMap();
00392     const Map<LocalOrdinal,GlobalOrdinal,Node> & target = *getTargetMap();
00393     ArrayView<const GlobalOrdinal> sourceGIDs = source.getNodeElementList();
00394     ArrayView<const GlobalOrdinal> targetGIDs = target.getNodeElementList();
00395 
00396     // -- compute numSameIDs_ ---
00397     // go through GID lists of source and target. if the ith GID on both is the same, 
00398     // increment numSameIDs_ and try the next. as soon as you come to a pair that don't
00399     // match, give up.
00400     typename ArrayView<const GlobalOrdinal>::iterator sourceIter = sourceGIDs.begin(),
00401                                                       targetIter = targetGIDs.begin();
00402     while( sourceIter != sourceGIDs.end() && targetIter != targetGIDs.end() && *sourceIter == *targetIter )
00403     {
00404       ++ImportData_->numSameIDs_;
00405       ++sourceIter;
00406       ++targetIter;
00407     }
00408     // targetIter should now point to the GID of the first non-same entry or the end of targetGIDs
00409 
00410     // -- compute numPermuteIDs and numRemoteIDs --
00411     // -- fill permuteToLIDs_, permuteFromLIDs_, remoteGIDs_, and remoteLIDs_ --
00412     // go through remaining entries in targetGIDs. if source owns that GID, 
00413     // increment numPermuteIDs_, and add entries to permuteToLIDs_ and permuteFromLIDs_.
00414     // otherwise increment numRemoteIDs_ and add entries to remoteLIDs_ and remoteGIDs_.
00415     remoteGIDs_ = rcp( new Array<GlobalOrdinal>() );
00416     for (; targetIter != targetGIDs.end(); ++targetIter) {
00417       if (source.isNodeGlobalElement(*targetIter)) {
00418         // both source and target list this GID (*targetIter)
00419         // determine the LIDs for this GID on both Maps and add them to the permutation lists
00420         ImportData_->permuteToLIDs_.push_back(target.getLocalElement(*targetIter));
00421         ImportData_->permuteFromLIDs_.push_back(source.getLocalElement(*targetIter));
00422       }
00423       else {
00424         // this GID is on another processor; store it, along with its destination LID on this processor
00425         remoteGIDs_->push_back(*targetIter);
00426         ImportData_->remoteLIDs_.push_back(target.getLocalElement(*targetIter));
00427       }
00428     }
00429 
00430     TPETRA_ABUSE_WARNING((getNumRemoteIDs() > 0) && !source.isDistributed(),std::runtime_error,
00431         "::setupSamePermuteRemote(): Target has remote LIDs but Source is not distributed globally."
00432         << std::endl << "Importing to a submap of the target map.");
00433   }
00434 
00435 
00436   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00437   void Import<LocalOrdinal,GlobalOrdinal,Node>::setupExport() {
00438     typedef typename Array<int>::difference_type size_type;
00439     const Map<LocalOrdinal,GlobalOrdinal,Node> & source = *getSourceMap();
00440 
00441     // For each entry remoteGIDs[i], remoteImageIDs[i] will contain
00442     // the process ID of the process that owns that GID.
00443     ArrayView<GlobalOrdinal> remoteGIDs = (*remoteGIDs_)();
00444     Array<int> remoteImageIDs(remoteGIDs.size());
00445     // lookup == IDNotPresent means that the source Map wasn't able to
00446     // figure out to which processes one or more of the GIDs in the
00447     // given list of remoteGIDs belong.
00448     //
00449     // The previous abuse warning said "The target Map has GIDs not
00450     // found in the source Map."  This statement could be confusing,
00451     // because it doesn't refer to ownership by the current process,
00452     // but rather to ownership by _any_ process participating in the
00453     // Map.  (It could not possibly refer to ownership by the current
00454     // process, since remoteGIDs is exactly the list of GIDs owned by
00455     // the target Map but not owned by the source Map.  It was
00456     // constructed that way by setupSamePermuteRemote().)  
00457     //
00458     // What this statement means is that the source and target Maps
00459     // don't contain the same set of GIDs globally (over all
00460     // processes).  That is, there is at least one GID owned by some
00461     // process in the target Map, which is not owned by _any_ process
00462     // in the source Map.
00463     const LookupStatus lookup = source.getRemoteIndexList(remoteGIDs, remoteImageIDs());
00464     TPETRA_ABUSE_WARNING( lookup == IDNotPresent, std::runtime_error, 
00465       "::setupExport(): the source Map wasn't able to figure out which process "
00466       "owns one or more of the GIDs in the list of remote GIDs.  This probably "
00467       "means that there is at least one GID owned by some process in the target"
00468       " Map which is not owned by any process in the source Map.  (That is, the"
00469       " source and target Maps do not contain the same set of GIDs globally.)");
00470 
00471     // Ignore remote GIDs that aren't owned by any process in the
00472     // source Map.  getRemoteIndexList() gives each of these a process
00473     // ID of -1.
00474     if ( lookup == IDNotPresent ) {
00475       const size_type numInvalidRemote = std::count_if( remoteImageIDs.begin(), remoteImageIDs.end(), std::bind1st(std::equal_to<int>(),-1) );
00476       // if all of them are invalid, we can delete the whole array
00477       const size_type totalNumRemote = getNumRemoteIDs();
00478       if (numInvalidRemote == totalNumRemote) {
00479         // all remotes are invalid; we have no remotes; we can delete the remotes
00480         remoteImageIDs.clear();
00481         (*remoteGIDs_).clear();
00482         ImportData_->remoteLIDs_.clear();
00483       }
00484       else {
00485         // Some remotes are valid; we need to keep the valid ones.
00486         // Pack and resize remoteImageIDs, remoteGIDs_, and
00487         // remoteLIDs_.
00488         size_type numValidRemote = 0;
00489         for (size_type r = 0; r < totalNumRemote; ++r) {
00490           if (remoteImageIDs[r] != -1) {
00491             remoteImageIDs[numValidRemote] = remoteImageIDs[r];
00492             (*remoteGIDs_)[numValidRemote] = (*remoteGIDs_)[r];
00493             ImportData_->remoteLIDs_[numValidRemote] = ImportData_->remoteLIDs_[r];
00494             ++numValidRemote;
00495           }
00496         }
00497         TEUCHOS_TEST_FOR_EXCEPTION( numValidRemote != totalNumRemote - numInvalidRemote, std::logic_error,
00498     "Tpetra::Import::setupExport(): After removing invalid remote GIDs and packing the valid remote "
00499           "GIDs, numValidRemote = " << numValidRemote << " != totalNumRemote - numInvalidRemote = "
00500           << totalNumRemote - numInvalidRemote << ".  Please report this bug to the Tpetra developers.");
00501 
00502         remoteImageIDs.resize(numValidRemote);
00503         (*remoteGIDs_).resize(numValidRemote);
00504         ImportData_->remoteLIDs_.resize(numValidRemote);
00505       }
00506       remoteGIDs = (*remoteGIDs_)();
00507     }
00508 
00509     // Sort remoteImageIDs in ascending order.  Apply the resulting
00510     // permutation to remoteGIDs_, so that remoteImageIDs[i] and
00511     // remoteGIDs_[i] refer to the same thing.
00512     sort2(remoteImageIDs.begin(), remoteImageIDs.end(), remoteGIDs.begin());
00513 
00514     // Call the Distributor's createFromRecvs() method to turn the
00515     // remote GIDs and their owning processes into a send-and-receive
00516     // communication plan.  remoteGIDs and remoteImageIDs_ are input;
00517     // exportGIDs and exportImageIDs_ are output arrays which are
00518     // allocated by createFromRecvs().
00519     ArrayRCP<GlobalOrdinal> exportGIDs;
00520     ImportData_->distributor_.createFromRecvs(remoteGIDs().getConst(), remoteImageIDs, exportGIDs, ImportData_->exportImageIDs_);
00521 
00522     // Find the LIDs corresponding to the GIDs in exportGIDs.  For
00523     // sparse matrix-vector multiply, this tells the calling process
00524     // how to index into the source vector to get the elements which
00525     // it needs to send.
00526     if (exportGIDs != null) {
00527       ImportData_->exportLIDs_ = arcp<LocalOrdinal>(exportGIDs.size());
00528     }
00529     typename ArrayRCP<LocalOrdinal>::iterator dst = ImportData_->exportLIDs_.begin();
00530     typename ArrayRCP<GlobalOrdinal>::const_iterator src = exportGIDs.begin();
00531     while (src != exportGIDs.end())
00532     {
00533       (*dst++) = source.getLocalElement(*src++);
00534     }
00535   }
00536 
00546   template <class LO, class GO, class Node> 
00547   RCP< const Import<LO,GO,Node> >
00548   createImport( const RCP<const Map<LO,GO,Node> > & src, 
00549                 const RCP<const Map<LO,GO,Node> > & tgt )
00550   {
00551     if (src == tgt) return null;
00552 #ifdef HAVE_TPETRA_DEBUG
00553     TEUCHOS_TEST_FOR_EXCEPTION(src == null || tgt == null, std::runtime_error,
00554         "Tpetra::createImport(): neither source nor target map may be null:\nsource: " << src << "\ntarget: " << tgt << "\n");
00555 #endif
00556     return rcp(new Import<LO,GO,Node>(src,tgt));
00557   }
00558 
00561   template <class LO, class GO, class Node> 
00562   RCP< const Import<LO,GO,Node> >
00563   TPETRA_DEPRECATED makeImport( const RCP<const Map<LO,GO,Node> > & src, 
00564                                 const RCP<const Map<LO,GO,Node> > & tgt )
00565   {
00566     return createImport<LO,GO,Node>(src,tgt);
00567   }
00568 
00569 } // namespace Tpetra
00570 
00571 #endif // TPETRA_IMPORT_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines