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