Epetra_BlockMap.cpp

Go to the documentation of this file.
00001 
00002 //@HEADER
00003 // ************************************************************************
00004 // 
00005 //               Epetra: Linear Algebra Services Package 
00006 //                 Copyright (2001) Sandia Corporation
00007 // 
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 // 
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //  
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //  
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA
00025 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00026 // 
00027 // ************************************************************************
00028 //@HEADER
00029 
00030 #include "Epetra_BlockMap.h"
00031 #include "Epetra_Comm.h"
00032 #include "Epetra_Directory.h"
00033 #include "Epetra_IntSerialDenseVector.h"
00034 #include "Epetra_HashTable.h"
00035 
00036 // Use the new LID hash table approach by default
00037 #define EPETRA_BLOCKMAP_NEW_LID
00038 
00039 //==============================================================================
00040 // Epetra_BlockMap constructor for a Epetra-defined uniform linear distribution of constant size elements.
00041 Epetra_BlockMap::Epetra_BlockMap(int NumGlobalElements, int ElementSize, int IndexBase, const Epetra_Comm& Comm)
00042   : Epetra_Object("Epetra::BlockMap"),
00043     BlockMapData_(0)
00044 {
00045   
00046   // Each processor gets roughly numGlobalPoints/p points
00047   // This routine automatically defines a linear partitioning of a
00048   // map with numGlobalPoints across the processors
00049   // specified in the given Epetra_Comm
00050   
00051   if (NumGlobalElements < 0) 
00052     throw ReportError("NumGlobalElements = " + toString(NumGlobalElements) + ".  Should be >= 0.", -1);
00053   if (ElementSize <= 0) 
00054     throw ReportError("ElementSize = " + toString(ElementSize) + ".  Should be > 0.", -2);
00055   
00056   BlockMapData_ = new Epetra_BlockMapData(NumGlobalElements, ElementSize, IndexBase, Comm);
00057   int NumProc = Comm.NumProc();
00058   BlockMapData_->ConstantElementSize_ = true;
00059   BlockMapData_->LinearMap_ = true;
00060 
00061   int MyPID = Comm.MyPID();
00062   BlockMapData_->NumMyElements_ = BlockMapData_->NumGlobalElements_ / NumProc;
00063   int remainder = BlockMapData_->NumGlobalElements_ % NumProc;
00064   int start_index = MyPID * (BlockMapData_->NumMyElements_ + 1);
00065 
00066   if (MyPID < remainder) 
00067     BlockMapData_->NumMyElements_++;
00068   else 
00069     start_index -= (MyPID - remainder);
00070 
00071   BlockMapData_->NumGlobalPoints_ = BlockMapData_->NumGlobalElements_ * BlockMapData_->ElementSize_;
00072   BlockMapData_->NumMyPoints_ = BlockMapData_->NumMyElements_ * BlockMapData_->ElementSize_;
00073 
00074   BlockMapData_->MinMyElementSize_ = BlockMapData_->ElementSize_;
00075   BlockMapData_->MaxMyElementSize_ = BlockMapData_->ElementSize_;
00076   BlockMapData_->MinElementSize_ = BlockMapData_->ElementSize_;
00077   BlockMapData_->MaxElementSize_ = BlockMapData_->ElementSize_;
00078 
00079   BlockMapData_->MinAllGID_ = BlockMapData_->IndexBase_;
00080   BlockMapData_->MaxAllGID_ = BlockMapData_->MinAllGID_ + BlockMapData_->NumGlobalElements_ - 1;
00081   BlockMapData_->MinMyGID_ = start_index + BlockMapData_->IndexBase_;
00082   BlockMapData_->MaxMyGID_ = BlockMapData_->MinMyGID_ + BlockMapData_->NumMyElements_ - 1;
00083   BlockMapData_->DistributedGlobal_ = IsDistributedGlobal(BlockMapData_->NumGlobalElements_, BlockMapData_->NumMyElements_);
00084   BlockMapData_->OneToOne_ = DetermineIsOneToOne();
00085 
00086   EndOfConstructorOps();
00087 }
00088 
00089 //==============================================================================
00090 // Epetra_BlockMap constructor for a user-defined linear distribution of constant size elements.
00091 Epetra_BlockMap::Epetra_BlockMap(int NumGlobalElements, int NumMyElements, 
00092          int ElementSize, int IndexBase, const Epetra_Comm& Comm)
00093   : Epetra_Object("Epetra::BlockMap"),
00094     BlockMapData_(0)
00095 {
00096   if (NumGlobalElements < -1) 
00097     throw ReportError("NumGlobalElements = " + toString(NumGlobalElements) + ".  Should be >= -1.", -1);
00098   if (NumMyElements < 0) 
00099     throw ReportError("NumMyElements = " + toString(NumMyElements) + ".  Should be >= 0.", -2);
00100   if (ElementSize <= 0) 
00101     throw ReportError("ElementSize = " + toString(ElementSize) + ". Should be > 0.", -3);
00102 
00103   BlockMapData_ = new Epetra_BlockMapData(NumGlobalElements, ElementSize, IndexBase, Comm);
00104   BlockMapData_->NumMyElements_ = NumMyElements;
00105   BlockMapData_->MinMyElementSize_ = BlockMapData_->ElementSize_;
00106   BlockMapData_->MaxMyElementSize_ = BlockMapData_->ElementSize_;
00107   BlockMapData_->MinElementSize_ = BlockMapData_->ElementSize_;
00108   BlockMapData_->MaxElementSize_ = BlockMapData_->ElementSize_;
00109   BlockMapData_->ConstantElementSize_ = true;
00110   BlockMapData_->LinearMap_ = true;
00111 
00112   // Each processor gets NumMyElements points
00113   
00114 
00115   // Get processor information
00116 
00117   int NumProc = Comm.NumProc();
00118 
00119   BlockMapData_->DistributedGlobal_ = IsDistributedGlobal(NumGlobalElements, NumMyElements);
00120 
00121   // Local Map and uniprocessor case:  Each processor gets a complete copy of all elements
00122   if (!BlockMapData_->DistributedGlobal_ || NumProc==1) {
00123     BlockMapData_->NumGlobalElements_ = BlockMapData_->NumMyElements_;
00124     CheckValidNGE(NumGlobalElements);
00125     
00126     BlockMapData_->NumGlobalPoints_ = BlockMapData_->NumGlobalElements_ * BlockMapData_->ElementSize_;
00127     BlockMapData_->NumMyPoints_ = BlockMapData_->NumMyElements_ * BlockMapData_->ElementSize_;
00128     
00129     BlockMapData_-> MinAllGID_ = BlockMapData_->IndexBase_;
00130     BlockMapData_->MaxAllGID_ = BlockMapData_->MinAllGID_ + BlockMapData_->NumGlobalElements_ - 1;
00131     BlockMapData_->MinMyGID_ = BlockMapData_->IndexBase_;
00132     BlockMapData_->MaxMyGID_ = BlockMapData_->MinMyGID_ + BlockMapData_->NumMyElements_ - 1;
00133   }
00134   else if (NumProc > 1) {
00135     // Sum up all local element counts to get global count
00136     BlockMapData_->Comm_->SumAll(&BlockMapData_->NumMyElements_, &BlockMapData_->NumGlobalElements_, 1);
00137     
00138     CheckValidNGE(NumGlobalElements);
00139     
00140     BlockMapData_->NumGlobalPoints_ = BlockMapData_->NumGlobalElements_ * BlockMapData_->ElementSize_;
00141     BlockMapData_->NumMyPoints_ = BlockMapData_->NumMyElements_ * BlockMapData_->ElementSize_;
00142     
00143     BlockMapData_->MinAllGID_ = BlockMapData_->IndexBase_;
00144     BlockMapData_->MaxAllGID_ = BlockMapData_->MinAllGID_ + BlockMapData_->NumGlobalElements_ - 1;
00145     
00146     // Use the ScanSum function to compute a prefix sum of the number of points
00147     BlockMapData_->Comm_->ScanSum(&BlockMapData_->NumMyElements_, &BlockMapData_->MaxMyGID_, 1);
00148     
00149     int start_index = BlockMapData_->MaxMyGID_ - BlockMapData_->NumMyElements_;
00150     BlockMapData_->MinMyGID_ = start_index + BlockMapData_->IndexBase_;
00151     BlockMapData_->MaxMyGID_ = BlockMapData_->MinMyGID_ + BlockMapData_->NumMyElements_ - 1;
00152   }
00153   else
00154     throw ReportError("Internal Error.  Report to Epetra developer", -99);
00155   
00156   BlockMapData_->OneToOne_ = DetermineIsOneToOne();
00157 
00158   EndOfConstructorOps();
00159 }
00160 
00161 //==============================================================================
00162 // Epetra_BlockMap constructor for a user-defined arbitrary distribution of constant size elements.
00163 Epetra_BlockMap::Epetra_BlockMap(int NumGlobalElements, int NumMyElements,
00164                                  const int * MyGlobalElements, 
00165          int ElementSize, int IndexBase,
00166                                  const Epetra_Comm& Comm)
00167   : Epetra_Object("Epetra::BlockMap"),
00168     BlockMapData_(0)
00169 {
00170   int i;
00171   // Each processor gets NumMyElements points
00172 
00173   if (NumGlobalElements < -1) 
00174     throw ReportError("NumGlobalElements = " + toString(NumGlobalElements) + ".  Should be >= -1.", -1);
00175   if (NumMyElements < 0) 
00176     throw ReportError("NumMyElements = " + toString(NumMyElements) + ".  Should be >= 0.", -2);
00177   if (ElementSize <= 0) 
00178     throw ReportError("ElementSize = " + toString(ElementSize) + ". Should be > 0.", -3);
00179 
00180   // Allocate storage for global index list information
00181 
00182   BlockMapData_ = new Epetra_BlockMapData(NumGlobalElements, ElementSize, IndexBase, Comm);
00183   if (NumMyElements > 0) {
00184     int errorcode = BlockMapData_->MyGlobalElements_.Size(NumMyElements);
00185     if(errorcode != 0)
00186       throw ReportError("Error with MyGlobalElements allocation.", -99);
00187   }
00188 
00189   BlockMapData_->NumMyElements_ = NumMyElements;
00190   BlockMapData_->MinMyElementSize_ = BlockMapData_->ElementSize_;
00191   BlockMapData_->MaxMyElementSize_ = BlockMapData_->ElementSize_;
00192   BlockMapData_->MinElementSize_ = BlockMapData_->ElementSize_;
00193   BlockMapData_->MaxElementSize_ = BlockMapData_->ElementSize_;
00194   BlockMapData_->ConstantElementSize_ = true;
00195   BlockMapData_->LinearMap_ = false;
00196   // Get processor information
00197 
00198   int NumProc = Comm.NumProc();
00199   if (NumMyElements > 0) {
00200     // Compute min/max GID on this processor
00201     BlockMapData_->MinMyGID_ = MyGlobalElements[0];
00202     BlockMapData_->MaxMyGID_ = MyGlobalElements[0];
00203     for (i = 0; i < NumMyElements; i++) {
00204       BlockMapData_->MyGlobalElements_[i] = MyGlobalElements[i];
00205       BlockMapData_->MinMyGID_ = EPETRA_MIN(BlockMapData_->MinMyGID_,MyGlobalElements[i]);
00206       BlockMapData_->MaxMyGID_ = EPETRA_MAX(BlockMapData_->MaxMyGID_,MyGlobalElements[i]);
00207     }
00208   }
00209   else {
00210     BlockMapData_->MinMyGID_ = BlockMapData_->IndexBase_;
00211     BlockMapData_->MaxMyGID_ = BlockMapData_->IndexBase_ - 1;
00212   }
00213   
00214   BlockMapData_->DistributedGlobal_ = IsDistributedGlobal(NumGlobalElements, NumMyElements);
00215 
00216   // Local Map and uniprocessor case:  Each processor gets a complete copy of all elements
00217   if (!BlockMapData_->DistributedGlobal_ || NumProc==1) {
00218     BlockMapData_->NumGlobalElements_ = BlockMapData_->NumMyElements_;
00219     CheckValidNGE(NumGlobalElements);
00220     BlockMapData_->NumGlobalPoints_ = BlockMapData_->NumGlobalElements_ * BlockMapData_->ElementSize_;
00221     BlockMapData_->NumMyPoints_ = BlockMapData_->NumMyElements_ * BlockMapData_->ElementSize_;
00222     
00223     BlockMapData_->MinAllGID_ = BlockMapData_->MinMyGID_;
00224     BlockMapData_->MaxAllGID_ = BlockMapData_->MaxMyGID_;
00225   }
00226   else if (NumProc > 1) {
00227     // Sum up all local element counts to get global count
00228     BlockMapData_->Comm_->SumAll(&BlockMapData_->NumMyElements_, &BlockMapData_->NumGlobalElements_, 1);
00229     CheckValidNGE(NumGlobalElements);
00230     
00231     BlockMapData_->NumGlobalPoints_ = BlockMapData_->NumGlobalElements_ * BlockMapData_->ElementSize_;
00232     BlockMapData_->NumMyPoints_ = BlockMapData_->NumMyElements_ * BlockMapData_->ElementSize_;
00233     
00234     // Use the Allreduce function to find min/max GID 
00235     int *tmp_send = new int[2];
00236     int *tmp_recv = new int[2];
00237     tmp_send[0] = - BlockMapData_->MinMyGID_; // Negative sign lets us do one reduction
00238     tmp_send[1] =   BlockMapData_->MaxMyGID_;
00239     BlockMapData_->Comm_->MaxAll(tmp_send, tmp_recv, 2);
00240     BlockMapData_->MinAllGID_ = - tmp_recv[0];
00241     BlockMapData_->MaxAllGID_ =   tmp_recv[1];
00242     delete [] tmp_send;
00243     delete [] tmp_recv;
00244     if (BlockMapData_->MinAllGID_ < BlockMapData_->IndexBase_)
00245       throw ReportError("Minimum global element index = " + toString(BlockMapData_->MinAllGID_) + 
00246       " is less than index base = " + toString(BlockMapData_->IndexBase_) +".", -5);
00247   }
00248   else
00249     throw ReportError("Internal Error.  Report to Epetra developer", -99);
00250   
00251   BlockMapData_->OneToOne_ = DetermineIsOneToOne();
00252 
00253   EndOfConstructorOps();
00254 }
00255 
00256 //==============================================================================
00257 // Epetra_BlockMap constructor for a user-defined arbitrary distribution of variable size elements.
00258 Epetra_BlockMap::Epetra_BlockMap(int NumGlobalElements, int NumMyElements,
00259                                  const int * MyGlobalElements, 
00260          const int *ElementSizeList, int IndexBase,
00261                                  const Epetra_Comm& Comm)
00262   : Epetra_Object("Epetra::BlockMap"),
00263     BlockMapData_(0)
00264 {
00265 
00266   int i;
00267   // Each processor gets NumMyElements points
00268 
00269   if (NumGlobalElements < -1) 
00270     throw ReportError("NumGlobalElements = " + toString(NumGlobalElements) + ".  Should be >= -1.", -1);
00271   if (NumMyElements < 0) 
00272     throw ReportError("NumMyElements = " + toString(NumMyElements) + ".  Should be >= 0.", -2);
00273   for (i = 0; i < NumMyElements; i++)
00274     if (ElementSizeList[i] <= 0) 
00275       throw ReportError("ElementSizeList["+toString(i)+"] = " + toString(ElementSizeList[i]) + ". Should be > 0.", -3);
00276   
00277   BlockMapData_ = new Epetra_BlockMapData(NumGlobalElements, 0, IndexBase, Comm);
00278   BlockMapData_->NumMyElements_ = NumMyElements;
00279   BlockMapData_->ConstantElementSize_ = false;
00280   BlockMapData_->LinearMap_ = false;
00281   // Allocate storage for global index list and element size information
00282 
00283   if (NumMyElements > 0) {
00284     int errorcode = BlockMapData_->MyGlobalElements_.Size(NumMyElements);
00285     if(errorcode != 0)
00286       throw ReportError("Error with MyGlobalElements allocation.", -99);
00287     errorcode = BlockMapData_->ElementSizeList_.Size(NumMyElements);
00288     if(errorcode != 0)
00289       throw ReportError("Error with ElementSizeList allocation.", -99);
00290   }
00291   // Get processor information
00292 
00293   int NumProc = Comm.NumProc();
00294   
00295   if (NumMyElements > 0) {
00296     // Compute min/max GID and element size, number of points on this processor
00297     BlockMapData_->MinMyGID_ = MyGlobalElements[0];
00298     BlockMapData_->MaxMyGID_ = MyGlobalElements[0];
00299     BlockMapData_->MinMyElementSize_ = ElementSizeList[0];
00300     BlockMapData_->MaxMyElementSize_ = ElementSizeList[0];
00301     BlockMapData_->NumMyPoints_ = 0;
00302     for (i = 0; i < NumMyElements; i++) {
00303       BlockMapData_->MyGlobalElements_[i] = MyGlobalElements[i];
00304       BlockMapData_->ElementSizeList_[i] = ElementSizeList[i];
00305       BlockMapData_->MinMyGID_ = EPETRA_MIN(BlockMapData_->MinMyGID_,MyGlobalElements[i]);
00306       BlockMapData_->MaxMyGID_ = EPETRA_MAX(BlockMapData_->MaxMyGID_,MyGlobalElements[i]);
00307       BlockMapData_->MinMyElementSize_ = EPETRA_MIN(BlockMapData_->MinMyElementSize_,ElementSizeList[i]);
00308       BlockMapData_->MaxMyElementSize_ = EPETRA_MAX(BlockMapData_->MaxMyElementSize_,ElementSizeList[i]);
00309       BlockMapData_->NumMyPoints_ += ElementSizeList[i];
00310     }
00311   }
00312   else {
00313     BlockMapData_->MinMyGID_ = BlockMapData_->IndexBase_;
00314     BlockMapData_->MaxMyGID_ = BlockMapData_->IndexBase_ - 1;
00315     BlockMapData_->MinMyElementSize_ = 1;
00316     BlockMapData_->MaxMyElementSize_ = 1;
00317     BlockMapData_->NumMyPoints_ = 0;
00318   }
00319 
00320   BlockMapData_->DistributedGlobal_ = IsDistributedGlobal(NumGlobalElements, NumMyElements);  
00321 
00322   // Local Map and uniprocessor case:  Each processor gets a complete copy of all elements
00323   if (!BlockMapData_->DistributedGlobal_ || NumProc == 1) {
00324     BlockMapData_->NumGlobalElements_ = BlockMapData_->NumMyElements_;
00325     CheckValidNGE(NumGlobalElements);
00326     BlockMapData_->NumGlobalPoints_ = BlockMapData_->NumMyPoints_;
00327     
00328     BlockMapData_->MinAllGID_ = BlockMapData_->MinMyGID_;
00329     BlockMapData_->MaxAllGID_ = BlockMapData_->MaxMyGID_;
00330     BlockMapData_->MinElementSize_ = BlockMapData_->MinMyElementSize_;
00331     BlockMapData_->MaxElementSize_ = BlockMapData_->MaxMyElementSize_;
00332   }
00333   else if (NumProc > 1) {
00334     // Sum up all local element and point counts to get global counts
00335     int *tmp_send = new int[4];
00336     int *tmp_recv = new int[4];
00337     tmp_send[0] = BlockMapData_->NumMyElements_;
00338     tmp_send[1] = BlockMapData_->NumMyPoints_;
00339     BlockMapData_->Comm_->SumAll(tmp_send, tmp_recv, 2);
00340     BlockMapData_->NumGlobalElements_ =  tmp_recv[0];
00341     BlockMapData_->NumGlobalPoints_ = tmp_recv[1];
00342     
00343     CheckValidNGE(NumGlobalElements);
00344     
00345     // Use the MaxAll function to find min/max GID 
00346     tmp_send[0] = - BlockMapData_->MinMyGID_; // Negative signs lets us do one reduction
00347     tmp_send[1] =   BlockMapData_->MaxMyGID_;
00348     tmp_send[2] = - BlockMapData_->MinMyElementSize_;
00349     if (BlockMapData_->NumMyElements_ == 0) 
00350       tmp_send[2] = - BlockMapData_->NumGlobalPoints_; // This processor has no elements, so should not sizes.
00351     tmp_send[3] =   BlockMapData_->MaxMyElementSize_;
00352     
00353     BlockMapData_->Comm_->MaxAll(tmp_send, tmp_recv, 4);
00354     
00355     BlockMapData_->MinAllGID_ =      - tmp_recv[0];
00356     BlockMapData_->MaxAllGID_ =        tmp_recv[1];
00357     BlockMapData_->MinElementSize_ = - tmp_recv[2];
00358     BlockMapData_->MaxElementSize_ =   tmp_recv[3];
00359     
00360     delete [] tmp_send;
00361     delete [] tmp_recv;
00362 
00363     // Check for constant element size
00364     if (BlockMapData_->MinElementSize_==BlockMapData_->MaxElementSize_) {
00365       BlockMapData_->ElementSize_ = BlockMapData_->MinElementSize_;
00366       BlockMapData_->ConstantElementSize_ = true;
00367     }
00368     
00369     if (BlockMapData_->MinAllGID_ < BlockMapData_->IndexBase_)
00370       throw ReportError("Minimum global element index = " + toString(BlockMapData_->MinAllGID_) + 
00371       " is less than index base = " + toString(BlockMapData_->IndexBase_) +".", -5);
00372   }
00373   else
00374     throw ReportError("Internal Error.  Report to Epetra developer", -99);
00375   
00376   BlockMapData_->OneToOne_ = DetermineIsOneToOne();
00377 
00378   EndOfConstructorOps();
00379 }
00380 
00381 //==============================================================================
00382 Epetra_BlockMap::Epetra_BlockMap(const Epetra_BlockMap& map)
00383   : Epetra_Object(map.Label()),
00384     BlockMapData_(map.BlockMapData_)
00385 {
00386   BlockMapData_->IncrementReferenceCount();
00387   
00388   // This call appears to be unnecessary overhead.  Removed 10-Aug-2004 maherou.
00389   // GlobalToLocalSetup(); // Setup any information for making global index to local index translation fast.
00390 }
00391 
00392 //==============================================================================
00393 bool Epetra_BlockMap::SameAs(const Epetra_BlockMap & Map) const {
00394 
00395   // Quickest test: See if both maps share an inner data class
00396   if (this->BlockMapData_ == Map.BlockMapData_) 
00397     return(true);
00398 
00399 
00400   // Next check other global properties that are easy global attributes
00401   if (BlockMapData_->MinAllGID_ != Map.MinAllGID() ||
00402       BlockMapData_->MaxAllGID_ != Map.MaxAllGID() ||
00403       BlockMapData_->NumGlobalElements_ != Map.NumGlobalElements() ||
00404       BlockMapData_->IndexBase_ != Map.IndexBase()) 
00405     return(false);
00406   
00407   // Last possible global check for constant element sizes
00408   if (BlockMapData_->ConstantElementSize_ && BlockMapData_->ElementSize_!=Map.ElementSize()) 
00409     return(false);
00410 
00411   // If we get this far, we need to check local properties and then check across
00412   // all processors to see if local properties are all true
00413  
00414   int numMyElements = BlockMapData_->NumMyElements_;
00415 
00416   int MySameMap = 1; // Assume not needed
00417   
00418   // First check if number of element is the same in each map
00419   if (numMyElements != Map.NumMyElements()) MySameMap = 0;
00420   
00421   if (MySameMap==1) // If numMyElements is the same, check to see that list of GIDs is the same
00422     for (int i = 0; i < numMyElements; i++)
00423       if (GID(i) != Map.GID(i)) MySameMap = 0;
00424 
00425   // If GIDs are the same, check to see element sizes are the same
00426   if (MySameMap==1 && !BlockMapData_->ConstantElementSize_) {
00427     int * sizeList1 = ElementSizeList();
00428     int * sizeList2 = Map.ElementSizeList();
00429     for (int i = 0; i < numMyElements; i++) if (sizeList1[i] != sizeList2[i]) MySameMap=0;
00430   }
00431   // Now get min of MySameMap across all processors
00432 
00433   int GlobalSameMap = 0;
00434   int err = Comm().MinAll(&MySameMap, &GlobalSameMap, 1);
00435   assert(err==0);
00436 
00437   return(GlobalSameMap==1);
00438 }
00439 
00440 //==============================================================================
00441 bool Epetra_BlockMap::PointSameAs(const Epetra_BlockMap & Map) const
00442 {
00443   if (this->BlockMapData_ == Map.BlockMapData_) 
00444     return(true);
00445   
00446   if (BlockMapData_->NumGlobalPoints_ != Map.NumGlobalPoints() ) 
00447     return(false);
00448   
00449   // If we get this far, we need to check local properties and then check across
00450   // all processors to see if local properties are all true
00451 
00452   int MySameMap = 1; // Assume not needed
00453   if (BlockMapData_->NumMyPoints_ != Map.NumMyPoints())
00454     MySameMap = 0;
00455 
00456   int GlobalSameMap = 0;
00457   int err = Comm().MinAll(&MySameMap, &GlobalSameMap, 1);
00458   assert( err == 0 );
00459 
00460   return(GlobalSameMap==1);
00461 }
00462 
00463 //==============================================================================
00464 int Epetra_BlockMap::MyGlobalElements(int * MyGlobalElements) const
00465 {
00466   // If the global element list is not create, then do so.  This can only happen when
00467   // a linear distribution has been specified.  Thus we can easily construct the update
00468   // list in this case.
00469 
00470   int i;
00471   int numMyElements = BlockMapData_->NumMyElements_;
00472   
00473   if (BlockMapData_->MyGlobalElements_.Length() == 0)
00474     for (i = 0; i < numMyElements; i++)
00475       MyGlobalElements[i] = BlockMapData_->MinMyGID_ + i;
00476   else
00477     for (i = 0; i < numMyElements; i++)
00478       MyGlobalElements[i] = BlockMapData_->MyGlobalElements_[i];
00479   return(0);
00480 }
00481 
00482 //==============================================================================
00483 int * Epetra_BlockMap::MyGlobalElements() const {
00484   int numMyElements = BlockMapData_->NumMyElements_;  
00485 
00486   // If ElementSizeList not built, do so
00487   if(BlockMapData_->MyGlobalElements_.Length() == 0 && numMyElements > 0) {
00488     int errorcode = BlockMapData_->MyGlobalElements_.Size(numMyElements + 1);
00489     if(errorcode != 0)
00490       throw ReportError("Error with MyGlobalElements allocation.", -99);
00491     
00492     // Build the array
00493     for (int i = 0; i < numMyElements; i++)
00494       BlockMapData_->MyGlobalElements_[i] = BlockMapData_->MinMyGID_ + i;
00495   }
00496   return(BlockMapData_->MyGlobalElements_.Values());
00497 }
00498 
00499 //==============================================================================
00500 int Epetra_BlockMap::FirstPointInElement(int LID) const
00501 {
00502   if (!MyLID(LID)) 
00503     EPETRA_CHK_ERR(-1);
00504   
00505   int entry;
00506 
00507   if (ConstantElementSize())
00508     entry = MaxElementSize() * LID; // convert to vector entry
00509   else {
00510     int * entrylist = FirstPointInElementList(); // get entry list
00511     entry = entrylist[LID];
00512   }
00513   return(entry);
00514 }
00515 
00516 //==============================================================================
00517 int Epetra_BlockMap::FirstPointInElementList(int * FirstPointInElementList) const
00518 {
00519   // If the first element entry list is not create, then do so.  
00520 
00521   // Note: This array is of length NumMyElement+1
00522 
00523   int i;
00524   int numMyElements = BlockMapData_->NumMyElements_;
00525 
00526   if (BlockMapData_->FirstPointInElementList_.Length() == 0) {
00527     FirstPointInElementList[0] = 0; // First element of first entry is always zero
00528     
00529     if (BlockMapData_->ConstantElementSize_)
00530       for (i = 0; i < numMyElements; i++)
00531   FirstPointInElementList[i+1] = FirstPointInElementList[i] + BlockMapData_->ElementSize_;
00532     else
00533       for (i = 0; i < numMyElements; i++)
00534   FirstPointInElementList[i+1] = FirstPointInElementList[i] + BlockMapData_->ElementSizeList_[i];
00535   }
00536   else 
00537     for (i = 0; i <= numMyElements; i++)
00538       FirstPointInElementList[i] = BlockMapData_->FirstPointInElementList_[i];
00539   return(0);
00540 }
00541 
00542 //==============================================================================
00543 int * Epetra_BlockMap::FirstPointInElementList() const {
00544   int numMyElements = BlockMapData_->NumMyElements_;
00545 
00546   // If ElementSizeList not built, do so
00547   if ((BlockMapData_->FirstPointInElementList_.Length() == 0) && (numMyElements > 0)) {
00548     BlockMapData_->FirstPointInElementList_.Size(BlockMapData_->NumMyElements_ + 1);
00549     BlockMapData_->FirstPointInElementList_[0] = 0; // First element of first entry is always zero
00550     if (BlockMapData_->ConstantElementSize_)
00551       for (int i = 0; i < numMyElements; i++)
00552   BlockMapData_->FirstPointInElementList_[i+1] = BlockMapData_->FirstPointInElementList_[i] + BlockMapData_->ElementSize_;
00553     else
00554       for (int i = 0; i < numMyElements; i++)
00555   BlockMapData_->FirstPointInElementList_[i+1] = BlockMapData_->FirstPointInElementList_[i] + BlockMapData_->ElementSizeList_[i];
00556   }
00557   return(BlockMapData_->FirstPointInElementList_.Values());
00558 }
00559 
00560 //==============================================================================
00561 int Epetra_BlockMap::ElementSizeList(int * ElementSizeList) const
00562 {
00563   // If the element size list is not create, then do so.  This can only happen when
00564   // a constant element size has been specified.  Thus we can easily construct the element size
00565   // list in this case.
00566 
00567   int i;
00568   int numMyElements = BlockMapData_->NumMyElements_;
00569 
00570   if (BlockMapData_->ElementSizeList_.Length() == 0)
00571     for (i = 0; i < numMyElements; i++)
00572       ElementSizeList[i] = BlockMapData_->ElementSize_;
00573   else
00574     for (i = 0; i < numMyElements; i++)
00575       ElementSizeList[i] = BlockMapData_->ElementSizeList_[i];
00576   
00577   return(0);
00578 }
00579 
00580 //==============================================================================
00581 int * Epetra_BlockMap::ElementSizeList() const {
00582   int numMyElements = BlockMapData_->NumMyElements_;
00583 
00584   // If ElementSizeList not built, do so
00585   if ((BlockMapData_->ElementSizeList_.Length() == 0) && (numMyElements > 0)) {
00586     BlockMapData_->ElementSizeList_.Size(numMyElements);
00587     for (int i = 0; i < numMyElements; i++)
00588       BlockMapData_->ElementSizeList_[i] = BlockMapData_->ElementSize_;
00589   }
00590   return(BlockMapData_->ElementSizeList_.Values());
00591 }
00592 
00593 //==============================================================================
00594 int Epetra_BlockMap::PointToElementList(int * PointToElementList) const {
00595   // Build an array such that the local element ID is stored for each point
00596 
00597   int i;
00598   if (BlockMapData_->PointToElementList_.Length() == 0) {
00599     int numMyElements = BlockMapData_->NumMyElements_;
00600     int * ptr = PointToElementList;
00601     for (i = 0; i < numMyElements; i++) {
00602       int Size = ElementSize(i);
00603       for (int j = 0; j < Size; j++) 
00604   *ptr++ = i;
00605     }
00606   }
00607   else {
00608     int numMyPoints = BlockMapData_->NumMyPoints_;
00609     for (i = 0; i < numMyPoints; i++)
00610       PointToElementList[i] = BlockMapData_->PointToElementList_[i];
00611   }
00612   return(0);
00613 }
00614 
00615 //==============================================================================
00616 int * Epetra_BlockMap::PointToElementList() const {
00617 
00618   // If PointToElementList not built, do so
00619   if ((BlockMapData_->PointToElementList_.Length() == 0) && (BlockMapData_->NumMyPoints_ > 0)) {
00620     BlockMapData_->PointToElementList_.Size(BlockMapData_->NumMyPoints_);
00621     int numMyElements = BlockMapData_->NumMyElements_;
00622     int * ptr = BlockMapData_->PointToElementList_.Values();
00623     for (int i = 0; i < numMyElements; i++) {
00624       int Size = ElementSize(i);
00625       for (int j = 0; j < Size; j++) 
00626   *ptr++ = i;
00627     }
00628   }
00629   return(BlockMapData_->PointToElementList_.Values());
00630 }
00631 
00632 //==============================================================================
00633 int Epetra_BlockMap::ElementSize(int LID) const {
00634 
00635   if (ConstantElementSize()) 
00636     return(BlockMapData_->ElementSize_);
00637   else
00638     return(BlockMapData_->ElementSizeList_[LID]);
00639 }
00640 
00641 //==============================================================================
00642 void Epetra_BlockMap::GlobalToLocalSetup()
00643 {
00644   int i;
00645   int numMyElements = BlockMapData_->NumMyElements_;
00646 
00647   if (BlockMapData_->NumGlobalElements_ == 0) {
00648     return; // Nothing to do
00649   }
00650 
00651   if (LinearMap() || numMyElements == 0) {
00652     return; // Nothing else to do
00653   }
00654 
00655   // Build LID_ vector to make look up of local index values fast
00656 
00657 #ifdef EPETRA_BLOCKMAP_NEW_LID
00658 
00659   //check for initial contiguous block
00660   int val = BlockMapData_->MyGlobalElements_[0];
00661   for( i = 0 ; i < numMyElements; ++i ) {
00662     if (val != BlockMapData_->MyGlobalElements_[i]) break;
00663     ++val;
00664   }
00665   BlockMapData_->LastContiguousGIDLoc_ = i - 1;
00666   if (BlockMapData_->LastContiguousGIDLoc_ < 0) {
00667     BlockMapData_->LastContiguousGID_ = BlockMapData_->MyGlobalElements_[0];
00668   }
00669   else {
00670     BlockMapData_->LastContiguousGID_ =
00671       BlockMapData_->MyGlobalElements_[BlockMapData_->LastContiguousGIDLoc_];
00672   }
00673 
00674   //Hash everything else
00675   if(i < numMyElements) {
00676     if (BlockMapData_->LIDHash_ != NULL) {
00677       delete BlockMapData_->LIDHash_;
00678     }
00679 
00680     BlockMapData_->LIDHash_ = new Epetra_HashTable(numMyElements - i + 1 );
00681     for(; i < numMyElements; ++i )
00682       BlockMapData_->LIDHash_->Add( BlockMapData_->MyGlobalElements_[i], i );
00683   }
00684     
00685 #else
00686     
00687   int SpanGID = BlockMapData_->MaxMyGID_ - BlockMapData_->MinMyGID_ + 1;
00688   BlockMapData_->LID_.Size(SpanGID);
00689     
00690   for (i = 0; i < SpanGID; i++) 
00691     BlockMapData_->LID_[i] = -1; // Fill all locations with -1
00692     
00693   for (i = 0; i < numMyElements; i++) {
00694     int tmp = BlockMapData_->MyGlobalElements_[i] - BlockMapData_->MinMyGID_;
00695     assert(tmp >= 0); 
00696     assert(tmp < SpanGID);
00697     BlockMapData_->LID_[BlockMapData_->MyGlobalElements_[i] - BlockMapData_->MinMyGID_] = i; // Spread local indices
00698   }
00699 
00700 #endif
00701 
00702 }
00703 
00704 //==============================================================================
00705 int Epetra_BlockMap::LID(int GID) const
00706 {
00707   if ((GID < BlockMapData_->MinMyGID_) || 
00708       (GID > BlockMapData_->MaxMyGID_)) {
00709     return(-1); // Out of range
00710   }
00711 
00712   if (BlockMapData_->LinearMap_) {
00713     return(GID - BlockMapData_->MinMyGID_); // Can compute with an offset
00714   }
00715 
00716   if( GID >= BlockMapData_->MyGlobalElements_[0] &&
00717       GID <= BlockMapData_->LastContiguousGID_ ) {
00718     return( GID - BlockMapData_->MyGlobalElements_[0] );
00719   }
00720 
00721 #ifdef EPETRA_BLOCKMAP_NEW_LID
00722   return BlockMapData_->LIDHash_->Get( GID );
00723 #else
00724   return(BlockMapData_->LID_[GID - BlockMapData_->MinMyGID_]); // Find it in LID array  
00725 #endif
00726 }
00727 
00728 //==============================================================================
00729 int Epetra_BlockMap::GID(int LID) const
00730 {
00731   if ((BlockMapData_->NumMyElements_==0) ||
00732       (LID < BlockMapData_->MinLID_) || 
00733       (LID > BlockMapData_->MaxLID_)) {
00734     return(BlockMapData_->IndexBase_ - 1); // Out of range
00735   }
00736 
00737   if (LinearMap()) {
00738     return(LID + BlockMapData_->MinMyGID_); // Can compute with an offset
00739   }
00740 
00741   return(BlockMapData_->MyGlobalElements_[LID]); // Find it in MyGlobalElements array
00742 }
00743 
00744 //==============================================================================
00745 int Epetra_BlockMap::FindLocalElementID(int PointID, int & ElementID, int & ElementOffset) const {
00746 
00747   if (PointID >= BlockMapData_->NumMyPoints_)
00748     return(-1); // Point is out of range
00749   
00750   if (ConstantElementSize()) {
00751     ElementID = PointID / BlockMapData_->MaxElementSize_;
00752     ElementOffset = PointID % BlockMapData_->MaxElementSize_;
00753     return(0);
00754   }
00755   else {
00756     int * tmpPointToElementList = PointToElementList();
00757     int * tmpFirstPointInElementList = FirstPointInElementList();
00758     ElementID = tmpPointToElementList[PointID];
00759     ElementOffset = PointID - tmpFirstPointInElementList[ElementID];
00760     return(0);
00761   }
00762 }
00763 
00764 //==============================================================================
00765 int Epetra_BlockMap::RemoteIDList(int NumIDs, const int * GIDList,
00766           int * PIDList, int * LIDList,
00767           int * SizeList) const
00768 {
00769   if (BlockMapData_->Directory_ == NULL) {
00770     BlockMapData_->Directory_ = Comm().CreateDirectory(*this);
00771   }
00772 
00773   Epetra_Directory* directory = BlockMapData_->Directory_;
00774   if (directory == NULL) {
00775     return(-1);
00776   }
00777 
00778   EPETRA_CHK_ERR( directory->GetDirectoryEntries(*this, NumIDs, GIDList,
00779              PIDList, LIDList, SizeList) );
00780 
00781   return(0);
00782 }
00783 
00784 //==============================================================================
00785 bool Epetra_BlockMap::DetermineIsOneToOne()
00786 {
00787   if (Comm().NumProc() < 2) {
00788     return(true);
00789   }
00790   
00791   if (BlockMapData_->Directory_ == NULL) {
00792     BlockMapData_->Directory_ = Comm().CreateDirectory(*this);
00793   }
00794 
00795   Epetra_Directory* directory = BlockMapData_->Directory_;
00796   if (directory == NULL) {
00797     throw ReportError("Epetra_BlockMap::IsOneToOne ERROR, CreateDirectory failed.",-1);
00798   }
00799 
00800   return(directory->GIDsAllUniquelyOwned());
00801 }
00802 
00803 //==============================================================================
00804 bool Epetra_BlockMap::IsDistributedGlobal(int NumGlobalElements, int NumMyElements) const {
00805 
00806   bool DistributedGlobal = false; // Assume map is not global distributed
00807   if (BlockMapData_->Comm_->NumProc() > 1) {
00808     int LocalReplicated = 0;
00809     int AllLocalReplicated;
00810     if (NumGlobalElements == NumMyElements) 
00811       LocalReplicated=1;
00812     BlockMapData_->Comm_->MinAll(&LocalReplicated, &AllLocalReplicated, 1);
00813     
00814     // If any PE has LocalReplicated=0, then map is distributed global
00815     if (AllLocalReplicated != 1) 
00816       DistributedGlobal = true;
00817   }
00818   return(DistributedGlobal);
00819 }
00820 
00821 //==============================================================================
00822 void Epetra_BlockMap::CheckValidNGE(int NumGlobalElements) {
00823   // Check to see if user's value for NumGlobalElements is either -1 
00824   // (in which case we use our computed value) or matches ours.
00825   if ((NumGlobalElements != -1) && (NumGlobalElements != BlockMapData_->NumGlobalElements_)) {
00826     int BmdNumGlobalElements = BlockMapData_->NumGlobalElements_;
00827     CleanupData();
00828     throw ReportError("Invalid NumGlobalElements.  NumGlobalElements = " + toString(NumGlobalElements) + 
00829           ".  Should equal " + toString(BmdNumGlobalElements) + 
00830           ", or be set to -1 to compute automatically", -4);
00831   }
00832 }
00833 
00834 //==============================================================================
00835 void Epetra_BlockMap::EndOfConstructorOps() {
00836   BlockMapData_->MinLID_ = 0;
00837   BlockMapData_->MaxLID_ = EPETRA_MAX(BlockMapData_->NumMyElements_ - 1, 0);
00838   
00839   GlobalToLocalSetup(); // Setup any information for making global index to local index translation fast.
00840 }
00841 
00842 //==============================================================================
00843 void Epetra_BlockMap::Print(ostream & os) const
00844 {
00845   int * MyGlobalElements1 = MyGlobalElements();
00846   int * FirstPointInElementList1 = 0;
00847   int * ElementSizeList1 = 0;
00848   if (!ConstantElementSize()) {
00849     FirstPointInElementList1 = FirstPointInElementList();
00850     ElementSizeList1 = ElementSizeList();
00851   }
00852   int MyPID = Comm().MyPID();
00853   int NumProc = Comm().NumProc();
00854   
00855   for (int iproc = 0; iproc < NumProc; iproc++) {
00856     if (MyPID == iproc) {
00857       if (MyPID == 0) {
00858   os <<  "\nNumber of Global Elements  = "; os << NumGlobalElements(); os << endl;
00859   os <<    "Number of Global Points    = "; os << NumGlobalPoints(); os << endl;
00860   os <<    "Maximum of all GIDs        = "; os << MaxAllGID(); os << endl;
00861   os <<    "Minimum of all GIDs        = "; os << MinAllGID(); os << endl;
00862   os <<    "Index Base                 = "; os << IndexBase(); os << endl;
00863   if (ConstantElementSize())
00864     os <<  "Constant Element Size      = "; os << ElementSize(); os << endl;
00865       }
00866       os << endl;
00867       
00868       os <<    "Number of Local Elements   = "; os << NumMyElements(); os << endl;
00869       os <<    "Number of Local Points     = "; os << NumMyPoints(); os << endl;
00870       os <<    "Maximum of my GIDs         = "; os << MaxMyGID(); os << endl;
00871       os <<    "Minimum of my GIDs         = "; os << MinMyGID(); os << endl;
00872       os << endl;
00873       
00874       os.width(14);
00875       os <<  "     MyPID"; os << "    ";
00876       os.width(14);
00877       os <<  "       Local Index "; os << " ";
00878       os.width(14);
00879       os <<  "      Global Index "; os << " ";
00880       if (!ConstantElementSize()) {
00881   os.width(14);
00882   os <<" FirstPointInElement "; os << " ";
00883   os.width(14);
00884   os <<"   ElementSize "; os << " ";
00885       }
00886       os << endl;
00887       
00888       for (int i = 0; i < NumMyElements(); i++) {
00889   os.width(14);
00890   os <<  MyPID; os << "    ";
00891   os.width(14);
00892   os <<  i; os << "    ";
00893   os.width(14);
00894   os <<  MyGlobalElements1[i]; os << "    ";
00895   if (!ConstantElementSize()) {   
00896     os.width(14);
00897     os << FirstPointInElementList1[i]; os << "    ";
00898     os.width(14);
00899     os << ElementSizeList1[i]; os << "    ";
00900   }
00901   os << endl;
00902       }
00903       
00904       os << flush;
00905       
00906     }
00907     // Do a few global ops to give I/O a chance to complete
00908     Comm().Barrier();
00909     Comm().Barrier();
00910     Comm().Barrier();
00911   }
00912   return;
00913 }
00914 
00915 //==============================================================================
00916 Epetra_BlockMap::~Epetra_BlockMap()  {
00917   CleanupData();
00918 }
00919 
00920 //==============================================================================
00921 void Epetra_BlockMap::CleanupData()
00922 {
00923   if(BlockMapData_ != 0) {
00924 
00925     BlockMapData_->DecrementReferenceCount();
00926     if(BlockMapData_->ReferenceCount() == 0) {
00927       delete BlockMapData_;
00928       BlockMapData_ = 0;
00929     }
00930   }
00931 }
00932 
00933 //=============================================================================
00934 Epetra_BlockMap & Epetra_BlockMap::operator= (const Epetra_BlockMap & map)
00935 {
00936   if((this != &map) && (BlockMapData_ != map.BlockMapData_)) {
00937     CleanupData();
00938     BlockMapData_ = map.BlockMapData_;
00939     BlockMapData_->IncrementReferenceCount();
00940   }
00941 
00942   return(*this);
00943 }

Generated on Wed May 12 21:41:04 2010 for Epetra Package Browser (Single Doxygen Collection) by  doxygen 1.4.7