Epetra Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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
379  bool high_rank_proc_owns_shared)
380 {
381 // FIXME long long
382 
383  //if usermap is already 1-to-1 then we'll just return a copy of it.
384  if (usermap.IsOneToOne()) {
385  Epetra_BlockMap newmap(usermap);
386  return(newmap);
387  }
388 
389  int myPID = usermap.Comm().MyPID();
390  Epetra_Directory* directory = usermap.Comm().CreateDirectory(usermap);
391 
392  int numMyElems = usermap.NumMyElements();
393  const int* myElems = usermap.MyGlobalElements();
394 
395  int* owner_procs = new int[numMyElems*2];
396  int* sizes = owner_procs+numMyElems;
397 
398  directory->GetDirectoryEntries(usermap, numMyElems, myElems, owner_procs,
399  0, sizes, high_rank_proc_owns_shared);
400 
401  //we'll fill a list of map-elements which belong on this processor
402 
403  int* myOwnedElems = new int[numMyElems*2];
404  int* ownedSizes = myOwnedElems+numMyElems;
405  int numMyOwnedElems = 0;
406 
407  for(int i=0; i<numMyElems; ++i) {
408  int GID = myElems[i];
409  int owner = owner_procs[i];
410 
411  if (myPID == owner) {
412  ownedSizes[numMyOwnedElems] = sizes[i];
413  myOwnedElems[numMyOwnedElems++] = GID;
414  }
415  }
416 
417  Epetra_BlockMap one_to_one_map(-1, numMyOwnedElems, myOwnedElems,
418  sizes, usermap.IndexBase(), usermap.Comm()); // CJ TODO FIXME long long
419 
420  delete [] myOwnedElems;
421  delete [] owner_procs;
422  delete directory;
423 
424  return(one_to_one_map);
425 }
426 #endif // EPETRA_NO_32BIT_GLOBAL_INDICES
427 
428 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
431  bool high_rank_proc_owns_shared)
432 {
433  //if usermap is already 1-to-1 then we'll just return a copy of it.
434  if (usermap.IsOneToOne()) {
435  Epetra_BlockMap newmap(usermap);
436  return(newmap);
437  }
438 
439  int myPID = usermap.Comm().MyPID();
440  Epetra_Directory* directory = usermap.Comm().CreateDirectory(usermap);
441 
442  int numMyElems = usermap.NumMyElements();
443  const long long* myElems = usermap.MyGlobalElements64();
444 
445  int* owner_procs = new int[numMyElems*2];
446  int* sizes = owner_procs+numMyElems;
447 
448  directory->GetDirectoryEntries(usermap, numMyElems, myElems, owner_procs,
449  0, sizes, high_rank_proc_owns_shared);
450 
451  //we'll fill a list of map-elements which belong on this processor
452 
453  long long* myOwnedElems = new long long[numMyElems*2];
454  long long* ownedSizes = myOwnedElems+numMyElems;
455  int numMyOwnedElems = 0;
456 
457  for(int i=0; i<numMyElems; ++i) {
458  long long GID = myElems[i];
459  int owner = owner_procs[i];
460 
461  if (myPID == owner) {
462  ownedSizes[numMyOwnedElems] = sizes[i];
463  myOwnedElems[numMyOwnedElems++] = GID;
464  }
465  }
466 
467  Epetra_BlockMap one_to_one_map((long long)-1, numMyOwnedElems, myOwnedElems,
468  sizes, usermap.IndexBase64(), usermap.Comm());
469 
470  delete [] myOwnedElems;
471  delete [] owner_procs;
472  delete directory;
473 
474  return(one_to_one_map);
475 }
476 #endif // EPETRA_NO_64BIT_GLOBAL_INDICES
477 
478 //----------------------------------------------------------------------------
479 int Epetra_Util::SortCrsEntries(int NumRows, const int *CRS_rowptr, int *CRS_colind, double *CRS_vals){
480  // For each row, sort column entries from smallest to largest.
481  // Use shell sort. Stable sort so it is fast if indices are already sorted.
482  // Code copied from Epetra_CrsMatrix::SortEntries()
483  int nnz = CRS_rowptr[NumRows];
484 
485  for(int i = 0; i < NumRows; i++){
486  int start=CRS_rowptr[i];
487  if(start >= nnz) continue;
488 
489  double* locValues = &CRS_vals[start];
490  int NumEntries = CRS_rowptr[i+1] - start;
491  int* locIndices = &CRS_colind[start];
492 
493  int n = NumEntries;
494  int m = n/2;
495 
496  while(m > 0) {
497  int max = n - m;
498  for(int j = 0; j < max; j++) {
499  for(int k = j; k >= 0; k-=m) {
500  if(locIndices[k+m] >= locIndices[k])
501  break;
502  double dtemp = locValues[k+m];
503  locValues[k+m] = locValues[k];
504  locValues[k] = dtemp;
505  int itemp = locIndices[k+m];
506  locIndices[k+m] = locIndices[k];
507  locIndices[k] = itemp;
508  }
509  }
510  m = m/2;
511  }
512  }
513  return(0);
514 }
515 
516 //----------------------------------------------------------------------------
517 int Epetra_Util::SortCrsEntries(int NumRows, const size_t *CRS_rowptr, int *CRS_colind, double *CRS_vals){
518  // For each row, sort column entries from smallest to largest.
519  // Use shell sort. Stable sort so it is fast if indices are already sorted.
520  // Code copied from Epetra_CrsMatrix::SortEntries()
521  size_t nnz = CRS_rowptr[NumRows];
522 
523  for(int i = 0; i < NumRows; i++){
524  size_t start = CRS_rowptr[i];
525  if(start >= nnz) continue;
526 
527  double* locValues = &CRS_vals[start];
528  int NumEntries = static_cast<int>(CRS_rowptr[i+1] - start);
529  int* locIndices = &CRS_colind[start];
530 
531  int n = NumEntries;
532  int m = n/2;
533 
534  while(m > 0) {
535  int max = n - m;
536  for(int j = 0; j < max; j++) {
537  for(int k = j; k >= 0; k-=m) {
538  if(locIndices[k+m] >= locIndices[k])
539  break;
540  double dtemp = locValues[k+m];
541  locValues[k+m] = locValues[k];
542  locValues[k] = dtemp;
543  int itemp = locIndices[k+m];
544  locIndices[k+m] = locIndices[k];
545  locIndices[k] = itemp;
546  }
547  }
548  m = m/2;
549  }
550  }
551  return(0);
552 }
553 
554 //----------------------------------------------------------------------------
555 int Epetra_Util::SortAndMergeCrsEntries(int NumRows, int *CRS_rowptr, int *CRS_colind, double *CRS_vals){
556  // For each row, sort column entries from smallest to largest, merging column ids that are identical by adding values.
557  // Use shell sort. Stable sort so it is fast if indices are already sorted.
558  // Code copied from Epetra_CrsMatrix::SortEntries()
559 
560  int nnz = CRS_rowptr[NumRows];
561  int new_curr=CRS_rowptr[0], old_curr=CRS_rowptr[0];
562 
563  for(int i = 0; i < NumRows; i++){
564  int start=CRS_rowptr[i];
565  if(start >= nnz) continue;
566 
567  double* locValues = &CRS_vals[start];
568  int NumEntries = CRS_rowptr[i+1] - start;
569  int* locIndices = &CRS_colind[start];
570 
571  // Sort phase
572  int n = NumEntries;
573  int m = n/2;
574 
575  while(m > 0) {
576  int max = n - m;
577  for(int j = 0; j < max; j++) {
578  for(int k = j; k >= 0; k-=m) {
579  if(locIndices[k+m] >= locIndices[k])
580  break;
581  double dtemp = locValues[k+m];
582  locValues[k+m] = locValues[k];
583  locValues[k] = dtemp;
584  int itemp = locIndices[k+m];
585  locIndices[k+m] = locIndices[k];
586  locIndices[k] = itemp;
587  }
588  }
589  m = m/2;
590  }
591 
592  // Merge & shrink
593  for(int j=CRS_rowptr[i]; j < CRS_rowptr[i+1]; j++) {
594  if(j > CRS_rowptr[i] && CRS_colind[j]==CRS_colind[new_curr-1]) {
595  CRS_vals[new_curr-1] += CRS_vals[j];
596  }
597  else if(new_curr==j) {
598  new_curr++;
599  }
600  else {
601  CRS_colind[new_curr] = CRS_colind[j];
602  CRS_vals[new_curr] = CRS_vals[j];
603  new_curr++;
604  }
605  }
606 
607  CRS_rowptr[i] = old_curr;
608  old_curr=new_curr;
609  }
610 
611  CRS_rowptr[NumRows] = new_curr;
612  return (0);
613 }
614 
615 //----------------------------------------------------------------------------
616 int Epetra_Util::SortAndMergeCrsEntries(int NumRows, size_t *CRS_rowptr, int *CRS_colind, double *CRS_vals){
617  // For each row, sort column entries from smallest to largest, merging column ids that are identical by adding values.
618  // Use shell sort. Stable sort so it is fast if indices are already sorted.
619  // Code copied from Epetra_CrsMatrix::SortEntries()
620 
621  size_t nnz = CRS_rowptr[NumRows];
622  size_t new_curr=CRS_rowptr[0], old_curr=CRS_rowptr[0];
623 
624  for(int i = 0; i < NumRows; i++){
625  size_t start=CRS_rowptr[i];
626  if(start >= nnz) continue;
627 
628  double* locValues = &CRS_vals[start];
629  int NumEntries = static_cast<int>(CRS_rowptr[i+1] - start);
630  int* locIndices = &CRS_colind[start];
631 
632  // Sort phase
633  int n = NumEntries;
634  int m = n/2;
635 
636  while(m > 0) {
637  int max = n - m;
638  for(int j = 0; j < max; j++) {
639  for(int k = j; k >= 0; k-=m) {
640  if(locIndices[k+m] >= locIndices[k])
641  break;
642  double dtemp = locValues[k+m];
643  locValues[k+m] = locValues[k];
644  locValues[k] = dtemp;
645  int itemp = locIndices[k+m];
646  locIndices[k+m] = locIndices[k];
647  locIndices[k] = itemp;
648  }
649  }
650  m = m/2;
651  }
652 
653  // Merge & shrink
654  for(size_t j=CRS_rowptr[i]; j < CRS_rowptr[i+1]; j++) {
655  if(j > CRS_rowptr[i] && CRS_colind[j]==CRS_colind[new_curr-1]) {
656  CRS_vals[new_curr-1] += CRS_vals[j];
657  }
658  else if(new_curr==j) {
659  new_curr++;
660  }
661  else {
662  CRS_colind[new_curr] = CRS_colind[j];
663  CRS_vals[new_curr] = CRS_vals[j];
664  new_curr++;
665  }
666  }
667 
668  CRS_rowptr[i] = old_curr;
669  old_curr=new_curr;
670  }
671 
672  CRS_rowptr[NumRows] = new_curr;
673  return (0);
674 }
675 
676 
677 //----------------------------------------------------------------------------
678 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
679 int Epetra_Util::GetPidGidPairs(const Epetra_Import & Importer,std::vector< std::pair<int,int> > & gpids, bool use_minus_one_for_local){
680  // 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.
681  // This only works if we have an MpiDistributor in our Importer. Otherwise return an error.
682 #ifdef HAVE_MPI
683  Epetra_MpiDistributor *D=dynamic_cast<Epetra_MpiDistributor*>(&Importer.Distributor());
684  if(!D) EPETRA_CHK_ERR(-2);
685 
686  int i,j,k;
687  int mypid=Importer.TargetMap().Comm().MyPID();
688  int N=Importer.TargetMap().NumMyElements();
689 
690  // Get the importer's data
691  const int *RemoteLIDs = Importer.RemoteLIDs();
692 
693  // Get the distributor's data
694  int NumReceives = D->NumReceives();
695  const int *ProcsFrom = D->ProcsFrom();
696  const int *LengthsFrom = D->LengthsFrom();
697 
698  // Resize the outgoing data structure
699  gpids.resize(N);
700 
701  // Start by claiming that I own all the data
702  if(use_minus_one_for_local)
703  for(i=0;i <N; i++) gpids[i]=std::make_pair(-1,Importer.TargetMap().GID(i));
704  else
705  for(i=0;i <N; i++) gpids[i]=std::make_pair(mypid,Importer.TargetMap().GID(i));
706 
707  // Now, for each remote ID, record who actually owns it. This loop follows the operation order in the
708  // MpiDistributor so it ought to duplicate that effect.
709  for(i=0,j=0;i<NumReceives;i++){
710  int pid=ProcsFrom[i];
711  for(k=0;k<LengthsFrom[i];k++){
712  if(pid!=mypid) gpids[RemoteLIDs[j]].first=pid;
713  j++;
714  }
715  }
716  return 0;
717 #else
718  EPETRA_CHK_ERR(-10);
719  // The above macro may not necessarily invoke 'return', or at least
720  // GCC 4.8.2 can't tell if it does. That's why I put an extra
721  // return statement here.
722  return 0;
723 #endif
724 }
725 #endif
726 
727 //----------------------------------------------------------------------------
728 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
729 int Epetra_Util::GetPidGidPairs(const Epetra_Import & Importer,std::vector< std::pair<int,long long> > & gpids, bool use_minus_one_for_local){
730  // 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.
731  // This only works if we have an MpiDistributor in our Importer. Otheriwise return an error.
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  gpids.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++) gpids[i]=std::make_pair(-1,Importer.TargetMap().GID64(i));
754  else
755  for(i=0;i <N; i++) gpids[i]=std::make_pair(mypid,Importer.TargetMap().GID64(i));
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) gpids[RemoteLIDs[j]].first=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 #endif
776 
777 
778 //----------------------------------------------------------------------------
779 int Epetra_Util::GetPids(const Epetra_Import & Importer, std::vector<int> &pids, bool use_minus_one_for_local){
780 #ifdef HAVE_MPI
781  Epetra_MpiDistributor *D=dynamic_cast<Epetra_MpiDistributor*>(&Importer.Distributor());
782  if(!D) EPETRA_CHK_ERR(-2);
783 
784  int i,j,k;
785  int mypid=Importer.TargetMap().Comm().MyPID();
786  int N=Importer.TargetMap().NumMyElements();
787 
788  // Get the importer's data
789  const int *RemoteLIDs = Importer.RemoteLIDs();
790 
791  // Get the distributor's data
792  int NumReceives = D->NumReceives();
793  const int *ProcsFrom = D->ProcsFrom();
794  const int *LengthsFrom = D->LengthsFrom();
795 
796  // Resize the outgoing data structure
797  pids.resize(N);
798 
799  // Start by claiming that I own all the data
800  if(use_minus_one_for_local)
801  for(i=0; i<N; i++) pids[i]=-1;
802  else
803  for(i=0; i<N; i++) pids[i]=mypid;
804 
805  // Now, for each remote ID, record who actually owns it. This loop follows the operation order in the
806  // MpiDistributor so it ought to duplicate that effect.
807  for(i=0,j=0;i<NumReceives;i++){
808  int pid=ProcsFrom[i];
809  for(k=0;k<LengthsFrom[i];k++){
810  if(pid!=mypid) pids[RemoteLIDs[j]]=pid;
811  j++;
812  }
813  }
814  return 0;
815 #else
816  EPETRA_CHK_ERR(-10);
817  // The above macro may not necessarily invoke 'return', or at least
818  // GCC 4.8.2 can't tell if it does. That's why I put an extra
819  // return statement here.
820  return 0;
821 #endif
822 }
823 
824 //----------------------------------------------------------------------------
825 int Epetra_Util::GetRemotePIDs(const Epetra_Import & Importer, std::vector<int> &RemotePIDs){
826 #ifdef HAVE_MPI
827  const Epetra_Distributor* Dptr = Importer.DistributorPtr();
828  if (Dptr == NULL) {
829  RemotePIDs.resize(0);
830  return 0;
831  }
832  const Epetra_MpiDistributor *D=dynamic_cast<const Epetra_MpiDistributor*>(Dptr);
833  if(!D) {
834  RemotePIDs.resize(0);
835  return 0;
836  }
837 
838  int i,j,k;
839 
840  // Get the distributor's data
841  int NumReceives = D->NumReceives();
842  const int *ProcsFrom = D->ProcsFrom();
843  const int *LengthsFrom = D->LengthsFrom();
844 
845  // Resize the outgoing data structure
846  RemotePIDs.resize(Importer.NumRemoteIDs());
847 
848  // Now, for each remote ID, record who actually owns it. This loop follows the operation order in the
849  // MpiDistributor so it ought to duplicate that effect.
850  for(i=0,j=0;i<NumReceives;i++){
851  int pid=ProcsFrom[i];
852  for(k=0;k<LengthsFrom[i];k++){
853  RemotePIDs[j]=pid;
854  j++;
855  }
856  }
857  return 0;
858 #else
859  RemotePIDs.resize(0);
860  return 0;
861 #endif
862 }
863 
864 //----------------------------------------------------------------------------
865 template<typename T>
867  const T* list,
868  int len,
869  int& insertPoint)
870 {
871  if (len < 1) {
872  insertPoint = 0;
873  return(-1);
874  }
875 
876  unsigned start = 0, end = len - 1;
877 
878  while(end - start > 1) {
879  unsigned mid = (start + end) >> 1;
880  if (list[mid] < item) start = mid;
881  else end = mid;
882  }
883 
884  if (list[start] == item) return(start);
885  if (list[end] == item) return(end);
886 
887  if (list[end] < item) {
888  insertPoint = end+1;
889  return(-1);
890  }
891 
892  if (list[start] < item) insertPoint = end;
893  else insertPoint = start;
894 
895  return(-1);
896 }
897 
899  const int* list,
900  int len,
901  int& insertPoint)
902 {
903  return Epetra_Util_binary_search<int>(item, list, len, insertPoint);
904 }
905 
906 //----------------------------------------------------------------------------
907 int Epetra_Util_binary_search(long long item,
908  const long long* list,
909  int len,
910  int& insertPoint)
911 {
912  return Epetra_Util_binary_search<long long>(item, list, len, insertPoint);
913 }
914 
915 //----------------------------------------------------------------------------
916 template<typename T>
918  const int* list,
919  const T* aux_list,
920  int len,
921  int& insertPoint)
922 {
923  if (len < 1) {
924  insertPoint = 0;
925  return(-1);
926  }
927 
928  unsigned start = 0, end = len - 1;
929 
930  while(end - start > 1) {
931  unsigned mid = (start + end) >> 1;
932  if (aux_list[list[mid]] < item) start = mid;
933  else end = mid;
934  }
935 
936  if (aux_list[list[start]] == item) return(start);
937  if (aux_list[list[end]] == item) return(end);
938 
939  if (aux_list[list[end]] < item) {
940  insertPoint = end+1;
941  return(-1);
942  }
943 
944  if (aux_list[list[start]] < item) insertPoint = end;
945  else insertPoint = start;
946 
947  return(-1);
948 }
949 
950 //----------------------------------------------------------------------------
952  const int* list,
953  const int* aux_list,
954  int len,
955  int& insertPoint)
956 {
957  return Epetra_Util_binary_search_aux<int>(item, list, aux_list, len, insertPoint);
958 }
959 
960 //----------------------------------------------------------------------------
961 int Epetra_Util_binary_search_aux(long long item,
962  const int* list,
963  const long long* aux_list,
964  int len,
965  int& insertPoint)
966 {
967  return Epetra_Util_binary_search_aux<long long>(item, list, aux_list, len, insertPoint);
968 }
969 
970 
971 
972 //=========================================================================
974  Epetra_MultiVector * RHS,
975  int & M, int & N, int & nz, int * & ptr,
976  int * & ind, double * & val, int & Nrhs,
977  double * & rhs, int & ldrhs,
978  double * & lhs, int & ldlhs) {
979 
980  int ierr = 0;
981  if (A==0) EPETRA_CHK_ERR(-1); // This matrix is defined
982  if (!A->IndicesAreContiguous()) { // Data must be contiguous for this to work
983  EPETRA_CHK_ERR(A->MakeDataContiguous()); // Call MakeDataContiguous() method on the matrix
984  ierr = 1; // Warn User that we changed the matrix
985  }
986 
987  M = A->NumMyRows();
988  N = A->NumMyCols();
989  nz = A->NumMyNonzeros();
990  val = (*A)[0]; // Dangerous, but cheap and effective way to access first element in
991 
992  const Epetra_CrsGraph & Graph = A->Graph();
993  ind = Graph[0]; // list of values and indices
994 
995  Nrhs = 0; // Assume no rhs, lhs
996 
997  if (RHS!=0) {
998  Nrhs = RHS->NumVectors();
999  if (Nrhs>1)
1000  if (!RHS->ConstantStride()) {EPETRA_CHK_ERR(-2)}; // Must have strided vectors
1001  ldrhs = RHS->Stride();
1002  rhs = (*RHS)[0]; // Dangerous but effective (again)
1003  }
1004  if (LHS!=0) {
1005  int Nlhs = LHS->NumVectors();
1006  if (Nlhs!=Nrhs) {EPETRA_CHK_ERR(-3)}; // Must have same number of rhs and lhs
1007  if (Nlhs>1)
1008  if (!LHS->ConstantStride()) {EPETRA_CHK_ERR(-4)}; // Must have strided vectors
1009  ldlhs = LHS->Stride();
1010  lhs = (*LHS)[0];
1011  }
1012 
1013  // Finally build ptr vector
1014 
1015  if (ptr==0) {
1016  ptr = new int[M+1];
1017  ptr[0] = 0;
1018  for (int i=0; i<M; i++) ptr[i+1] = ptr[i] + Graph.NumMyIndices(i);
1019  }
1020  EPETRA_CHK_ERR(ierr);
1021  return(0);
1022 }
1023 
1024 // Explicitly instantiate these two even though a call to them is made.
1025 // Possible fix for a bug reported by Kenneth Belcourt.
1026 template void Epetra_Util::Sort<int> (bool, int, int *, int, double **, int, int **, int, long long **);
1027 template void Epetra_Util::Sort<long long>(bool, int, long long *, int, double **, int, int **, int, long long **);
1028 
1029 
const int * LengthsFrom() const
Number of values we&#39;re receiving from each proc.
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
long long MaxAllGID64() const
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:119
unsigned int Seed_
Definition: Epetra_Util.h:265
int NumRemoteIDs() const
Returns the number of elements that are not on the calling processor.
bool GlobalIndicesLongLong() const
Returns true if map create with long long NumGlobalElements.
unsigned int Seed() const
Get seed from Random function.
Epetra_Distributor: The Epetra Gather/Scatter Setup Base Class.
static Epetra_Map Create_OneToOne_Map(const Epetra_Map &usermap, bool high_rank_proc_owns_shared=false)
Epetra_Util Create_OneToOne_Map function.
int NumReceives() const
The number of procs from which we will receive data.
long long IndexBase64() const
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
long long NumGlobalElements64() const
bool UniqueGIDs() const
Returns true if map GIDs are 1-to-1.
virtual Epetra_Directory * CreateDirectory(const Epetra_BlockMap &Map) const =0
Create a directory object for the given Epetra_BlockMap.
static const double chopVal_
Definition: Epetra_Util.h:262
#define EPETRA_CHK_ERR(a)
Epetra_Distributor & Distributor() const
bool GlobalIndicesInt() const
Returns true if map create with int NumGlobalElements.
MPI implementation of Epetra_Distributor.
virtual int MinAll(double *PartialMins, double *GlobalMins, int Count) const =0
Epetra_Comm Global Min function.
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
Epetra_Import: This class builds an import object for efficient importing of off-processor elements...
Definition: Epetra_Import.h:63
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.
int NumVectors() const
Returns the number of vectors in the multi-vector.
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
int IndexBase() const
Index base for this map.
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...
int NumMyCols() const
Returns the number of entries in the set of column-indices that appear on this processor.
static int GetPidGidPairs(const Epetra_Import &Importer, std::vector< std::pair< int, int > > &gpids, bool use_minus_one_for_local)
Epetra_Util GetPidGidPairs function.
long long GID64(int LID) const
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
const Epetra_BlockMap & TargetMap() const
Returns the TargetMap used to construct this importer.
static Epetra_Map Create_Root_Map(const Epetra_Map &usermap, int root=0)
Epetra_Util Create_Root_Map function.
int GID(int LID) const
Returns global ID of local ID, return IndexBase-1 if not found on this processor. ...
static Epetra_BlockMap Create_OneToOne_BlockMap64(const Epetra_BlockMap &usermap, bool high_rank_proc_owns_shared=false)
Epetra_Util Create_OneToOne_Map function for long long ordinals.
static Epetra_Map TCreate_Root_Map(const Epetra_Map &usermap, int root)
int NumMyRows() const
Returns the number of matrix rows owned by the calling processor.
Epetra_GIDTypeVector: A class for constructing and using dense &quot;int&quot; and &quot;long long&quot; vectors on a par...
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
bool IndicesAreContiguous() const
If matrix indices are packed into single array (done in OptimizeStorage()) return true...
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.
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.
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...
int Stride() const
Returns the stride between vectors in the multi-vector (only meaningful if ConstantStride() is true)...
bool ConstantStride() const
Returns true if this multi-vector has constant stride between vectors.
int NumMyIndices(int Row) const
Returns the current number of nonzero entries in specified local row on this processor.
static int GetRemotePIDs(const Epetra_Import &Importer, std::vector< int > &RemotePIDs)
Epetra_Util GetRemotePIDs.
const int * ProcsFrom() const
A list of procs sending values to this proc.
virtual int NumProc() const =0
Returns total number of processes.
int * RemoteLIDs() const
List of elements in the target map that are coming from other processors.
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...
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)
virtual int ReportError(const std::string Message, int ErrorCode) const
Error reporting method.
const Epetra_Distributor * DistributorPtr() const
const Epetra_CrsGraph & Graph() const
Returns a reference to the Epetra_CrsGraph object associated with this matrix.
long long * MyGlobalElements64() const
int n
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.
int NumMyNonzeros() const
Returns the number of nonzero entries in the calling processor&#39;s portion of the matrix.