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_ConfigDefs.h"
00044 #include "Epetra_BlockMap.h"
00045 #include "Epetra_Comm.h"
00046 #include "Epetra_Directory.h"
00047 #include "Epetra_IntSerialDenseVector.h"
00048 #include "Epetra_HashTable.h"
00049 #include <limits>
00050 
00051 #ifdef HAVE_MPI
00052 #include "Epetra_MpiComm.h"
00053 #endif
00054 
00055 // Use the new LID hash table approach by default
00056 #define EPETRA_BLOCKMAP_NEW_LID
00057 
00058 //==============================================================================
00059 // Epetra_BlockMap constructor function for a Epetra-defined uniform linear distribution of constant size elements.
00060 void Epetra_BlockMap::ConstructAutoUniform(long long NumGlobal_Elements, int Element_Size, long long Index_Base, const Epetra_Comm& comm, bool IsLongLong)
00061 {
00062 
00063   // Each processor gets roughly numGlobalPoints/p points
00064   // This routine automatically defines a linear partitioning of a
00065   // map with numGlobalPoints across the processors
00066   // specified in the given Epetra_Comm
00067 
00068   if (NumGlobal_Elements < 0)
00069     throw ReportError("NumGlobal_Elements = " + toString(NumGlobal_Elements) + ".  Should be >= 0.", -1);
00070   if (Element_Size <= 0)
00071     throw ReportError("ElementSize = " + toString(Element_Size) + ".  Should be > 0.", -2);
00072 
00073   BlockMapData_ = new Epetra_BlockMapData(NumGlobal_Elements, Element_Size, Index_Base, comm, IsLongLong);
00074   int NumProc = comm.NumProc();
00075   BlockMapData_->ConstantElementSize_ = true;
00076   BlockMapData_->LinearMap_ = true;
00077 
00078   int MyPID = comm.MyPID();
00079 
00080   if(BlockMapData_->NumGlobalElements_ / NumProc > (long long) std::numeric_limits<int>::max())
00081     throw ReportError("Epetra_BlockMap::ConstructAutoUniform: Error. Not enough space for elements on each processor", -99);
00082 
00083   BlockMapData_->NumMyElements_ = (int) (BlockMapData_->NumGlobalElements_ / NumProc);
00084   int remainder = (int) (BlockMapData_->NumGlobalElements_ % NumProc); // remainder will fit int
00085   long long start_index = ((long long)(MyPID)) * (BlockMapData_->NumMyElements_ + 1);
00086 
00087   if (MyPID < remainder)
00088     BlockMapData_->NumMyElements_++;
00089   else
00090     start_index -= (MyPID - remainder);
00091 
00092   BlockMapData_->NumGlobalPoints_ = BlockMapData_->NumGlobalElements_ * BlockMapData_->ElementSize_;
00093   BlockMapData_->NumMyPoints_ = BlockMapData_->NumMyElements_ * BlockMapData_->ElementSize_;
00094 
00095   BlockMapData_->MinMyElementSize_ = BlockMapData_->ElementSize_;
00096   BlockMapData_->MaxMyElementSize_ = BlockMapData_->ElementSize_;
00097   BlockMapData_->MinElementSize_ = BlockMapData_->ElementSize_;
00098   BlockMapData_->MaxElementSize_ = BlockMapData_->ElementSize_;
00099 
00100   BlockMapData_->MinAllGID_ = BlockMapData_->IndexBase_;
00101   BlockMapData_->MaxAllGID_ = BlockMapData_->MinAllGID_ + BlockMapData_->NumGlobalElements_ - 1;
00102   BlockMapData_->MinMyGID_ = start_index + BlockMapData_->IndexBase_;
00103   BlockMapData_->MaxMyGID_ = BlockMapData_->MinMyGID_ + BlockMapData_->NumMyElements_ - 1;
00104   BlockMapData_->DistributedGlobal_ = IsDistributedGlobal(BlockMapData_->NumGlobalElements_, BlockMapData_->NumMyElements_);
00105 
00106   EndOfConstructorOps();
00107 }
00108 
00109 //==============================================================================
00110 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00111 Epetra_BlockMap::Epetra_BlockMap(long long NumGlobal_Elements, int Element_Size, int Index_Base, const Epetra_Comm& comm)
00112   : Epetra_Object("Epetra::BlockMap"),
00113     BlockMapData_(0)
00114 {
00115   const bool IsLongLong = true;
00116   ConstructAutoUniform(NumGlobal_Elements, Element_Size, static_cast<long long>(Index_Base), comm, IsLongLong);
00117 }
00118 
00119 Epetra_BlockMap::Epetra_BlockMap(long long NumGlobal_Elements, int Element_Size, long long Index_Base, const Epetra_Comm& comm)
00120   : Epetra_Object("Epetra::BlockMap"),
00121     BlockMapData_(0)
00122 {
00123   const bool IsLongLong = true;
00124   ConstructAutoUniform(NumGlobal_Elements, Element_Size, Index_Base, comm, IsLongLong);
00125 }
00126 #endif
00127 
00128 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00129 Epetra_BlockMap::Epetra_BlockMap(int NumGlobal_Elements, int Element_Size, int Index_Base, const Epetra_Comm& comm)
00130   : Epetra_Object("Epetra::BlockMap"),
00131     BlockMapData_(0)
00132 {
00133   const bool IsLongLong = false;
00134   ConstructAutoUniform((long long)NumGlobal_Elements, Element_Size, Index_Base, comm, IsLongLong);
00135 }
00136 #endif
00137 //==============================================================================
00138 
00139 // Epetra_BlockMap constructor function for a user-defined linear distribution of constant size elements.
00140 void Epetra_BlockMap::ConstructUserLinear(
00141     long long NumGlobal_Elements, int NumMy_Elements,
00142     int Element_Size, long long Index_Base, const Epetra_Comm& comm, bool IsLongLong)
00143 {
00144   if (NumGlobal_Elements < -1)
00145     throw ReportError("NumGlobal_Elements = " + toString(NumGlobal_Elements) + ".  Should be >= -1.", -1);
00146   if (NumMy_Elements < 0)
00147     throw ReportError("NumMy_Elements = " + toString(NumMy_Elements) + ".  Should be >= 0.", -2);
00148   if (Element_Size <= 0)
00149     throw ReportError("ElementSize = " + toString(Element_Size) + ". Should be > 0.", -3);
00150 
00151   BlockMapData_ = new Epetra_BlockMapData(NumGlobal_Elements, Element_Size, Index_Base, comm, IsLongLong);
00152   BlockMapData_->NumMyElements_ = NumMy_Elements;
00153   BlockMapData_->MinMyElementSize_ = BlockMapData_->ElementSize_;
00154   BlockMapData_->MaxMyElementSize_ = BlockMapData_->ElementSize_;
00155   BlockMapData_->MinElementSize_ = BlockMapData_->ElementSize_;
00156   BlockMapData_->MaxElementSize_ = BlockMapData_->ElementSize_;
00157   BlockMapData_->ConstantElementSize_ = true;
00158   BlockMapData_->LinearMap_ = true;
00159 
00160   // Each processor gets NumMyElements points
00161 
00162 
00163   // Get processor information
00164 
00165   int NumProc = comm.NumProc();
00166 
00167   BlockMapData_->DistributedGlobal_ = IsDistributedGlobal(NumGlobal_Elements, NumMy_Elements);
00168 
00169   // Local Map and uniprocessor case:  Each processor gets a complete copy of all elements
00170   if (!BlockMapData_->DistributedGlobal_ || NumProc==1) {
00171     BlockMapData_->NumGlobalElements_ = BlockMapData_->NumMyElements_;
00172     CheckValidNGE(NumGlobal_Elements);
00173 
00174     BlockMapData_->NumGlobalPoints_ = BlockMapData_->NumGlobalElements_ * BlockMapData_->ElementSize_;
00175     BlockMapData_->NumMyPoints_ = BlockMapData_->NumMyElements_ * BlockMapData_->ElementSize_;
00176 
00177     BlockMapData_->  MinAllGID_ = BlockMapData_->IndexBase_;
00178     BlockMapData_->MaxAllGID_ = BlockMapData_->MinAllGID_ + BlockMapData_->NumGlobalElements_ - 1;
00179     BlockMapData_->MinMyGID_ = BlockMapData_->IndexBase_;
00180     BlockMapData_->MaxMyGID_ = BlockMapData_->MinMyGID_ + BlockMapData_->NumMyElements_ - 1;
00181   }
00182   else if (NumProc > 1) {
00183     // Sum up all local element counts to get global count
00184     long long tmp_NumMyElements = BlockMapData_->NumMyElements_;
00185     BlockMapData_->Comm_->SumAll(&tmp_NumMyElements, &BlockMapData_->NumGlobalElements_, 1);
00186 
00187     CheckValidNGE(NumGlobal_Elements);
00188 
00189     BlockMapData_->NumGlobalPoints_ = BlockMapData_->NumGlobalElements_ * BlockMapData_->ElementSize_;
00190     BlockMapData_->NumMyPoints_ = BlockMapData_->NumMyElements_ * BlockMapData_->ElementSize_;
00191 
00192     BlockMapData_->MinAllGID_ = BlockMapData_->IndexBase_;
00193     BlockMapData_->MaxAllGID_ = BlockMapData_->MinAllGID_ + BlockMapData_->NumGlobalElements_ - 1;
00194 
00195     // Use the ScanSum function to compute a prefix sum of the number of points
00196     long long tmp2_NumMyElements = BlockMapData_->NumMyElements_;
00197     BlockMapData_->Comm_->ScanSum(&tmp2_NumMyElements, &BlockMapData_->MaxMyGID_, 1);
00198 
00199     long long start_index = BlockMapData_->MaxMyGID_ - BlockMapData_->NumMyElements_;
00200     BlockMapData_->MinMyGID_ = start_index + BlockMapData_->IndexBase_;
00201     BlockMapData_->MaxMyGID_ = BlockMapData_->MinMyGID_ + BlockMapData_->NumMyElements_ - 1;
00202   }
00203   else
00204     throw ReportError("Internal Error.  Report to Epetra developer", -99);
00205 
00206 
00207   EndOfConstructorOps();
00208 }
00209 
00210 //==============================================================================
00211 
00212 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00213 Epetra_BlockMap::Epetra_BlockMap(long long NumGlobal_Elements, int NumMy_Elements,
00214       int Element_Size, int Index_Base, const Epetra_Comm& comm)
00215   : Epetra_Object("Epetra::BlockMap"),
00216     BlockMapData_(0)
00217 {
00218   const bool IsLongLong = true;
00219   ConstructUserLinear(NumGlobal_Elements, NumMy_Elements, Element_Size,static_cast<long long>(Index_Base), comm, IsLongLong);
00220 }
00221 
00222 Epetra_BlockMap::Epetra_BlockMap(long long NumGlobal_Elements, int NumMy_Elements,
00223       int Element_Size, long long Index_Base, const Epetra_Comm& comm)
00224   : Epetra_Object("Epetra::BlockMap"),
00225     BlockMapData_(0)
00226 {
00227   const bool IsLongLong = true;
00228   ConstructUserLinear(NumGlobal_Elements, NumMy_Elements, Element_Size,Index_Base, comm, IsLongLong);
00229 }
00230 #endif
00231 
00232 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00233 Epetra_BlockMap::Epetra_BlockMap(int NumGlobal_Elements, int NumMy_Elements,
00234       int Element_Size, int Index_Base, const Epetra_Comm& comm)
00235   : Epetra_Object("Epetra::BlockMap"),
00236     BlockMapData_(0)
00237 {
00238   const bool IsLongLong = false;
00239   ConstructUserLinear((long long)NumGlobal_Elements, NumMy_Elements, Element_Size,Index_Base, comm, IsLongLong);
00240 }
00241 #endif
00242 
00243 // Epetra_BlockMap constructor for a user-defined arbitrary distribution of constant size elements.
00244 template<typename int_type>
00245 void Epetra_BlockMap::ConstructUserConstant(int_type NumGlobal_Elements, int NumMy_Elements,
00246                                  const int_type * myGlobalElements,
00247          int Element_Size, int_type indexBase,
00248                                  const Epetra_Comm& comm, bool IsLongLong)
00249 {
00250   int i;
00251   // Each processor gets NumMyElements points
00252 
00253   if (NumGlobal_Elements < -1)
00254     throw ReportError("NumGlobal_Elements = " + toString(NumGlobal_Elements) + ".  Should be >= -1.", -1);
00255   if (NumMy_Elements < 0)
00256     throw ReportError("NumMy_Elements = " + toString(NumMy_Elements) + ".  Should be >= 0.", -2);
00257   if (Element_Size <= 0)
00258     throw ReportError("ElementSize = " + toString(Element_Size) + ". Should be > 0.", -3);
00259 
00260   // Allocate storage for global index list information
00261 
00262   BlockMapData_ = new Epetra_BlockMapData(NumGlobal_Elements, Element_Size, indexBase, comm, IsLongLong);
00263   if (NumMy_Elements > 0) {
00264     int errorcode = SizeMyGlobalElement<int_type>(NumMy_Elements);
00265     if(errorcode != 0)
00266       throw ReportError("Error with MyGlobalElements allocation.", -99);
00267   }
00268 
00269   BlockMapData_->NumMyElements_ = NumMy_Elements;
00270   BlockMapData_->MinMyElementSize_ = BlockMapData_->ElementSize_;
00271   BlockMapData_->MaxMyElementSize_ = BlockMapData_->ElementSize_;
00272   BlockMapData_->MinElementSize_ = BlockMapData_->ElementSize_;
00273   BlockMapData_->MaxElementSize_ = BlockMapData_->ElementSize_;
00274   BlockMapData_->ConstantElementSize_ = true;
00275   BlockMapData_->LinearMap_ = false;
00276   // Get processor information
00277 
00278   int NumProc = comm.NumProc();
00279   if (NumMy_Elements > 0) {
00280     // Compute min/max GID on this processor
00281     BlockMapData_->MinMyGID_ = myGlobalElements[0];
00282     BlockMapData_->MaxMyGID_ = myGlobalElements[0];
00283     for (i = 0; i < NumMy_Elements; i++) {
00284       MyGlobalElementVal<int_type>(i) = myGlobalElements[i];
00285       BlockMapData_->MinMyGID_ = EPETRA_MIN(BlockMapData_->MinMyGID_, (long long) myGlobalElements[i]);
00286       BlockMapData_->MaxMyGID_ = EPETRA_MAX(BlockMapData_->MaxMyGID_, (long long) myGlobalElements[i]);
00287     }
00288   }
00289   else {
00290     BlockMapData_->MinMyGID_ = BlockMapData_->IndexBase_;
00291     BlockMapData_->MaxMyGID_ = BlockMapData_->IndexBase_ - 1;
00292   }
00293 
00294   BlockMapData_->DistributedGlobal_ = IsDistributedGlobal(NumGlobal_Elements, NumMy_Elements);
00295 
00296   // Local Map and uniprocessor case:  Each processor gets a complete copy of all elements
00297   if (!BlockMapData_->DistributedGlobal_ || NumProc==1) {
00298     BlockMapData_->NumGlobalElements_ = BlockMapData_->NumMyElements_;
00299     CheckValidNGE(NumGlobal_Elements);
00300     BlockMapData_->NumGlobalPoints_ = BlockMapData_->NumGlobalElements_ * BlockMapData_->ElementSize_;
00301     BlockMapData_->NumMyPoints_ = BlockMapData_->NumMyElements_ * BlockMapData_->ElementSize_;
00302 
00303     BlockMapData_->MinAllGID_ = BlockMapData_->MinMyGID_;
00304     BlockMapData_->MaxAllGID_ = BlockMapData_->MaxMyGID_;
00305   }
00306   else if (NumProc > 1) {
00307     // Sum up all local element counts to get global count
00308     long long tmp_NumMyElements = BlockMapData_->NumMyElements_;
00309     BlockMapData_->Comm_->SumAll(&tmp_NumMyElements, &BlockMapData_->NumGlobalElements_, 1);
00310     CheckValidNGE(NumGlobal_Elements);
00311 
00312     BlockMapData_->NumGlobalPoints_ = BlockMapData_->NumGlobalElements_ * BlockMapData_->ElementSize_;
00313     BlockMapData_->NumMyPoints_ = BlockMapData_->NumMyElements_ * BlockMapData_->ElementSize_;
00314 
00315     // Use the Allreduce function to find min/max GID
00316     long long *tmp_send = new long long[2];
00317     long long *tmp_recv = new long long[2];
00318     if (BlockMapData_->NumMyElements_ > 0 ||
00319         BlockMapData_->NumGlobalElements_ == 0) // This test restores behavior for empty map
00320       tmp_send[0] = - BlockMapData_->MinMyGID_; // Negative sign lets us do one reduction
00321     else
00322       tmp_send[0] = -(std::numeric_limits<int_type>::max())-1;
00323     tmp_send[1] =   BlockMapData_->MaxMyGID_;
00324     BlockMapData_->Comm_->MaxAll(tmp_send, tmp_recv, 2);
00325     BlockMapData_->MinAllGID_ = - tmp_recv[0];
00326     BlockMapData_->MaxAllGID_ =   tmp_recv[1];
00327     delete [] tmp_send;
00328     delete [] tmp_recv;
00329     if (BlockMapData_->MinAllGID_ < BlockMapData_->IndexBase_)
00330       throw ReportError("Minimum global element index = " + toString(BlockMapData_->MinAllGID_) +
00331       " is less than index base = " + toString(BlockMapData_->IndexBase_) +".", -5);
00332   }
00333   else
00334     throw ReportError("Internal Error.  Report to Epetra developer", -99);
00335 
00336 
00337   EndOfConstructorOps();
00338 }
00339 
00340 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00341 Epetra_BlockMap::Epetra_BlockMap(long long NumGlobal_Elements, int NumMy_Elements,
00342                                  const long long * myGlobalElements,
00343          int Element_Size, int indexBase,
00344                                  const Epetra_Comm& comm)
00345   : Epetra_Object("Epetra::BlockMap"),
00346     BlockMapData_(0)
00347 {
00348   const bool IsLongLong = true;
00349   ConstructUserConstant(NumGlobal_Elements, NumMy_Elements, myGlobalElements,
00350     Element_Size, static_cast<long long>(indexBase), comm, IsLongLong);
00351 }
00352 
00353 Epetra_BlockMap::Epetra_BlockMap(long long NumGlobal_Elements, int NumMy_Elements,
00354                                  const long long * myGlobalElements,
00355          int Element_Size, long long indexBase,
00356                                  const Epetra_Comm& comm)
00357   : Epetra_Object("Epetra::BlockMap"),
00358     BlockMapData_(0)
00359 {
00360   const bool IsLongLong = true;
00361   ConstructUserConstant(NumGlobal_Elements, NumMy_Elements, myGlobalElements,
00362     Element_Size, indexBase, comm, IsLongLong);
00363 }
00364 #endif
00365 
00366 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00367 Epetra_BlockMap::Epetra_BlockMap(int NumGlobal_Elements, int NumMy_Elements,
00368                                  const int * myGlobalElements,
00369          int Element_Size, int indexBase,
00370                                  const Epetra_Comm& comm)
00371   : Epetra_Object("Epetra::BlockMap"),
00372     BlockMapData_(0)
00373 {
00374   const bool IsLongLong = false;
00375   ConstructUserConstant(NumGlobal_Elements, NumMy_Elements, myGlobalElements,
00376     Element_Size, indexBase, comm, IsLongLong);
00377 }
00378 #endif
00379 
00380 
00381 //==============================================================================
00382 // Epetra_BlockMap constructor function for a user-defined arbitrary distribution of variable size elements.
00383 template<typename int_type>
00384 void Epetra_BlockMap::ConstructUserVariable(int_type NumGlobal_Elements, int NumMy_Elements,
00385                                  const int_type * myGlobalElements,
00386          const int *elementSizeList, int_type indexBase,
00387                                  const Epetra_Comm& comm, bool IsLongLong)
00388 {
00389 
00390   int i;
00391   // Each processor gets NumMyElements points
00392 
00393   if (NumGlobal_Elements < -1)
00394     throw ReportError("NumGlobal_Elements = " + toString(NumGlobal_Elements) + ".  Should be >= -1.", -1);
00395   if (NumMy_Elements < 0)
00396     throw ReportError("NumMy_Elements = " + toString(NumMy_Elements) + ".  Should be >= 0.", -2);
00397   for (i = 0; i < NumMy_Elements; i++)
00398     if (elementSizeList[i] <= 0)
00399       throw ReportError("elementSizeList["+toString(i)+"] = " + toString(elementSizeList[i]) + ". Should be > 0.", -3);
00400 
00401   BlockMapData_ = new Epetra_BlockMapData(NumGlobal_Elements, 0, indexBase, comm, IsLongLong);
00402   BlockMapData_->NumMyElements_ = NumMy_Elements;
00403   BlockMapData_->ConstantElementSize_ = false;
00404   BlockMapData_->LinearMap_ = false;
00405   // Allocate storage for global index list and element size information
00406 
00407   if (NumMy_Elements > 0) {
00408     int errorcode = SizeMyGlobalElement<int_type>(NumMy_Elements);
00409     if(errorcode != 0)
00410       throw ReportError("Error with MyGlobalElements allocation.", -99);
00411     errorcode = BlockMapData_->ElementSizeList_.Size(NumMy_Elements);
00412     if(errorcode != 0)
00413       throw ReportError("Error with ElementSizeList allocation.", -99);
00414   }
00415   // Get processor information
00416 
00417   int NumProc = comm.NumProc();
00418 
00419   if (NumMy_Elements > 0) {
00420     // Compute min/max GID and element size, number of points on this processor
00421     BlockMapData_->MinMyGID_ = myGlobalElements[0];
00422     BlockMapData_->MaxMyGID_ = myGlobalElements[0];
00423     BlockMapData_->MinMyElementSize_ = elementSizeList[0];
00424     BlockMapData_->MaxMyElementSize_ = elementSizeList[0];
00425     BlockMapData_->NumMyPoints_ = 0;
00426     for (i = 0; i < NumMy_Elements; i++) {
00427       MyGlobalElementVal<int_type>(i) = myGlobalElements[i];
00428       BlockMapData_->ElementSizeList_[i] = elementSizeList[i];
00429       BlockMapData_->MinMyGID_ = EPETRA_MIN(BlockMapData_->MinMyGID_,(long long) myGlobalElements[i]);
00430       BlockMapData_->MaxMyGID_ = EPETRA_MAX(BlockMapData_->MaxMyGID_,(long long) myGlobalElements[i]);
00431       BlockMapData_->MinMyElementSize_ = EPETRA_MIN(BlockMapData_->MinMyElementSize_,elementSizeList[i]);
00432       BlockMapData_->MaxMyElementSize_ = EPETRA_MAX(BlockMapData_->MaxMyElementSize_,elementSizeList[i]);
00433       BlockMapData_->NumMyPoints_ += elementSizeList[i];
00434     }
00435   }
00436   else {
00437     BlockMapData_->MinMyGID_ = BlockMapData_->IndexBase_;
00438     BlockMapData_->MaxMyGID_ = BlockMapData_->IndexBase_ - 1;
00439     BlockMapData_->MinMyElementSize_ = 1;
00440     BlockMapData_->MaxMyElementSize_ = 1;
00441     BlockMapData_->NumMyPoints_ = 0;
00442   }
00443 
00444   BlockMapData_->DistributedGlobal_ = IsDistributedGlobal(NumGlobal_Elements, NumMy_Elements);
00445 
00446   // Local Map and uniprocessor case:  Each processor gets a complete copy of all elements
00447   if (!BlockMapData_->DistributedGlobal_ || NumProc == 1) {
00448     BlockMapData_->NumGlobalElements_ = BlockMapData_->NumMyElements_;
00449     CheckValidNGE(NumGlobal_Elements);
00450     BlockMapData_->NumGlobalPoints_ = BlockMapData_->NumMyPoints_;
00451 
00452     BlockMapData_->MinAllGID_ = BlockMapData_->MinMyGID_;
00453     BlockMapData_->MaxAllGID_ = BlockMapData_->MaxMyGID_;
00454     BlockMapData_->MinElementSize_ = BlockMapData_->MinMyElementSize_;
00455     BlockMapData_->MaxElementSize_ = BlockMapData_->MaxMyElementSize_;
00456   }
00457   else if (NumProc > 1) {
00458     // Sum up all local element and point counts to get global counts
00459     int_type *tmp_send = new int_type[4];
00460     int_type *tmp_recv = new int_type[4];
00461     tmp_send[0] = BlockMapData_->NumMyElements_;
00462     tmp_send[1] = BlockMapData_->NumMyPoints_;
00463     BlockMapData_->Comm_->SumAll(tmp_send, tmp_recv, 2);
00464     BlockMapData_->NumGlobalElements_ =  tmp_recv[0];
00465     BlockMapData_->NumGlobalPoints_ = tmp_recv[1];
00466 
00467     CheckValidNGE(NumGlobal_Elements);
00468 
00469     // Use the MaxAll function to find min/max GID
00470     if (BlockMapData_->NumMyElements_ > 0 ||
00471         BlockMapData_->NumGlobalElements_ == 0) // This test restores behavior for empty map
00472       tmp_send[0] = - static_cast<int_type>(BlockMapData_->MinMyGID_); // Negative sign lets us do one reduction
00473     else
00474       tmp_send[0] = -(std::numeric_limits<int_type>::max())-1;
00475     tmp_send[1] =   static_cast<int_type>(BlockMapData_->MaxMyGID_);
00476     tmp_send[2] = - BlockMapData_->MinMyElementSize_;
00477     if (BlockMapData_->NumMyElements_ == 0)
00478       tmp_send[2] = - static_cast<int_type>(BlockMapData_->NumGlobalPoints_); // This processor has no elements, so should not sizes.
00479     tmp_send[3] =   BlockMapData_->MaxMyElementSize_;
00480 
00481     BlockMapData_->Comm_->MaxAll(tmp_send, tmp_recv, 4);
00482 
00483     BlockMapData_->MinAllGID_ =      - tmp_recv[0];
00484     BlockMapData_->MaxAllGID_ =        tmp_recv[1];
00485     BlockMapData_->MinElementSize_ = - (int) tmp_recv[2];
00486     BlockMapData_->MaxElementSize_ =   (int) tmp_recv[3];
00487 
00488     delete [] tmp_send;
00489     delete [] tmp_recv;
00490 
00491     // Check for constant element size
00492     if (BlockMapData_->MinElementSize_==BlockMapData_->MaxElementSize_) {
00493       BlockMapData_->ElementSize_ = BlockMapData_->MinElementSize_;
00494       BlockMapData_->ConstantElementSize_ = true;
00495     }
00496 
00497     if (BlockMapData_->MinAllGID_ < BlockMapData_->IndexBase_)
00498       throw ReportError("Minimum global element index = " + toString(BlockMapData_->MinAllGID_) +
00499       " is less than index base = " + toString(BlockMapData_->IndexBase_) +".", -5);
00500   }
00501   else
00502     throw ReportError("Internal Error.  Report to Epetra developer", -99);
00503 
00504 
00505   EndOfConstructorOps();
00506 }
00507 
00508 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00509 Epetra_BlockMap::Epetra_BlockMap(long long NumGlobal_Elements, int NumMy_Elements,
00510                                  const long long * myGlobalElements,
00511          const int *elementSizeList, int indexBase,
00512                                  const Epetra_Comm& comm)
00513   : Epetra_Object("Epetra::BlockMap"),
00514     BlockMapData_(0)
00515 {
00516   const bool IsLongLong = true;
00517   ConstructUserVariable(NumGlobal_Elements, NumMy_Elements, myGlobalElements,
00518     elementSizeList, static_cast<long long>(indexBase), comm, IsLongLong);
00519 }
00520 
00521 Epetra_BlockMap::Epetra_BlockMap(long long NumGlobal_Elements, int NumMy_Elements,
00522                                  const long long * myGlobalElements,
00523          const int *elementSizeList, long long indexBase,
00524                                  const Epetra_Comm& comm)
00525   : Epetra_Object("Epetra::BlockMap"),
00526     BlockMapData_(0)
00527 {
00528   const bool IsLongLong = true;
00529   ConstructUserVariable(NumGlobal_Elements, NumMy_Elements, myGlobalElements,
00530     elementSizeList, indexBase, comm, IsLongLong);
00531 }
00532 #endif
00533 
00534 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00535 Epetra_BlockMap::Epetra_BlockMap(int NumGlobal_Elements, int NumMy_Elements,
00536                                  const int * myGlobalElements,
00537          const int *elementSizeList, int indexBase,
00538                                  const Epetra_Comm& comm)
00539   : Epetra_Object("Epetra::BlockMap"),
00540     BlockMapData_(0)
00541 {
00542   const bool IsLongLong = false;
00543   ConstructUserVariable(NumGlobal_Elements, NumMy_Elements, myGlobalElements,
00544     elementSizeList, indexBase, comm, IsLongLong);
00545 }
00546 #endif
00547 
00548 
00549 //==============================================================================
00550 // Epetra_BlockMap constructor function for a user-defined arbitrary distribution of variable size elements,
00551 // with all the information on globals provided by the user.
00552 template<typename int_type>
00553 void Epetra_BlockMap::ConstructUserConstantNoComm(int_type NumGlobal_Elements, int NumMy_Elements,
00554                                                   const int_type * myGlobalElements,
00555                                                   int Element_Size, int_type indexBase,
00556                                                   const Epetra_Comm& comm, bool IsLongLong,
00557                                                   bool UserIsDistributedGlobal,
00558                                                   int_type UserMinAllGID, int_type UserMaxAllGID)
00559 {
00560 
00561 
00562   int i;
00563   // Each processor gets NumMyElements points
00564 
00565   if (NumGlobal_Elements < -1)
00566     throw ReportError("NumGlobal_Elements = " + toString(NumGlobal_Elements) + ".  Should be >= -1.", -1);
00567   if (NumMy_Elements < 0)
00568     throw ReportError("NumMy_Elements = " + toString(NumMy_Elements) + ".  Should be >= 0.", -2);
00569   if (Element_Size <= 0)
00570     throw ReportError("ElementSize = " + toString(Element_Size) + ". Should be > 0.", -3);
00571 
00572   // Allocate storage for global index list information
00573 
00574   BlockMapData_ = new Epetra_BlockMapData(NumGlobal_Elements, Element_Size, indexBase, comm, IsLongLong);
00575   if (NumMy_Elements > 0) {
00576     int errorcode = SizeMyGlobalElement<int_type>(NumMy_Elements);
00577     if(errorcode != 0)
00578       throw ReportError("Error with MyGlobalElements allocation.", -99);
00579   }
00580 
00581   BlockMapData_->NumMyElements_ = NumMy_Elements;
00582   BlockMapData_->MinMyElementSize_ = BlockMapData_->ElementSize_;
00583   BlockMapData_->MaxMyElementSize_ = BlockMapData_->ElementSize_;
00584   BlockMapData_->MinElementSize_ = BlockMapData_->ElementSize_;
00585   BlockMapData_->MaxElementSize_ = BlockMapData_->ElementSize_;
00586   BlockMapData_->ConstantElementSize_ = true;
00587   BlockMapData_->LinearMap_ = false;
00588   // Get processor information
00589 
00590   int NumProc = comm.NumProc();
00591   if (NumMy_Elements > 0) {
00592     // Compute min/max GID on this processor
00593     BlockMapData_->MinMyGID_ = myGlobalElements[0];
00594     BlockMapData_->MaxMyGID_ = myGlobalElements[0];
00595     for (i = 0; i < NumMy_Elements; i++) {
00596       MyGlobalElementVal<int_type>(i) = myGlobalElements[i];
00597       BlockMapData_->MinMyGID_ = EPETRA_MIN(BlockMapData_->MinMyGID_, (long long) myGlobalElements[i]);
00598       BlockMapData_->MaxMyGID_ = EPETRA_MAX(BlockMapData_->MaxMyGID_, (long long) myGlobalElements[i]);
00599     }
00600   }
00601   else {
00602     BlockMapData_->MinMyGID_ = BlockMapData_->IndexBase_;
00603     BlockMapData_->MaxMyGID_ = BlockMapData_->IndexBase_ - 1;
00604   }
00605 
00606   BlockMapData_->DistributedGlobal_ = UserIsDistributedGlobal;
00607 
00608   // Local Map and uniprocessor case:  Each processor gets a complete copy of all elements
00609   if (!BlockMapData_->DistributedGlobal_ || NumProc==1) {
00610     BlockMapData_->NumGlobalElements_ = BlockMapData_->NumMyElements_;
00611     CheckValidNGE(NumGlobal_Elements);
00612     BlockMapData_->NumGlobalPoints_ = BlockMapData_->NumGlobalElements_ * BlockMapData_->ElementSize_;
00613     BlockMapData_->NumMyPoints_ = BlockMapData_->NumMyElements_ * BlockMapData_->ElementSize_;
00614 
00615     BlockMapData_->MinAllGID_ = BlockMapData_->MinMyGID_;
00616     BlockMapData_->MaxAllGID_ = BlockMapData_->MaxMyGID_;
00617   }
00618   else if (NumProc > 1) {
00619     if(NumGlobal_Elements==-1) {
00620       // Sum up all local element counts to get global count
00621       long long tmp_NumMyElements = BlockMapData_->NumMyElements_;
00622       BlockMapData_->Comm_->SumAll(&tmp_NumMyElements, &BlockMapData_->NumGlobalElements_, 1);
00623     }
00624     else {
00625       // User provides this information
00626       BlockMapData_->NumGlobalElements_ = NumGlobal_Elements;
00627     }
00628     CheckValidNGE(BlockMapData_->NumGlobalElements_);
00629 
00630     BlockMapData_->NumGlobalPoints_ = BlockMapData_->NumGlobalElements_ * BlockMapData_->ElementSize_;
00631     BlockMapData_->NumMyPoints_     = BlockMapData_->NumMyElements_ * BlockMapData_->ElementSize_;
00632 
00633     BlockMapData_->MinAllGID_       = UserMinAllGID;
00634     BlockMapData_->MaxAllGID_       = UserMaxAllGID;
00635     if (BlockMapData_->MinAllGID_ < BlockMapData_->IndexBase_)
00636       throw ReportError("Minimum global element index = " + toString(BlockMapData_->MinAllGID_) +
00637       " is less than index base = " + toString(BlockMapData_->IndexBase_) +".", -5);
00638   }
00639   else
00640     throw ReportError("Internal Error.  Report to Epetra developer", -99);
00641 
00642 
00643   EndOfConstructorOps();
00644 }
00645 
00646 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00647 Epetra_BlockMap::Epetra_BlockMap(long long NumGlobal_Elements, int NumMy_Elements,
00648                                  const long long * myGlobalElements,
00649                                  int theElementSize, int indexBase,
00650                                  const Epetra_Comm& comm,
00651                                  bool UserIsDistributedGlobal,
00652                                  long long UserMinAllGID, long long UserMaxAllGID)
00653   : Epetra_Object("Epetra::BlockMap"),
00654     BlockMapData_(0)
00655 {
00656   const bool IsLongLong = true;
00657   ConstructUserConstantNoComm(NumGlobal_Elements, NumMy_Elements, myGlobalElements,
00658                               theElementSize, (long long) indexBase, comm, IsLongLong,
00659                               UserIsDistributedGlobal, UserMinAllGID, UserMaxAllGID);
00660 }
00661 Epetra_BlockMap::Epetra_BlockMap(long long NumGlobal_Elements, int NumMy_Elements,
00662                                  const long long * myGlobalElements,
00663                                  int theElementSize, long long indexBase,
00664                                  const Epetra_Comm& comm,
00665                                  bool UserIsDistributedGlobal,
00666                                  long long UserMinAllGID, long long UserMaxAllGID)
00667   : Epetra_Object("Epetra::BlockMap"),
00668     BlockMapData_(0)
00669 {
00670   const bool IsLongLong = true;
00671   ConstructUserConstantNoComm(NumGlobal_Elements, NumMy_Elements, myGlobalElements,
00672                               theElementSize, indexBase, comm, IsLongLong,
00673                               UserIsDistributedGlobal, UserMinAllGID, UserMaxAllGID);
00674 }
00675 #endif
00676 
00677 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00678 Epetra_BlockMap::Epetra_BlockMap(int NumGlobal_Elements, int NumMy_Elements,
00679                                  const int * myGlobalElements,
00680                                  int theElementSize, int indexBase,
00681                                  const Epetra_Comm& comm,
00682                                  bool UserIsDistributedGlobal,
00683                                  int UserMinAllGID, int UserMaxAllGID)
00684   : Epetra_Object("Epetra::BlockMap"),
00685     BlockMapData_(0)
00686 {
00687   const bool IsLongLong = false;
00688   ConstructUserConstantNoComm(NumGlobal_Elements, NumMy_Elements, myGlobalElements,
00689                              theElementSize, indexBase, comm, IsLongLong,
00690                              UserIsDistributedGlobal, UserMinAllGID, UserMaxAllGID);
00691 }
00692 #endif
00693 
00694 
00695 
00696 
00697 //==============================================================================
00698 Epetra_BlockMap::Epetra_BlockMap(const Epetra_BlockMap& map)
00699   : Epetra_Object(map.Label()),
00700     BlockMapData_(map.BlockMapData_)
00701 {
00702   BlockMapData_->IncrementReferenceCount();
00703 
00704   // This call appears to be unnecessary overhead.  Removed 10-Aug-2004 maherou.
00705   // GlobalToLocalSetup(); // Setup any information for making global index to local index translation fast.
00706 }
00707 
00708 //==============================================================================
00709 bool Epetra_BlockMap::SameAs(const Epetra_BlockMap & Map) const {
00710 
00711   // Quickest test: See if both maps share an inner data class
00712   if (this->BlockMapData_ == Map.BlockMapData_)
00713     return(true);
00714 
00715   if(!GlobalIndicesTypeMatch(Map))
00716     return(false);
00717 
00718   // Next check other global properties that are easy global attributes
00719   if (BlockMapData_->MinAllGID_ != Map.MinAllGID64() ||
00720       BlockMapData_->MaxAllGID_ != Map.MaxAllGID64() ||
00721       BlockMapData_->NumGlobalElements_ != Map.NumGlobalElements64() ||
00722       BlockMapData_->IndexBase_ != Map.IndexBase64())
00723     return(false);
00724 
00725   // Last possible global check for constant element sizes
00726   if (BlockMapData_->ConstantElementSize_ && BlockMapData_->ElementSize_!=Map.ElementSize())
00727     return(false);
00728 
00729   // If we get this far, we need to check local properties and then check across
00730   // all processors to see if local properties are all true
00731 
00732   int numMyElements = BlockMapData_->NumMyElements_;
00733 
00734   int MySameMap = 1; // Assume not needed
00735 
00736   // First check if number of element is the same in each map
00737   if (numMyElements != Map.NumMyElements()) MySameMap = 0;
00738 
00739   // If numMyElements is the same, check to see that list of GIDs is the same
00740   if (MySameMap==1) {
00741     if (LinearMap() && Map.LinearMap() ) {
00742       // For linear maps, just need to check whether lower bound is the same
00743       if (MinMyGID64() != Map.MinMyGID64() )
00744         MySameMap = 0;
00745     }
00746     else {
00747       for (int i = 0; i < numMyElements; i++) {
00748         if (GID64(i) != Map.GID64(i)) {
00749           MySameMap = 0;
00750           break;
00751         }
00752       }
00753     }
00754   }
00755 //    for (int i = 0; i < numMyElements; i++)
00756 //      if (GID64(i) != Map.GID64(i)) MySameMap = 0;
00757 
00758   // If GIDs are the same, check to see element sizes are the same
00759   if (MySameMap==1 && !BlockMapData_->ConstantElementSize_) {
00760     int * sizeList1 = ElementSizeList();
00761     int * sizeList2 = Map.ElementSizeList();
00762     for (int i = 0; i < numMyElements; i++) if (sizeList1[i] != sizeList2[i]) MySameMap=0;
00763   }
00764   // Now get min of MySameMap across all processors
00765 
00766   int GlobalSameMap = 0;
00767   int err = Comm().MinAll(&MySameMap, &GlobalSameMap, 1);
00768   assert(err==0);
00769 
00770   return(GlobalSameMap==1);
00771 }
00772 
00773 //==============================================================================
00774 bool Epetra_BlockMap::PointSameAs(const Epetra_BlockMap & Map) const
00775 {
00776   if (this->BlockMapData_ == Map.BlockMapData_)
00777     return(true);
00778 
00779   if(!GlobalIndicesTypeMatch(Map))
00780     return(false);
00781 
00782   if (BlockMapData_->NumGlobalPoints_ != Map.NumGlobalPoints64() )
00783     return(false);
00784 
00785   // If we get this far, we need to check local properties and then check across
00786   // all processors to see if local properties are all true
00787 
00788   int MySameMap = 1; // Assume not needed
00789   if (BlockMapData_->NumMyPoints_ != Map.NumMyPoints())
00790     MySameMap = 0;
00791 
00792   int GlobalSameMap = 0;
00793   int err = Comm().MinAll(&MySameMap, &GlobalSameMap, 1);
00794   assert( err == 0 );
00795 
00796   return(GlobalSameMap==1);
00797 }
00798 
00799 //==============================================================================
00800 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00801 int Epetra_BlockMap::MyGlobalElements(long long * myGlobalElements) const
00802 {
00803   // Although one can populate long long data from int data, we don't
00804   // allow it to maintain int/long long symmetry.
00805   if(!BlockMapData_->GlobalIndicesLongLong_)
00806     throw ReportError("Epetra_BlockMap::MyGlobalElements(long long *) ERROR, Can't call for non long long* map.",-1);
00807 
00808   // If the global element list is not create, then do so.  This can only happen when
00809   // a linear distribution has been specified.  Thus we can easily construct the update
00810   // list in this case.
00811 
00812   int i;
00813   int numMyElements = BlockMapData_->NumMyElements_;
00814 
00815   if (BlockMapData_->MyGlobalElements_LL_.Length() == 0)
00816     for (i = 0; i < numMyElements; i++)
00817       myGlobalElements[i] = BlockMapData_->MinMyGID_ + i;
00818   else
00819     for (i = 0; i < numMyElements; i++)
00820       myGlobalElements[i] = BlockMapData_->MyGlobalElements_LL_[i];
00821   return(0);
00822 }
00823 #endif
00824 
00825 //==============================================================================
00826 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00827 int Epetra_BlockMap::MyGlobalElements(int * myGlobalElements) const
00828 {
00829   if(!BlockMapData_->GlobalIndicesInt_)
00830     throw ReportError("Epetra_BlockMap::MyGlobalElements(int *) ERROR, Can't call for non int* map.",-1);
00831 
00832   // If the global element list is not create, then do so.  This can only happen when
00833   // a linear distribution has been specified.  Thus we can easily construct the update
00834   // list in this case.
00835 
00836   int i;
00837   int numMyElements = BlockMapData_->NumMyElements_;
00838 
00839   if (BlockMapData_->MyGlobalElements_int_.Length() == 0)
00840     for (i = 0; i < numMyElements; i++)
00841       myGlobalElements[i] = (int) BlockMapData_->MinMyGID_ + i;
00842   else
00843     for (i = 0; i < numMyElements; i++)
00844       myGlobalElements[i] = (int) BlockMapData_->MyGlobalElements_int_[i];
00845   return(0);
00846 }
00847 #endif
00848 
00849 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00850 int Epetra_BlockMap::MyGlobalElementsPtr(long long *& MyGlobalElementList) const
00851 {
00852   MyGlobalElementList = MyGlobalElements64();
00853   return(0);
00854 }
00855 #endif
00856 
00857 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00858 int Epetra_BlockMap::MyGlobalElementsPtr(int *& MyGlobalElementList) const
00859 {
00860   MyGlobalElementList = MyGlobalElements();
00861   return(0);
00862 }
00863 #endif
00864 //==============================================================================
00865 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00866 int * Epetra_BlockMap::MyGlobalElements() const {
00867   if(!BlockMapData_->GlobalIndicesInt_)
00868     throw ReportError("Epetra_BlockMap::MyGlobalElements() ERROR, Can't call for non int* map.",-1);
00869 
00870   int numMyElements = BlockMapData_->NumMyElements_;
00871 
00872   // If ElementSizeList not built, do so
00873   if(BlockMapData_->MyGlobalElements_int_.Length() == 0 && numMyElements > 0) {
00874     int errorcode = BlockMapData_->MyGlobalElements_int_.Size(numMyElements + 1);
00875     if(errorcode != 0)
00876       throw ReportError("Error with MyGlobalElements allocation.", -99);
00877 
00878     // Build the array
00879     for (int i = 0; i < numMyElements; i++)
00880       BlockMapData_->MyGlobalElements_int_[i] = (int) BlockMapData_->MinMyGID_ + i;
00881   }
00882   return(BlockMapData_->MyGlobalElements_int_.Values());
00883 }
00884 #endif
00885 //==============================================================================
00886 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00887 long long * Epetra_BlockMap::MyGlobalElements64() const {
00888   if(!BlockMapData_->GlobalIndicesLongLong_)
00889     throw ReportError("Epetra_BlockMap::MyGlobalElements64 ERROR, Can't call for non long long* map.",-1);
00890 
00891   int numMyElements = BlockMapData_->NumMyElements_;
00892 
00893   // If ElementSizeList not built, do so
00894   if(BlockMapData_->MyGlobalElements_LL_.Length() == 0 && numMyElements > 0) {
00895     int errorcode = BlockMapData_->MyGlobalElements_LL_.Size(numMyElements + 1);
00896     if(errorcode != 0)
00897       throw ReportError("Error with MyGlobalElements allocation.", -99);
00898 
00899     // Build the array
00900     for (int i = 0; i < numMyElements; i++)
00901       BlockMapData_->MyGlobalElements_LL_[i] = BlockMapData_->MinMyGID_ + i;
00902   }
00903   return(BlockMapData_->MyGlobalElements_LL_.Values());
00904 }
00905 #endif
00906 //==============================================================================
00907 int Epetra_BlockMap::FirstPointInElement(int lid) const
00908 {
00909   if (!MyLID(lid))
00910     EPETRA_CHK_ERR(-1);
00911 
00912   int entry;
00913 
00914   if (ConstantElementSize())
00915     entry = MaxElementSize() * lid; // convert to vector entry
00916   else {
00917     int * entrylist = FirstPointInElementList(); // get entry list
00918     entry = entrylist[lid];
00919   }
00920   return(entry);
00921 }
00922 
00923 //==============================================================================
00924 int Epetra_BlockMap::FirstPointInElementList(int * firstPointInElementList) const
00925 {
00926   // If the first element entry list is not create, then do so.
00927 
00928   // Note: This array is of length NumMyElement+1
00929 
00930   int i;
00931   int numMyElements = BlockMapData_->NumMyElements_;
00932 
00933   if (BlockMapData_->FirstPointInElementList_.Length() == 0) {
00934     firstPointInElementList[0] = 0; // First element of first entry is always zero
00935 
00936     if (BlockMapData_->ConstantElementSize_)
00937       for (i = 0; i < numMyElements; i++)
00938   firstPointInElementList[i+1] = firstPointInElementList[i] + BlockMapData_->ElementSize_;
00939     else
00940       for (i = 0; i < numMyElements; i++)
00941   firstPointInElementList[i+1] = firstPointInElementList[i] + BlockMapData_->ElementSizeList_[i];
00942   }
00943   else
00944     for (i = 0; i <= numMyElements; i++)
00945       firstPointInElementList[i] = BlockMapData_->FirstPointInElementList_[i];
00946   return(0);
00947 }
00948 
00949 //==============================================================================
00950 int * Epetra_BlockMap::FirstPointInElementList() const {
00951   int numMyElements = BlockMapData_->NumMyElements_;
00952 
00953   // If ElementSizeList not built, do so
00954   if ((BlockMapData_->FirstPointInElementList_.Length() == 0) && (numMyElements > 0)) {
00955     BlockMapData_->FirstPointInElementList_.Size(BlockMapData_->NumMyElements_ + 1);
00956     BlockMapData_->FirstPointInElementList_[0] = 0; // First element of first entry is always zero
00957     if (BlockMapData_->ConstantElementSize_)
00958       for (int i = 0; i < numMyElements; i++)
00959   BlockMapData_->FirstPointInElementList_[i+1] = BlockMapData_->FirstPointInElementList_[i] + BlockMapData_->ElementSize_;
00960     else
00961       for (int i = 0; i < numMyElements; i++)
00962   BlockMapData_->FirstPointInElementList_[i+1] = BlockMapData_->FirstPointInElementList_[i] + BlockMapData_->ElementSizeList_[i];
00963   }
00964   return(BlockMapData_->FirstPointInElementList_.Values());
00965 }
00966 
00967 //==============================================================================
00968 int Epetra_BlockMap::ElementSizeList(int * elementSizeList) const
00969 {
00970   // If the element size list is not create, then do so.  This can only happen when
00971   // a constant element size has been specified.  Thus we can easily construct the element size
00972   // list in this case.
00973 
00974   int i;
00975   int numMyElements = BlockMapData_->NumMyElements_;
00976 
00977   if (BlockMapData_->ElementSizeList_.Length() == 0)
00978     for (i = 0; i < numMyElements; i++)
00979       elementSizeList[i] = BlockMapData_->ElementSize_;
00980   else
00981     for (i = 0; i < numMyElements; i++)
00982       elementSizeList[i] = BlockMapData_->ElementSizeList_[i];
00983 
00984   return(0);
00985 }
00986 
00987 //==============================================================================
00988 int * Epetra_BlockMap::ElementSizeList() const {
00989   int numMyElements = BlockMapData_->NumMyElements_;
00990 
00991   // If ElementSizeList not built, do so
00992   if ((BlockMapData_->ElementSizeList_.Length() == 0) && (numMyElements > 0)) {
00993     BlockMapData_->ElementSizeList_.Size(numMyElements);
00994     for (int i = 0; i < numMyElements; i++)
00995       BlockMapData_->ElementSizeList_[i] = BlockMapData_->ElementSize_;
00996   }
00997   return(BlockMapData_->ElementSizeList_.Values());
00998 }
00999 
01000 //==============================================================================
01001 int Epetra_BlockMap::PointToElementList(int * pointToElementList) const {
01002   // Build an array such that the local element ID is stored for each point
01003 
01004   int i;
01005   if (BlockMapData_->PointToElementList_.Length() == 0) {
01006     int numMyElements = BlockMapData_->NumMyElements_;
01007     int * ptr = pointToElementList;
01008     for (i = 0; i < numMyElements; i++) {
01009       int Size = ElementSize(i);
01010       for (int j = 0; j < Size; j++)
01011   *ptr++ = i;
01012     }
01013   }
01014   else {
01015     int numMyPoints = BlockMapData_->NumMyPoints_;
01016     for (i = 0; i < numMyPoints; i++)
01017       pointToElementList[i] = BlockMapData_->PointToElementList_[i];
01018   }
01019   return(0);
01020 }
01021 
01022 //==============================================================================
01023 int * Epetra_BlockMap::PointToElementList() const {
01024 
01025   // If PointToElementList not built, do so
01026   if ((BlockMapData_->PointToElementList_.Length() == 0) && (BlockMapData_->NumMyPoints_ > 0)) {
01027     BlockMapData_->PointToElementList_.Size(BlockMapData_->NumMyPoints_);
01028     int numMyElements = BlockMapData_->NumMyElements_;
01029     int * ptr = BlockMapData_->PointToElementList_.Values();
01030     for (int i = 0; i < numMyElements; i++) {
01031       int Size = ElementSize(i);
01032       for (int j = 0; j < Size; j++)
01033   *ptr++ = i;
01034     }
01035   }
01036   return(BlockMapData_->PointToElementList_.Values());
01037 }
01038 
01039 //==============================================================================
01040 int Epetra_BlockMap::ElementSize(int lid) const {
01041 
01042   if (ConstantElementSize())
01043     return(BlockMapData_->ElementSize_);
01044   else
01045     return(BlockMapData_->ElementSizeList_[lid]);
01046 }
01047 
01048 //==============================================================================
01049 bool Epetra_BlockMap::IsOneToOne() const {
01050   if(!BlockMapData_->OneToOneIsDetermined_){
01051     BlockMapData_->OneToOne_ = DetermineIsOneToOne();
01052     BlockMapData_->OneToOneIsDetermined_ = true;
01053   }
01054   return(BlockMapData_->OneToOne_);
01055 }
01056 
01057 //==============================================================================
01058 template<typename int_type>
01059 void Epetra_BlockMap::TGlobalToLocalSetup()
01060 {
01061   int i;
01062   int numMyElements = BlockMapData_->NumMyElements_;
01063 
01064   if (BlockMapData_->NumGlobalElements_ == 0) {
01065     return; // Nothing to do
01066   }
01067 
01068   if (LinearMap() || numMyElements == 0) {
01069     return; // Nothing else to do
01070   }
01071 
01072   // Build LID_ vector to make look up of local index values fast
01073 
01074 #ifdef EPETRA_BLOCKMAP_NEW_LID
01075 
01076   // Here follows an optimization that checks for an initial block of
01077   // contiguous GIDs, and stores the GID->LID mapping for those in a
01078   // more efficient way.  This supports a common use case for
01079   // overlapping Maps, in which owned entries of a vector are ordered
01080   // before halo entries.
01081   //
01082   // Epetra defines EPETRA_BLOCKMAP_NEW_LID by default (see the top of
01083   // this file).
01084 
01085   //check for initial contiguous block
01086   int_type val = MyGlobalElementValGet<int_type>(0);
01087   for( i = 0 ; i < numMyElements; ++i ) {
01088     if (val != MyGlobalElementValGet<int_type>(i)) break;
01089     ++val;
01090   }
01091   BlockMapData_->LastContiguousGIDLoc_ = i - 1;
01092   if (BlockMapData_->LastContiguousGIDLoc_ < 0) {
01093     BlockMapData_->LastContiguousGID_ = MyGlobalElementValGet<int_type>(0);
01094   }
01095   else {
01096     BlockMapData_->LastContiguousGID_ =
01097       MyGlobalElementValGet<int_type>(BlockMapData_->LastContiguousGIDLoc_);
01098   }
01099 
01100   //Hash everything else
01101   if(i < numMyElements) {
01102     if (BlockMapData_->LIDHash_ != NULL) {
01103       delete BlockMapData_->LIDHash_;
01104     }
01105 
01106     BlockMapData_->LIDHash_ = new Epetra_HashTable<int>(numMyElements - i + 1 );
01107     for(; i < numMyElements; ++i )
01108       BlockMapData_->LIDHash_->Add( MyGlobalElementValGet<int_type>(i), i );
01109   }
01110 
01111 #else
01112 
01113   int SpanGID = BlockMapData_->MaxMyGID_ - BlockMapData_->MinMyGID_ + 1;
01114   BlockMapData_->LID_.Size(SpanGID);
01115 
01116   for (i = 0; i < SpanGID; i++)
01117     BlockMapData_->LID_[i] = -1; // Fill all locations with -1
01118 
01119   for (i = 0; i < numMyElements; i++) {
01120     int tmp = MyGlobalElementValGet<int_type>(i) - BlockMapData_->MinMyGID_;
01121     assert(tmp >= 0);
01122     assert(tmp < SpanGID);
01123     BlockMapData_->LID_[MyGlobalElementValGet<int_type>(i) - BlockMapData_->MinMyGID_] = i; // Spread local indices
01124   }
01125 
01126 #endif
01127 
01128 }
01129 
01130 void Epetra_BlockMap::GlobalToLocalSetup()
01131 {
01132   if(BlockMapData_->GlobalIndicesInt_)
01133   {
01134 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
01135     TGlobalToLocalSetup<int>();
01136 #else
01137   throw ReportError("Epetra_BlockMap::GlobalToLocalSetup ERROR, GlobalIndices int but no API for it.",-1);
01138 #endif
01139   }
01140   else if(BlockMapData_->GlobalIndicesLongLong_)
01141   {
01142 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
01143     TGlobalToLocalSetup<long long>();
01144 #else
01145   throw ReportError("Epetra_BlockMap::GlobalToLocalSetup ERROR, GlobalIndices long long but no API for it.",-1);
01146 #endif
01147   }
01148   else
01149   {
01150     throw ReportError("Epetra_BlockMap::GlobalToLocalSetup ERROR, GlobalIndices type unknown.",-1);
01151   }
01152 }
01153 
01154 //==============================================================================
01155 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
01156 int Epetra_BlockMap::LID(long long gid) const
01157 {
01158   if ((gid < BlockMapData_->MinMyGID_) ||
01159     (gid > BlockMapData_->MaxMyGID_)) {
01160   return(-1); // Out of range
01161   }
01162 
01163   if (BlockMapData_->LinearMap_) {
01164   return (int) (gid - BlockMapData_->MinMyGID_); // Can compute with an offset
01165   }
01166 
01167 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
01168   if(BlockMapData_->GlobalIndicesInt_) {
01169     if( (int) gid >= BlockMapData_->MyGlobalElements_int_[0] &&
01170       (int) gid <= BlockMapData_->LastContiguousGID_ ) {
01171     return (int) gid - BlockMapData_->MyGlobalElements_int_[0];
01172     }
01173   }
01174   else
01175 #endif
01176 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
01177   if(BlockMapData_->GlobalIndicesLongLong_) {
01178     if( gid >= BlockMapData_->MyGlobalElements_LL_[0] &&
01179       gid <= BlockMapData_->LastContiguousGID_ ) {
01180     return (int) ( gid - BlockMapData_->MyGlobalElements_LL_[0] );
01181     }
01182   }
01183   else
01184 #endif
01185     throw ReportError("Epetra_BlockMap::LID ERROR, GlobalIndices type unknown.",-1);
01186 
01187 #ifdef EPETRA_BLOCKMAP_NEW_LID
01188   return BlockMapData_->LIDHash_->Get( gid );
01189 #else
01190   return(BlockMapData_->LID_[gid - BlockMapData_->MinMyGID_]); // Find it in LID array
01191 #endif
01192 }
01193 #endif
01194 
01195 //==============================================================================
01196 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
01197 int Epetra_BlockMap::LID(int gid) const
01198 {
01199   if ((gid < (int) BlockMapData_->MinMyGID_) ||
01200     (gid > (int) BlockMapData_->MaxMyGID_)) {
01201   return(-1); // Out of range
01202   }
01203 
01204   if (BlockMapData_->LinearMap_) {
01205   return (int) (gid - (int) BlockMapData_->MinMyGID_); // Can compute with an offset
01206   }
01207 
01208   if(BlockMapData_->GlobalIndicesInt_) {
01209     if( gid >= BlockMapData_->MyGlobalElements_int_[0] &&
01210       gid <= (int) BlockMapData_->LastContiguousGID_ ) {
01211     return (int) ( gid - BlockMapData_->MyGlobalElements_int_[0] );
01212     }
01213   }
01214   else if(BlockMapData_->GlobalIndicesLongLong_) {
01215   throw ReportError("Epetra_BlockMap::LID ERROR, int version called for long long map.",-1);
01216   }
01217   else {
01218   throw ReportError("Epetra_BlockMap::LID ERROR, GlobalIndices type unknown.",-1);
01219   }
01220 
01221 #ifdef EPETRA_BLOCKMAP_NEW_LID
01222   return BlockMapData_->LIDHash_->Get( gid );
01223 #else
01224   return(BlockMapData_->LID_[gid - BlockMapData_->MinMyGID_]); // Find it in LID array
01225 #endif
01226 }
01227 #endif
01228 
01229 //==============================================================================
01230 
01231 long long Epetra_BlockMap::GID64(int lid) const
01232 {
01233   if ((BlockMapData_->NumMyElements_==0) ||
01234       (lid < BlockMapData_->MinLID_) ||
01235       (lid > BlockMapData_->MaxLID_)) {
01236     return(BlockMapData_->IndexBase_ - 1); // Out of range
01237   }
01238 
01239   if (LinearMap()) {
01240     return(lid + BlockMapData_->MinMyGID_); // Can compute with an offset
01241   }
01242 
01243 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
01244   if(BlockMapData_->GlobalIndicesInt_)
01245   {
01246     return(BlockMapData_->MyGlobalElements_int_[lid]); // Find it in MyGlobalElements array
01247   }
01248   else
01249 #endif
01250 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
01251   if(BlockMapData_->GlobalIndicesLongLong_)
01252   {
01253     return(BlockMapData_->MyGlobalElements_LL_[lid]); // Find it in MyGlobalElements array
01254   }
01255   else
01256 #endif
01257   throw ReportError("Epetra_BlockMap::GID64 ERROR, GlobalIndices type unknown.",-1);
01258 }
01259 
01260 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
01261 int Epetra_BlockMap::GID(int lid) const
01262 {
01263   if(BlockMapData_->GlobalIndicesInt_)
01264   {
01265     if ((BlockMapData_->NumMyElements_==0) ||
01266         (lid < BlockMapData_->MinLID_) ||
01267         (lid > BlockMapData_->MaxLID_)) {
01268       return((int) BlockMapData_->IndexBase_ - 1); // Out of range
01269     }
01270 
01271     if (LinearMap()) {
01272       return(lid + (int) BlockMapData_->MinMyGID_); // Can compute with an offset
01273     }
01274 
01275     return(BlockMapData_->MyGlobalElements_int_[lid]); // Find it in MyGlobalElements array
01276   }
01277 
01278   throw ReportError("Epetra_BlockMap::GID ERROR, GlobalIndices type unknown or long long.",-1);
01279 }
01280 #endif
01281 
01282 //==============================================================================
01283 int Epetra_BlockMap::FindLocalElementID(int PointID, int & ElementID, int & ElementOffset) const {
01284 
01285   if (PointID >= BlockMapData_->NumMyPoints_)
01286     return(-1); // Point is out of range
01287 
01288   if (ConstantElementSize()) {
01289     ElementID = PointID / BlockMapData_->MaxElementSize_;
01290     ElementOffset = PointID % BlockMapData_->MaxElementSize_;
01291     return(0);
01292   }
01293   else {
01294     int * tmpPointToElementList = PointToElementList();
01295     int * tmpFirstPointInElementList = FirstPointInElementList();
01296     ElementID = tmpPointToElementList[PointID];
01297     ElementOffset = PointID - tmpFirstPointInElementList[ElementID];
01298     return(0);
01299   }
01300 }
01301 
01302 //==============================================================================
01303 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
01304 int Epetra_BlockMap::RemoteIDList(int NumIDs, const int * GIDList,
01305           int * PIDList, int * LIDList,
01306           int * SizeList) const
01307 {
01308   if(!BlockMapData_->GlobalIndicesInt_)
01309     throw ReportError("Epetra_BlockMap::RemoteIDList ERROR, Can't call int* version for non int* map.",-1);
01310 
01311   if (BlockMapData_->Directory_ == NULL) {
01312     BlockMapData_->Directory_ = Comm().CreateDirectory(*this);
01313   }
01314 
01315   Epetra_Directory* directory = BlockMapData_->Directory_;
01316   if (directory == NULL) {
01317     return(-1);
01318   }
01319 
01320   EPETRA_CHK_ERR( directory->GetDirectoryEntries(*this, NumIDs, GIDList,
01321              PIDList, LIDList, SizeList) );
01322 
01323   return(0);
01324 }
01325 #endif
01326 //==============================================================================
01327 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
01328 int Epetra_BlockMap::RemoteIDList(int NumIDs, const long long * GIDList,
01329           int * PIDList, int * LIDList,
01330           int * SizeList) const
01331 {
01332   if(!BlockMapData_->GlobalIndicesLongLong_)
01333     throw ReportError("Epetra_BlockMap::RemoteIDList ERROR, Can't call long long* version for non long long* map.",-1);
01334 
01335   if (BlockMapData_->Directory_ == NULL) {
01336     BlockMapData_->Directory_ = Comm().CreateDirectory(*this);
01337   }
01338 
01339   Epetra_Directory* directory = BlockMapData_->Directory_;
01340   if (directory == NULL) {
01341     return(-1);
01342   }
01343 
01344   EPETRA_CHK_ERR( directory->GetDirectoryEntries(*this, NumIDs, GIDList,
01345              PIDList, LIDList, SizeList) );
01346 
01347   return(0);
01348 }
01349 #endif
01350 
01351 //==============================================================================
01352 bool Epetra_BlockMap::DetermineIsOneToOne() const
01353 {
01354   if (Comm().NumProc() < 2) {
01355     return(true);
01356   }
01357 
01358   if (BlockMapData_->Directory_ == NULL) {
01359     BlockMapData_->Directory_ = Comm().CreateDirectory(*this);
01360   }
01361 
01362   Epetra_Directory* directory = BlockMapData_->Directory_;
01363   if (directory == NULL) {
01364     throw ReportError("Epetra_BlockMap::IsOneToOne ERROR, CreateDirectory failed.",-1);
01365   }
01366 
01367   return(directory->GIDsAllUniquelyOwned());
01368 }
01369 
01370 //==============================================================================
01371 bool Epetra_BlockMap::IsDistributedGlobal(long long numGlobalElements, int numMyElements) const {
01372 
01373   bool isDistributedGlobal = false; // Assume map is not global distributed
01374   if (BlockMapData_->Comm_->NumProc() > 1) {
01375     int LocalReplicated = 0;
01376     int AllLocalReplicated;
01377     if (numGlobalElements == numMyElements)
01378       LocalReplicated=1;
01379     BlockMapData_->Comm_->MinAll(&LocalReplicated, &AllLocalReplicated, 1);
01380 
01381     // If any PE has LocalReplicated=0, then map is distributed global
01382     if (AllLocalReplicated != 1)
01383       isDistributedGlobal = true;
01384   }
01385   return(isDistributedGlobal);
01386 }
01387 
01388 //==============================================================================
01389 void Epetra_BlockMap::CheckValidNGE(long long numGlobalElements) {
01390   // Check to see if user's value for numGlobalElements is either -1
01391   // (in which case we use our computed value) or matches ours.
01392   if ((numGlobalElements != -1) && (numGlobalElements != BlockMapData_->NumGlobalElements_)) {
01393     long long BmdNumGlobalElements = BlockMapData_->NumGlobalElements_;
01394     CleanupData();
01395     throw ReportError("Invalid NumGlobalElements.  NumGlobalElements = " + toString(numGlobalElements) +
01396           ".  Should equal " + toString(BmdNumGlobalElements) +
01397           ", or be set to -1 to compute automatically", -4);
01398   }
01399 }
01400 
01401 //==============================================================================
01402 void Epetra_BlockMap::EndOfConstructorOps() {
01403   BlockMapData_->MinLID_ = 0;
01404   BlockMapData_->MaxLID_ = EPETRA_MAX(BlockMapData_->NumMyElements_ - 1, 0);
01405 
01406   GlobalToLocalSetup(); // Setup any information for making global index to local index translation fast.
01407 }
01408 
01409 //==============================================================================
01410 void Epetra_BlockMap::Print(std::ostream & os) const
01411 {
01412   int * FirstPointInElementList1 = 0;
01413   int * ElementSizeList1 = 0;
01414   if (!ConstantElementSize()) {
01415     FirstPointInElementList1 = FirstPointInElementList();
01416     ElementSizeList1 = ElementSizeList();
01417   }
01418   int MyPID = Comm().MyPID();
01419   int NumProc = Comm().NumProc();
01420 
01421   for (int iproc = 0; iproc < NumProc; iproc++) {
01422     if (MyPID == iproc) {
01423       if (MyPID == 0) {
01424   os <<  "\nNumber of Global Elements  = "; os << NumGlobalElements64(); os << std::endl;
01425   os <<    "Number of Global Points    = "; os << NumGlobalPoints64(); os << std::endl;
01426   os <<    "Maximum of all GIDs        = "; os << MaxAllGID64(); os << std::endl;
01427   os <<    "Minimum of all GIDs        = "; os << MinAllGID64(); os << std::endl;
01428   os <<    "Index Base                 = "; os << IndexBase64(); os << std::endl;
01429   if (ConstantElementSize())
01430     os <<  "Constant Element Size      = "; os << ElementSize(); os << std::endl;
01431       }
01432       os << std::endl;
01433 
01434       os <<    "Number of Local Elements   = "; os << NumMyElements(); os << std::endl;
01435       os <<    "Number of Local Points     = "; os << NumMyPoints(); os << std::endl;
01436       os <<    "Maximum of my GIDs         = "; os << MaxMyGID64(); os << std::endl;
01437       os <<    "Minimum of my GIDs         = "; os << MinMyGID64(); os << std::endl;
01438       os << std::endl;
01439 
01440       os.width(14);
01441       os <<  "     MyPID"; os << "    ";
01442       os.width(14);
01443       os <<  "       Local Index "; os << " ";
01444       os.width(14);
01445       os <<  "      Global Index "; os << " ";
01446       if (!ConstantElementSize()) {
01447   os.width(14);
01448   os <<" FirstPointInElement "; os << " ";
01449   os.width(14);
01450   os <<"   ElementSize "; os << " ";
01451       }
01452       os << std::endl;
01453 
01454       for (int i = 0; i < NumMyElements(); i++) {
01455   os.width(14);
01456   os <<  MyPID; os << "    ";
01457   os.width(14);
01458   os <<  i; os << "    ";
01459   os.width(14);
01460 
01461   if(BlockMapData_->GlobalIndicesLongLong_)
01462   {
01463 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
01464     long long * MyGlobalElements1 = MyGlobalElements64();
01465     os <<  MyGlobalElements1[i]; os << "    ";
01466 #else
01467     throw ReportError("Epetra_BlockMap::Print: ERROR, GlobalIndicesLongLong but no API for it.",-1);
01468 #endif
01469   }
01470   else if(BlockMapData_->GlobalIndicesInt_)
01471   {
01472 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
01473     int * MyGlobalElements1 = MyGlobalElements();
01474     os <<  MyGlobalElements1[i]; os << "    ";
01475 #else
01476     throw ReportError("Epetra_BlockMap::Print: ERROR, no GlobalIndicesLongLong but no API for it.",-1);
01477 #endif
01478   }
01479 
01480   if (!ConstantElementSize()) {
01481     os.width(14);
01482     os << FirstPointInElementList1[i]; os << "    ";
01483     os.width(14);
01484     os << ElementSizeList1[i]; os << "    ";
01485   }
01486   os << std::endl;
01487       }
01488 
01489       os << std::flush;
01490 
01491     }
01492     // Do a few global ops to give I/O a chance to complete
01493     Comm().Barrier();
01494     Comm().Barrier();
01495     Comm().Barrier();
01496   }
01497   return;
01498 }
01499 
01500 //==============================================================================
01501 Epetra_BlockMap::~Epetra_BlockMap()  {
01502   CleanupData();
01503 }
01504 
01505 //==============================================================================
01506 void Epetra_BlockMap::CleanupData()
01507 {
01508   if(BlockMapData_ != 0) {
01509 
01510     BlockMapData_->DecrementReferenceCount();
01511     if(BlockMapData_->ReferenceCount() == 0) {
01512       delete BlockMapData_;
01513       BlockMapData_ = 0;
01514     }
01515   }
01516 }
01517 
01518 //=============================================================================
01519 Epetra_BlockMap & Epetra_BlockMap::operator= (const Epetra_BlockMap & map)
01520 {
01521   if((this != &map) && (BlockMapData_ != map.BlockMapData_)) {
01522     CleanupData();
01523     BlockMapData_ = map.BlockMapData_;
01524     BlockMapData_->IncrementReferenceCount();
01525   }
01526 
01527   return(*this);
01528 }
01529 
01530 //=============================================================================
01531 Epetra_BlockMap * Epetra_BlockMap::RemoveEmptyProcesses() const
01532 {
01533 #ifdef HAVE_MPI
01534   const Epetra_MpiComm * MpiComm = dynamic_cast<const Epetra_MpiComm*>(&Comm());
01535 
01536   // If the Comm isn't MPI, just treat this as a copy constructor
01537   if(!MpiComm) return new Epetra_BlockMap(*this);
01538 
01539   MPI_Comm NewComm,MyMPIComm = MpiComm->Comm();
01540 
01541   // Create the new communicator.  MPI_Comm_split returns a valid
01542   // communicator on all processes.  On processes where color == MPI_UNDEFINED,
01543   // ignore the result.  Passing key == 0 tells MPI to order the
01544   // processes in the new communicator by their rank in the old
01545   // communicator.
01546   const int color = (NumMyElements() == 0) ? MPI_UNDEFINED : 1;
01547 
01548   // MPI_Comm_split must be called collectively over the original
01549   // communicator.  We can't just call it on processes with color
01550   // one, even though we will ignore its result on processes with
01551   // color zero.
01552   int rv = MPI_Comm_split(MyMPIComm,color,0,&NewComm);
01553   if(rv!=MPI_SUCCESS) throw ReportError("Epetra_BlockMap::RemoveEmptyProcesses: MPI_Comm_split failed.",-1);
01554 
01555   if(color == MPI_UNDEFINED)
01556     return 0; // We're not in the new map
01557   else {
01558     Epetra_MpiComm * NewEpetraComm = new Epetra_MpiComm(NewComm);
01559 
01560     // Use the copy constructor for a new map, but basically because it does nothing useful
01561     Epetra_BlockMap * NewMap = new Epetra_BlockMap(*this);
01562 
01563     // Get rid of the old BlockMapData, now make a new one from scratch...
01564     NewMap->CleanupData();
01565     if(GlobalIndicesInt()) {
01566 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
01567       NewMap->BlockMapData_ = new Epetra_BlockMapData(NumGlobalElements(),0,IndexBase(),*NewEpetraComm,false);
01568 #endif
01569     }
01570     else {
01571 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
01572       NewMap->BlockMapData_ = new Epetra_BlockMapData(NumGlobalElements64(),0,IndexBase64(),*NewEpetraComm,true);
01573 #endif
01574     }
01575     // Now copy all of the relevent bits of BlockMapData...
01576     //    NewMap->BlockMapData_->Comm_                    = NewEpetraComm;
01577     NewMap->BlockMapData_->LID_                     = BlockMapData_->LID_;
01578 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
01579     NewMap->BlockMapData_->MyGlobalElements_int_    = BlockMapData_->MyGlobalElements_int_;
01580 #endif
01581 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
01582     NewMap->BlockMapData_->MyGlobalElements_LL_     = BlockMapData_->MyGlobalElements_LL_;
01583 #endif
01584     NewMap->BlockMapData_->FirstPointInElementList_ = BlockMapData_->FirstPointInElementList_;
01585     NewMap->BlockMapData_->ElementSizeList_         = BlockMapData_->ElementSizeList_;
01586     NewMap->BlockMapData_->PointToElementList_      = BlockMapData_->PointToElementList_;
01587 
01588     NewMap->BlockMapData_->NumGlobalElements_       = BlockMapData_->NumGlobalElements_;
01589     NewMap->BlockMapData_->NumMyElements_           = BlockMapData_->NumMyElements_;
01590     NewMap->BlockMapData_->IndexBase_               = BlockMapData_->IndexBase_;
01591     NewMap->BlockMapData_->ElementSize_             = BlockMapData_->ElementSize_;
01592     NewMap->BlockMapData_->MinMyElementSize_        = BlockMapData_->MinMyElementSize_;
01593     NewMap->BlockMapData_->MaxMyElementSize_        = BlockMapData_->MaxMyElementSize_;
01594     NewMap->BlockMapData_->MinElementSize_          = BlockMapData_->MinElementSize_;
01595     NewMap->BlockMapData_->MaxElementSize_          = BlockMapData_->MaxElementSize_;
01596     NewMap->BlockMapData_->MinAllGID_               = BlockMapData_->MinAllGID_;
01597     NewMap->BlockMapData_->MaxAllGID_               = BlockMapData_->MaxAllGID_;
01598     NewMap->BlockMapData_->MinMyGID_                = BlockMapData_->MinMyGID_;
01599     NewMap->BlockMapData_->MaxMyGID_                = BlockMapData_->MaxMyGID_;
01600     NewMap->BlockMapData_->MinLID_                  = BlockMapData_->MinLID_;
01601     NewMap->BlockMapData_->MaxLID_                  = BlockMapData_->MaxLID_;
01602     NewMap->BlockMapData_->NumGlobalPoints_         = BlockMapData_->NumGlobalPoints_;
01603     NewMap->BlockMapData_->NumMyPoints_             = BlockMapData_->NumMyPoints_;
01604     NewMap->BlockMapData_->ConstantElementSize_     = BlockMapData_->ConstantElementSize_;
01605     NewMap->BlockMapData_->LinearMap_               = BlockMapData_->LinearMap_;
01606     NewMap->BlockMapData_->DistributedGlobal_       = NewEpetraComm->NumProc()==1 ? false : BlockMapData_->DistributedGlobal_;
01607     NewMap->BlockMapData_->OneToOneIsDetermined_    = BlockMapData_->OneToOneIsDetermined_;
01608     NewMap->BlockMapData_->OneToOne_                = BlockMapData_->OneToOne_;
01609     NewMap->BlockMapData_->GlobalIndicesInt_        = BlockMapData_->GlobalIndicesInt_;
01610     NewMap->BlockMapData_->GlobalIndicesLongLong_   = BlockMapData_->GlobalIndicesLongLong_;
01611     NewMap->BlockMapData_->LastContiguousGID_       = BlockMapData_->LastContiguousGID_;
01612     NewMap->BlockMapData_->LastContiguousGIDLoc_    = BlockMapData_->LastContiguousGIDLoc_;
01613     NewMap->BlockMapData_->LIDHash_                 = BlockMapData_->LIDHash_ ? new Epetra_HashTable<int>(*BlockMapData_->LIDHash_) : 0;
01614 
01615     // Delay directory construction
01616     NewMap->BlockMapData_->Directory_               = 0;
01617 
01618     // Cleanup
01619     delete NewEpetraComm;
01620 
01621     return NewMap;
01622   }
01623 #else
01624     // MPI isn't compiled, so just treat this as a copy constructor
01625     return new Epetra_BlockMap(*this);
01626 #endif
01627 }
01628 
01629 //=============================================================================
01630 Epetra_BlockMap* Epetra_BlockMap::ReplaceCommWithSubset(const Epetra_Comm * theComm) const
01631 {
01632   // mfh 26 Mar 2013: The lazy way to do this is simply to recreate
01633   // the Map by calling its ordinary public constructor, using the
01634   // original Map's data.  This only involves O(1) all-reduces over
01635   // the new communicator, which in the common case only includes a
01636   // small number of processes.
01637   Epetra_BlockMap * NewMap=0;
01638 
01639   // Create the Map to return (unless theComm is zero, in which case we return zero).
01640   if (theComm) {
01641     // Map requires that the index base equal the global min GID.
01642     // Figuring out the global min GID requires a reduction over all
01643     // processes in the new communicator.  It could be that some (or
01644     // even all) of these processes contain zero entries.  (Recall
01645     // that this method, unlike removeEmptyProcesses(), may remove
01646     // an arbitrary subset of processes.)  We deal with this by
01647     // doing a min over the min GID on each process if the process
01648     // has more than zero entries, or the global max GID, if that
01649     // process has zero entries.  If no processes have any entries,
01650     // then the index base doesn't matter anyway.
01651 
01652 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
01653     if(GlobalIndicesInt()) {
01654       int MyMin, theIndexBase;
01655       MyMin  = NumMyElements() > 0 ? MinMyGID() : MaxAllGID();
01656       theComm->MinAll(&MyMin,&theIndexBase,1);
01657       NewMap = new Epetra_BlockMap(-1,NumMyElements(),MyGlobalElements(),ElementSizeList(),theIndexBase,*theComm);
01658     }
01659     else
01660 #endif
01661 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
01662     if(GlobalIndicesLongLong()) {
01663       long long MyMin, theIndexBase;
01664       MyMin = NumMyElements() > 0 ? MinMyGID64() : MaxAllGID64();
01665       theComm->MinAll(&MyMin,&theIndexBase,1);
01666       NewMap = new Epetra_BlockMap(-1,NumMyElements(),MyGlobalElements64(),ElementSizeList(),theIndexBase,*theComm);
01667     }
01668     else
01669 #endif
01670     throw ReportError("Epetra_BlockMap::ReplaceCommWithSubset ERROR, GlobalIndices type unknown.",-1);
01671   }
01672   return NewMap;
01673 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines