Epetra Package Browser (Single Doxygen Collection)  Development
Epetra_Util.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_Util.h"
45 #include "Epetra_Object.h"
46 #include "Epetra_Comm.h"
47 #include "Epetra_Directory.h"
48 #include "Epetra_Map.h"
49 #include "Epetra_LocalMap.h"
50 #include "Epetra_CrsGraph.h"
51 #include "Epetra_CrsMatrix.h"
52 #include "Epetra_MultiVector.h"
53 #include "Epetra_IntVector.h"
54 #include "Epetra_GIDTypeVector.h"
55 #include "Epetra_Import.h"
56 #include <limits>
57 
58 #ifdef HAVE_MPI
59 #include "Epetra_MpiDistributor.h"
60 #endif
61 
62 const double Epetra_Util::chopVal_ = 1.0e-15;
63 
64 //=========================================================================
65 double Epetra_Util::Chop(const double & Value){
66  if (std::abs(Value) < chopVal_) return 0;
67  return Value;
68 }
69 
70 //=========================================================================
71 unsigned int Epetra_Util::RandomInt() {
72 
73  const int a = 16807;
74  const int m = 2147483647;
75  const int q = 127773;
76  const int r = 2836;
77 
78  int hi = Seed_ / q;
79  int lo = Seed_ % q;
80  int test = a * lo - r * hi;
81  if (test > 0)
82  Seed_ = test;
83  else
84  Seed_ = test + m;
85 
86  return(Seed_);
87 }
88 
89 //=========================================================================
91  const double Modulus = 2147483647.0;
92  const double DbleOne = 1.0;
93  const double DbleTwo = 2.0;
94 
95  double randdouble = RandomInt(); // implicit conversion from int to double
96  randdouble = DbleTwo * (randdouble / Modulus) - DbleOne; // scale to (-1.0, 1.0)
97 
98  return(randdouble);
99 }
100 
101 //=========================================================================
102 unsigned int Epetra_Util::Seed() const {
103  return(Seed_);
104 }
105 
106 //=========================================================================
107 int Epetra_Util::SetSeed(unsigned int Seed_in) {
108  Seed_ = Seed_in;
109  return(0);
110 }
111 
112 //=============================================================================
113 template<typename T>
114 void Epetra_Util::Sort(bool SortAscending, int NumKeys, T * Keys,
115  int NumDoubleCompanions,double ** DoubleCompanions,
116  int NumIntCompanions, int ** IntCompanions,
117  int NumLongLongCompanions, long long ** LongLongCompanions)
118 {
119  int i;
120 
121  int n = NumKeys;
122  T * const list = Keys;
123  int m = n/2;
124 
125  while (m > 0) {
126  int max = n - m;
127  for (int j=0; j<max; j++)
128  {
129  for (int k=j; k>=0; k-=m)
130  {
131  if ((SortAscending && list[k+m] >= list[k]) ||
132  ( !SortAscending && list[k+m] <= list[k]))
133  break;
134  T temp = list[k+m];
135  list[k+m] = list[k];
136  list[k] = temp;
137  for (i=0; i<NumDoubleCompanions; i++) {
138  double dtemp = DoubleCompanions[i][k+m];
139  DoubleCompanions[i][k+m] = DoubleCompanions[i][k];
140  DoubleCompanions[i][k] = dtemp;
141  }
142  for (i=0; i<NumIntCompanions; i++) {
143  int itemp = IntCompanions[i][k+m];
144  IntCompanions[i][k+m] = IntCompanions[i][k];
145  IntCompanions[i][k] = itemp;
146  }
147  for (i=0; i<NumLongLongCompanions; i++) {
148  long long LLtemp = LongLongCompanions[i][k+m];
149  LongLongCompanions[i][k+m] = LongLongCompanions[i][k];
150  LongLongCompanions[i][k] = LLtemp;
151  }
152  }
153  }
154  m = m/2;
155  }
156 }
157 
158 void Epetra_Util::Sort(bool SortAscending, int NumKeys, int * Keys,
159  int NumDoubleCompanions,double ** DoubleCompanions,
160  int NumIntCompanions, int ** IntCompanions,
161  int NumLongLongCompanions, long long ** LongLongCompanions)
162 {
163  Sort<int>(SortAscending, NumKeys, Keys,
164  NumDoubleCompanions, DoubleCompanions,
165  NumIntCompanions, IntCompanions,
166  NumLongLongCompanions, LongLongCompanions);
167 }
168 
169 void Epetra_Util::Sort(bool SortAscending, int NumKeys, long long * Keys,
170  int NumDoubleCompanions,double ** DoubleCompanions,
171  int NumIntCompanions, int ** IntCompanions,
172  int NumLongLongCompanions, long long ** LongLongCompanions)
173 {
174  Sort<long long>(SortAscending, NumKeys, Keys,
175  NumDoubleCompanions, DoubleCompanions,
176  NumIntCompanions, IntCompanions,
177  NumLongLongCompanions, LongLongCompanions);
178 }
179 
180 void Epetra_Util::Sort(bool SortAscending, int NumKeys, double * Keys,
181  int NumDoubleCompanions,double ** DoubleCompanions,
182  int NumIntCompanions, int ** IntCompanions,
183  int NumLongLongCompanions, long long ** LongLongCompanions)
184 {
185  Sort<double>(SortAscending, NumKeys, Keys,
186  NumDoubleCompanions, DoubleCompanions,
187  NumIntCompanions, IntCompanions,
188  NumLongLongCompanions, LongLongCompanions);
189 }
190 
191 
192 void Epetra_Util::Sort(bool SortAscending, int NumKeys, int * Keys,
193  int NumDoubleCompanions,double ** DoubleCompanions,
194  int NumIntCompanions, int ** IntCompanions)
195 {
196  int i;
197 
198  int n = NumKeys;
199  int * const list = Keys;
200  int m = n/2;
201 
202  while (m > 0) {
203  int max = n - m;
204  for (int j=0; j<max; j++)
205  {
206  for (int k=j; k>=0; k-=m)
207  {
208  if ((SortAscending && list[k+m] >= list[k]) ||
209  ( !SortAscending && list[k+m] <= list[k]))
210  break;
211  int temp = list[k+m];
212  list[k+m] = list[k];
213  list[k] = temp;
214  for (i=0; i<NumDoubleCompanions; i++) {
215  double dtemp = DoubleCompanions[i][k+m];
216  DoubleCompanions[i][k+m] = DoubleCompanions[i][k];
217  DoubleCompanions[i][k] = dtemp;
218  }
219  for (i=0; i<NumIntCompanions; i++) {
220  int itemp = IntCompanions[i][k+m];
221  IntCompanions[i][k+m] = IntCompanions[i][k];
222  IntCompanions[i][k] = itemp;
223  }
224  }
225  }
226  m = m/2;
227  }
228 }
229 
230 //----------------------------------------------------------------------------
231 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES // FIXME
232 // FIXME long long
235  bool high_rank_proc_owns_shared)
236 {
237  //if usermap is already 1-to-1 then we'll just return a copy of it.
238  if (usermap.IsOneToOne()) {
239  Epetra_Map newmap(usermap);
240  return(newmap);
241  }
242 
243  int myPID = usermap.Comm().MyPID();
244  Epetra_Directory* directory = usermap.Comm().CreateDirectory(usermap);
245 
246  int numMyElems = usermap.NumMyElements();
247  const int* myElems = usermap.MyGlobalElements();
248 
249  int* owner_procs = new int[numMyElems];
250 
251  directory->GetDirectoryEntries(usermap, numMyElems, myElems, owner_procs,
252  0, 0, high_rank_proc_owns_shared);
253 
254  //we'll fill a list of map-elements which belong on this processor
255 
256  int* myOwnedElems = new int[numMyElems];
257  int numMyOwnedElems = 0;
258 
259  for(int i=0; i<numMyElems; ++i) {
260  int GID = myElems[i];
261  int owner = owner_procs[i];
262 
263  if (myPID == owner) {
264  myOwnedElems[numMyOwnedElems++] = GID;
265  }
266  }
267 
268  Epetra_Map one_to_one_map(-1, numMyOwnedElems, myOwnedElems,
269  usermap.IndexBase(), usermap.Comm()); // CJ TODO FIXME long long
270 
271  delete [] myOwnedElems;
272  delete [] owner_procs;
273  delete directory;
274 
275  return(one_to_one_map);
276 }
277 #endif // EPETRA_NO_32BIT_GLOBAL_INDICES
278 //----------------------------------------------------------------------------
279 
280 template<typename int_type>
281 static Epetra_Map TCreate_Root_Map(const Epetra_Map& usermap,
282  int root)
283 {
284  int numProc = usermap.Comm().NumProc();
285  if (numProc==1) {
286  Epetra_Map newmap(usermap);
287  return(newmap);
288  }
289 
290  const Epetra_Comm & comm = usermap.Comm();
291  bool isRoot = usermap.Comm().MyPID()==root;
292 
293  //if usermap is already completely owned by root then we'll just return a copy of it.
294  int quickreturn = 0;
295  int globalquickreturn = 0;
296 
297  if (isRoot) {
298  if (usermap.NumMyElements()==usermap.NumGlobalElements64()) quickreturn = 1;
299  }
300  else {
301  if (usermap.NumMyElements()==0) quickreturn = 1;
302  }
303  usermap.Comm().MinAll(&quickreturn, &globalquickreturn, 1);
304 
305  if (globalquickreturn==1) {
306  Epetra_Map newmap(usermap);
307  return(newmap);
308  }
309 
310  // Linear map: Simple case, just put all GIDs linearly on root processor
311  if (usermap.LinearMap() && root!=-1) {
312  int numMyElements = 0;
313  if(usermap.MaxAllGID64()+1 > std::numeric_limits<int>::max())
314  throw "Epetra_Util::Create_Root_Map: cannot fit all gids in int";
315  if (isRoot) numMyElements = (int)(usermap.MaxAllGID64()+1);
316  Epetra_Map newmap((int_type) -1, numMyElements, (int_type)usermap.IndexBase64(), comm);
317  return(newmap);
318  }
319 
320  if (!usermap.UniqueGIDs())
321  throw usermap.ReportError("usermap must have unique GIDs",-1);
322 
323  // General map
324 
325  // Build IntVector of the GIDs, then ship them to root processor
326  int numMyElements = usermap.NumMyElements();
327  Epetra_Map allGidsMap((int_type) -1, numMyElements, (int_type) 0, comm);
328  typename Epetra_GIDTypeVector<int_type>::impl allGids(allGidsMap);
329  for (int i=0; i<numMyElements; i++) allGids[i] = (int_type) usermap.GID64(i);
330 
331  if(usermap.MaxAllGID64() > std::numeric_limits<int>::max())
332  throw "Epetra_Util::Create_Root_Map: cannot fit all gids in int";
333  int numGlobalElements = (int) usermap.NumGlobalElements64();
334  if (root!=-1) {
335  int n1 = 0; if (isRoot) n1 = numGlobalElements;
336  Epetra_Map allGidsOnRootMap((int_type) -1, n1, (int_type) 0, comm);
337  Epetra_Import importer(allGidsOnRootMap, allGidsMap);
338  typename Epetra_GIDTypeVector<int_type>::impl allGidsOnRoot(allGidsOnRootMap);
339  allGidsOnRoot.Import(allGids, importer, Insert);
340 
341  Epetra_Map rootMap((int_type)-1, allGidsOnRoot.MyLength(), allGidsOnRoot.Values(), (int_type)usermap.IndexBase64(), comm);
342  return(rootMap);
343  }
344  else {
345  int n1 = numGlobalElements;
346  Epetra_LocalMap allGidsOnRootMap((int_type) n1, (int_type) 0, comm);
347  Epetra_Import importer(allGidsOnRootMap, allGidsMap);
348  typename Epetra_GIDTypeVector<int_type>::impl allGidsOnRoot(allGidsOnRootMap);
349  allGidsOnRoot.Import(allGids, importer, Insert);
350 
351  Epetra_Map rootMap((int_type) -1, allGidsOnRoot.MyLength(), allGidsOnRoot.Values(), (int_type)usermap.IndexBase64(), comm);
352 
353  return(rootMap);
354  }
355 }
356 
358  int root)
359 {
360 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
361  if(usermap.GlobalIndicesInt()) {
362  return TCreate_Root_Map<int>(usermap, root);
363  }
364  else
365 #endif
366 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
367  if(usermap.GlobalIndicesLongLong()) {
368  return TCreate_Root_Map<long long>(usermap, root);
369  }
370  else
371 #endif
372  throw "Epetra_Util::Create_Root_Map: GlobalIndices type unknown";
373 }
374 
375 //----------------------------------------------------------------------------
376 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES // FIXME
377 // FIXME long long
380  bool high_rank_proc_owns_shared)
381 {
382 // FIXME long long
383 
384  //if usermap is already 1-to-1 then we'll just return a copy of it.
385  if (usermap.IsOneToOne()) {
386  Epetra_BlockMap newmap(usermap);
387  return(newmap);
388  }
389 
390  int myPID = usermap.Comm().MyPID();
391  Epetra_Directory* directory = usermap.Comm().CreateDirectory(usermap);
392 
393  int numMyElems = usermap.NumMyElements();
394  const int* myElems = usermap.MyGlobalElements();
395 
396  int* owner_procs = new int[numMyElems*2];
397  int* sizes = owner_procs+numMyElems;
398 
399  directory->GetDirectoryEntries(usermap, numMyElems, myElems, owner_procs,
400  0, sizes, high_rank_proc_owns_shared);
401 
402  //we'll fill a list of map-elements which belong on this processor
403 
404  int* myOwnedElems = new int[numMyElems*2];
405  int* ownedSizes = myOwnedElems+numMyElems;
406  int numMyOwnedElems = 0;
407 
408  for(int i=0; i<numMyElems; ++i) {
409  int GID = myElems[i];
410  int owner = owner_procs[i];
411 
412  if (myPID == owner) {
413  ownedSizes[numMyOwnedElems] = sizes[i];
414  myOwnedElems[numMyOwnedElems++] = GID;
415  }
416  }
417 
418  Epetra_BlockMap one_to_one_map(-1, numMyOwnedElems, myOwnedElems,
419  sizes, usermap.IndexBase(), usermap.Comm()); // CJ TODO FIXME long long
420 
421  delete [] myOwnedElems;
422  delete [] owner_procs;
423  delete directory;
424 
425  return(one_to_one_map);
426 }
427 #endif // EPETRA_NO_32BIT_GLOBAL_INDICES
428 
429 
430 //----------------------------------------------------------------------------
431 int Epetra_Util::SortCrsEntries(int NumRows, const int *CRS_rowptr, int *CRS_colind, double *CRS_vals){
432  // For each row, sort column entries from smallest to largest.
433  // Use shell sort. Stable sort so it is fast if indices are already sorted.
434  // Code copied from Epetra_CrsMatrix::SortEntries()
435  int nnz = CRS_rowptr[NumRows];
436 
437  for(int i = 0; i < NumRows; i++){
438  int start=CRS_rowptr[i];
439  if(start >= nnz) continue;
440 
441  double* locValues = &CRS_vals[start];
442  int NumEntries = CRS_rowptr[i+1] - start;
443  int* locIndices = &CRS_colind[start];
444 
445  int n = NumEntries;
446  int m = n/2;
447 
448  while(m > 0) {
449  int max = n - m;
450  for(int j = 0; j < max; j++) {
451  for(int k = j; k >= 0; k-=m) {
452  if(locIndices[k+m] >= locIndices[k])
453  break;
454  double dtemp = locValues[k+m];
455  locValues[k+m] = locValues[k];
456  locValues[k] = dtemp;
457  int itemp = locIndices[k+m];
458  locIndices[k+m] = locIndices[k];
459  locIndices[k] = itemp;
460  }
461  }
462  m = m/2;
463  }
464  }
465  return(0);
466 }
467 
468 //----------------------------------------------------------------------------
469 int Epetra_Util::SortCrsEntries(int NumRows, const size_t *CRS_rowptr, int *CRS_colind, double *CRS_vals){
470  // For each row, sort column entries from smallest to largest.
471  // Use shell sort. Stable sort so it is fast if indices are already sorted.
472  // Code copied from Epetra_CrsMatrix::SortEntries()
473  size_t nnz = CRS_rowptr[NumRows];
474 
475  for(int i = 0; i < NumRows; i++){
476  size_t start = CRS_rowptr[i];
477  if(start >= nnz) continue;
478 
479  double* locValues = &CRS_vals[start];
480  int NumEntries = static_cast<int>(CRS_rowptr[i+1] - start);
481  int* locIndices = &CRS_colind[start];
482 
483  int n = NumEntries;
484  int m = n/2;
485 
486  while(m > 0) {
487  int max = n - m;
488  for(int j = 0; j < max; j++) {
489  for(int k = j; k >= 0; k-=m) {
490  if(locIndices[k+m] >= locIndices[k])
491  break;
492  double dtemp = locValues[k+m];
493  locValues[k+m] = locValues[k];
494  locValues[k] = dtemp;
495  int itemp = locIndices[k+m];
496  locIndices[k+m] = locIndices[k];
497  locIndices[k] = itemp;
498  }
499  }
500  m = m/2;
501  }
502  }
503  return(0);
504 }
505 
506 //----------------------------------------------------------------------------
507 int Epetra_Util::SortAndMergeCrsEntries(int NumRows, int *CRS_rowptr, int *CRS_colind, double *CRS_vals){
508  // For each row, sort column entries from smallest to largest, merging column ids that are identical by adding values.
509  // Use shell sort. Stable sort so it is fast if indices are already sorted.
510  // Code copied from Epetra_CrsMatrix::SortEntries()
511 
512  int nnz = CRS_rowptr[NumRows];
513  int new_curr=CRS_rowptr[0], old_curr=CRS_rowptr[0];
514 
515  for(int i = 0; i < NumRows; i++){
516  int start=CRS_rowptr[i];
517  if(start >= nnz) continue;
518 
519  double* locValues = &CRS_vals[start];
520  int NumEntries = CRS_rowptr[i+1] - start;
521  int* locIndices = &CRS_colind[start];
522 
523  // Sort phase
524  int n = NumEntries;
525  int m = n/2;
526 
527  while(m > 0) {
528  int max = n - m;
529  for(int j = 0; j < max; j++) {
530  for(int k = j; k >= 0; k-=m) {
531  if(locIndices[k+m] >= locIndices[k])
532  break;
533  double dtemp = locValues[k+m];
534  locValues[k+m] = locValues[k];
535  locValues[k] = dtemp;
536  int itemp = locIndices[k+m];
537  locIndices[k+m] = locIndices[k];
538  locIndices[k] = itemp;
539  }
540  }
541  m = m/2;
542  }
543 
544  // Merge & shrink
545  for(int j=CRS_rowptr[i]; j < CRS_rowptr[i+1]; j++) {
546  if(j > CRS_rowptr[i] && CRS_colind[j]==CRS_colind[new_curr-1]) {
547  CRS_vals[new_curr-1] += CRS_vals[j];
548  }
549  else if(new_curr==j) {
550  new_curr++;
551  }
552  else {
553  CRS_colind[new_curr] = CRS_colind[j];
554  CRS_vals[new_curr] = CRS_vals[j];
555  new_curr++;
556  }
557  }
558 
559  CRS_rowptr[i] = old_curr;
560  old_curr=new_curr;
561  }
562 
563  CRS_rowptr[NumRows] = new_curr;
564  return (0);
565 }
566 
567 //----------------------------------------------------------------------------
568 int Epetra_Util::SortAndMergeCrsEntries(int NumRows, size_t *CRS_rowptr, int *CRS_colind, double *CRS_vals){
569  // For each row, sort column entries from smallest to largest, merging column ids that are identical by adding values.
570  // Use shell sort. Stable sort so it is fast if indices are already sorted.
571  // Code copied from Epetra_CrsMatrix::SortEntries()
572 
573  size_t nnz = CRS_rowptr[NumRows];
574  size_t new_curr=CRS_rowptr[0], old_curr=CRS_rowptr[0];
575 
576  for(int i = 0; i < NumRows; i++){
577  size_t start=CRS_rowptr[i];
578  if(start >= nnz) continue;
579 
580  double* locValues = &CRS_vals[start];
581  int NumEntries = static_cast<int>(CRS_rowptr[i+1] - start);
582  int* locIndices = &CRS_colind[start];
583 
584  // Sort phase
585  int n = NumEntries;
586  int m = n/2;
587 
588  while(m > 0) {
589  int max = n - m;
590  for(int j = 0; j < max; j++) {
591  for(int k = j; k >= 0; k-=m) {
592  if(locIndices[k+m] >= locIndices[k])
593  break;
594  double dtemp = locValues[k+m];
595  locValues[k+m] = locValues[k];
596  locValues[k] = dtemp;
597  int itemp = locIndices[k+m];
598  locIndices[k+m] = locIndices[k];
599  locIndices[k] = itemp;
600  }
601  }
602  m = m/2;
603  }
604 
605  // Merge & shrink
606  for(size_t j=CRS_rowptr[i]; j < CRS_rowptr[i+1]; j++) {
607  if(j > CRS_rowptr[i] && CRS_colind[j]==CRS_colind[new_curr-1]) {
608  CRS_vals[new_curr-1] += CRS_vals[j];
609  }
610  else if(new_curr==j) {
611  new_curr++;
612  }
613  else {
614  CRS_colind[new_curr] = CRS_colind[j];
615  CRS_vals[new_curr] = CRS_vals[j];
616  new_curr++;
617  }
618  }
619 
620  CRS_rowptr[i] = old_curr;
621  old_curr=new_curr;
622  }
623 
624  CRS_rowptr[NumRows] = new_curr;
625  return (0);
626 }
627 
628 
629 //----------------------------------------------------------------------------
630 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
631 int Epetra_Util::GetPidGidPairs(const Epetra_Import & Importer,std::vector< std::pair<int,int> > & gpids, bool use_minus_one_for_local){
632  // Put the (PID,GID) pair in member of Importer.TargetMap() in gpids. If use_minus_one_for_local==true, put in -1 instead of MyPID.
633  // This only works if we have an MpiDistributor in our Importer. Otherwise return an error.
634 #ifdef HAVE_MPI
635  Epetra_MpiDistributor *D=dynamic_cast<Epetra_MpiDistributor*>(&Importer.Distributor());
636  if(!D) EPETRA_CHK_ERR(-2);
637 
638  int i,j,k;
639  int mypid=Importer.TargetMap().Comm().MyPID();
640  int N=Importer.TargetMap().NumMyElements();
641 
642  // Get the importer's data
643  const int *RemoteLIDs = Importer.RemoteLIDs();
644 
645  // Get the distributor's data
646  int NumReceives = D->NumReceives();
647  const int *ProcsFrom = D->ProcsFrom();
648  const int *LengthsFrom = D->LengthsFrom();
649 
650  // Resize the outgoing data structure
651  gpids.resize(N);
652 
653  // Start by claiming that I own all the data
654  if(use_minus_one_for_local)
655  for(i=0;i <N; i++) gpids[i]=std::make_pair(-1,Importer.TargetMap().GID(i));
656  else
657  for(i=0;i <N; i++) gpids[i]=std::make_pair(mypid,Importer.TargetMap().GID(i));
658 
659  // Now, for each remote ID, record who actually owns it. This loop follows the operation order in the
660  // MpiDistributor so it ought to duplicate that effect.
661  for(i=0,j=0;i<NumReceives;i++){
662  int pid=ProcsFrom[i];
663  for(k=0;k<LengthsFrom[i];k++){
664  if(pid!=mypid) gpids[RemoteLIDs[j]].first=pid;
665  j++;
666  }
667  }
668  return 0;
669 #else
670  EPETRA_CHK_ERR(-10);
671  // The above macro may not necessarily invoke 'return', or at least
672  // GCC 4.8.2 can't tell if it does. That's why I put an extra
673  // return statement here.
674  return 0;
675 #endif
676 }
677 #endif
678 
679 //----------------------------------------------------------------------------
680 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
681 int Epetra_Util::GetPidGidPairs(const Epetra_Import & Importer,std::vector< std::pair<int,long long> > & gpids, bool use_minus_one_for_local){
682  // Put the (PID,GID) pair in member of Importer.TargetMap() in gpids. If use_minus_one_for_local==true, put in -1 instead of MyPID.
683  // This only works if we have an MpiDistributor in our Importer. Otheriwise return an error.
684 #ifdef HAVE_MPI
685  Epetra_MpiDistributor *D=dynamic_cast<Epetra_MpiDistributor*>(&Importer.Distributor());
686  if(!D) EPETRA_CHK_ERR(-2);
687 
688  int i,j,k;
689  int mypid=Importer.TargetMap().Comm().MyPID();
690  int N=Importer.TargetMap().NumMyElements();
691 
692  // Get the importer's data
693  const int *RemoteLIDs = Importer.RemoteLIDs();
694 
695  // Get the distributor's data
696  int NumReceives = D->NumReceives();
697  const int *ProcsFrom = D->ProcsFrom();
698  const int *LengthsFrom = D->LengthsFrom();
699 
700  // Resize the outgoing data structure
701  gpids.resize(N);
702 
703  // Start by claiming that I own all the data
704  if(use_minus_one_for_local)
705  for(i=0;i <N; i++) gpids[i]=std::make_pair(-1,Importer.TargetMap().GID64(i));
706  else
707  for(i=0;i <N; i++) gpids[i]=std::make_pair(mypid,Importer.TargetMap().GID64(i));
708 
709  // Now, for each remote ID, record who actually owns it. This loop follows the operation order in the
710  // MpiDistributor so it ought to duplicate that effect.
711  for(i=0,j=0;i<NumReceives;i++){
712  int pid=ProcsFrom[i];
713  for(k=0;k<LengthsFrom[i];k++){
714  if(pid!=mypid) gpids[RemoteLIDs[j]].first=pid;
715  j++;
716  }
717  }
718  return 0;
719 #else
720  EPETRA_CHK_ERR(-10);
721  // The above macro may not necessarily invoke 'return', or at least
722  // GCC 4.8.2 can't tell if it does. That's why I put an extra
723  // return statement here.
724  return 0;
725 #endif
726 }
727 #endif
728 
729 
730 //----------------------------------------------------------------------------
731 int Epetra_Util::GetPids(const Epetra_Import & Importer, std::vector<int> &pids, bool use_minus_one_for_local){
732 #ifdef HAVE_MPI
733  Epetra_MpiDistributor *D=dynamic_cast<Epetra_MpiDistributor*>(&Importer.Distributor());
734  if(!D) EPETRA_CHK_ERR(-2);
735 
736  int i,j,k;
737  int mypid=Importer.TargetMap().Comm().MyPID();
738  int N=Importer.TargetMap().NumMyElements();
739 
740  // Get the importer's data
741  const int *RemoteLIDs = Importer.RemoteLIDs();
742 
743  // Get the distributor's data
744  int NumReceives = D->NumReceives();
745  const int *ProcsFrom = D->ProcsFrom();
746  const int *LengthsFrom = D->LengthsFrom();
747 
748  // Resize the outgoing data structure
749  pids.resize(N);
750 
751  // Start by claiming that I own all the data
752  if(use_minus_one_for_local)
753  for(i=0; i<N; i++) pids[i]=-1;
754  else
755  for(i=0; i<N; i++) pids[i]=mypid;
756 
757  // Now, for each remote ID, record who actually owns it. This loop follows the operation order in the
758  // MpiDistributor so it ought to duplicate that effect.
759  for(i=0,j=0;i<NumReceives;i++){
760  int pid=ProcsFrom[i];
761  for(k=0;k<LengthsFrom[i];k++){
762  if(pid!=mypid) pids[RemoteLIDs[j]]=pid;
763  j++;
764  }
765  }
766  return 0;
767 #else
768  EPETRA_CHK_ERR(-10);
769  // The above macro may not necessarily invoke 'return', or at least
770  // GCC 4.8.2 can't tell if it does. That's why I put an extra
771  // return statement here.
772  return 0;
773 #endif
774 }
775 
776 //----------------------------------------------------------------------------
777 int Epetra_Util::GetRemotePIDs(const Epetra_Import & Importer, std::vector<int> &RemotePIDs){
778 #ifdef HAVE_MPI
779  Epetra_MpiDistributor *D=dynamic_cast<Epetra_MpiDistributor*>(&Importer.Distributor());
780  if(!D) {
781  RemotePIDs.resize(0);
782  return 0;
783  }
784 
785  int i,j,k;
786 
787  // Get the distributor's data
788  int NumReceives = D->NumReceives();
789  const int *ProcsFrom = D->ProcsFrom();
790  const int *LengthsFrom = D->LengthsFrom();
791 
792  // Resize the outgoing data structure
793  RemotePIDs.resize(Importer.NumRemoteIDs());
794 
795  // Now, for each remote ID, record who actually owns it. This loop follows the operation order in the
796  // MpiDistributor so it ought to duplicate that effect.
797  for(i=0,j=0;i<NumReceives;i++){
798  int pid=ProcsFrom[i];
799  for(k=0;k<LengthsFrom[i];k++){
800  RemotePIDs[j]=pid;
801  j++;
802  }
803  }
804  return 0;
805 #else
806  RemotePIDs.resize(0);
807  return 0;
808 #endif
809 }
810 
811 //----------------------------------------------------------------------------
812 template<typename T>
814  const T* list,
815  int len,
816  int& insertPoint)
817 {
818  if (len < 1) {
819  insertPoint = 0;
820  return(-1);
821  }
822 
823  unsigned start = 0, end = len - 1;
824 
825  while(end - start > 1) {
826  unsigned mid = (start + end) >> 1;
827  if (list[mid] < item) start = mid;
828  else end = mid;
829  }
830 
831  if (list[start] == item) return(start);
832  if (list[end] == item) return(end);
833 
834  if (list[end] < item) {
835  insertPoint = end+1;
836  return(-1);
837  }
838 
839  if (list[start] < item) insertPoint = end;
840  else insertPoint = start;
841 
842  return(-1);
843 }
844 
846  const int* list,
847  int len,
848  int& insertPoint)
849 {
850  return Epetra_Util_binary_search<int>(item, list, len, insertPoint);
851 }
852 
853 //----------------------------------------------------------------------------
854 int Epetra_Util_binary_search(long long item,
855  const long long* list,
856  int len,
857  int& insertPoint)
858 {
859  return Epetra_Util_binary_search<long long>(item, list, len, insertPoint);
860 }
861 
862 //----------------------------------------------------------------------------
863 template<typename T>
865  const int* list,
866  const T* aux_list,
867  int len,
868  int& insertPoint)
869 {
870  if (len < 1) {
871  insertPoint = 0;
872  return(-1);
873  }
874 
875  unsigned start = 0, end = len - 1;
876 
877  while(end - start > 1) {
878  unsigned mid = (start + end) >> 1;
879  if (aux_list[list[mid]] < item) start = mid;
880  else end = mid;
881  }
882 
883  if (aux_list[list[start]] == item) return(start);
884  if (aux_list[list[end]] == item) return(end);
885 
886  if (aux_list[list[end]] < item) {
887  insertPoint = end+1;
888  return(-1);
889  }
890 
891  if (aux_list[list[start]] < item) insertPoint = end;
892  else insertPoint = start;
893 
894  return(-1);
895 }
896 
897 //----------------------------------------------------------------------------
899  const int* list,
900  const int* aux_list,
901  int len,
902  int& insertPoint)
903 {
904  return Epetra_Util_binary_search_aux<int>(item, list, aux_list, len, insertPoint);
905 }
906 
907 //----------------------------------------------------------------------------
908 int Epetra_Util_binary_search_aux(long long item,
909  const int* list,
910  const long long* aux_list,
911  int len,
912  int& insertPoint)
913 {
914  return Epetra_Util_binary_search_aux<long long>(item, list, aux_list, len, insertPoint);
915 }
916 
917 
918 
919 //=========================================================================
921  Epetra_MultiVector * RHS,
922  int & M, int & N, int & nz, int * & ptr,
923  int * & ind, double * & val, int & Nrhs,
924  double * & rhs, int & ldrhs,
925  double * & lhs, int & ldlhs) {
926 
927  int ierr = 0;
928  if (A==0) EPETRA_CHK_ERR(-1); // This matrix is defined
929  if (!A->IndicesAreContiguous()) { // Data must be contiguous for this to work
930  EPETRA_CHK_ERR(A->MakeDataContiguous()); // Call MakeDataContiguous() method on the matrix
931  ierr = 1; // Warn User that we changed the matrix
932  }
933 
934  M = A->NumMyRows();
935  N = A->NumMyCols();
936  nz = A->NumMyNonzeros();
937  val = (*A)[0]; // Dangerous, but cheap and effective way to access first element in
938 
939  const Epetra_CrsGraph & Graph = A->Graph();
940  ind = Graph[0]; // list of values and indices
941 
942  Nrhs = 0; // Assume no rhs, lhs
943 
944  if (RHS!=0) {
945  Nrhs = RHS->NumVectors();
946  if (Nrhs>1)
947  if (!RHS->ConstantStride()) {EPETRA_CHK_ERR(-2)}; // Must have strided vectors
948  ldrhs = RHS->Stride();
949  rhs = (*RHS)[0]; // Dangerous but effective (again)
950  }
951  if (LHS!=0) {
952  int Nlhs = LHS->NumVectors();
953  if (Nlhs!=Nrhs) {EPETRA_CHK_ERR(-3)}; // Must have same number of rhs and lhs
954  if (Nlhs>1)
955  if (!LHS->ConstantStride()) {EPETRA_CHK_ERR(-4)}; // Must have strided vectors
956  ldlhs = LHS->Stride();
957  lhs = (*LHS)[0];
958  }
959 
960  // Finally build ptr vector
961 
962  if (ptr==0) {
963  ptr = new int[M+1];
964  ptr[0] = 0;
965  for (int i=0; i<M; i++) ptr[i+1] = ptr[i] + Graph.NumMyIndices(i);
966  }
967  EPETRA_CHK_ERR(ierr);
968  return(0);
969 }
970 
971 // Explicitly instantiate these two even though a call to them is made.
972 // Possible fix for a bug reported by Kenneth Belcourt.
973 template void Epetra_Util::Sort<int> (bool, int, int *, int, double **, int, int **, int, long long **);
974 template void Epetra_Util::Sort<long long>(bool, int, long long *, int, double **, int, int **, int, long long **);
975 
976 
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:119
unsigned int Seed_
Definition: Epetra_Util.h:256
bool UniqueGIDs() const
Returns true if map GIDs are 1-to-1.
bool ConstantStride() const
Returns true if this multi-vector has constant stride between vectors.
const Epetra_CrsGraph & Graph() const
Returns a reference to the Epetra_CrsGraph object associated with this matrix.
static Epetra_Map Create_OneToOne_Map(const Epetra_Map &usermap, bool high_rank_proc_owns_shared=false)
Epetra_Util Create_OneToOne_Map function.
int NumMyRows() const
Returns the number of matrix rows owned by the calling processor.
int NumMyCols() const
Returns the number of entries in the set of column-indices that appear on this processor.
virtual Epetra_Directory * CreateDirectory(const Epetra_BlockMap &Map) const =0
Create a directory object for the given Epetra_BlockMap.
int NumRemoteIDs() const
Returns the number of elements that are not on the calling processor.
static const double chopVal_
Definition: Epetra_Util.h:253
#define EPETRA_CHK_ERR(a)
int Stride() const
Returns the stride between vectors in the multi-vector (only meaningful if ConstantStride() is true)...
int IndexBase() const
Index base for this map.
const Epetra_BlockMap & TargetMap() const
Returns the TargetMap used to construct this importer.
unsigned int Seed() const
Get seed from Random function.
MPI implementation of Epetra_Distributor.
virtual int ReportError(const std::string Message, int ErrorCode) const
Error reporting method.
virtual int MinAll(double *PartialMins, double *GlobalMins, int Count) const =0
Epetra_Comm Global Min function.
long long NumGlobalElements64() const
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.
Epetra_Import: This class builds an import object for efficient importing of off-processor elements...
Definition: Epetra_Import.h:63
int * RemoteLIDs() const
List of elements in the target map that are coming from other processors.
static int SortAndMergeCrsEntries(int NumRows, int *CRS_rowptr, int *CRS_colind, double *CRS_vals)
Epetra_Util SortAndMergeCrsEntries function.
static int GetPids(const Epetra_Import &Importer, std::vector< int > &pids, bool use_minus_one_for_local)
Epetra_Util GetPids function.
virtual int MyPID() const =0
Return my process ID.
int MakeDataContiguous()
Eliminates memory that is used for construction. Make consecutive row index sections contiguous...
unsigned int RandomInt()
Returns a random integer on the interval (0, 2^31-1)
Definition: Epetra_Util.cpp:71
long long MaxAllGID64() const
bool IndicesAreContiguous() const
If matrix indices are packed into single array (done in OptimizeStorage()) return true...
Epetra_Directory: This class is a pure virtual class whose interface allows Epetra_Map and Epetr_Bloc...
static int GetPidGidPairs(const Epetra_Import &Importer, std::vector< std::pair< int, int > > &gpids, bool use_minus_one_for_local)
Epetra_Util GetPidGidPairs function.
static int SortCrsEntries(int NumRows, const int *CRS_rowptr, int *CRS_colind, double *CRS_vals)
Epetra_Util SortCrsEntries function.
int SetSeed(unsigned int Seed_in)
Set seed for Random function.
static Epetra_BlockMap Create_OneToOne_BlockMap(const Epetra_BlockMap &usermap, bool high_rank_proc_owns_shared=false)
Epetra_Util Create_OneToOne_Map function.
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition: Epetra_Comm.h:73
long long GID64(int LID) const
bool IsOneToOne() const
static Epetra_Map Create_Root_Map(const Epetra_Map &usermap, int root=0)
Epetra_Util Create_Root_Map function.
int NumVectors() const
Returns the number of vectors in the multi-vector.
bool GlobalIndicesInt() const
Returns true if map create with int NumGlobalElements.
int NumMyElements() const
Number of elements on the calling processor.
static Epetra_Map TCreate_Root_Map(const Epetra_Map &usermap, int root)
Epetra_GIDTypeVector: A class for constructing and using dense "int" and "long long" vectors on a par...
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
int Epetra_Util_binary_search_aux(T item, const int *list, const T *aux_list, int len, int &insertPoint)
Utility function to perform a binary-search on a list of data.
Epetra_Distributor & Distributor() const
int Epetra_Util_binary_search(T item, const T *list, int len, int &insertPoint)
Utility function to perform a binary-search on a list of data.
int NumMyNonzeros() const
Returns the number of nonzero entries in the calling processor&#39;s portion of the matrix.
long long IndexBase64() const
static int GetRemotePIDs(const Epetra_Import &Importer, std::vector< int > &RemotePIDs)
Epetra_Util GetRemotePIDs.
virtual int NumProc() const =0
Returns total number of processes.
const int * LengthsFrom() const
Number of values we&#39;re receiving from each proc.
double RandomDouble()
Returns a random double on the interval (-1.0,1.0)
Definition: Epetra_Util.cpp:90
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
int NumMyIndices(int Row) const
Returns the current number of nonzero entries in specified local row on this processor.
static double Chop(const double &Value)
Epetra_Util Chop method. Return zero if input Value is less than ChopValue.
Definition: Epetra_Util.cpp:65
static void Sort(bool SortAscending, int NumKeys, T *Keys, int NumDoubleCompanions, double **DoubleCompanions, int NumIntCompanions, int **IntCompanions, int NumLongLongCompanions, long long **LongLongCompanions)
Epetra_Util Sort Routine (Shell sort)
int NumReceives() const
The number of procs from which we will receive data.
bool LinearMap() const
Returns true if the global ID space is contiguously divided (but not necessarily uniformly) across al...
const int * ProcsFrom() const
A list of procs sending values to this proc.
bool GlobalIndicesLongLong() const
Returns true if map create with long long NumGlobalElements.
int GID(int LID) const
Returns global ID of local ID, return IndexBase-1 if not found on this processor. ...
Epetra_CrsGraph: A class for constructing and using sparse compressed row graphs. ...
Epetra_LocalMap: A class for replicating vectors and matrices across multiple processors.
int Epetra_Util_ExtractHbData(Epetra_CrsMatrix *A, Epetra_MultiVector *LHS, Epetra_MultiVector *RHS, int &M, int &N, int &nz, int *&ptr, int *&ind, double *&val, int &Nrhs, double *&rhs, int &ldrhs, double *&lhs, int &ldlhs)
Harwell-Boeing data extraction routine.