Epetra Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Epetra_BlockMap.cpp
Go to the documentation of this file.
1 
2 //@HEADER
3 // ************************************************************************
4 //
5 // Epetra: Linear Algebra Services Package
6 // Copyright 2011 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ************************************************************************
41 //@HEADER
42 
43 #include "Epetra_ConfigDefs.h"
44 #include "Epetra_BlockMap.h"
45 #include "Epetra_Comm.h"
46 #include "Epetra_Directory.h"
48 #include "Epetra_HashTable.h"
49 #include <limits>
50 
51 #ifdef HAVE_MPI
52 #include "Epetra_MpiComm.h"
53 #endif
54 
55 // Use the new LID hash table approach by default
56 #define EPETRA_BLOCKMAP_NEW_LID
57 
58 //==============================================================================
59 // Epetra_BlockMap constructor function for a Epetra-defined uniform linear distribution of constant size elements.
60 void Epetra_BlockMap::ConstructAutoUniform(long long NumGlobal_Elements, int Element_Size, long long Index_Base, const Epetra_Comm& comm, bool IsLongLong)
61 {
62 
63  // Each processor gets roughly numGlobalPoints/p points
64  // This routine automatically defines a linear partitioning of a
65  // map with numGlobalPoints across the processors
66  // specified in the given Epetra_Comm
67 
68  if (NumGlobal_Elements < 0)
69  throw ReportError("NumGlobal_Elements = " + toString(NumGlobal_Elements) + ". Should be >= 0.", -1);
70  if (Element_Size <= 0)
71  throw ReportError("ElementSize = " + toString(Element_Size) + ". Should be > 0.", -2);
72 
73  BlockMapData_ = new Epetra_BlockMapData(NumGlobal_Elements, Element_Size, Index_Base, comm, IsLongLong);
74  int NumProc = comm.NumProc();
76  BlockMapData_->LinearMap_ = true;
77 
78  int MyPID = comm.MyPID();
79 
80  if(BlockMapData_->NumGlobalElements_ / NumProc > (long long) std::numeric_limits<int>::max())
81  throw ReportError("Epetra_BlockMap::ConstructAutoUniform: Error. Not enough space for elements on each processor", -99);
82 
84  int remainder = (int) (BlockMapData_->NumGlobalElements_ % NumProc); // remainder will fit int
85  long long start_index = ((long long)(MyPID)) * (BlockMapData_->NumMyElements_ + 1);
86 
87  if (MyPID < remainder)
89  else
90  start_index -= (MyPID - remainder);
91 
94 
99 
105 
107 }
108 
109 //==============================================================================
110 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
111 Epetra_BlockMap::Epetra_BlockMap(long long NumGlobal_Elements, int Element_Size, int Index_Base, const Epetra_Comm& comm)
112  : Epetra_Object("Epetra::BlockMap"),
113  BlockMapData_(0)
114 {
115  const bool IsLongLong = true;
116  ConstructAutoUniform(NumGlobal_Elements, Element_Size, static_cast<long long>(Index_Base), comm, IsLongLong);
117 }
118 
119 Epetra_BlockMap::Epetra_BlockMap(long long NumGlobal_Elements, int Element_Size, long long Index_Base, const Epetra_Comm& comm)
120  : Epetra_Object("Epetra::BlockMap"),
121  BlockMapData_(0)
122 {
123  const bool IsLongLong = true;
124  ConstructAutoUniform(NumGlobal_Elements, Element_Size, Index_Base, comm, IsLongLong);
125 }
126 #endif
127 
128 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
129 Epetra_BlockMap::Epetra_BlockMap(int NumGlobal_Elements, int Element_Size, int Index_Base, const Epetra_Comm& comm)
130  : Epetra_Object("Epetra::BlockMap"),
131  BlockMapData_(0)
132 {
133  const bool IsLongLong = false;
134  ConstructAutoUniform((long long)NumGlobal_Elements, Element_Size, Index_Base, comm, IsLongLong);
135 }
136 #endif
137 //==============================================================================
138 
139 // Epetra_BlockMap constructor function for a user-defined linear distribution of constant size elements.
141  long long NumGlobal_Elements, int NumMy_Elements,
142  int Element_Size, long long Index_Base, const Epetra_Comm& comm, bool IsLongLong)
143 {
144  if (NumGlobal_Elements < -1)
145  throw ReportError("NumGlobal_Elements = " + toString(NumGlobal_Elements) + ". Should be >= -1.", -1);
146  if (NumMy_Elements < 0)
147  throw ReportError("NumMy_Elements = " + toString(NumMy_Elements) + ". Should be >= 0.", -2);
148  if (Element_Size <= 0)
149  throw ReportError("ElementSize = " + toString(Element_Size) + ". Should be > 0.", -3);
150 
151  BlockMapData_ = new Epetra_BlockMapData(NumGlobal_Elements, Element_Size, Index_Base, comm, IsLongLong);
152  BlockMapData_->NumMyElements_ = NumMy_Elements;
158  BlockMapData_->LinearMap_ = true;
159 
160  // Each processor gets NumMyElements points
161 
162 
163  // Get processor information
164 
165  int NumProc = comm.NumProc();
166 
167  BlockMapData_->DistributedGlobal_ = IsDistributedGlobal(NumGlobal_Elements, NumMy_Elements);
168 
169  // Local Map and uniprocessor case: Each processor gets a complete copy of all elements
170  if (!BlockMapData_->DistributedGlobal_ || NumProc==1) {
172  CheckValidNGE(NumGlobal_Elements);
173 
176 
177  BlockMapData_-> MinAllGID_ = BlockMapData_->IndexBase_;
181  }
182  else if (NumProc > 1) {
183  // Sum up all local element counts to get global count
184  long long tmp_NumMyElements = BlockMapData_->NumMyElements_;
185  BlockMapData_->Comm_->SumAll(&tmp_NumMyElements, &BlockMapData_->NumGlobalElements_, 1);
186 
187  CheckValidNGE(NumGlobal_Elements);
188 
191 
194 
195  // Use the ScanSum function to compute a prefix sum of the number of points
196  long long tmp2_NumMyElements = BlockMapData_->NumMyElements_;
197  BlockMapData_->Comm_->ScanSum(&tmp2_NumMyElements, &BlockMapData_->MaxMyGID_, 1);
198 
199  long long start_index = BlockMapData_->MaxMyGID_ - BlockMapData_->NumMyElements_;
202  }
203  else
204  throw ReportError("Internal Error. Report to Epetra developer", -99);
205 
206 
208 }
209 
210 //==============================================================================
211 
212 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
213 Epetra_BlockMap::Epetra_BlockMap(long long NumGlobal_Elements, int NumMy_Elements,
214  int Element_Size, int Index_Base, const Epetra_Comm& comm)
215  : Epetra_Object("Epetra::BlockMap"),
216  BlockMapData_(0)
217 {
218  const bool IsLongLong = true;
219  ConstructUserLinear(NumGlobal_Elements, NumMy_Elements, Element_Size,static_cast<long long>(Index_Base), comm, IsLongLong);
220 }
221 
222 Epetra_BlockMap::Epetra_BlockMap(long long NumGlobal_Elements, int NumMy_Elements,
223  int Element_Size, long long Index_Base, const Epetra_Comm& comm)
224  : Epetra_Object("Epetra::BlockMap"),
225  BlockMapData_(0)
226 {
227  const bool IsLongLong = true;
228  ConstructUserLinear(NumGlobal_Elements, NumMy_Elements, Element_Size,Index_Base, comm, IsLongLong);
229 }
230 #endif
231 
232 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
233 Epetra_BlockMap::Epetra_BlockMap(int NumGlobal_Elements, int NumMy_Elements,
234  int Element_Size, int Index_Base, const Epetra_Comm& comm)
235  : Epetra_Object("Epetra::BlockMap"),
236  BlockMapData_(0)
237 {
238  const bool IsLongLong = false;
239  ConstructUserLinear((long long)NumGlobal_Elements, NumMy_Elements, Element_Size,Index_Base, comm, IsLongLong);
240 }
241 #endif
242 
243 // Epetra_BlockMap constructor for a user-defined arbitrary distribution of constant size elements.
244 template<typename int_type>
245 void Epetra_BlockMap::ConstructUserConstant(int_type NumGlobal_Elements, int NumMy_Elements,
246  const int_type * myGlobalElements,
247  int Element_Size, int_type indexBase,
248  const Epetra_Comm& comm, bool IsLongLong)
249 {
250  int i;
251  // Each processor gets NumMyElements points
252 
253  if (NumGlobal_Elements < -1)
254  throw ReportError("NumGlobal_Elements = " + toString(NumGlobal_Elements) + ". Should be >= -1.", -1);
255  if (NumMy_Elements < 0)
256  throw ReportError("NumMy_Elements = " + toString(NumMy_Elements) + ". Should be >= 0.", -2);
257  if (Element_Size <= 0)
258  throw ReportError("ElementSize = " + toString(Element_Size) + ". Should be > 0.", -3);
259 
260  // Allocate storage for global index list information
261 
262  BlockMapData_ = new Epetra_BlockMapData(NumGlobal_Elements, Element_Size, indexBase, comm, IsLongLong);
263  if (NumMy_Elements > 0) {
264  int errorcode = SizeMyGlobalElement<int_type>(NumMy_Elements);
265  if(errorcode != 0)
266  throw ReportError("Error with MyGlobalElements allocation.", -99);
267  }
268 
269  BlockMapData_->NumMyElements_ = NumMy_Elements;
275  BlockMapData_->LinearMap_ = false;
276  // Get processor information
277 
278  int NumProc = comm.NumProc();
279  if (NumMy_Elements > 0) {
280  // Compute min/max GID on this processor
281  BlockMapData_->MinMyGID_ = myGlobalElements[0];
282  BlockMapData_->MaxMyGID_ = myGlobalElements[0];
283  for (i = 0; i < NumMy_Elements; i++) {
284  MyGlobalElementVal<int_type>(i) = myGlobalElements[i];
285  BlockMapData_->MinMyGID_ = EPETRA_MIN(BlockMapData_->MinMyGID_, (long long) myGlobalElements[i]);
286  BlockMapData_->MaxMyGID_ = EPETRA_MAX(BlockMapData_->MaxMyGID_, (long long) myGlobalElements[i]);
287  }
288  }
289  else {
292  }
293 
294  BlockMapData_->DistributedGlobal_ = IsDistributedGlobal(NumGlobal_Elements, NumMy_Elements);
295 
296  // Local Map and uniprocessor case: Each processor gets a complete copy of all elements
297  if (!BlockMapData_->DistributedGlobal_ || NumProc==1) {
299  CheckValidNGE(NumGlobal_Elements);
302 
305  }
306  else if (NumProc > 1) {
307  // Sum up all local element counts to get global count
308  long long tmp_NumMyElements = BlockMapData_->NumMyElements_;
309  BlockMapData_->Comm_->SumAll(&tmp_NumMyElements, &BlockMapData_->NumGlobalElements_, 1);
310  CheckValidNGE(NumGlobal_Elements);
311 
314 
315  // Use the Allreduce function to find min/max GID
316  long long *tmp_send = new long long[2];
317  long long *tmp_recv = new long long[2];
318  if (BlockMapData_->NumMyElements_ > 0 ||
319  BlockMapData_->NumGlobalElements_ == 0) // This test restores behavior for empty map
320  tmp_send[0] = - BlockMapData_->MinMyGID_; // Negative sign lets us do one reduction
321  else
322  tmp_send[0] = -(std::numeric_limits<int_type>::max())-1;
323  tmp_send[1] = BlockMapData_->MaxMyGID_;
324  BlockMapData_->Comm_->MaxAll(tmp_send, tmp_recv, 2);
325  BlockMapData_->MinAllGID_ = - tmp_recv[0];
326  BlockMapData_->MaxAllGID_ = tmp_recv[1];
327  delete [] tmp_send;
328  delete [] tmp_recv;
330  throw ReportError("Minimum global element index = " + toString(BlockMapData_->MinAllGID_) +
331  " is less than index base = " + toString(BlockMapData_->IndexBase_) +".", -5);
332  }
333  else
334  throw ReportError("Internal Error. Report to Epetra developer", -99);
335 
336 
338 }
339 
340 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
341 Epetra_BlockMap::Epetra_BlockMap(long long NumGlobal_Elements, int NumMy_Elements,
342  const long long * myGlobalElements,
343  int Element_Size, int indexBase,
344  const Epetra_Comm& comm)
345  : Epetra_Object("Epetra::BlockMap"),
346  BlockMapData_(0)
347 {
348  const bool IsLongLong = true;
349  ConstructUserConstant(NumGlobal_Elements, NumMy_Elements, myGlobalElements,
350  Element_Size, static_cast<long long>(indexBase), comm, IsLongLong);
351 }
352 
353 Epetra_BlockMap::Epetra_BlockMap(long long NumGlobal_Elements, int NumMy_Elements,
354  const long long * myGlobalElements,
355  int Element_Size, long long indexBase,
356  const Epetra_Comm& comm)
357  : Epetra_Object("Epetra::BlockMap"),
358  BlockMapData_(0)
359 {
360  const bool IsLongLong = true;
361  ConstructUserConstant(NumGlobal_Elements, NumMy_Elements, myGlobalElements,
362  Element_Size, indexBase, comm, IsLongLong);
363 }
364 #endif
365 
366 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
367 Epetra_BlockMap::Epetra_BlockMap(int NumGlobal_Elements, int NumMy_Elements,
368  const int * myGlobalElements,
369  int Element_Size, int indexBase,
370  const Epetra_Comm& comm)
371  : Epetra_Object("Epetra::BlockMap"),
372  BlockMapData_(0)
373 {
374  const bool IsLongLong = false;
375  ConstructUserConstant(NumGlobal_Elements, NumMy_Elements, myGlobalElements,
376  Element_Size, indexBase, comm, IsLongLong);
377 }
378 #endif
379 
380 
381 //==============================================================================
382 // Epetra_BlockMap constructor function for a user-defined arbitrary distribution of variable size elements.
383 template<typename int_type>
384 void Epetra_BlockMap::ConstructUserVariable(int_type NumGlobal_Elements, int NumMy_Elements,
385  const int_type * myGlobalElements,
386  const int *elementSizeList, int_type indexBase,
387  const Epetra_Comm& comm, bool IsLongLong)
388 {
389 
390  int i;
391  // Each processor gets NumMyElements points
392 
393  if (NumGlobal_Elements < -1)
394  throw ReportError("NumGlobal_Elements = " + toString(NumGlobal_Elements) + ". Should be >= -1.", -1);
395  if (NumMy_Elements < 0)
396  throw ReportError("NumMy_Elements = " + toString(NumMy_Elements) + ". Should be >= 0.", -2);
397  for (i = 0; i < NumMy_Elements; i++)
398  if (elementSizeList[i] <= 0)
399  throw ReportError("elementSizeList["+toString(i)+"] = " + toString(elementSizeList[i]) + ". Should be > 0.", -3);
400 
401  BlockMapData_ = new Epetra_BlockMapData(NumGlobal_Elements, 0, indexBase, comm, IsLongLong);
402  BlockMapData_->NumMyElements_ = NumMy_Elements;
404  BlockMapData_->LinearMap_ = false;
405  // Allocate storage for global index list and element size information
406 
407  if (NumMy_Elements > 0) {
408  int errorcode = SizeMyGlobalElement<int_type>(NumMy_Elements);
409  if(errorcode != 0)
410  throw ReportError("Error with MyGlobalElements allocation.", -99);
411  errorcode = BlockMapData_->ElementSizeList_.Size(NumMy_Elements);
412  if(errorcode != 0)
413  throw ReportError("Error with ElementSizeList allocation.", -99);
414  }
415  // Get processor information
416 
417  int NumProc = comm.NumProc();
418 
419  if (NumMy_Elements > 0) {
420  // Compute min/max GID and element size, number of points on this processor
421  BlockMapData_->MinMyGID_ = myGlobalElements[0];
422  BlockMapData_->MaxMyGID_ = myGlobalElements[0];
423  BlockMapData_->MinMyElementSize_ = elementSizeList[0];
424  BlockMapData_->MaxMyElementSize_ = elementSizeList[0];
426  for (i = 0; i < NumMy_Elements; i++) {
427  MyGlobalElementVal<int_type>(i) = myGlobalElements[i];
428  BlockMapData_->ElementSizeList_[i] = elementSizeList[i];
429  BlockMapData_->MinMyGID_ = EPETRA_MIN(BlockMapData_->MinMyGID_,(long long) myGlobalElements[i]);
430  BlockMapData_->MaxMyGID_ = EPETRA_MAX(BlockMapData_->MaxMyGID_,(long long) myGlobalElements[i]);
433  BlockMapData_->NumMyPoints_ += elementSizeList[i];
434  }
435  }
436  else {
442  }
443 
444  BlockMapData_->DistributedGlobal_ = IsDistributedGlobal(NumGlobal_Elements, NumMy_Elements);
445 
446  // Local Map and uniprocessor case: Each processor gets a complete copy of all elements
447  if (!BlockMapData_->DistributedGlobal_ || NumProc == 1) {
449  CheckValidNGE(NumGlobal_Elements);
451 
456  }
457  else if (NumProc > 1) {
458  // Sum up all local element and point counts to get global counts
459  int_type *tmp_send = new int_type[4];
460  int_type *tmp_recv = new int_type[4];
461  tmp_send[0] = BlockMapData_->NumMyElements_;
462  tmp_send[1] = BlockMapData_->NumMyPoints_;
463  BlockMapData_->Comm_->SumAll(tmp_send, tmp_recv, 2);
464  BlockMapData_->NumGlobalElements_ = tmp_recv[0];
465  BlockMapData_->NumGlobalPoints_ = tmp_recv[1];
466 
467  CheckValidNGE(NumGlobal_Elements);
468 
469  // Use the MaxAll function to find min/max GID
470  if (BlockMapData_->NumMyElements_ > 0 ||
471  BlockMapData_->NumGlobalElements_ == 0) // This test restores behavior for empty map
472  tmp_send[0] = - static_cast<int_type>(BlockMapData_->MinMyGID_); // Negative sign lets us do one reduction
473  else
474  tmp_send[0] = -(std::numeric_limits<int_type>::max())-1;
475  tmp_send[1] = static_cast<int_type>(BlockMapData_->MaxMyGID_);
476  tmp_send[2] = - BlockMapData_->MinMyElementSize_;
477  if (BlockMapData_->NumMyElements_ == 0)
478  tmp_send[2] = - static_cast<int_type>(BlockMapData_->NumGlobalPoints_); // This processor has no elements, so should not sizes.
479  tmp_send[3] = BlockMapData_->MaxMyElementSize_;
480 
481  BlockMapData_->Comm_->MaxAll(tmp_send, tmp_recv, 4);
482 
483  BlockMapData_->MinAllGID_ = - tmp_recv[0];
484  BlockMapData_->MaxAllGID_ = tmp_recv[1];
485  BlockMapData_->MinElementSize_ = - (int) tmp_recv[2];
486  BlockMapData_->MaxElementSize_ = (int) tmp_recv[3];
487 
488  delete [] tmp_send;
489  delete [] tmp_recv;
490 
491  // Check for constant element size
495  }
496 
498  throw ReportError("Minimum global element index = " + toString(BlockMapData_->MinAllGID_) +
499  " is less than index base = " + toString(BlockMapData_->IndexBase_) +".", -5);
500  }
501  else
502  throw ReportError("Internal Error. Report to Epetra developer", -99);
503 
504 
506 }
507 
508 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
509 Epetra_BlockMap::Epetra_BlockMap(long long NumGlobal_Elements, int NumMy_Elements,
510  const long long * myGlobalElements,
511  const int *elementSizeList, int indexBase,
512  const Epetra_Comm& comm)
513  : Epetra_Object("Epetra::BlockMap"),
514  BlockMapData_(0)
515 {
516  const bool IsLongLong = true;
517  ConstructUserVariable(NumGlobal_Elements, NumMy_Elements, myGlobalElements,
518  elementSizeList, static_cast<long long>(indexBase), comm, IsLongLong);
519 }
520 
521 Epetra_BlockMap::Epetra_BlockMap(long long NumGlobal_Elements, int NumMy_Elements,
522  const long long * myGlobalElements,
523  const int *elementSizeList, long long indexBase,
524  const Epetra_Comm& comm)
525  : Epetra_Object("Epetra::BlockMap"),
526  BlockMapData_(0)
527 {
528  const bool IsLongLong = true;
529  ConstructUserVariable(NumGlobal_Elements, NumMy_Elements, myGlobalElements,
530  elementSizeList, indexBase, comm, IsLongLong);
531 }
532 #endif
533 
534 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
535 Epetra_BlockMap::Epetra_BlockMap(int NumGlobal_Elements, int NumMy_Elements,
536  const int * myGlobalElements,
537  const int *elementSizeList, int indexBase,
538  const Epetra_Comm& comm)
539  : Epetra_Object("Epetra::BlockMap"),
540  BlockMapData_(0)
541 {
542  const bool IsLongLong = false;
543  ConstructUserVariable(NumGlobal_Elements, NumMy_Elements, myGlobalElements,
544  elementSizeList, indexBase, comm, IsLongLong);
545 }
546 #endif
547 
548 
549 //==============================================================================
550 // Epetra_BlockMap constructor function for a user-defined arbitrary distribution of variable size elements,
551 // with all the information on globals provided by the user.
552 template<typename int_type>
553 void Epetra_BlockMap::ConstructUserConstantNoComm(int_type NumGlobal_Elements, int NumMy_Elements,
554  const int_type * myGlobalElements,
555  int Element_Size, int_type indexBase,
556  const Epetra_Comm& comm, bool IsLongLong,
557  bool UserIsDistributedGlobal,
558  int_type UserMinAllGID, int_type UserMaxAllGID)
559 {
560 
561 
562  int i;
563  // Each processor gets NumMyElements points
564 
565  if (NumGlobal_Elements < -1)
566  throw ReportError("NumGlobal_Elements = " + toString(NumGlobal_Elements) + ". Should be >= -1.", -1);
567  if (NumMy_Elements < 0)
568  throw ReportError("NumMy_Elements = " + toString(NumMy_Elements) + ". Should be >= 0.", -2);
569  if (Element_Size <= 0)
570  throw ReportError("ElementSize = " + toString(Element_Size) + ". Should be > 0.", -3);
571 
572  // Allocate storage for global index list information
573 
574  BlockMapData_ = new Epetra_BlockMapData(NumGlobal_Elements, Element_Size, indexBase, comm, IsLongLong);
575  if (NumMy_Elements > 0) {
576  int errorcode = SizeMyGlobalElement<int_type>(NumMy_Elements);
577  if(errorcode != 0)
578  throw ReportError("Error with MyGlobalElements allocation.", -99);
579  }
580 
581  BlockMapData_->NumMyElements_ = NumMy_Elements;
587  BlockMapData_->LinearMap_ = false;
588  // Get processor information
589 
590  int NumProc = comm.NumProc();
591  if (NumMy_Elements > 0) {
592  // Compute min/max GID on this processor
593  BlockMapData_->MinMyGID_ = myGlobalElements[0];
594  BlockMapData_->MaxMyGID_ = myGlobalElements[0];
595  for (i = 0; i < NumMy_Elements; i++) {
596  MyGlobalElementVal<int_type>(i) = myGlobalElements[i];
597  BlockMapData_->MinMyGID_ = EPETRA_MIN(BlockMapData_->MinMyGID_, (long long) myGlobalElements[i]);
598  BlockMapData_->MaxMyGID_ = EPETRA_MAX(BlockMapData_->MaxMyGID_, (long long) myGlobalElements[i]);
599  }
600  }
601  else {
604  }
605 
606  BlockMapData_->DistributedGlobal_ = UserIsDistributedGlobal;
607 
608  // Local Map and uniprocessor case: Each processor gets a complete copy of all elements
609  if (!BlockMapData_->DistributedGlobal_ || NumProc==1) {
611  CheckValidNGE(NumGlobal_Elements);
614 
617  }
618  else if (NumProc > 1) {
619  if(NumGlobal_Elements==-1) {
620  // Sum up all local element counts to get global count
621  long long tmp_NumMyElements = BlockMapData_->NumMyElements_;
622  BlockMapData_->Comm_->SumAll(&tmp_NumMyElements, &BlockMapData_->NumGlobalElements_, 1);
623  }
624  else {
625  // User provides this information
626  BlockMapData_->NumGlobalElements_ = NumGlobal_Elements;
627  }
629 
632 
633  BlockMapData_->MinAllGID_ = UserMinAllGID;
634  BlockMapData_->MaxAllGID_ = UserMaxAllGID;
636  throw ReportError("Minimum global element index = " + toString(BlockMapData_->MinAllGID_) +
637  " is less than index base = " + toString(BlockMapData_->IndexBase_) +".", -5);
638  }
639  else
640  throw ReportError("Internal Error. Report to Epetra developer", -99);
641 
642 
644 }
645 
646 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
647 Epetra_BlockMap::Epetra_BlockMap(long long NumGlobal_Elements, int NumMy_Elements,
648  const long long * myGlobalElements,
649  int theElementSize, int indexBase,
650  const Epetra_Comm& comm,
651  bool UserIsDistributedGlobal,
652  long long UserMinAllGID, long long UserMaxAllGID)
653  : Epetra_Object("Epetra::BlockMap"),
654  BlockMapData_(0)
655 {
656  const bool IsLongLong = true;
657  ConstructUserConstantNoComm(NumGlobal_Elements, NumMy_Elements, myGlobalElements,
658  theElementSize, (long long) indexBase, comm, IsLongLong,
659  UserIsDistributedGlobal, UserMinAllGID, UserMaxAllGID);
660 }
661 Epetra_BlockMap::Epetra_BlockMap(long long NumGlobal_Elements, int NumMy_Elements,
662  const long long * myGlobalElements,
663  int theElementSize, long long indexBase,
664  const Epetra_Comm& comm,
665  bool UserIsDistributedGlobal,
666  long long UserMinAllGID, long long UserMaxAllGID)
667  : Epetra_Object("Epetra::BlockMap"),
668  BlockMapData_(0)
669 {
670  const bool IsLongLong = true;
671  ConstructUserConstantNoComm(NumGlobal_Elements, NumMy_Elements, myGlobalElements,
672  theElementSize, indexBase, comm, IsLongLong,
673  UserIsDistributedGlobal, UserMinAllGID, UserMaxAllGID);
674 }
675 #endif
676 
677 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
678 Epetra_BlockMap::Epetra_BlockMap(int NumGlobal_Elements, int NumMy_Elements,
679  const int * myGlobalElements,
680  int theElementSize, int indexBase,
681  const Epetra_Comm& comm,
682  bool UserIsDistributedGlobal,
683  int UserMinAllGID, int UserMaxAllGID)
684  : Epetra_Object("Epetra::BlockMap"),
685  BlockMapData_(0)
686 {
687  const bool IsLongLong = false;
688  ConstructUserConstantNoComm(NumGlobal_Elements, NumMy_Elements, myGlobalElements,
689  theElementSize, indexBase, comm, IsLongLong,
690  UserIsDistributedGlobal, UserMinAllGID, UserMaxAllGID);
691 }
692 #endif
693 
694 
695 
696 
697 //==============================================================================
699  : Epetra_Object(map.Label()),
700  BlockMapData_(map.BlockMapData_)
701 {
703 
704  // This call appears to be unnecessary overhead. Removed 10-Aug-2004 maherou.
705  // GlobalToLocalSetup(); // Setup any information for making global index to local index translation fast.
706 }
707 
708 //==============================================================================
710 
711  // Quickest test: See if both maps share an inner data class
712  if (this->BlockMapData_ == Map.BlockMapData_)
713  return(true);
714  return false;
715 }
716 
717 //==============================================================================
718 bool Epetra_BlockMap::SameAs(const Epetra_BlockMap & Map) const {
719 
720  // Quickest test: See if both maps share an inner data class
721  if (this->BlockMapData_ == Map.BlockMapData_)
722  return(true);
723 
724  if(!GlobalIndicesTypeMatch(Map))
725  return(false);
726 
727  // Next check other global properties that are easy global attributes
728  if (BlockMapData_->MinAllGID_ != Map.MinAllGID64() ||
732  return(false);
733 
734  // Last possible global check for constant element sizes
736  return(false);
737 
738  // If we get this far, we need to check local properties and then check across
739  // all processors to see if local properties are all true
740 
741  int numMyElements = BlockMapData_->NumMyElements_;
742 
743  int MySameMap = 1; // Assume not needed
744 
745  // First check if number of element is the same in each map
746  if (numMyElements != Map.NumMyElements()) MySameMap = 0;
747 
748  // If numMyElements is the same, check to see that list of GIDs is the same
749  if (MySameMap==1) {
750  if (LinearMap() && Map.LinearMap() ) {
751  // For linear maps, just need to check whether lower bound is the same
752  if (MinMyGID64() != Map.MinMyGID64() )
753  MySameMap = 0;
754  }
755  else {
756  for (int i = 0; i < numMyElements; i++) {
757  if (GID64(i) != Map.GID64(i)) {
758  MySameMap = 0;
759  break;
760  }
761  }
762  }
763  }
764 // for (int i = 0; i < numMyElements; i++)
765 // if (GID64(i) != Map.GID64(i)) MySameMap = 0;
766 
767  // If GIDs are the same, check to see element sizes are the same
768  if (MySameMap==1 && !BlockMapData_->ConstantElementSize_) {
769  int * sizeList1 = ElementSizeList();
770  int * sizeList2 = Map.ElementSizeList();
771  for (int i = 0; i < numMyElements; i++) if (sizeList1[i] != sizeList2[i]) MySameMap=0;
772  }
773  // Now get min of MySameMap across all processors
774 
775  int GlobalSameMap = 0;
776 #ifdef NDEBUG
777  (void) Comm().MinAll(&MySameMap, &GlobalSameMap, 1);
778 #else
779  // Don't declare 'err' unless we actually use it (in the assert(),
780  // which gets defined away in a release build).
781  int err = Comm().MinAll(&MySameMap, &GlobalSameMap, 1);
782  assert(err==0);
783 #endif // NDEBUG
784  return(GlobalSameMap==1);
785 }
786 
787 //==============================================================================
789 {
790  if (this->BlockMapData_ == Map.BlockMapData_)
791  return(true);
792 
793  if(!GlobalIndicesTypeMatch(Map))
794  return(false);
795 
797  return(false);
798 
799  // If we get this far, we need to check local properties and then check across
800  // all processors to see if local properties are all true
801 
802  int MySameMap = 1; // Assume not needed
803  if (BlockMapData_->NumMyPoints_ != Map.NumMyPoints())
804  MySameMap = 0;
805 
806  int GlobalSameMap = 0;
807 
808 #ifdef NDEBUG
809  (void) Comm().MinAll(&MySameMap, &GlobalSameMap, 1);
810 #else
811  // Don't declare 'err' unless we actually use it (in the assert(),
812  // which gets defined away in a release build).
813  int err = Comm().MinAll(&MySameMap, &GlobalSameMap, 1);
814  assert( err == 0 );
815 #endif // NDEBUG
816 
817  return(GlobalSameMap==1);
818 }
819 
820 //==============================================================================
821 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
822 int Epetra_BlockMap::MyGlobalElements(long long * myGlobalElements) const
823 {
824  // Although one can populate long long data from int data, we don't
825  // allow it to maintain int/long long symmetry.
827  throw ReportError("Epetra_BlockMap::MyGlobalElements(long long *) ERROR, Can't call for non long long* map.",-1);
828 
829  // If the global element list is not create, then do so. This can only happen when
830  // a linear distribution has been specified. Thus we can easily construct the update
831  // list in this case.
832 
833  int i;
834  int numMyElements = BlockMapData_->NumMyElements_;
835 
837  for (i = 0; i < numMyElements; i++)
838  myGlobalElements[i] = BlockMapData_->MinMyGID_ + i;
839  else
840  for (i = 0; i < numMyElements; i++)
841  myGlobalElements[i] = BlockMapData_->MyGlobalElements_LL_[i];
842  return(0);
843 }
844 #endif
845 
846 //==============================================================================
847 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
848 int Epetra_BlockMap::MyGlobalElements(int * myGlobalElements) const
849 {
851  throw ReportError("Epetra_BlockMap::MyGlobalElements(int *) ERROR, Can't call for non int* map.",-1);
852 
853  // If the global element list is not create, then do so. This can only happen when
854  // a linear distribution has been specified. Thus we can easily construct the update
855  // list in this case.
856 
857  int i;
858  int numMyElements = BlockMapData_->NumMyElements_;
859 
861  for (i = 0; i < numMyElements; i++)
862  myGlobalElements[i] = (int) BlockMapData_->MinMyGID_ + i;
863  else
864  for (i = 0; i < numMyElements; i++)
865  myGlobalElements[i] = (int) BlockMapData_->MyGlobalElements_int_[i];
866  return(0);
867 }
868 #endif
869 
870 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
871 int Epetra_BlockMap::MyGlobalElementsPtr(long long *& MyGlobalElementList) const
872 {
873  MyGlobalElementList = MyGlobalElements64();
874  return(0);
875 }
876 #endif
877 
878 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
879 int Epetra_BlockMap::MyGlobalElementsPtr(int *& MyGlobalElementList) const
880 {
881  MyGlobalElementList = MyGlobalElements();
882  return(0);
883 }
884 #endif
885 //==============================================================================
886 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
889  throw ReportError("Epetra_BlockMap::MyGlobalElements() ERROR, Can't call for non int* map.",-1);
890 
891  int numMyElements = BlockMapData_->NumMyElements_;
892 
893  // If ElementSizeList not built, do so
894  if(BlockMapData_->MyGlobalElements_int_.Length() == 0 && numMyElements > 0) {
895  int errorcode = BlockMapData_->MyGlobalElements_int_.Size(numMyElements + 1);
896  if(errorcode != 0)
897  throw ReportError("Error with MyGlobalElements allocation.", -99);
898 
899  // Build the array
900  for (int i = 0; i < numMyElements; i++)
902  }
904 }
905 #endif
906 //==============================================================================
907 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
910  throw ReportError("Epetra_BlockMap::MyGlobalElements64 ERROR, Can't call for non long long* map.",-1);
911 
912  int numMyElements = BlockMapData_->NumMyElements_;
913 
914  // If ElementSizeList not built, do so
915  if(BlockMapData_->MyGlobalElements_LL_.Length() == 0 && numMyElements > 0) {
916  int errorcode = BlockMapData_->MyGlobalElements_LL_.Size(numMyElements + 1);
917  if(errorcode != 0)
918  throw ReportError("Error with MyGlobalElements allocation.", -99);
919 
920  // Build the array
921  for (int i = 0; i < numMyElements; i++)
923  }
925 }
926 #endif
927 //==============================================================================
929 {
930  if (!MyLID(lid))
931  EPETRA_CHK_ERR(-1);
932 
933  int entry;
934 
935  if (ConstantElementSize())
936  entry = MaxElementSize() * lid; // convert to vector entry
937  else {
938  int * entrylist = FirstPointInElementList(); // get entry list
939  entry = entrylist[lid];
940  }
941  return(entry);
942 }
943 
944 //==============================================================================
945 int Epetra_BlockMap::FirstPointInElementList(int * firstPointInElementList) const
946 {
947  // If the first element entry list is not create, then do so.
948 
949  // Note: This array is of length NumMyElement+1
950 
951  int i;
952  int numMyElements = BlockMapData_->NumMyElements_;
953 
955  firstPointInElementList[0] = 0; // First element of first entry is always zero
956 
958  for (i = 0; i < numMyElements; i++)
959  firstPointInElementList[i+1] = firstPointInElementList[i] + BlockMapData_->ElementSize_;
960  else
961  for (i = 0; i < numMyElements; i++)
962  firstPointInElementList[i+1] = firstPointInElementList[i] + BlockMapData_->ElementSizeList_[i];
963  }
964  else
965  for (i = 0; i <= numMyElements; i++)
966  firstPointInElementList[i] = BlockMapData_->FirstPointInElementList_[i];
967  return(0);
968 }
969 
970 //==============================================================================
972  int numMyElements = BlockMapData_->NumMyElements_;
973 
974  // If ElementSizeList not built, do so
975  if ((BlockMapData_->FirstPointInElementList_.Length() == 0) && (numMyElements > 0)) {
977  BlockMapData_->FirstPointInElementList_[0] = 0; // First element of first entry is always zero
979  for (int i = 0; i < numMyElements; i++)
981  else
982  for (int i = 0; i < numMyElements; i++)
984  }
986 }
987 
988 //==============================================================================
989 int Epetra_BlockMap::ElementSizeList(int * elementSizeList) const
990 {
991  // If the element size list is not create, then do so. This can only happen when
992  // a constant element size has been specified. Thus we can easily construct the element size
993  // list in this case.
994 
995  int i;
996  int numMyElements = BlockMapData_->NumMyElements_;
997 
999  for (i = 0; i < numMyElements; i++)
1000  elementSizeList[i] = BlockMapData_->ElementSize_;
1001  else
1002  for (i = 0; i < numMyElements; i++)
1003  elementSizeList[i] = BlockMapData_->ElementSizeList_[i];
1004 
1005  return(0);
1006 }
1007 
1008 //==============================================================================
1010  int numMyElements = BlockMapData_->NumMyElements_;
1011 
1012  // If ElementSizeList not built, do so
1013  if ((BlockMapData_->ElementSizeList_.Length() == 0) && (numMyElements > 0)) {
1014  BlockMapData_->ElementSizeList_.Size(numMyElements);
1015  for (int i = 0; i < numMyElements; i++)
1017  }
1019 }
1020 
1021 //==============================================================================
1022 int Epetra_BlockMap::PointToElementList(int * pointToElementList) const {
1023  // Build an array such that the local element ID is stored for each point
1024 
1025  int i;
1027  int numMyElements = BlockMapData_->NumMyElements_;
1028  int * ptr = pointToElementList;
1029  for (i = 0; i < numMyElements; i++) {
1030  int Size = ElementSize(i);
1031  for (int j = 0; j < Size; j++)
1032  *ptr++ = i;
1033  }
1034  }
1035  else {
1036  int numMyPoints = BlockMapData_->NumMyPoints_;
1037  for (i = 0; i < numMyPoints; i++)
1038  pointToElementList[i] = BlockMapData_->PointToElementList_[i];
1039  }
1040  return(0);
1041 }
1042 
1043 //==============================================================================
1045 
1046  // If PointToElementList not built, do so
1049  int numMyElements = BlockMapData_->NumMyElements_;
1050  int * ptr = BlockMapData_->PointToElementList_.Values();
1051  for (int i = 0; i < numMyElements; i++) {
1052  int Size = ElementSize(i);
1053  for (int j = 0; j < Size; j++)
1054  *ptr++ = i;
1055  }
1056  }
1058 }
1059 
1060 //==============================================================================
1061 int Epetra_BlockMap::ElementSize(int lid) const {
1062 
1063  if (ConstantElementSize())
1064  return(BlockMapData_->ElementSize_);
1065  else
1066  return(BlockMapData_->ElementSizeList_[lid]);
1067 }
1068 
1069 //==============================================================================
1074  }
1075  return(BlockMapData_->OneToOne_);
1076 }
1077 
1078 //==============================================================================
1079 template<typename int_type>
1081 {
1082  int i;
1083  int numMyElements = BlockMapData_->NumMyElements_;
1084 
1085  if (BlockMapData_->NumGlobalElements_ == 0) {
1086  return; // Nothing to do
1087  }
1088 
1089  if (LinearMap() || numMyElements == 0) {
1090  return; // Nothing else to do
1091  }
1092 
1093  // Build LID_ vector to make look up of local index values fast
1094 
1095 #ifdef EPETRA_BLOCKMAP_NEW_LID
1096 
1097  // Here follows an optimization that checks for an initial block of
1098  // contiguous GIDs, and stores the GID->LID mapping for those in a
1099  // more efficient way. This supports a common use case for
1100  // overlapping Maps, in which owned entries of a vector are ordered
1101  // before halo entries.
1102  //
1103  // Epetra defines EPETRA_BLOCKMAP_NEW_LID by default (see the top of
1104  // this file).
1105 
1106  //check for initial contiguous block
1107  int_type val = MyGlobalElementValGet<int_type>(0);
1108  for( i = 0 ; i < numMyElements; ++i ) {
1109  if (val != MyGlobalElementValGet<int_type>(i)) break;
1110  ++val;
1111  }
1114  BlockMapData_->LastContiguousGID_ = MyGlobalElementValGet<int_type>(0);
1115  }
1116  else {
1118  MyGlobalElementValGet<int_type>(BlockMapData_->LastContiguousGIDLoc_);
1119  }
1120 
1121  //Hash everything else
1122  if(i < numMyElements) {
1123  if (BlockMapData_->LIDHash_ != NULL) {
1124  delete BlockMapData_->LIDHash_;
1125  }
1126 
1127  BlockMapData_->LIDHash_ = new Epetra_HashTable<int>(numMyElements - i + 1 );
1128  for(; i < numMyElements; ++i )
1129  BlockMapData_->LIDHash_->Add( MyGlobalElementValGet<int_type>(i), i );
1130  }
1131 
1132 #else
1133 
1134  int SpanGID = BlockMapData_->MaxMyGID_ - BlockMapData_->MinMyGID_ + 1;
1135  BlockMapData_->LID_.Size(SpanGID);
1136 
1137  for (i = 0; i < SpanGID; i++)
1138  BlockMapData_->LID_[i] = -1; // Fill all locations with -1
1139 
1140  for (i = 0; i < numMyElements; i++) {
1141  int tmp = MyGlobalElementValGet<int_type>(i) - BlockMapData_->MinMyGID_;
1142  assert(tmp >= 0);
1143  assert(tmp < SpanGID);
1144  BlockMapData_->LID_[MyGlobalElementValGet<int_type>(i) - BlockMapData_->MinMyGID_] = i; // Spread local indices
1145  }
1146 
1147 #endif
1148 
1149 }
1150 
1152 {
1154  {
1155 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
1156  TGlobalToLocalSetup<int>();
1157 #else
1158  throw ReportError("Epetra_BlockMap::GlobalToLocalSetup ERROR, GlobalIndices int but no API for it.",-1);
1159 #endif
1160  }
1162  {
1163 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1164  TGlobalToLocalSetup<long long>();
1165 #else
1166  throw ReportError("Epetra_BlockMap::GlobalToLocalSetup ERROR, GlobalIndices long long but no API for it.",-1);
1167 #endif
1168  }
1169  else
1170  {
1171  throw ReportError("Epetra_BlockMap::GlobalToLocalSetup ERROR, GlobalIndices type unknown.",-1);
1172  }
1173 }
1174 
1175 //==============================================================================
1176 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1177 int Epetra_BlockMap::LID(long long gid) const
1178 {
1179  if ((gid < BlockMapData_->MinMyGID_) ||
1180  (gid > BlockMapData_->MaxMyGID_)) {
1181  return(-1); // Out of range
1182  }
1183 
1184  if (BlockMapData_->LinearMap_) {
1185  return (int) (gid - BlockMapData_->MinMyGID_); // Can compute with an offset
1186  }
1187 
1188 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
1190  if( (int) gid >= BlockMapData_->MyGlobalElements_int_[0] &&
1191  (int) gid <= BlockMapData_->LastContiguousGID_ ) {
1192  return (int) gid - BlockMapData_->MyGlobalElements_int_[0];
1193  }
1194  }
1195  else
1196 #endif
1197 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1199  if( gid >= BlockMapData_->MyGlobalElements_LL_[0] &&
1200  gid <= BlockMapData_->LastContiguousGID_ ) {
1201  return (int) ( gid - BlockMapData_->MyGlobalElements_LL_[0] );
1202  }
1203  }
1204  else
1205 #endif
1206  throw ReportError("Epetra_BlockMap::LID ERROR, GlobalIndices type unknown.",-1);
1207 
1208 #ifdef EPETRA_BLOCKMAP_NEW_LID
1209  return BlockMapData_->LIDHash_->Get( gid );
1210 #else
1211  return(BlockMapData_->LID_[gid - BlockMapData_->MinMyGID_]); // Find it in LID array
1212 #endif
1213 }
1214 #endif
1215 
1216 //==============================================================================
1217 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
1218 int Epetra_BlockMap::LID(int gid) const
1219 {
1220  if ((gid < (int) BlockMapData_->MinMyGID_) ||
1221  (gid > (int) BlockMapData_->MaxMyGID_)) {
1222  return(-1); // Out of range
1223  }
1224 
1225  if (BlockMapData_->LinearMap_) {
1226  return (int) (gid - (int) BlockMapData_->MinMyGID_); // Can compute with an offset
1227  }
1228 
1230  if( gid >= BlockMapData_->MyGlobalElements_int_[0] &&
1231  gid <= (int) BlockMapData_->LastContiguousGID_ ) {
1232  return (int) ( gid - BlockMapData_->MyGlobalElements_int_[0] );
1233  }
1234  }
1236  throw ReportError("Epetra_BlockMap::LID ERROR, int version called for long long map.",-1);
1237  }
1238  else {
1239  throw ReportError("Epetra_BlockMap::LID ERROR, GlobalIndices type unknown.",-1);
1240  }
1241 
1242 #ifdef EPETRA_BLOCKMAP_NEW_LID
1243  return BlockMapData_->LIDHash_->Get( gid );
1244 #else
1245  return(BlockMapData_->LID_[gid - BlockMapData_->MinMyGID_]); // Find it in LID array
1246 #endif
1247 }
1248 #endif
1249 
1250 //==============================================================================
1251 
1252 long long Epetra_BlockMap::GID64(int lid) const
1253 {
1254  if ((BlockMapData_->NumMyElements_==0) ||
1255  (lid < BlockMapData_->MinLID_) ||
1256  (lid > BlockMapData_->MaxLID_)) {
1257  return(BlockMapData_->IndexBase_ - 1); // Out of range
1258  }
1259 
1260  if (LinearMap()) {
1261  return(lid + BlockMapData_->MinMyGID_); // Can compute with an offset
1262  }
1263 
1264 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
1266  {
1267  return(BlockMapData_->MyGlobalElements_int_[lid]); // Find it in MyGlobalElements array
1268  }
1269  else
1270 #endif
1271 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1273  {
1274  return(BlockMapData_->MyGlobalElements_LL_[lid]); // Find it in MyGlobalElements array
1275  }
1276  else
1277 #endif
1278  throw ReportError("Epetra_BlockMap::GID64 ERROR, GlobalIndices type unknown.",-1);
1279 }
1280 
1281 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
1282 int Epetra_BlockMap::GID(int lid) const
1283 {
1285  {
1286  if ((BlockMapData_->NumMyElements_==0) ||
1287  (lid < BlockMapData_->MinLID_) ||
1288  (lid > BlockMapData_->MaxLID_)) {
1289  return((int) BlockMapData_->IndexBase_ - 1); // Out of range
1290  }
1291 
1292  if (LinearMap()) {
1293  return(lid + (int) BlockMapData_->MinMyGID_); // Can compute with an offset
1294  }
1295 
1296  return(BlockMapData_->MyGlobalElements_int_[lid]); // Find it in MyGlobalElements array
1297  }
1298 
1299  throw ReportError("Epetra_BlockMap::GID ERROR, GlobalIndices type unknown or long long.",-1);
1300 }
1301 #endif
1302 
1303 //==============================================================================
1304 int Epetra_BlockMap::FindLocalElementID(int PointID, int & ElementID, int & ElementOffset) const {
1305 
1306  if (PointID >= BlockMapData_->NumMyPoints_)
1307  return(-1); // Point is out of range
1308 
1309  if (ConstantElementSize()) {
1310  ElementID = PointID / BlockMapData_->MaxElementSize_;
1311  ElementOffset = PointID % BlockMapData_->MaxElementSize_;
1312  return(0);
1313  }
1314  else {
1315  int * tmpPointToElementList = PointToElementList();
1316  int * tmpFirstPointInElementList = FirstPointInElementList();
1317  ElementID = tmpPointToElementList[PointID];
1318  ElementOffset = PointID - tmpFirstPointInElementList[ElementID];
1319  return(0);
1320  }
1321 }
1322 
1323 //==============================================================================
1324 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
1325 int Epetra_BlockMap::RemoteIDList(int NumIDs, const int * GIDList,
1326  int * PIDList, int * LIDList,
1327  int * SizeList) const
1328 {
1330  throw ReportError("Epetra_BlockMap::RemoteIDList ERROR, Can't call int* version for non int* map.",-1);
1331 
1332  if (BlockMapData_->Directory_ == NULL) {
1334  }
1335 
1337  if (directory == NULL) {
1338  return(-1);
1339  }
1340 
1341  EPETRA_CHK_ERR( directory->GetDirectoryEntries(*this, NumIDs, GIDList,
1342  PIDList, LIDList, SizeList) );
1343 
1344  return(0);
1345 }
1346 #endif
1347 //==============================================================================
1348 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1349 int Epetra_BlockMap::RemoteIDList(int NumIDs, const long long * GIDList,
1350  int * PIDList, int * LIDList,
1351  int * SizeList) const
1352 {
1354  throw ReportError("Epetra_BlockMap::RemoteIDList ERROR, Can't call long long* version for non long long* map.",-1);
1355 
1356  if (BlockMapData_->Directory_ == NULL) {
1358  }
1359 
1361  if (directory == NULL) {
1362  return(-1);
1363  }
1364 
1365  EPETRA_CHK_ERR( directory->GetDirectoryEntries(*this, NumIDs, GIDList,
1366  PIDList, LIDList, SizeList) );
1367 
1368  return(0);
1369 }
1370 #endif
1371 
1372 //==============================================================================
1374 {
1375  if (Comm().NumProc() < 2) {
1376  return(true);
1377  }
1378 
1379  if (BlockMapData_->Directory_ == NULL) {
1381  }
1382 
1384  if (directory == NULL) {
1385  throw ReportError("Epetra_BlockMap::IsOneToOne ERROR, CreateDirectory failed.",-1);
1386  }
1387 
1388  return(directory->GIDsAllUniquelyOwned());
1389 }
1390 
1391 //==============================================================================
1392 bool Epetra_BlockMap::IsDistributedGlobal(long long numGlobalElements, int numMyElements) const {
1393 
1394  bool isDistributedGlobal = false; // Assume map is not global distributed
1395  if (BlockMapData_->Comm_->NumProc() > 1) {
1396  int LocalReplicated = 0;
1397  int AllLocalReplicated;
1398  if (numGlobalElements == numMyElements)
1399  LocalReplicated=1;
1400  BlockMapData_->Comm_->MinAll(&LocalReplicated, &AllLocalReplicated, 1);
1401 
1402  // If any PE has LocalReplicated=0, then map is distributed global
1403  if (AllLocalReplicated != 1)
1404  isDistributedGlobal = true;
1405  }
1406  return(isDistributedGlobal);
1407 }
1408 
1409 //==============================================================================
1410 void Epetra_BlockMap::CheckValidNGE(long long numGlobalElements) {
1411  // Check to see if user's value for numGlobalElements is either -1
1412  // (in which case we use our computed value) or matches ours.
1413  if ((numGlobalElements != -1) && (numGlobalElements != BlockMapData_->NumGlobalElements_)) {
1414  long long BmdNumGlobalElements = BlockMapData_->NumGlobalElements_;
1415  CleanupData();
1416  throw ReportError("Invalid NumGlobalElements. NumGlobalElements = " + toString(numGlobalElements) +
1417  ". Should equal " + toString(BmdNumGlobalElements) +
1418  ", or be set to -1 to compute automatically", -4);
1419  }
1420 }
1421 
1422 //==============================================================================
1424  BlockMapData_->MinLID_ = 0;
1426 
1427  GlobalToLocalSetup(); // Setup any information for making global index to local index translation fast.
1428 }
1429 
1430 //==============================================================================
1431 void Epetra_BlockMap::Print(std::ostream & os) const
1432 {
1433  int * FirstPointInElementList1 = 0;
1434  int * ElementSizeList1 = 0;
1435  if (!ConstantElementSize()) {
1436  FirstPointInElementList1 = FirstPointInElementList();
1437  ElementSizeList1 = ElementSizeList();
1438  }
1439  int MyPID = Comm().MyPID();
1440  int NumProc = Comm().NumProc();
1441 
1442  for (int iproc = 0; iproc < NumProc; iproc++) {
1443  if (MyPID == iproc) {
1444  if (MyPID == 0) {
1445  os << "\nNumber of Global Elements = "; os << NumGlobalElements64(); os << std::endl;
1446  os << "Number of Global Points = "; os << NumGlobalPoints64(); os << std::endl;
1447  os << "Maximum of all GIDs = "; os << MaxAllGID64(); os << std::endl;
1448  os << "Minimum of all GIDs = "; os << MinAllGID64(); os << std::endl;
1449  os << "Index Base = "; os << IndexBase64(); os << std::endl;
1450  if (ConstantElementSize()) {
1451  os << "Constant Element Size = "; os << ElementSize(); os << std::endl;
1452  }
1453  }
1454  os << std::endl;
1455 
1456  os << "Number of Local Elements = "; os << NumMyElements(); os << std::endl;
1457  os << "Number of Local Points = "; os << NumMyPoints(); os << std::endl;
1458  os << "Maximum of my GIDs = "; os << MaxMyGID64(); os << std::endl;
1459  os << "Minimum of my GIDs = "; os << MinMyGID64(); os << std::endl;
1460  os << std::endl;
1461 
1462  os.width(14);
1463  os << " MyPID"; os << " ";
1464  os.width(14);
1465  os << " Local Index "; os << " ";
1466  os.width(14);
1467  os << " Global Index "; os << " ";
1468  if (!ConstantElementSize()) {
1469  os.width(14);
1470  os <<" FirstPointInElement "; os << " ";
1471  os.width(14);
1472  os <<" ElementSize "; os << " ";
1473  }
1474  os << std::endl;
1475 
1476  for (int i = 0; i < NumMyElements(); i++) {
1477  os.width(14);
1478  os << MyPID; os << " ";
1479  os.width(14);
1480  os << i; os << " ";
1481  os.width(14);
1482 
1484  {
1485 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1486  long long * MyGlobalElements1 = MyGlobalElements64();
1487  os << MyGlobalElements1[i]; os << " ";
1488 #else
1489  throw ReportError("Epetra_BlockMap::Print: ERROR, GlobalIndicesLongLong but no API for it.",-1);
1490 #endif
1491  }
1493  {
1494 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
1495  int * MyGlobalElements1 = MyGlobalElements();
1496  os << MyGlobalElements1[i]; os << " ";
1497 #else
1498  throw ReportError("Epetra_BlockMap::Print: ERROR, no GlobalIndicesLongLong but no API for it.",-1);
1499 #endif
1500  }
1501 
1502  if (!ConstantElementSize()) {
1503  os.width(14);
1504  os << FirstPointInElementList1[i]; os << " ";
1505  os.width(14);
1506  os << ElementSizeList1[i]; os << " ";
1507  }
1508  os << std::endl;
1509  }
1510 
1511  os << std::flush;
1512 
1513  }
1514  // Do a few global ops to give I/O a chance to complete
1515  Comm().Barrier();
1516  Comm().Barrier();
1517  Comm().Barrier();
1518  }
1519  return;
1520 }
1521 
1522 //==============================================================================
1524  CleanupData();
1525 }
1526 
1527 //==============================================================================
1529 {
1530  if(BlockMapData_ != 0) {
1531 
1533  if(BlockMapData_->ReferenceCount() == 0) {
1534  delete BlockMapData_;
1535  BlockMapData_ = 0;
1536  }
1537  }
1538 }
1539 
1540 //=============================================================================
1542 {
1543  if((this != &map) && (BlockMapData_ != map.BlockMapData_)) {
1544  CleanupData();
1547  }
1548 
1549  return(*this);
1550 }
1551 
1552 //=============================================================================
1554 {
1555 #ifdef HAVE_MPI
1556  const Epetra_MpiComm * MpiComm = dynamic_cast<const Epetra_MpiComm*>(&Comm());
1557 
1558  // If the Comm isn't MPI, just treat this as a copy constructor
1559  if(!MpiComm) return new Epetra_BlockMap(*this);
1560 
1561  MPI_Comm NewComm,MyMPIComm = MpiComm->Comm();
1562 
1563  // Create the new communicator. MPI_Comm_split returns a valid
1564  // communicator on all processes. On processes where color == MPI_UNDEFINED,
1565  // ignore the result. Passing key == 0 tells MPI to order the
1566  // processes in the new communicator by their rank in the old
1567  // communicator.
1568  const int color = (NumMyElements() == 0) ? MPI_UNDEFINED : 1;
1569 
1570  // MPI_Comm_split must be called collectively over the original
1571  // communicator. We can't just call it on processes with color
1572  // one, even though we will ignore its result on processes with
1573  // color zero.
1574  int rv = MPI_Comm_split(MyMPIComm,color,0,&NewComm);
1575  if(rv!=MPI_SUCCESS) throw ReportError("Epetra_BlockMap::RemoveEmptyProcesses: MPI_Comm_split failed.",-1);
1576 
1577  if(color == MPI_UNDEFINED)
1578  return 0; // We're not in the new map
1579  else {
1580  Epetra_MpiComm * NewEpetraComm = new Epetra_MpiComm(NewComm);
1581 
1582  // Use the copy constructor for a new map, but basically because it does nothing useful
1583  Epetra_BlockMap * NewMap = new Epetra_BlockMap(*this);
1584 
1585  // Get rid of the old BlockMapData, now make a new one from scratch...
1586  NewMap->CleanupData();
1587  if(GlobalIndicesInt()) {
1588 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
1589  NewMap->BlockMapData_ = new Epetra_BlockMapData(NumGlobalElements(),0,IndexBase(),*NewEpetraComm,false);
1590 #endif
1591  }
1592  else {
1593 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1594  NewMap->BlockMapData_ = new Epetra_BlockMapData(NumGlobalElements64(),0,IndexBase64(),*NewEpetraComm,true);
1595 #endif
1596  }
1597  // Now copy all of the relevent bits of BlockMapData...
1598  // NewMap->BlockMapData_->Comm_ = NewEpetraComm;
1599  NewMap->BlockMapData_->LID_ = BlockMapData_->LID_;
1600 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
1602 #endif
1603 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1605 #endif
1609 
1628  NewMap->BlockMapData_->DistributedGlobal_ = NewEpetraComm->NumProc()==1 ? false : BlockMapData_->DistributedGlobal_;
1636 
1637  // Delay directory construction
1638  NewMap->BlockMapData_->Directory_ = 0;
1639 
1640  // Cleanup
1641  delete NewEpetraComm;
1642 
1643  return NewMap;
1644  }
1645 #else
1646  // MPI isn't compiled, so just treat this as a copy constructor
1647  return new Epetra_BlockMap(*this);
1648 #endif
1649 }
1650 
1651 //=============================================================================
1653 {
1654  // mfh 26 Mar 2013: The lazy way to do this is simply to recreate
1655  // the Map by calling its ordinary public constructor, using the
1656  // original Map's data. This only involves O(1) all-reduces over
1657  // the new communicator, which in the common case only includes a
1658  // small number of processes.
1659  Epetra_BlockMap * NewMap=0;
1660 
1661  // Create the Map to return (unless theComm is zero, in which case we return zero).
1662  if (theComm) {
1663  // Map requires that the index base equal the global min GID.
1664  // Figuring out the global min GID requires a reduction over all
1665  // processes in the new communicator. It could be that some (or
1666  // even all) of these processes contain zero entries. (Recall
1667  // that this method, unlike removeEmptyProcesses(), may remove
1668  // an arbitrary subset of processes.) We deal with this by
1669  // doing a min over the min GID on each process if the process
1670  // has more than zero entries, or the global max GID, if that
1671  // process has zero entries. If no processes have any entries,
1672  // then the index base doesn't matter anyway.
1673 
1674 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
1675  if(GlobalIndicesInt()) {
1676  int MyMin, theIndexBase;
1677  MyMin = NumMyElements() > 0 ? MinMyGID() : MaxAllGID();
1678  theComm->MinAll(&MyMin,&theIndexBase,1);
1679  NewMap = new Epetra_BlockMap(-1,NumMyElements(),MyGlobalElements(),ElementSizeList(),theIndexBase,*theComm);
1680  }
1681  else
1682 #endif
1683 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1684  if(GlobalIndicesLongLong()) {
1685  long long MyMin, theIndexBase;
1686  MyMin = NumMyElements() > 0 ? MinMyGID64() : MaxAllGID64();
1687  theComm->MinAll(&MyMin,&theIndexBase,1);
1688  NewMap = new Epetra_BlockMap(-1,NumMyElements(),MyGlobalElements64(),ElementSizeList(),theIndexBase,*theComm);
1689  }
1690  else
1691 #endif
1692  throw ReportError("Epetra_BlockMap::ReplaceCommWithSubset ERROR, GlobalIndices type unknown.",-1);
1693  }
1694  return NewMap;
1695 }
long long MinMyGID64() const
Epetra_IntSerialDenseVector PointToElementList_
bool PointSameAs(const Epetra_BlockMap &Map) const
Returns true if this and Map have identical point-wise structure.
int NumGlobalElements() const
Number of elements across all processors.
Epetra_BlockMapData: The Epetra BlockMap Data Class.
int * FirstPointInElementList() const
Pointer to internal array containing a mapping between the local elements and the first local point n...
long long MaxAllGID64() const
Epetra_IntSerialDenseVector ElementSizeList_
int ElementSize() const
Returns the size of elements in the map; only valid if map has constant element size.
void ConstructUserConstantNoComm(int_type NumGlobal_Elements, int NumMy_Elements, const int_type *myGlobalElements, int ElementSize, int_type indexBase, const Epetra_Comm &comm, bool IsLongLong, bool UserIsDistributedGlobal, int_type UserMinAllGID, int_type UserMaxAllGID)
const Epetra_Comm * Comm_
int Length() const
Returns length of vector.
Epetra_BlockMapData * BlockMapData_
bool GlobalIndicesLongLong() const
Returns true if map create with long long NumGlobalElements.
bool SameAs(const Epetra_BlockMap &Map) const
Returns true if this and Map are identical maps.
long long * Values()
Returns pointer to the values in vector.
long long IndexBase64() const
void DecrementReferenceCount()
Decrement reference count.
Definition: Epetra_Data.cpp:66
virtual void Print(std::ostream &os) const
Print object to an output stream.
long long NumGlobalElements64() const
value_type Get(const long long key)
int Size(int Length_in)
Set length of a Epetra_IntSerialDenseVector object; init values to zero.
virtual Epetra_Directory * CreateDirectory(const Epetra_BlockMap &Map) const =0
Create a directory object for the given Epetra_BlockMap.
int NumProc() const
Returns total number of processes.
bool ConstantElementSize() const
Returns true if map has constant element size.
virtual bool GIDsAllUniquelyOwned() const =0
GIDsAllUniquelyOwned: returns true if all GIDs appear on just one processor.
Epetra_HashTable< int > * LIDHash_
#define EPETRA_CHK_ERR(a)
int * ElementSizeList() const
List of the element sizes corresponding to the array MyGlobalElements().
Epetra_IntSerialDenseVector MyGlobalElements_int_
Epetra_BlockMap * RemoveEmptyProcesses() const
Return a new BlockMap with processes with zero elements removed.
#define EPETRA_MIN(x, y)
bool GlobalIndicesInt() const
Returns true if map create with int NumGlobalElements.
long long NumGlobalPoints64() const
virtual int MinAll(double *PartialMins, double *GlobalMins, int Count) const =0
Epetra_Comm Global Min function.
Epetra_MpiComm: The Epetra MPI Communication Class.
virtual int GetDirectoryEntries(const Epetra_BlockMap &Map, const int NumEntries, const int *GlobalEntries, int *Procs, int *LocalEntries, int *EntrySizes, bool high_rank_sharing_procs=false) const =0
GetDirectoryEntries : Returns proc and local id info for non-local map entries.
bool IsOneToOne() const
virtual void Barrier() const =0
Epetra_Comm Barrier function.
virtual int MyPID() const =0
Return my process ID.
int Length() const
Returns length of vector.
int IndexBase() const
Index base for this map.
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const =0
Epetra_Comm Global Max function.
int Size(int Length_in)
Set length of a Epetra_LongLongSerialDenseVector object; init values to zero.
void CheckValidNGE(long long NumGlobalElements)
bool SameBlockMapDataAs(const Epetra_BlockMap &Map) const
Returns true if maps share same block map data underneath.
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
Epetra_Comm Global Sum function.
int NumMyElements() const
Number of elements on the calling processor.
Epetra_Directory: This class is a pure virtual class whose interface allows Epetra_Map and Epetr_Bloc...
long long GID64(int LID) const
std::string toString(const int &x) const
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition: Epetra_Comm.h:81
bool DetermineIsOneToOne() const
long long MinAllGID64() const
Epetra_Directory * Directory_
int FirstPointInElement(int LID) const
Returns the requested entry in the FirstPointInElementList; see FirstPointInElementList() for details...
Epetra_Object: The base Epetra class.
Definition: Epetra_Object.h:65
bool IsDistributedGlobal(long long NumGlobalElements, int NumMyElements) const
int GID(int LID) const
Returns global ID of local ID, return IndexBase-1 if not found on this processor. ...
Epetra_BlockMap & operator=(const Epetra_BlockMap &map)
Assignment Operator.
int ReferenceCount() const
Get reference count.
Definition: Epetra_Data.cpp:71
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
MPI_Comm Comm() const
Extract MPI Communicator from a Epetra_MpiComm object.
int MyGlobalElementsPtr(int *&MyGlobalElementList) const
Epetra_IntSerialDenseVector LID_
int FindLocalElementID(int PointID, int &ElementID, int &ElementOffset) const
Returns the LID of the element that contains the given local PointID, and the Offset of the point in ...
int LID(int GID) const
Returns local ID of global ID, return -1 if not found on this processor.
Epetra_LongLongSerialDenseVector MyGlobalElements_LL_
Epetra_BlockMap * ReplaceCommWithSubset(const Epetra_Comm *Comm) const
Replace this BlockMap&#39;s communicator with a subset communicator.
void IncrementReferenceCount()
Increment reference count.
Definition: Epetra_Data.cpp:61
int MinMyGID() const
Returns the minimum global ID owned by this processor.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
bool LinearMap() const
Returns true if the global ID space is contiguously divided (but not necessarily uniformly) across al...
void ConstructUserLinear(long long NumGlobal_Elements, int NumMy_Elements, int Element_Size, long long Index_Base, const Epetra_Comm &comm, bool IsLongLong)
int * PointToElementList() const
For each local point, indicates the local element ID that the point belongs to.
virtual int NumProc() const =0
Returns total number of processes.
void ConstructUserVariable(int_type NumGlobal_Elements, int NumMy_Elements, const int_type *myGlobalElements, const int *elementSizeList, int_type indexBase, const Epetra_Comm &comm, bool IsLongLong)
Epetra_BlockMap(int NumGlobalElements, int ElementSize, int IndexBase, const Epetra_Comm &Comm)
Epetra_BlockMap constructor for a Epetra-defined uniform linear distribution of constant size element...
virtual ~Epetra_BlockMap(void)
Epetra_BlockMap destructor.
void ConstructAutoUniform(long long NumGlobal_Elements, int Element_Size, long long Index_Base, const Epetra_Comm &comm, bool IsLongLong)
int MaxElementSize() const
Maximum element size across all processors.
Epetra_IntSerialDenseVector FirstPointInElementList_
int * MyGlobalElements() const
Pointer to internal array containing list of global IDs assigned to the calling processor.
virtual int ReportError(const std::string Message, int ErrorCode) const
Error reporting method.
bool MyLID(int lid) const
Returns true if the LID passed in belongs to the calling processor in this map, otherwise returns fal...
int NumMyPoints() const
Number of local points for this map; equals the sum of all element sizes on the calling processor...
long long * MyGlobalElements64() const
int * Values()
Returns pointer to the values in vector.
virtual int ScanSum(double *MyVals, double *ScanSums, int Count) const =0
Epetra_Comm Scan Sum function.
void Add(const long long key, const value_type value)
long long MaxMyGID64() const
#define EPETRA_MAX(x, y)
void ConstructUserConstant(int_type NumGlobal_Elements, int NumMy_Elements, const int_type *myGlobalElements, int Element_Size, int_type indexBase, const Epetra_Comm &comm, bool IsLongLong)
int RemoteIDList(int NumIDs, const int *GIDList, int *PIDList, int *LIDList) const
Returns the processor IDs and corresponding local index value for a given list of global indices...
bool GlobalIndicesTypeMatch(const Epetra_BlockMap &other) const
int MaxAllGID() const
Returns the maximum global ID across the entire map.