Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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"
55#include "Epetra_Import.h"
56#include <limits>
57
58#ifdef HAVE_MPI
60#endif
61
62const double Epetra_Util::chopVal_ = 1.0e-15;
63
64//=========================================================================
65double Epetra_Util::Chop(const double & Value){
66 if (std::abs(Value) < chopVal_) return 0;
67 return Value;
68}
69
70//=========================================================================
71unsigned 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//=========================================================================
102unsigned int Epetra_Util::Seed() const {
103 return(Seed_);
104}
105
106//=========================================================================
107int Epetra_Util::SetSeed(unsigned int Seed_in) {
108 Seed_ = Seed_in;
109 return(0);
110}
111
112//=============================================================================
113template<typename T>
114void 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
158void 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
169void 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
180void 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
192void 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
280template<typename int_type>
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//----------------------------------------------------------------------------
479int 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//----------------------------------------------------------------------------
517int 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//----------------------------------------------------------------------------
555int 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//----------------------------------------------------------------------------
616int 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
679int 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
729int 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//----------------------------------------------------------------------------
779int 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//----------------------------------------------------------------------------
825int 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//----------------------------------------------------------------------------
865template<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//----------------------------------------------------------------------------
907int 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//----------------------------------------------------------------------------
916template<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//----------------------------------------------------------------------------
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.
1026template void Epetra_Util::Sort<int> (bool, int, int *, int, double **, int, int **, int, long long **);
1027template void Epetra_Util::Sort<long long>(bool, int, long long *, int, double **, int, int **, int, long long **);
1028
1029
@ Insert
#define EPETRA_CHK_ERR(a)
static Epetra_Map TCreate_Root_Map(const Epetra_Map &usermap, int root)
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.
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.
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
long long IndexBase64() const
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
long long * MyGlobalElements64() const
long long MaxAllGID64() const
int GID(int LID) const
Returns global ID of local ID, return IndexBase-1 if not found on this processor.
bool LinearMap() const
Returns true if the global ID space is contiguously divided (but not necessarily uniformly) across al...
int IndexBase() const
Index base for this map.
long long GID64(int LID) const
bool UniqueGIDs() const
Returns true if map GIDs are 1-to-1.
long long NumGlobalElements64() const
bool GlobalIndicesInt() const
Returns true if map create with int NumGlobalElements.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
int NumMyElements() const
Number of elements on the calling processor.
bool IsOneToOne() const
bool GlobalIndicesLongLong() const
Returns true if map create with long long NumGlobalElements.
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition: Epetra_Comm.h:73
virtual int NumProc() const =0
Returns total number of processes.
virtual Epetra_Directory * CreateDirectory(const Epetra_BlockMap &Map) const =0
Create a directory object for the given Epetra_BlockMap.
virtual int MinAll(double *PartialMins, double *GlobalMins, int Count) const =0
Epetra_Comm Global Min function.
virtual int MyPID() const =0
Return my process ID.
Epetra_CrsGraph: A class for constructing and using sparse compressed row graphs.
int NumMyIndices(int Row) const
Returns the current number of nonzero entries in specified local row on this processor.
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
int NumMyCols() const
Returns the number of entries in the set of column-indices that appear on this processor.
int MakeDataContiguous()
Eliminates memory that is used for construction. Make consecutive row index sections contiguous.
int NumMyNonzeros() const
Returns the number of nonzero entries in the calling processor's portion of the matrix.
const Epetra_CrsGraph & Graph() const
Returns a reference to the Epetra_CrsGraph object associated with this matrix.
int NumMyRows() const
Returns the number of matrix rows owned by the calling processor.
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...
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_Distributor: The Epetra Gather/Scatter Setup Base Class.
Epetra_GIDTypeVector: A class for constructing and using dense "int" and "long long" vectors on a par...
Epetra_Import: This class builds an import object for efficient importing of off-processor elements.
Definition: Epetra_Import.h:63
Epetra_Distributor & Distributor() const
int * RemoteLIDs() const
List of elements in the target map that are coming from other processors.
int NumRemoteIDs() const
Returns the number of elements that are not on the calling processor.
const Epetra_Distributor * DistributorPtr() const
const Epetra_BlockMap & TargetMap() const
Returns the TargetMap used to construct this importer.
Epetra_LocalMap: A class for replicating vectors and matrices across multiple processors.
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:119
MPI implementation of Epetra_Distributor.
const int * ProcsFrom() const
A list of procs sending values to this proc.
const int * LengthsFrom() const
Number of values we're receiving from each proc.
int NumReceives() const
The number of procs from which we will receive data.
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
int NumVectors() const
Returns the number of vectors in the multi-vector.
bool ConstantStride() const
Returns true if this multi-vector has constant stride between vectors.
int Stride() const
Returns the stride between vectors in the multi-vector (only meaningful if ConstantStride() is true).
virtual int ReportError(const std::string Message, int ErrorCode) const
Error reporting method.
int SetSeed(unsigned int Seed_in)
Set seed for Random function.
static int SortCrsEntries(int NumRows, const int *CRS_rowptr, int *CRS_colind, double *CRS_vals)
Epetra_Util SortCrsEntries function.
unsigned int Seed_
Definition: Epetra_Util.h:265
static Epetra_BlockMap Create_OneToOne_BlockMap(const Epetra_BlockMap &usermap, bool high_rank_proc_owns_shared=false)
Epetra_Util Create_OneToOne_Map function.
static const double chopVal_
Definition: Epetra_Util.h:262
unsigned int Seed() const
Get seed from Random function.
static Epetra_Map Create_Root_Map(const Epetra_Map &usermap, int root=0)
Epetra_Util Create_Root_Map function.
static int GetPids(const Epetra_Import &Importer, std::vector< int > &pids, bool use_minus_one_for_local)
Epetra_Util GetPids function.
static int GetRemotePIDs(const Epetra_Import &Importer, std::vector< int > &RemotePIDs)
Epetra_Util GetRemotePIDs.
static int SortAndMergeCrsEntries(int NumRows, int *CRS_rowptr, int *CRS_colind, double *CRS_vals)
Epetra_Util SortAndMergeCrsEntries function.
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 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 double Chop(const double &Value)
Epetra_Util Chop method. Return zero if input Value is less than ChopValue.
Definition: Epetra_Util.cpp:65
double RandomDouble()
Returns a random double on the interval (-1.0,1.0)
Definition: Epetra_Util.cpp:90
unsigned int RandomInt()
Returns a random integer on the interval (0, 2^31-1)
Definition: Epetra_Util.cpp:71
static void EPETRA_LIB_DLL_EXPORT 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)
static Epetra_Map Create_OneToOne_Map(const Epetra_Map &usermap, bool high_rank_proc_owns_shared=false)
Epetra_Util Create_OneToOne_Map function.