Epetra Package Browser (Single Doxygen Collection) Development
Epetra_Import.cpp
Go to the documentation of this file.
00001 /*
00002 //@HEADER
00003 // ************************************************************************
00004 // 
00005 //               Epetra: Linear Algebra Services Package 
00006 //                 Copyright 2011 Sandia Corporation
00007 // 
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00039 // 
00040 // ************************************************************************
00041 //@HEADER
00042 */
00043 
00044 #include "Epetra_Import.h"
00045 #include "Epetra_BlockMap.h"
00046 #include "Epetra_Distributor.h"
00047 #include "Epetra_Comm.h"
00048 #include "Epetra_Util.h"
00049 
00050 #include <algorithm>
00051 #include <vector>
00052 
00053 //==============================================================================
00054 // Epetra_Import constructor for a Epetra_BlockMap object
00055 Epetra_Import::Epetra_Import( const Epetra_BlockMap &  targetMap, const Epetra_BlockMap & sourceMap)
00056   : Epetra_Object("Epetra::Import"),
00057     TargetMap_(targetMap),
00058     SourceMap_(sourceMap),
00059     NumSameIDs_(0),
00060     NumPermuteIDs_(0),
00061     PermuteToLIDs_(0),
00062     PermuteFromLIDs_(0),
00063     NumRemoteIDs_(0),
00064     RemoteLIDs_(0),
00065     NumExportIDs_(0),
00066     ExportLIDs_(0),
00067     ExportPIDs_(0),
00068     NumSend_(0),
00069     NumRecv_(0),
00070     Distor_(0)
00071 {
00072 
00073   int i;
00074   
00075   // Build three ID lists:
00076   // NumSameIDs - Number of IDs in TargetMap and SourceMap that are identical, up to the first
00077   //              nonidentical ID.
00078   // NumPermuteIDs - Number of IDs in SourceMap that must be indirectly loaded but are on this processor.
00079   // NumRemoteIDs - Number of IDs that are in SourceMap but not in TargetMap, and thus must be imported.
00080   
00081   int NumSourceIDs = sourceMap.NumMyElements();
00082   int NumTargetIDs = targetMap.NumMyElements();
00083   
00084   int *TargetGIDs = 0;
00085   if (NumTargetIDs>0) {
00086     TargetGIDs = new int[NumTargetIDs];
00087     targetMap.MyGlobalElements(TargetGIDs);
00088   }
00089   
00090   int * SourceGIDs = 0;
00091   if (NumSourceIDs>0) {
00092     SourceGIDs = new int[NumSourceIDs];
00093     sourceMap.MyGlobalElements(SourceGIDs);
00094   }
00095   
00096   int MinIDs = EPETRA_MIN(NumSourceIDs, NumTargetIDs);
00097   
00098   
00099   NumSameIDs_ = 0;
00100   for (i=0; i< MinIDs; i++) if (TargetGIDs[i]==SourceGIDs[i]) NumSameIDs_++; else break;
00101   
00102   
00103   // Find count of Target IDs that are truly remote and those that are local but permuted
00104 
00105   NumPermuteIDs_ = 0;
00106   NumRemoteIDs_ = 0;
00107   for (i=NumSameIDs_; i< NumTargetIDs; i++) 
00108     if (sourceMap.MyGID(TargetGIDs[i])) NumPermuteIDs_++; // Check if Target GID is a local Source GID
00109     else NumRemoteIDs_++; // If not, then it is remote
00110   
00111   
00112   
00113   // Define remote and permutation lists
00114   
00115   int * RemoteGIDs=0;
00116   RemoteLIDs_ = 0;
00117   if (NumRemoteIDs_>0) {
00118     RemoteLIDs_ = new int[NumRemoteIDs_];
00119     RemoteGIDs = new int[NumRemoteIDs_];
00120   }
00121   if (NumPermuteIDs_>0)  {
00122     PermuteToLIDs_ = new int[NumPermuteIDs_];
00123     PermuteFromLIDs_ = new int[NumPermuteIDs_];
00124   }
00125   
00126   NumPermuteIDs_ = 0;
00127   NumRemoteIDs_ = 0;
00128   for (i=NumSameIDs_; i< NumTargetIDs; i++) {
00129     if (sourceMap.MyGID(TargetGIDs[i])) {
00130       PermuteToLIDs_[NumPermuteIDs_] = i;
00131       PermuteFromLIDs_[NumPermuteIDs_++] = sourceMap.LID(TargetGIDs[i]);
00132     }
00133     else {
00134       //NumRecv_ +=TargetMap.ElementSize(i); // Count total number of entries to receive
00135       NumRecv_ +=targetMap.MaxElementSize(); // Count total number of entries to receive (currently need max)
00136       RemoteGIDs[NumRemoteIDs_] = TargetGIDs[i];
00137       RemoteLIDs_[NumRemoteIDs_++] = i;
00138     }
00139   }
00140 
00141   if( NumRemoteIDs_>0 && !sourceMap.DistributedGlobal() )
00142     ReportError("Warning in Epetra_Import: Serial Import has remote IDs. (Importing to Subset of Target Map)", 1);
00143   
00144   // Test for distributed cases
00145   
00146   int * RemotePIDs = 0;
00147 
00148   if (sourceMap.DistributedGlobal()) {
00149     
00150     if (NumRemoteIDs_>0)  RemotePIDs = new int[NumRemoteIDs_];
00151     int ierr = sourceMap.RemoteIDList(NumRemoteIDs_, RemoteGIDs, RemotePIDs, 0); // Get remote PIDs
00152     if (ierr) throw ReportError("Error in sourceMap.RemoteIDList call", ierr);
00153 
00154     //Get rid of IDs that don't exist in SourceMap
00155     if(NumRemoteIDs_>0) {
00156       int cnt = 0;
00157       for( i = 0; i < NumRemoteIDs_; ++i )
00158         if( RemotePIDs[i] == -1 ) ++cnt;
00159       if( cnt ) {
00160         if( NumRemoteIDs_-cnt ) {
00161           int * NewRemoteGIDs = new int[NumRemoteIDs_-cnt];
00162           int * NewRemotePIDs = new int[NumRemoteIDs_-cnt];
00163           int * NewRemoteLIDs = new int[NumRemoteIDs_-cnt];
00164           cnt = 0;
00165           for( i = 0; i < NumRemoteIDs_; ++i )
00166             if( RemotePIDs[i] != -1 ) {
00167               NewRemoteGIDs[cnt] = RemoteGIDs[i];
00168               NewRemotePIDs[cnt] = RemotePIDs[i];
00169               NewRemoteLIDs[cnt] = targetMap.LID(RemoteGIDs[i]);
00170               ++cnt;
00171             }
00172           NumRemoteIDs_ = cnt;
00173           delete [] RemoteGIDs;
00174           delete [] RemotePIDs;
00175           delete [] RemoteLIDs_;
00176           RemoteGIDs = NewRemoteGIDs;
00177           RemotePIDs = NewRemotePIDs;
00178           RemoteLIDs_ = NewRemoteLIDs;
00179           ReportError("Warning in Epetra_Import: Target IDs not found in Source Map (Do you want to import to subset of Target Map?)", 1);
00180         }
00181         else { //valid RemoteIDs empty
00182           NumRemoteIDs_ = 0;
00183           delete [] RemoteGIDs;
00184           RemoteGIDs = 0;
00185           delete [] RemotePIDs;
00186           RemotePIDs = 0;
00187         }
00188       }
00189     }
00190 
00191     //Sort Remote IDs by processor so DoReverses will work
00192     Epetra_Util util;
00193     int * tmpPtr[2];
00194     tmpPtr[0] = RemoteLIDs_, tmpPtr[1] = RemoteGIDs;
00195     util.Sort(true,NumRemoteIDs_,RemotePIDs,0,0,2,tmpPtr);
00196 
00197     Distor_ = sourceMap.Comm().CreateDistributor();
00198     
00199     // Construct list of exports that calling processor needs to send as a result
00200     // of everyone asking for what it needs to receive.
00201     
00202     bool Deterministic = true;
00203     ierr = Distor_->CreateFromRecvs( NumRemoteIDs_, RemoteGIDs, RemotePIDs,
00204                        Deterministic, NumExportIDs_, ExportLIDs_, ExportPIDs_ );
00205     if (ierr!=0) throw ReportError("Error in Epetra_Distributor.CreateFromRecvs()", ierr);
00206 
00207     // Export IDs come in as GIDs, convert to LIDs
00208     for (i=0; i< NumExportIDs_; i++) {
00209       if (ExportPIDs_[i] < 0) throw ReportError("targetMap requested a GID that is not in the sourceMap.", -1);
00210       ExportLIDs_[i] = sourceMap.LID(ExportLIDs_[i]);
00211       NumSend_ += sourceMap.MaxElementSize(); // Count total number of entries to send (currently need max)
00212     }
00213   }
00214 
00215   if( NumRemoteIDs_>0 ) delete [] RemoteGIDs;
00216   if( NumRemoteIDs_>0 ) delete [] RemotePIDs;
00217 
00218   if (NumTargetIDs>0) delete [] TargetGIDs;
00219   if (NumSourceIDs>0) delete [] SourceGIDs;
00220   
00221   return;
00222 }
00223 
00224 //==============================================================================
00225 // Epetra_Import copy constructor 
00226 Epetra_Import::Epetra_Import(const Epetra_Import & Importer)
00227   : Epetra_Object(Importer),
00228     TargetMap_(Importer.TargetMap_),
00229     SourceMap_(Importer.SourceMap_),
00230     NumSameIDs_(Importer.NumSameIDs_),
00231     NumPermuteIDs_(Importer.NumPermuteIDs_),
00232     PermuteToLIDs_(0),
00233     PermuteFromLIDs_(0),
00234     NumRemoteIDs_(Importer.NumRemoteIDs_),
00235     RemoteLIDs_(0),
00236     NumExportIDs_(Importer.NumExportIDs_),
00237     ExportLIDs_(0),
00238     ExportPIDs_(0),
00239     NumSend_(Importer.NumSend_),
00240     NumRecv_(Importer.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] = Importer.PermuteToLIDs_[i];
00249       PermuteFromLIDs_[i] = Importer.PermuteFromLIDs_[i];
00250     }
00251   }
00252 
00253   if (NumRemoteIDs_>0) {
00254     RemoteLIDs_ = new int[NumRemoteIDs_];
00255     for (i=0; i< NumRemoteIDs_; i++) RemoteLIDs_[i] = Importer.RemoteLIDs_[i];
00256   }
00257 
00258   if (NumExportIDs_>0) {
00259     ExportLIDs_ = new int[NumExportIDs_];
00260     ExportPIDs_ = new int[NumExportIDs_];
00261     for (i=0; i< NumExportIDs_; i++) {
00262       ExportLIDs_[i] = Importer.ExportLIDs_[i];
00263       ExportPIDs_[i] = Importer.ExportPIDs_[i];
00264     }
00265   }
00266 
00267   if (Importer.Distor_!=0) Distor_ = Importer.Distor_->Clone();
00268 
00269 }
00270 
00271 //==============================================================================
00272 // Epetra_Import destructor 
00273 Epetra_Import::~Epetra_Import()
00274 {
00275   if( Distor_ != 0 ) delete Distor_;
00276   if (RemoteLIDs_ != 0) delete [] RemoteLIDs_;
00277   if (PermuteToLIDs_ != 0) delete [] PermuteToLIDs_;
00278   if (PermuteFromLIDs_ != 0) delete [] PermuteFromLIDs_;
00279 
00280   if( ExportPIDs_ != 0 ) delete [] ExportPIDs_; // These were created by Distor_
00281   if( ExportLIDs_ != 0 ) delete [] ExportLIDs_;
00282 }
00283 //=============================================================================
00284 void Epetra_Import::Print(ostream & os) const
00285 {
00286   // mfh 14 Dec 2011: The implementation of Print() I found here
00287   // previously didn't print much at all, and it included a message
00288   // saying that it wasn't finished ("Epetra_Import::Print needs
00289   // attention!!!").  What you see below is a port of
00290   // Tpetra::Import::print, which does have a full implementation.
00291   // This should allow a side-by-side comparison of Epetra_Import with
00292   // Tpetra::Import.
00293 
00294   // If true, then copy the array data and sort it before printing.
00295   // Otherwise, leave the data in its original order.  
00296   //
00297   // NOTE: Do NOT sort the arrays in place!  Only sort in the copy.
00298   // Epetra depends on the order being preserved, and some arrays'
00299   // orders are coupled.
00300   const bool sortIDs = false;
00301 
00302   const Epetra_Comm& comm = SourceMap_.Comm();
00303   const int myRank = comm.MyPID();
00304   const int numProcs = comm.NumProc();
00305   
00306   if (myRank == 0) {
00307     os << "Import Data Members:" << endl;
00308   }
00309   // We don't need a barrier before this for loop, because Proc 0 is
00310   // the first one to do anything in the for loop anyway.
00311   for (int p = 0; p < numProcs; ++p) {
00312     if (myRank == p) {
00313       os << "Image ID       : " << myRank << endl;
00314 
00315       os << "permuteFromLIDs:";
00316       if (PermuteFromLIDs_ == NULL) {
00317   os << " NULL";
00318       } else {
00319   std::vector<int> permuteFromLIDs (NumPermuteIDs_);
00320   std::copy (PermuteFromLIDs_, PermuteFromLIDs_ + NumPermuteIDs_, 
00321        permuteFromLIDs.begin());
00322   if (sortIDs) {
00323     std::sort (permuteFromLIDs.begin(), permuteFromLIDs.end());
00324   }
00325   os << " {";
00326   for (int i = 0; i < NumPermuteIDs_; ++i) {
00327     os << permuteFromLIDs[i];
00328     if (i < NumPermuteIDs_ - 1) {
00329       os << ", ";
00330     }
00331   }
00332   os << "}";
00333       }
00334       os << endl;
00335 
00336       os << "permuteToLIDs  :";
00337       if (PermuteToLIDs_ == NULL) {
00338   os << " NULL";
00339       } else {
00340   std::vector<int> permuteToLIDs (NumPermuteIDs_);
00341   std::copy (PermuteToLIDs_, PermuteToLIDs_ + NumPermuteIDs_, 
00342        permuteToLIDs.begin());
00343   if (sortIDs) {
00344     std::sort (permuteToLIDs.begin(), permuteToLIDs.end());
00345   }
00346   os << " {";
00347   for (int i = 0; i < NumPermuteIDs_; ++i) {
00348     os << permuteToLIDs[i];
00349     if (i < NumPermuteIDs_ - 1) {
00350       os << ", ";
00351     }
00352   }
00353   os << "}";
00354       }
00355       os << endl;
00356 
00357       os << "remoteLIDs     :";
00358       if (RemoteLIDs_ == NULL) {
00359   os << " NULL";
00360       } else {
00361   std::vector<int> remoteLIDs (NumRemoteIDs_);
00362   std::copy (RemoteLIDs_, RemoteLIDs_ + NumRemoteIDs_, 
00363        remoteLIDs.begin());
00364   if (sortIDs) {
00365     std::sort (remoteLIDs.begin(), remoteLIDs.end());
00366   }
00367   os << " {";
00368   for (int i = 0; i < NumRemoteIDs_; ++i) {
00369     os << remoteLIDs[i];
00370     if (i < NumRemoteIDs_ - 1) {
00371       os << ", ";
00372     }
00373   }
00374   os << "}";
00375       }
00376       os << endl;
00377 
00378       // If sorting for output, the export LIDs and export PIDs have
00379       // to be sorted together.  We can use Epetra_Util::Sort, using
00380       // the PIDs as the keys to match Tpetra::Import.
00381       std::vector<int> exportLIDs (NumExportIDs_);
00382       std::vector<int> exportPIDs (NumExportIDs_);
00383       if (ExportLIDs_ != NULL) {
00384   std::copy (ExportLIDs_, ExportLIDs_ + NumExportIDs_, exportLIDs.begin());
00385   std::copy (ExportPIDs_, ExportPIDs_ + NumExportIDs_, exportPIDs.begin());
00386 
00387   if (sortIDs && NumExportIDs_ > 0) {
00388     int* intCompanions[1]; // Input for Epetra_Util::Sort().
00389     intCompanions[0] = &exportLIDs[0];
00390     Epetra_Util::Sort (true, NumExportIDs_, &exportPIDs[0], 
00391            0, (double**) NULL, 1, intCompanions);
00392   }
00393       }
00394 
00395       os << "exportLIDs     :";
00396       if (ExportLIDs_ == NULL) {
00397   os << " NULL";
00398       } else {
00399   os << " {";
00400   for (int i = 0; i < NumExportIDs_; ++i) {
00401     os << exportLIDs[i];
00402     if (i < NumExportIDs_ - 1) {
00403       os << ", ";
00404     }
00405   }
00406   os << "}";
00407       }
00408       os << endl;
00409 
00410       os << "exportImageIDs :";
00411       if (ExportPIDs_ == NULL) {
00412   os << " NULL";
00413       } else {
00414   os << " {";
00415   for (int i = 0; i < NumExportIDs_; ++i) {
00416     os << exportPIDs[i];
00417     if (i < NumExportIDs_ - 1) {
00418       os << ", ";
00419     }
00420   }
00421   os << "}";
00422       }
00423       os << endl;
00424 
00425       os << "numSameIDs     : " << NumSameIDs_ << endl;
00426       os << "numPermuteIDs  : " << NumPermuteIDs_ << endl;
00427       os << "numRemoteIDs   : " << NumRemoteIDs_ << endl;
00428       os << "numExportIDs   : " << NumExportIDs_ << endl;
00429 
00430       // Epetra keeps NumSend_ and NumRecv_, whereas in Tpetra, these
00431       // are stored in the Distributor object.  This is why we print
00432       // them here.
00433       os << "Number of sends: " << NumSend_ << endl;
00434       os << "Number of recvs: " << NumRecv_ << endl;
00435     } // if my rank is p
00436 
00437     // A few global barriers give I/O a chance to complete.
00438     comm.Barrier();
00439     comm.Barrier();
00440     comm.Barrier();
00441   } // for each rank p
00442 
00443   const bool printMaps = false;
00444   if (printMaps) {
00445     // The original implementation printed the Maps first.  We moved
00446     // printing the Maps to the end, for easy comparison with the
00447     // output of Tpetra::Import::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 
00463   if (myRank == 0) {
00464     os << endl << endl << "Distributor:" << endl << std::flush;
00465   }
00466   comm.Barrier();
00467   if (Distor_ == NULL) {
00468     if (myRank == 0) {
00469       os << " is NULL." << endl;
00470     }
00471   } else {
00472     Distor_->Print(os); // Printing the Distributor is itself distributed.
00473   }
00474   comm.Barrier();
00475 }
00476 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines