Epetra_BlockMap.cpp

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

Generated on Thu Sep 18 12:37:57 2008 for Epetra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1