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::SortAndMergeCrsEntries(int NumRows, int *CRS_rowptr, int *CRS_colind, double *CRS_vals){
470  // For each row, sort column entries from smallest to largest, merging column ids that are identical by adding values.
471  // Use shell sort. Stable sort so it is fast if indices are already sorted.
472  // Code copied from Epetra_CrsMatrix::SortEntries()
473 
474  int nnz = CRS_rowptr[NumRows];
475  int new_curr=CRS_rowptr[0], old_curr=CRS_rowptr[0];
476 
477  for(int i = 0; i < NumRows; i++){
478  int start=CRS_rowptr[i];
479  if(start >= nnz) continue;
480 
481  double* locValues = &CRS_vals[start];
482  int NumEntries = CRS_rowptr[i+1] - start;
483  int* locIndices = &CRS_colind[start];
484 
485  // Sort phase
486  int n = NumEntries;
487  int m = n/2;
488 
489  while(m > 0) {
490  int max = n - m;
491  for(int j = 0; j < max; j++) {
492  for(int k = j; k >= 0; k-=m) {
493  if(locIndices[k+m] >= locIndices[k])
494  break;
495  double dtemp = locValues[k+m];
496  locValues[k+m] = locValues[k];
497  locValues[k] = dtemp;
498  int itemp = locIndices[k+m];
499  locIndices[k+m] = locIndices[k];
500  locIndices[k] = itemp;
501  }
502  }
503  m = m/2;
504  }
505 
506  // Merge & shrink
507  for(int j=CRS_rowptr[i]; j < CRS_rowptr[i+1]; j++) {
508  if(j > CRS_rowptr[i] && CRS_colind[j]==CRS_colind[new_curr-1]) {
509  CRS_vals[new_curr-1] += CRS_vals[j];
510  }
511  else if(new_curr==j) {
512  new_curr++;
513  }
514  else {
515  CRS_colind[new_curr] = CRS_colind[j];
516  CRS_vals[new_curr] = CRS_vals[j];
517  new_curr++;
518  }
519  }
520 
521  CRS_rowptr[i] = old_curr;
522  old_curr=new_curr;
523  }
524 
525  CRS_rowptr[NumRows] = new_curr;
526  return (0);
527 }
528 
529 
530 //----------------------------------------------------------------------------
531 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
532 int Epetra_Util::GetPidGidPairs(const Epetra_Import & Importer,std::vector< std::pair<int,int> > & gpids, bool use_minus_one_for_local){
533  // 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.
534  // This only works if we have an MpiDistributor in our Importer. Otherwise return an error.
535 #ifdef HAVE_MPI
536  Epetra_MpiDistributor *D=dynamic_cast<Epetra_MpiDistributor*>(&Importer.Distributor());
537  if(!D) EPETRA_CHK_ERR(-2);
538 
539  int i,j,k;
540  int mypid=Importer.TargetMap().Comm().MyPID();
541  int N=Importer.TargetMap().NumMyElements();
542 
543  // Get the importer's data
544  const int *RemoteLIDs = Importer.RemoteLIDs();
545 
546  // Get the distributor's data
547  int NumReceives = D->NumReceives();
548  const int *ProcsFrom = D->ProcsFrom();
549  const int *LengthsFrom = D->LengthsFrom();
550 
551  // Resize the outgoing data structure
552  gpids.resize(N);
553 
554  // Start by claiming that I own all the data
555  if(use_minus_one_for_local)
556  for(i=0;i <N; i++) gpids[i]=std::make_pair(-1,Importer.TargetMap().GID(i));
557  else
558  for(i=0;i <N; i++) gpids[i]=std::make_pair(mypid,Importer.TargetMap().GID(i));
559 
560  // Now, for each remote ID, record who actually owns it. This loop follows the operation order in the
561  // MpiDistributor so it ought to duplicate that effect.
562  for(i=0,j=0;i<NumReceives;i++){
563  int pid=ProcsFrom[i];
564  for(k=0;k<LengthsFrom[i];k++){
565  if(pid!=mypid) gpids[RemoteLIDs[j]].first=pid;
566  j++;
567  }
568  }
569  return 0;
570 #else
571  EPETRA_CHK_ERR(-10);
572  // The above macro may not necessarily invoke 'return', or at least
573  // GCC 4.8.2 can't tell if it does. That's why I put an extra
574  // return statement here.
575  return 0;
576 #endif
577 }
578 #endif
579 
580 //----------------------------------------------------------------------------
581 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
582 int Epetra_Util::GetPidGidPairs(const Epetra_Import & Importer,std::vector< std::pair<int,long long> > & gpids, bool use_minus_one_for_local){
583  // 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.
584  // This only works if we have an MpiDistributor in our Importer. Otheriwise return an error.
585 #ifdef HAVE_MPI
586  Epetra_MpiDistributor *D=dynamic_cast<Epetra_MpiDistributor*>(&Importer.Distributor());
587  if(!D) EPETRA_CHK_ERR(-2);
588 
589  int i,j,k;
590  int mypid=Importer.TargetMap().Comm().MyPID();
591  int N=Importer.TargetMap().NumMyElements();
592 
593  // Get the importer's data
594  const int *RemoteLIDs = Importer.RemoteLIDs();
595 
596  // Get the distributor's data
597  int NumReceives = D->NumReceives();
598  const int *ProcsFrom = D->ProcsFrom();
599  const int *LengthsFrom = D->LengthsFrom();
600 
601  // Resize the outgoing data structure
602  gpids.resize(N);
603 
604  // Start by claiming that I own all the data
605  if(use_minus_one_for_local)
606  for(i=0;i <N; i++) gpids[i]=std::make_pair(-1,Importer.TargetMap().GID64(i));
607  else
608  for(i=0;i <N; i++) gpids[i]=std::make_pair(mypid,Importer.TargetMap().GID64(i));
609 
610  // Now, for each remote ID, record who actually owns it. This loop follows the operation order in the
611  // MpiDistributor so it ought to duplicate that effect.
612  for(i=0,j=0;i<NumReceives;i++){
613  int pid=ProcsFrom[i];
614  for(k=0;k<LengthsFrom[i];k++){
615  if(pid!=mypid) gpids[RemoteLIDs[j]].first=pid;
616  j++;
617  }
618  }
619  return 0;
620 #else
621  EPETRA_CHK_ERR(-10);
622  // The above macro may not necessarily invoke 'return', or at least
623  // GCC 4.8.2 can't tell if it does. That's why I put an extra
624  // return statement here.
625  return 0;
626 #endif
627 }
628 #endif
629 
630 
631 //----------------------------------------------------------------------------
632 int Epetra_Util::GetPids(const Epetra_Import & Importer, std::vector<int> &pids, bool use_minus_one_for_local){
633 #ifdef HAVE_MPI
634  Epetra_MpiDistributor *D=dynamic_cast<Epetra_MpiDistributor*>(&Importer.Distributor());
635  if(!D) EPETRA_CHK_ERR(-2);
636 
637  int i,j,k;
638  int mypid=Importer.TargetMap().Comm().MyPID();
639  int N=Importer.TargetMap().NumMyElements();
640 
641  // Get the importer's data
642  const int *RemoteLIDs = Importer.RemoteLIDs();
643 
644  // Get the distributor's data
645  int NumReceives = D->NumReceives();
646  const int *ProcsFrom = D->ProcsFrom();
647  const int *LengthsFrom = D->LengthsFrom();
648 
649  // Resize the outgoing data structure
650  pids.resize(N);
651 
652  // Start by claiming that I own all the data
653  if(use_minus_one_for_local)
654  for(i=0; i<N; i++) pids[i]=-1;
655  else
656  for(i=0; i<N; i++) pids[i]=mypid;
657 
658  // Now, for each remote ID, record who actually owns it. This loop follows the operation order in the
659  // MpiDistributor so it ought to duplicate that effect.
660  for(i=0,j=0;i<NumReceives;i++){
661  int pid=ProcsFrom[i];
662  for(k=0;k<LengthsFrom[i];k++){
663  if(pid!=mypid) pids[RemoteLIDs[j]]=pid;
664  j++;
665  }
666  }
667  return 0;
668 #else
669  EPETRA_CHK_ERR(-10);
670  // The above macro may not necessarily invoke 'return', or at least
671  // GCC 4.8.2 can't tell if it does. That's why I put an extra
672  // return statement here.
673  return 0;
674 #endif
675 }
676 
677 //----------------------------------------------------------------------------
678 int Epetra_Util::GetRemotePIDs(const Epetra_Import & Importer, std::vector<int> &RemotePIDs){
679 #ifdef HAVE_MPI
680  Epetra_MpiDistributor *D=dynamic_cast<Epetra_MpiDistributor*>(&Importer.Distributor());
681  if(!D) {
682  RemotePIDs.resize(0);
683  return 0;
684  }
685 
686  int i,j,k;
687 
688  // Get the distributor's data
689  int NumReceives = D->NumReceives();
690  const int *ProcsFrom = D->ProcsFrom();
691  const int *LengthsFrom = D->LengthsFrom();
692 
693  // Resize the outgoing data structure
694  RemotePIDs.resize(Importer.NumRemoteIDs());
695 
696  // Now, for each remote ID, record who actually owns it. This loop follows the operation order in the
697  // MpiDistributor so it ought to duplicate that effect.
698  for(i=0,j=0;i<NumReceives;i++){
699  int pid=ProcsFrom[i];
700  for(k=0;k<LengthsFrom[i];k++){
701  RemotePIDs[j]=pid;
702  j++;
703  }
704  }
705  return 0;
706 #else
707  RemotePIDs.resize(0);
708  return 0;
709 #endif
710 }
711 
712 //----------------------------------------------------------------------------
713 template<typename T>
715  const T* list,
716  int len,
717  int& insertPoint)
718 {
719  if (len < 1) {
720  insertPoint = 0;
721  return(-1);
722  }
723 
724  unsigned start = 0, end = len - 1;
725 
726  while(end - start > 1) {
727  unsigned mid = (start + end) >> 1;
728  if (list[mid] < item) start = mid;
729  else end = mid;
730  }
731 
732  if (list[start] == item) return(start);
733  if (list[end] == item) return(end);
734 
735  if (list[end] < item) {
736  insertPoint = end+1;
737  return(-1);
738  }
739 
740  if (list[start] < item) insertPoint = end;
741  else insertPoint = start;
742 
743  return(-1);
744 }
745 
747  const int* list,
748  int len,
749  int& insertPoint)
750 {
751  return Epetra_Util_binary_search<int>(item, list, len, insertPoint);
752 }
753 
754 //----------------------------------------------------------------------------
755 int Epetra_Util_binary_search(long long item,
756  const long long* list,
757  int len,
758  int& insertPoint)
759 {
760  return Epetra_Util_binary_search<long long>(item, list, len, insertPoint);
761 }
762 
763 //----------------------------------------------------------------------------
764 template<typename T>
766  const int* list,
767  const T* aux_list,
768  int len,
769  int& insertPoint)
770 {
771  if (len < 1) {
772  insertPoint = 0;
773  return(-1);
774  }
775 
776  unsigned start = 0, end = len - 1;
777 
778  while(end - start > 1) {
779  unsigned mid = (start + end) >> 1;
780  if (aux_list[list[mid]] < item) start = mid;
781  else end = mid;
782  }
783 
784  if (aux_list[list[start]] == item) return(start);
785  if (aux_list[list[end]] == item) return(end);
786 
787  if (aux_list[list[end]] < item) {
788  insertPoint = end+1;
789  return(-1);
790  }
791 
792  if (aux_list[list[start]] < item) insertPoint = end;
793  else insertPoint = start;
794 
795  return(-1);
796 }
797 
798 //----------------------------------------------------------------------------
800  const int* list,
801  const int* aux_list,
802  int len,
803  int& insertPoint)
804 {
805  return Epetra_Util_binary_search_aux<int>(item, list, aux_list, len, insertPoint);
806 }
807 
808 //----------------------------------------------------------------------------
809 int Epetra_Util_binary_search_aux(long long item,
810  const int* list,
811  const long long* aux_list,
812  int len,
813  int& insertPoint)
814 {
815  return Epetra_Util_binary_search_aux<long long>(item, list, aux_list, len, insertPoint);
816 }
817 
818 
819 
820 //=========================================================================
822  Epetra_MultiVector * RHS,
823  int & M, int & N, int & nz, int * & ptr,
824  int * & ind, double * & val, int & Nrhs,
825  double * & rhs, int & ldrhs,
826  double * & lhs, int & ldlhs) {
827 
828  int ierr = 0;
829  if (A==0) EPETRA_CHK_ERR(-1); // This matrix is defined
830  if (!A->IndicesAreContiguous()) { // Data must be contiguous for this to work
831  EPETRA_CHK_ERR(A->MakeDataContiguous()); // Call MakeDataContiguous() method on the matrix
832  ierr = 1; // Warn User that we changed the matrix
833  }
834 
835  M = A->NumMyRows();
836  N = A->NumMyCols();
837  nz = A->NumMyNonzeros();
838  val = (*A)[0]; // Dangerous, but cheap and effective way to access first element in
839 
840  const Epetra_CrsGraph & Graph = A->Graph();
841  ind = Graph[0]; // list of values and indices
842 
843  Nrhs = 0; // Assume no rhs, lhs
844 
845  if (RHS!=0) {
846  Nrhs = RHS->NumVectors();
847  if (Nrhs>1)
848  if (!RHS->ConstantStride()) {EPETRA_CHK_ERR(-2)}; // Must have strided vectors
849  ldrhs = RHS->Stride();
850  rhs = (*RHS)[0]; // Dangerous but effective (again)
851  }
852  if (LHS!=0) {
853  int Nlhs = LHS->NumVectors();
854  if (Nlhs!=Nrhs) {EPETRA_CHK_ERR(-3)}; // Must have same number of rhs and lhs
855  if (Nlhs>1)
856  if (!LHS->ConstantStride()) {EPETRA_CHK_ERR(-4)}; // Must have strided vectors
857  ldlhs = LHS->Stride();
858  lhs = (*LHS)[0];
859  }
860 
861  // Finally build ptr vector
862 
863  if (ptr==0) {
864  ptr = new int[M+1];
865  ptr[0] = 0;
866  for (int i=0; i<M; i++) ptr[i+1] = ptr[i] + Graph.NumMyIndices(i);
867  }
868  EPETRA_CHK_ERR(ierr);
869  return(0);
870 }
871 
872 // Explicitly instantiate these two even though a call to them is made.
873 // Possible fix for a bug reported by Kenneth Belcourt.
874 template void Epetra_Util::Sort<int> (bool, int, int *, int, double **, int, int **, int, long long **);
875 template void Epetra_Util::Sort<long long>(bool, int, long long *, int, double **, int, int **, int, long long **);
876 
877 
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:246
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:243
#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.