Epetra Package Browser (Single Doxygen Collection) Development
Epetra_Export.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 // 
00004 //               Epetra: Linear Algebra Services Package 
00005 //                 Copyright 2011 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 #include "Epetra_Export.h"
00043 #include "Epetra_BlockMap.h"
00044 #include "Epetra_Distributor.h"
00045 #include "Epetra_Comm.h"
00046 #include "Epetra_Util.h"
00047 #include <vector>
00048 
00049 //==============================================================================
00050 // Epetra_Export constructor for a Epetra_BlockMap object
00051 Epetra_Export::Epetra_Export( const Epetra_BlockMap &  sourceMap, const Epetra_BlockMap & targetMap)
00052   : Epetra_Object("Epetra::Export"), 
00053     TargetMap_(targetMap),
00054     SourceMap_(sourceMap),
00055     NumSameIDs_(0),
00056     NumPermuteIDs_(0),
00057     PermuteToLIDs_(0),
00058     PermuteFromLIDs_(0),
00059     NumRemoteIDs_(0),
00060     RemoteLIDs_(0),
00061     NumExportIDs_(0),
00062     ExportLIDs_(0),
00063     ExportPIDs_(0),
00064     NumSend_(0),
00065     NumRecv_(0),
00066     Distor_(0)
00067 {
00068 
00069   int i;
00070 
00071   // Build three ID lists:
00072   // NumSameIDs - Number of IDs in TargetMap and SourceMap that are identical, up to the first
00073   //              nonidentical ID.
00074   // NumPermuteIDs - Number of IDs in SourceMap that must be indirectly loaded but are on this processor.
00075   // NumExportIDs - Number of IDs that are in SourceMap but not in TargetMap, and thus must be exported.
00076 
00077   int NumSourceIDs = sourceMap.NumMyElements();
00078   int NumTargetIDs = targetMap.NumMyElements();
00079 
00080   int *TargetGIDs = 0;
00081   if (NumTargetIDs>0) {
00082     TargetGIDs = new int[NumTargetIDs];
00083     targetMap.MyGlobalElements(TargetGIDs);
00084   }
00085 
00086   int * SourceGIDs = 0;
00087   if (NumSourceIDs>0) {
00088     SourceGIDs = new int[NumSourceIDs];
00089     sourceMap.MyGlobalElements(SourceGIDs);
00090   }
00091 
00092   int MinIDs = EPETRA_MIN(NumSourceIDs, NumTargetIDs);
00093 
00094   NumSameIDs_ = 0;
00095   for (i=0; i< MinIDs; i++) if (TargetGIDs[i]==SourceGIDs[i]) NumSameIDs_++; else break;
00096 
00097   // Find count of Source IDs that are truly remote and those that are local but permuted
00098 
00099   NumPermuteIDs_ = 0;
00100   NumExportIDs_ = 0;
00101   for (i=NumSameIDs_; i< NumSourceIDs; i++) 
00102     if (targetMap.MyGID(SourceGIDs[i])) NumPermuteIDs_++; // Check if Source GID is a local Target GID
00103     else NumExportIDs_++; // If not, then it is remote
00104 
00105   // Define remote and permutation lists
00106 
00107   int * ExportGIDs = 0;
00108   if (NumExportIDs_>0) {
00109     ExportLIDs_ = new int[NumExportIDs_];
00110     ExportGIDs = new int[NumExportIDs_];
00111   }
00112   if (NumPermuteIDs_>0)  {
00113     PermuteToLIDs_ = new int[NumPermuteIDs_];
00114     PermuteFromLIDs_ = new int[NumPermuteIDs_];
00115   }
00116 
00117   NumPermuteIDs_ = 0;
00118   NumExportIDs_ = 0;
00119   for (i=NumSameIDs_; i< NumSourceIDs; i++) {
00120     if (targetMap.MyGID(SourceGIDs[i])) {
00121       PermuteFromLIDs_[NumPermuteIDs_] = i;
00122       PermuteToLIDs_[NumPermuteIDs_++] = targetMap.LID(SourceGIDs[i]);
00123     }
00124     else {
00125       //NumSend_ +=sourceMap.ElementSize(i); // Count total number of entries to send
00126       NumSend_ +=sourceMap.MaxElementSize(); // Count total number of entries to send (currently need max)
00127       ExportGIDs[NumExportIDs_] = SourceGIDs[i];
00128       ExportLIDs_[NumExportIDs_++] = i;
00129     }
00130   }
00131      
00132   if ( NumExportIDs_>0 && !sourceMap.DistributedGlobal()) 
00133     ReportError("Warning in Epetra_Export: Serial Export has remote IDs. (Exporting from Subset of Source Map)", 1);
00134 
00135   // Test for distributed cases
00136   int ierr = 0;
00137 
00138   if (sourceMap.DistributedGlobal()) {
00139 
00140     if (NumExportIDs_>0) ExportPIDs_ = new int[NumExportIDs_];
00141     ierr = targetMap.RemoteIDList(NumExportIDs_, ExportGIDs, ExportPIDs_, 0); // Get remote PIDs
00142     if( ierr ) throw ReportError("Error in Epetra_BlockMap::RemoteIDList", ierr);
00143 
00144     //Get rid of IDs not in Target Map
00145     if(NumExportIDs_>0) {
00146       int cnt = 0;
00147       for( i = 0; i < NumExportIDs_; ++i )
00148   if( ExportPIDs_[i] == -1 ) ++cnt;
00149       if( cnt ) {
00150   int * NewExportGIDs = 0;
00151   int * NewExportPIDs = 0;
00152   int * NewExportLIDs = 0;
00153   int cnt1 = NumExportIDs_-cnt;
00154   if (cnt1) {
00155     NewExportGIDs = new int[cnt1];
00156     NewExportPIDs = new int[cnt1];
00157     NewExportLIDs = new int[cnt1];
00158   }
00159   cnt = 0;
00160   for( i = 0; i < NumExportIDs_; ++i )
00161     if( ExportPIDs_[i] != -1 ) {
00162       NewExportGIDs[cnt] = ExportGIDs[i];
00163       NewExportPIDs[cnt] = ExportPIDs_[i];
00164       NewExportLIDs[cnt] = ExportLIDs_[i];
00165       ++cnt;
00166           }
00167   assert(cnt==cnt1); // Sanity test
00168   NumExportIDs_ = cnt;
00169   delete [] ExportGIDs;
00170   delete [] ExportPIDs_;
00171   delete [] ExportLIDs_;
00172   ExportGIDs = NewExportGIDs;
00173   ExportPIDs_ = NewExportPIDs;
00174   ExportLIDs_ = NewExportLIDs;
00175   ReportError("Warning in Epetra_Export: Source IDs not found in Target Map (Do you want to export from subset of Source Map?)", 1 );
00176       }
00177     }
00178     
00179     //Make sure Export IDs are ordered by processor
00180     Epetra_Util util;
00181     int * tmpPtr[2];
00182     tmpPtr[0] = ExportLIDs_, tmpPtr[1] = ExportGIDs;
00183     util.Sort(true,NumExportIDs_,ExportPIDs_,0,0,2,tmpPtr);
00184 
00185     Distor_ = sourceMap.Comm().CreateDistributor();
00186     
00187     // Construct list of exports that calling processor needs to send as a result
00188     // of everyone asking for what it needs to receive.
00189     
00190     ierr = Distor_->CreateFromSends( NumExportIDs_, ExportPIDs_, true, NumRemoteIDs_);
00191     if (ierr!=0) throw ReportError("Error in Epetra_Distributor.CreateFromSends()", ierr);
00192     
00193     // Use comm plan with ExportGIDs to find out who is sending to us and
00194     // get proper ordering of GIDs for remote entries 
00195     // (that we will convert to LIDs when done).
00196     
00197     if (NumRemoteIDs_>0) RemoteLIDs_ = new int[NumRemoteIDs_]; // Allocate space for LIDs in target that are
00198     // going to get something from off-processor.
00199     char * cRemoteGIDs = 0; //Do will alloc memory for this object
00200     int LenCRemoteGIDs = 0;
00201     ierr = Distor_->Do(reinterpret_cast<char *> (ExportGIDs), 
00202     sizeof( int ),
00203     LenCRemoteGIDs,
00204     cRemoteGIDs);
00205     if (ierr) throw ReportError("Error in Epetra_Distributor.Do()", ierr);
00206     int * RemoteGIDs = reinterpret_cast<int*>(cRemoteGIDs);
00207 
00208     // Remote IDs come in as GIDs, convert to LIDs
00209     for (i=0; i< NumRemoteIDs_; i++) {
00210       RemoteLIDs_[i] = targetMap.LID(RemoteGIDs[i]);
00211       //NumRecv_ += targetMap.ElementSize(RemoteLIDs_[i]); // Count total number of entries to receive
00212       NumRecv_ += targetMap.MaxElementSize(); // Count total number of entries to receive (currently need max)
00213     }
00214 
00215     if (LenCRemoteGIDs>0) delete [] cRemoteGIDs;
00216   }
00217   if (NumExportIDs_>0) delete [] ExportGIDs;
00218   if (NumTargetIDs>0) delete [] TargetGIDs;
00219   if (NumSourceIDs>0) delete [] SourceGIDs;
00220   
00221   return;
00222 }
00223 
00224 //==============================================================================
00225 // Epetra_Export copy constructor 
00226 Epetra_Export::Epetra_Export(const Epetra_Export & Exporter)
00227   : Epetra_Object(Exporter), 
00228      TargetMap_(Exporter.TargetMap_),
00229     SourceMap_(Exporter.SourceMap_),
00230     NumSameIDs_(Exporter.NumSameIDs_),
00231     NumPermuteIDs_(Exporter.NumPermuteIDs_),
00232     PermuteToLIDs_(0),
00233     PermuteFromLIDs_(0),
00234     NumRemoteIDs_(Exporter.NumRemoteIDs_),
00235     RemoteLIDs_(0),
00236     NumExportIDs_(Exporter.NumExportIDs_),
00237     ExportLIDs_(0),
00238     ExportPIDs_(0),
00239     NumSend_(Exporter.NumSend_),
00240     NumRecv_(Exporter.NumRecv_),
00241     Distor_(0)
00242 {
00243   int i;
00244   if (NumPermuteIDs_>0) {
00245     PermuteToLIDs_ = new int[NumPermuteIDs_];
00246     PermuteFromLIDs_ = new int[NumPermuteIDs_];
00247     for (i=0; i< NumPermuteIDs_; i++) {
00248       PermuteToLIDs_[i] = Exporter.PermuteToLIDs_[i];
00249       PermuteFromLIDs_[i] = Exporter.PermuteFromLIDs_[i];
00250     }
00251   }
00252 
00253   if (NumRemoteIDs_>0) {
00254     RemoteLIDs_ = new int[NumRemoteIDs_];
00255     for (i=0; i< NumRemoteIDs_; i++) RemoteLIDs_[i] = Exporter.RemoteLIDs_[i];
00256   }
00257 
00258   TargetMap().Comm().Barrier();
00259   if (NumExportIDs_>0) {
00260     ExportLIDs_ = new int[NumExportIDs_];
00261     ExportPIDs_ = new int[NumExportIDs_];
00262     for (i=0; i< NumExportIDs_; i++) {
00263       ExportLIDs_[i] = Exporter.ExportLIDs_[i];
00264       ExportPIDs_[i] = Exporter.ExportPIDs_[i];
00265     }
00266   }
00267 
00268   if (Exporter.Distor_!=0) Distor_ = Exporter.Distor_->Clone();
00269 
00270 }
00271 
00272 //==============================================================================
00273 // Epetra_Export destructor 
00274 Epetra_Export::~Epetra_Export()
00275 {
00276   if( Distor_ != 0 ) delete Distor_;
00277 
00278   if (RemoteLIDs_ != 0) delete [] RemoteLIDs_;
00279   if (PermuteToLIDs_ != 0) delete [] PermuteToLIDs_;
00280   if (PermuteFromLIDs_ != 0) delete [] PermuteFromLIDs_;
00281 
00282   if( ExportPIDs_ != 0 ) delete [] ExportPIDs_; // These were created by GSPlan
00283   if( ExportLIDs_ != 0 ) delete [] ExportLIDs_;
00284 }
00285 //=============================================================================
00286 void Epetra_Export::Print(ostream & os) const
00287 {
00288   // mfh 05 Jan 2012: The implementation of Print() I found here
00289   // previously didn't print much at all, and it included a message
00290   // saying that it wasn't finished ("Epetra_Export Print needs
00291   // attention!!!!").  What you see below is a port of
00292   // Tpetra::Export::print, which does have a full implementation.
00293   // This should allow a side-by-side comparison of Epetra_Export with
00294   // Tpetra::Export.
00295 
00296   // If true, then copy the array data and sort it before printing.
00297   // Otherwise, leave the data in its original order.  
00298   //
00299   // NOTE: Do NOT sort the arrays in place!  Only sort in the copy.
00300   // Epetra depends on the order being preserved, and some arrays'
00301   // orders are coupled.
00302   const bool sortIDs = true;
00303 
00304   const Epetra_Comm& comm = SourceMap_.Comm();
00305   const int myRank = comm.MyPID();
00306   const int numProcs = comm.NumProc();
00307 
00308   if (myRank == 0) {
00309     os << "Export Data Members:" << endl;
00310   }
00311   // We don't need a barrier before this for loop, because Proc 0 is
00312   // the first one to do anything in the for loop anyway.
00313   for (int p = 0; p < numProcs; ++p) {
00314     if (myRank == p) {
00315       os << "Image ID       : " << myRank << endl;
00316 
00317       os << "permuteFromLIDs:";
00318       if (PermuteFromLIDs_ == NULL) {
00319   os << " NULL";
00320       } else {
00321   std::vector<int> permuteFromLIDs (NumPermuteIDs_);
00322   std::copy (PermuteFromLIDs_, PermuteFromLIDs_ + NumPermuteIDs_, 
00323        permuteFromLIDs.begin());
00324   if (sortIDs) {
00325     std::sort (permuteFromLIDs.begin(), permuteFromLIDs.end());
00326   }
00327   os << " {";
00328   for (int i = 0; i < NumPermuteIDs_; ++i) {
00329     os << permuteFromLIDs[i];
00330     if (i < NumPermuteIDs_ - 1) {
00331       os << " ";
00332     }
00333   }
00334   os << "}";
00335       }
00336       os << endl;
00337 
00338       os << "permuteToLIDs  :";
00339       if (PermuteToLIDs_ == NULL) {
00340   os << " NULL";
00341       } else {
00342   std::vector<int> permuteToLIDs (NumPermuteIDs_);
00343   std::copy (PermuteToLIDs_, PermuteToLIDs_ + NumPermuteIDs_, 
00344        permuteToLIDs.begin());
00345   if (sortIDs) {
00346     std::sort (permuteToLIDs.begin(), permuteToLIDs.end());
00347   }
00348   os << " {";
00349   for (int i = 0; i < NumPermuteIDs_; ++i) {
00350     os << permuteToLIDs[i];
00351     if (i < NumPermuteIDs_ - 1) {
00352       os << " ";
00353     }
00354   }
00355   os << "}";
00356       }
00357       os << endl;
00358 
00359       os << "remoteLIDs     :";
00360       if (RemoteLIDs_ == NULL) {
00361   os << " NULL";
00362       } else {
00363   std::vector<int> remoteLIDs (NumRemoteIDs_);
00364   std::copy (RemoteLIDs_, RemoteLIDs_ + NumRemoteIDs_, 
00365        remoteLIDs.begin());
00366   if (sortIDs) {
00367     std::sort (remoteLIDs.begin(), remoteLIDs.end());
00368   }
00369   os << " {";
00370   for (int i = 0; i < NumRemoteIDs_; ++i) {
00371     os << remoteLIDs[i];
00372     if (i < NumRemoteIDs_ - 1) {
00373       os << " ";
00374     }
00375   }
00376   os << "}";
00377       }
00378       os << endl;
00379 
00380       // If sorting for output, the export LIDs and export PIDs have
00381       // to be sorted together.  We can use Epetra_Util::Sort, using
00382       // the PIDs as the keys to match Tpetra::Export.
00383       std::vector<int> exportLIDs (NumExportIDs_);
00384       std::vector<int> exportPIDs (NumExportIDs_);
00385       if (ExportLIDs_ != NULL) {
00386   std::copy (ExportLIDs_, ExportLIDs_ + NumExportIDs_, exportLIDs.begin());
00387   std::copy (ExportPIDs_, ExportPIDs_ + NumExportIDs_, exportPIDs.begin());
00388 
00389   if (sortIDs && NumExportIDs_ > 0) {
00390     int* intCompanions[1]; // Input for Epetra_Util::Sort().
00391     intCompanions[0] = &exportLIDs[0];
00392     Epetra_Util::Sort (true, NumExportIDs_, &exportPIDs[0], 
00393            0, (double**) NULL, 1, intCompanions);
00394   }
00395       }
00396 
00397       os << "exportLIDs     :";
00398       if (ExportLIDs_ == NULL) {
00399   os << " NULL";
00400       } else {
00401   os << " {";
00402   for (int i = 0; i < NumExportIDs_; ++i) {
00403     os << exportLIDs[i];
00404     if (i < NumExportIDs_ - 1) {
00405       os << " ";
00406     }
00407   }
00408   os << "}";
00409       }
00410       os << endl;
00411 
00412       os << "exportImageIDs :";
00413       if (ExportPIDs_ == NULL) {
00414   os << " NULL";
00415       } else {
00416   os << " {";
00417   for (int i = 0; i < NumExportIDs_; ++i) {
00418     os << exportPIDs[i];
00419     if (i < NumExportIDs_ - 1) {
00420       os << " ";
00421     }
00422   }
00423   os << "}";
00424       }
00425       os << endl;
00426 
00427       os << "numSameIDs     : " << NumSameIDs_ << endl;
00428       os << "numPermuteIDs  : " << NumPermuteIDs_ << endl;
00429       os << "numRemoteIDs   : " << NumRemoteIDs_ << endl;
00430       os << "numExportIDs   : " << NumExportIDs_ << endl;
00431 
00432       // Epetra keeps NumSend_ and NumRecv_, whereas in Tpetra, these
00433       // are stored in the Distributor object.  This is why we print
00434       // them here.
00435       os << "Number of sends: " << NumSend_ << endl;
00436       os << "Number of recvs: " << NumRecv_ << endl;
00437     } // if my rank is p
00438 
00439     // A few global barriers give I/O a chance to complete.
00440     comm.Barrier();
00441     comm.Barrier();
00442     comm.Barrier();
00443   } // for each rank p
00444 
00445   // The original implementation printed the Maps first.  We moved
00446   // printing the Maps to the end, for easy comparison with the output
00447   // of Tpetra::Export::print().
00448   if (myRank == 0) {
00449     os << endl << endl << "Source Map:" << endl << std::flush;
00450   }
00451   comm.Barrier();
00452   SourceMap_.Print(os);
00453   comm.Barrier();
00454   
00455   if (myRank == 0) {
00456     os << endl << endl << "Target Map:" << endl << std::flush;
00457   }
00458   comm.Barrier();
00459   TargetMap_.Print(os);
00460   comm.Barrier();
00461 
00462   if (myRank == 0) {
00463     os << endl << endl << "Distributor:" << endl << std::flush;
00464   }
00465   comm.Barrier();
00466   if (Distor_ == NULL) {
00467     if (myRank == 0) {
00468       os << " is NULL." << endl;
00469     }
00470   } else {
00471     Distor_->Print(os); // Printing the Distributor is itself distributed.
00472   }
00473   comm.Barrier();
00474 }
00475 
00476 //----------------------------------------------------------------------------
00477 Epetra_Export& Epetra_Export::operator=(const Epetra_Export& src)
00478 {
00479   (void)src;
00480   //not currently supported
00481   bool throw_err = true;
00482   if (throw_err) {
00483     throw ReportError("Epetra_Export::operator= not supported.",-1);
00484   }
00485   return(*this);
00486 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines