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