Epetra Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Epetra_IntMultiVector.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_IntMultiVector.h"
45 #include "Epetra_IntVector.h"
46 #include "Epetra_Comm.h"
47 #ifdef EPETRA_MPI
48 #include "Epetra_MpiComm.h"
49 #endif
50 #include "Epetra_BlockMap.h"
51 #include "Epetra_Map.h"
52 #include "Epetra_Import.h"
53 #include "Epetra_Export.h"
54 #include "Epetra_Distributor.h"
55 
56 //=============================================================================
57 
58 // Epetra_BlockMap Constructor
59 
60 Epetra_IntMultiVector::Epetra_IntMultiVector(const Epetra_BlockMap& map, int numVectors, bool zeroOut)
61  : Epetra_DistObject(map, "Epetra::MultiVector"),
63  Values_(0),
64  Pointers_(0),
65  MyLength_(map.NumMyPoints()),
66  GlobalLength_(map.NumGlobalPoints64()),
67  NumVectors_(numVectors),
68  UserAllocated_(false),
69  ConstantStride_(true),
70  Stride_(map.NumMyPoints()),
71  Allocated_(false)
72 {
73  Util_.SetSeed(1);
74 
76 
77  for (int i = 0; i < NumVectors_; i++) Pointers_[i] = Values_+i*Stride_;
78 
79  if(zeroOut) PutScalar(0); // Fill all vectors with zero.
80 }
81 //==========================================================================
82 
83 // Copy Constructor
84 
86  : Epetra_DistObject(Source),
87  Epetra_CompObject(Source),
88  Values_(0),
89  Pointers_(0),
90  MyLength_(Source.MyLength_),
91  GlobalLength_(Source.GlobalLength_),
92  NumVectors_(Source.NumVectors_),
93  UserAllocated_(false),
94  ConstantStride_(true),
95  Stride_(0),
96  Allocated_(false),
97  Util_(Source.Util_)
98 {
100 
101  int ** Source_Pointers = Source.Pointers();
102  for (int i = 0; i< NumVectors_; i++) Pointers_[i] = Source_Pointers[i];
103 
104  DoCopy();
105 
106 }
107 //==========================================================================
108 
109 // This constructor copies in or makes view of a standard Fortran array
110 
112  int *A, int MyLDA, int numVectors)
113  : Epetra_DistObject(map, "Epetra::MultiVector"),
115  Values_(0),
116  Pointers_(0),
117  MyLength_(map.NumMyPoints()),
118  GlobalLength_(map.NumGlobalPoints64()),
119  NumVectors_(numVectors),
120  UserAllocated_(false),
121  ConstantStride_(true),
122  Stride_(map.NumMyPoints()),
123  Allocated_(false)
124 {
125  Util_.SetSeed(1);
126 
127  if (CV==Copy) AllocateForCopy();
128  else AllocateForView();
129 
130  for (int i = 0; i< NumVectors_; i++) Pointers_[i] = A + i*MyLDA;
131 
132  if (CV==Copy) DoCopy();
133  else DoView();
134 
135 }
136 
137 //==========================================================================
138 
139 // This constructor copies in or makes view of a C/C++ array of pointer
140 
142  int **ArrayOfPointers, int numVectors)
143  : Epetra_DistObject(map, "Epetra::MultiVector"),
145  Values_(0),
146  Pointers_(0),
147  MyLength_(map.NumMyPoints()),
148  GlobalLength_(map.NumGlobalPoints64()),
149  NumVectors_(numVectors),
150  UserAllocated_(false),
151  ConstantStride_(true),
152  Stride_(map.NumMyPoints()),
153  Allocated_(false)
154 {
155  Util_.SetSeed(1);
156 
157  if (CV==Copy) AllocateForCopy();
158  else AllocateForView();
159 
160  for (int i = 0; i< NumVectors_; i++) Pointers_[i] = ArrayOfPointers[i];
161 
162  if (CV==Copy) DoCopy();
163  else DoView();
164 
165 }
166 
167 //==========================================================================
168 
169 // This constructor copies or makes view of selected vectors, specified in Indices,
170 // from an existing MultiVector
171 
173  int *Indices, int numVectors)
174  : Epetra_DistObject(Source.Map(), "Epetra::MultiVector"),
176  Values_(0),
177  Pointers_(0),
178  MyLength_(Source.MyLength_),
179  GlobalLength_(Source.GlobalLength_),
180  NumVectors_(numVectors),
181  UserAllocated_(false),
182  ConstantStride_(true),
183  Stride_(0),
184  Allocated_(false)
185 {
186  Util_.SetSeed(1);
187 
188  if (CV==Copy) AllocateForCopy();
189  else AllocateForView();
190 
191  int ** Source_Pointers = Source.Pointers();
192  for (int i = 0; i< NumVectors_; i++) Pointers_[i] = Source_Pointers[Indices[i]];
193 
194  if (CV==Copy) DoCopy();
195  else DoView();
196 
197 }
198 
199 //==========================================================================
200 
201 // This interface copies or makes view of a range of vectors from an existing IntMultiVector
202 
204  int StartIndex, int numVectors)
205  : Epetra_DistObject(Source.Map(), "Epetra::MultiVector"),
207  Values_(0),
208  Pointers_(0),
209  MyLength_(Source.MyLength_),
210  GlobalLength_(Source.GlobalLength_),
211  NumVectors_(numVectors),
212  UserAllocated_(false),
213  ConstantStride_(true),
214  Stride_(0),
215  Allocated_(false)
216 {
217  Util_.SetSeed(1);
218 
219  if (CV==Copy) AllocateForCopy();
220  else AllocateForView();
221 
222  int ** Source_Pointers = Source.Pointers();
223  for (int i = 0; i< NumVectors_; i++) Pointers_[i] = Source_Pointers[StartIndex+i];
224 
225  if (CV==Copy) DoCopy();
226  else DoView();
227 }
228 
229 //=========================================================================
231 
232  if (!Allocated_) return;
233 
234  delete [] Pointers_;
235  if (!UserAllocated_ && Values_!=0) delete [] Values_;
236 
237  if (IntVectors_!=0) {
238  for (int i=0; i<NumVectors_; i++) if (IntVectors_[i]!=0) delete IntVectors_[i];
239  delete [] IntVectors_;
240  }
241 
242 
243  if (OrdinalTemp_!=0) delete [] OrdinalTemp_;
244 
245 }
246 
247 //=========================================================================
249 {
250 
251  if (Allocated_) return(0);
252 
253  if (NumVectors_<=0)
254  throw ReportError("Number of Vectors = " + toString(NumVectors_) + ", but must be greater than zero", -1);
255 
257  if (Stride_>0) Values_ = new int[Stride_ * NumVectors_];
258  Pointers_ = new int *[NumVectors_];
259 
260  OrdinalTemp_ = 0;
261  IntVectors_ = 0;
262 
263  int randval = rand(); // Use POSIX standard random function
264  if (DistributedGlobal())
265  Util_.SetSeed(2*Comm_->MyPID() + randval);
266  else {
267  int locrandval = randval;
268  Comm_->MaxAll(&locrandval, &randval, 1);
269  Util_.SetSeed(randval); // Replicated Local MultiVectors must have same seeds
270  }
271 
272  Allocated_ = true;
273  UserAllocated_ = false;
274  return(0);
275 }
276 
277 //=========================================================================
279 {
280  // On entry Pointers_ contains pointers to the incoming vectors. These
281  // pointers are the only unique piece of information for each of the
282  // constructors.
283 
284  // \internal { Optimization of this function can impact performance since it
285  // involves a fair amount of memory traffic.}
286 
287  // On exit, Pointers_ is redefined to point to its own MultiVector vectors.
288 
289  for (int i = 0; i< NumVectors_; i++)
290  {
291  int * from = Pointers_[i];
292  int * to = Values_+i*Stride_;
293  Pointers_[i] = to;
294  const int myLength = MyLength_;
295 #ifdef EPETRA_HAVE_OMP
296 #pragma omp parallel for default(none) shared(to,from)
297  for (int j=0; j<myLength; j++) to[j] = from[j];
298 #else
299  memcpy(to, from, myLength*sizeof(int));
300 #endif
301  }
302 
303  return(0);
304 }
305 //=========================================================================
307 {
308 
309  if (NumVectors_<=0)
310  throw ReportError("Number of Vectors = " + toString(NumVectors_) + ", but must be greater than zero", -1);
311 
312  Pointers_ = new int *[NumVectors_];
313 
314  OrdinalTemp_ = 0;
315  IntVectors_ = 0;
316 
317  int randval = rand(); // Use POSIX standard random function
318  if (DistributedGlobal())
319  Util_.SetSeed(2*Comm_->MyPID() + randval);
320  else {
321  int locrandval = randval;
322  Comm_->MaxAll(&locrandval, &randval, 1);
323  Util_.SetSeed(randval); // Replicated Local MultiVectors must have same seeds
324  }
325 
326  Allocated_ = true;
327  UserAllocated_ = true;
328 
329  return(0);
330 }
331 
332 //=========================================================================
334 {
335  // On entry Pointers_ contains pointers to the incoming vectors. These
336  // pointers are the only unique piece of information for each of the
337  // constructors.
338 
339 
340  Values_ = Pointers_[0];
341 
342  if (NumVectors_ == 1) {
344  ConstantStride_ = true;
345  return(0);
346  }
347 
348  // Remainder of code checks if MultiVector has regular stride
349 
350  Stride_ = (int)(Pointers_[1] - Pointers_[0]);
351  ConstantStride_ = false;
352 
353  for (int i = 1; i < NumVectors_-1; i++) {
354  if (Pointers_[i+1] - Pointers_[i] != Stride_) return(0);
355  }
356 
357  ConstantStride_ = true;
358 
359  return(0);
360 }
361 //=========================================================================
362 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
363 int Epetra_IntMultiVector::ReplaceGlobalValue(int GlobalRow, int VectorIndex, int OrdinalValue) {
364 
365  // Use the more general method below
366  EPETRA_CHK_ERR(ChangeGlobalValue<int>(GlobalRow, 0, VectorIndex, OrdinalValue, false));
367  return(0);
368 }
369 #endif
370 //=========================================================================
371 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
372 int Epetra_IntMultiVector::ReplaceGlobalValue(long long GlobalRow, int VectorIndex, int OrdinalValue) {
373 
374  // Use the more general method below
375  EPETRA_CHK_ERR(ChangeGlobalValue<long long>(GlobalRow, 0, VectorIndex, OrdinalValue, false));
376  return(0);
377 }
378 #endif
379 //=========================================================================
380 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
381 int Epetra_IntMultiVector::ReplaceGlobalValue(int GlobalBlockRow, int BlockRowOffset,
382  int VectorIndex, int OrdinalValue) {
383  // Use the more general method below
384  EPETRA_CHK_ERR(ChangeGlobalValue<int>(GlobalBlockRow, BlockRowOffset, VectorIndex, OrdinalValue, false));
385  return(0);
386 }
387 #endif
388 //=========================================================================
389 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
390 int Epetra_IntMultiVector::ReplaceGlobalValue(long long GlobalBlockRow, int BlockRowOffset,
391  int VectorIndex, int OrdinalValue) {
392  // Use the more general method below
393  EPETRA_CHK_ERR(ChangeGlobalValue<long long>(GlobalBlockRow, BlockRowOffset, VectorIndex, OrdinalValue, false));
394  return(0);
395 }
396 #endif
397 //=========================================================================
398 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
399 int Epetra_IntMultiVector::SumIntoGlobalValue(int GlobalRow, int VectorIndex, int OrdinalValue) {
400 
401  // Use the more general method below
402  EPETRA_CHK_ERR(ChangeGlobalValue<int>(GlobalRow, 0, VectorIndex, OrdinalValue, true));
403  return(0);
404 }
405 #endif
406 //=========================================================================
407 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
408 int Epetra_IntMultiVector::SumIntoGlobalValue(long long GlobalRow, int VectorIndex, int OrdinalValue) {
409 
410  // Use the more general method below
411  EPETRA_CHK_ERR(ChangeGlobalValue<long long>(GlobalRow, 0, VectorIndex, OrdinalValue, true));
412  return(0);
413 }
414 #endif
415 //=========================================================================
416 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
417 int Epetra_IntMultiVector::SumIntoGlobalValue(int GlobalBlockRow, int BlockRowOffset,
418  int VectorIndex, int OrdinalValue) {
419  // Use the more general method below
420  EPETRA_CHK_ERR(ChangeGlobalValue<int>(GlobalBlockRow, BlockRowOffset, VectorIndex, OrdinalValue, true));
421  return(0);
422 }
423 #endif
424 //=========================================================================
425 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
426 int Epetra_IntMultiVector::SumIntoGlobalValue(long long GlobalBlockRow, int BlockRowOffset,
427  int VectorIndex, int OrdinalValue) {
428  // Use the more general method below
429  EPETRA_CHK_ERR(ChangeGlobalValue<long long>(GlobalBlockRow, BlockRowOffset, VectorIndex, OrdinalValue, true));
430  return(0);
431 }
432 #endif
433 //=========================================================================
434 int Epetra_IntMultiVector::ReplaceMyValue(int MyRow, int VectorIndex, int OrdinalValue) {
435 
436  // Use the more general method below
437  EPETRA_CHK_ERR(ChangeMyValue(MyRow, 0, VectorIndex, OrdinalValue, false));
438  return(0);
439 }
440 //=========================================================================
441 int Epetra_IntMultiVector::ReplaceMyValue(int MyBlockRow, int BlockRowOffset,
442  int VectorIndex, int OrdinalValue) {
443  // Use the more general method below
444  EPETRA_CHK_ERR(ChangeMyValue(MyBlockRow, BlockRowOffset, VectorIndex, OrdinalValue, false));
445  return(0);
446 }
447 //=========================================================================
448 int Epetra_IntMultiVector::SumIntoMyValue(int MyRow, int VectorIndex, int OrdinalValue) {
449  // Use the more general method below
450  EPETRA_CHK_ERR(ChangeMyValue(MyRow, 0, VectorIndex, OrdinalValue, true));
451  return(0);
452 }
453 //=========================================================================
454 int Epetra_IntMultiVector::SumIntoMyValue(int MyBlockRow, int BlockRowOffset,
455  int VectorIndex, int OrdinalValue) {
456  // Use the more general method below
457  EPETRA_CHK_ERR(ChangeMyValue(MyBlockRow, BlockRowOffset, VectorIndex, OrdinalValue, true));
458  return(0);
459 }
460 //=========================================================================
461 template<typename int_type>
462 int Epetra_IntMultiVector::ChangeGlobalValue(int_type GlobalBlockRow, int BlockRowOffset,
463  int VectorIndex, int OrdinalValue, bool SumInto) {
464 
465  if(!Map().template GlobalIndicesIsType<int_type>())
466  throw ReportError("Epetra_IntMultiVector::ChangeGlobalValues mismatch between argument types (int/long long) and map type.", -1);
467 
468  // Convert GID to LID and call LID version
469  EPETRA_CHK_ERR(ChangeMyValue(Map().LID(GlobalBlockRow), BlockRowOffset, VectorIndex, OrdinalValue, SumInto));
470  return(0);
471 }
472 //=========================================================================
473 int Epetra_IntMultiVector::ChangeMyValue(int MyBlockRow, int BlockRowOffset,
474  int VectorIndex, int OrdinalValue, bool SumInto) {
475 
476  if (!Map().MyLID(MyBlockRow)) EPETRA_CHK_ERR(1); // I don't own this one, return a warning flag
477  if (VectorIndex>= NumVectors_) EPETRA_CHK_ERR(-1); // Consider this a real error
478  if (BlockRowOffset<0 || BlockRowOffset>=Map().ElementSize(MyBlockRow)) EPETRA_CHK_ERR(-2); // Offset is out-of-range
479 
480  int entry = Map().FirstPointInElement(MyBlockRow);
481 
482  if (SumInto)
483  Pointers_[VectorIndex][entry+BlockRowOffset] += OrdinalValue;
484  else
485  Pointers_[VectorIndex][entry+BlockRowOffset] = OrdinalValue;
486 
487  return(0);
488 }
489 
490 //=========================================================================
491 
492 // Extract a copy of a Epetra_IntMultiVector. Put in a user's Fortran-style array
493 
494 int Epetra_IntMultiVector::ExtractCopy(int *A, int MyLDA) const {
495  if (NumVectors_>1 && Stride_ > MyLDA) EPETRA_CHK_ERR(-1); // LDA not big enough
496 
497  const int myLength = MyLength_;
498  for (int i=0; i< NumVectors_; i++)
499  {
500  int * from = Pointers_[i];
501  int * to = A + i*MyLDA;
502  for (int j=0; j<myLength; j++) *to++ = *from++;
503  }
504 
505  return(0);
506 }
507 
508 //=========================================================================
510 {
511  // mfh 28 Mar 2013: We can't check for compatibility across the
512  // whole communicator, unless we know that the current and new
513  // Maps are nonnull on _all_ participating processes.
514 
515  // So, we'll check to make sure that the maps are the same size on this processor and then
516  // just go with it.
517  if(Map().NumMyElements() == map.NumMyElements() && Map().NumMyPoints() == map.NumMyPoints()) {
519  return(0);
520  }
521 
522  return(-1);
523 }
524 
525 //=========================================================================
526 
527 // Extract a copy of a Epetra_IntMultiVector. Put in a user's array of pointers
528 
529 int Epetra_IntMultiVector::ExtractCopy(int **ArrayOfPointers) const {
530  const int myLength = MyLength_;
531  for (int i=0; i< NumVectors_; i++)
532  {
533  int * from = Pointers_[i];
534  int * to = ArrayOfPointers[i];
535  memcpy(to, from, myLength*sizeof(int));
536  }
537 
538  return(0);
539 }
540 
541 
542 
543 //=========================================================================
544 
545 // Extract a view of a Epetra_IntMultiVector. Set up a user's Fortran-style array
546 
547 int Epetra_IntMultiVector::ExtractView(int **A, int *MyLDA) const {
548  if (!ConstantStride_) EPETRA_CHK_ERR(-1); // Can't make a Fortran-style view if not constant stride
549  *MyLDA = Stride_; // Set user's LDA
550  *A = Values_; // Set user's value pointer
551  return(0);
552 }
553 
554 
555 
556 //=========================================================================
557 
558 // Extract a view of a Epetra_IntMultiVector. Put in a user's array of pointers
559 
560 int Epetra_IntMultiVector::ExtractView(int ***ArrayOfPointers) const {
561  *ArrayOfPointers = Pointers_;
562 
563  return(0);
564 }
565 
566 
567 //=========================================================================
568 int Epetra_IntMultiVector::PutScalar(int ScalarConstant) {
569 
570  // Fills MultiVector with the value ScalarConstant **/
571 
572  const int myLength = MyLength_;
573  for (int i = 0; i < NumVectors_; i++) {
574  int * const to = Pointers_[i];
575 #ifdef EPETRA_HAVE_OMP
576 #pragma omp parallel for default(none) shared(ScalarConstant)
577 #endif
578  for (int j=0; j<myLength; j++) to[j] = ScalarConstant;
579  }
580  return(0);
581 }
582 //=========================================================================
584 
585  const Epetra_IntMultiVector & A = dynamic_cast<const Epetra_IntMultiVector &>(Source);
586 
587  if (NumVectors()!=A.NumVectors()) {EPETRA_CHK_ERR(-3)};
588  return(0);
589 }
590 
591 //=========================================================================
593  int NumSameIDs,
594  int NumPermuteIDs,
595  int * PermuteToLIDs,
596  int *PermuteFromLIDs,
597  const Epetra_OffsetIndex * Indexor,
598  Epetra_CombineMode CombineMode)
599 {
600  (void)Indexor;
601 
602  const Epetra_IntMultiVector & A = dynamic_cast<const Epetra_IntMultiVector &>(Source);
603 
604  int **From = A.Pointers();
605  int **To = Pointers_;
606  int numVectors = NumVectors_;
607 
608  int * ToFirstPointInElementList = 0;
609  int * FromFirstPointInElementList = 0;
610  int * FromElementSizeList = 0;
611  int MaxElementSize = Map().MaxElementSize();
612  bool ConstantElementSize = Map().ConstantElementSize();
613 
614  if (!ConstantElementSize) {
615  ToFirstPointInElementList = Map().FirstPointInElementList();
616  FromFirstPointInElementList = A.Map().FirstPointInElementList();
617  FromElementSizeList = A.Map().ElementSizeList();
618  }
619  int jj, jjj, k;
620 
621  int NumSameEntries;
622 
623  bool Case1 = false;
624  bool Case2 = false;
625  // bool Case3 = false;
626 
627  if (MaxElementSize==1) {
628  Case1 = true;
629  NumSameEntries = NumSameIDs;
630  }
631  else if (ConstantElementSize) {
632  Case2 = true;
633  NumSameEntries = NumSameIDs * MaxElementSize;
634  }
635  else {
636  // Case3 = true;
637  NumSameEntries = FromFirstPointInElementList[NumSameIDs];
638  }
639 
640  // Short circuit for the case where the source and target vector is the same.
641  if (To==From) NumSameEntries = 0;
642 
643  // Do copy first
644  if (NumSameIDs>0)
645  if (To!=From) {
646  for (int i=0; i < numVectors; i++) {
647  if (CombineMode==Epetra_AddLocalAlso) for (int j=0; j<NumSameEntries; j++) To[i][j] += From[i][j]; // Add to existing value
648  else for (int j=0; j<NumSameEntries; j++) To[i][j] = From[i][j];
649  }
650  }
651  // Do local permutation next
652  if (NumPermuteIDs>0) {
653 
654  // Point entry case
655  if (Case1) {
656 
657  if (numVectors==1) {
658  if (CombineMode==Epetra_AddLocalAlso) for (int j=0; j<NumPermuteIDs; j++) To[0][PermuteToLIDs[j]] += From[0][PermuteFromLIDs[j]]; // Add to existing value
659  else for (int j=0; j<NumPermuteIDs; j++) To[0][PermuteToLIDs[j]] = From[0][PermuteFromLIDs[j]];
660  }
661  else {
662  for (int j=0; j<NumPermuteIDs; j++) {
663  jj = PermuteToLIDs[j];
664  jjj = PermuteFromLIDs[j];
665  if (CombineMode==Epetra_AddLocalAlso) for (int i=0; i<numVectors; i++) To[i][jj] += From[i][jjj]; // Add to existing value
666  else for (int i=0; i<numVectors; i++) To[i][jj] = From[i][jjj];
667  }
668  }
669  }
670  // constant element size case
671  else if (Case2) {
672 
673  for (int j=0; j<NumPermuteIDs; j++) {
674  jj = MaxElementSize*PermuteToLIDs[j];
675  jjj = MaxElementSize*PermuteFromLIDs[j];
676  if (CombineMode==Epetra_AddLocalAlso) for (int i=0; i<numVectors; i++) for (k=0; k<MaxElementSize; k++) To[i][jj+k] += From[i][jjj+k]; // Add to existing value
677  else for(int i=0; i<numVectors; i++) for (k=0; k<MaxElementSize; k++) To[i][jj+k] = From[i][jjj+k];
678  }
679  }
680 
681  // variable element size case
682  else {
683 
684  for (int j=0; j<NumPermuteIDs; j++) {
685  jj = ToFirstPointInElementList[PermuteToLIDs[j]];
686  jjj = FromFirstPointInElementList[PermuteFromLIDs[j]];
687  int ElementSize = FromElementSizeList[PermuteFromLIDs[j]];
688  if (CombineMode==Epetra_AddLocalAlso) for (int i=0; i<numVectors; i++) for (k=0; k<ElementSize; k++) To[i][jj+k] += From[i][jjj+k]; // Add to existing value
689  else for (int i=0; i<numVectors; i++) for (k=0; k<ElementSize; k++) To[i][jj+k] = From[i][jjj+k];
690  }
691  }
692  }
693  return(0);
694 }
695 
696 //=========================================================================
698  int NumExportIDs,
699  int * ExportLIDs,
700  int & LenExports,
701  char * & Exports,
702  int & SizeOfPacket,
703  int * Sizes,
704  bool & VarSizes,
705  Epetra_Distributor & Distor)
706 {
707  (void)Sizes;
708  (void)VarSizes;
709  (void)Distor;
710 
711  const Epetra_IntMultiVector & A = dynamic_cast<const Epetra_IntMultiVector &>(Source);
712  int jj, k;
713 
714  int **From = A.Pointers();
715  int MaxElementSize = Map().MaxElementSize();
716  int numVectors = NumVectors_;
717  bool ConstantElementSize = Map().ConstantElementSize();
718 
719  int * FromFirstPointInElementList = 0;
720  int * FromElementSizeList = 0;
721 
722  if (!ConstantElementSize) {
723  FromFirstPointInElementList = A.Map().FirstPointInElementList();
724  FromElementSizeList = A.Map().ElementSizeList();
725  }
726 
727  int * OrdinalExports = 0;
728 
729  SizeOfPacket = numVectors*MaxElementSize*(int)sizeof(int);
730 
731  if(SizeOfPacket*NumExportIDs>LenExports) {
732  if (LenExports>0) delete [] Exports;
733  LenExports = SizeOfPacket*NumExportIDs;
734  OrdinalExports = new int[numVectors*MaxElementSize*NumExportIDs];
735  Exports = (char *) OrdinalExports;
736  }
737 
738  int * ptr;
739 
740  if (NumExportIDs>0) {
741  ptr = (int *) Exports;
742 
743  // Point entry case
744  if (MaxElementSize==1) {
745 
746  if (numVectors==1)
747  for (int j=0; j<NumExportIDs; j++)
748  *ptr++ = From[0][ExportLIDs[j]];
749 
750  else {
751  for (int j=0; j<NumExportIDs; j++) {
752  jj = ExportLIDs[j];
753  for (int i=0; i<numVectors; i++)
754  *ptr++ = From[i][jj];
755  }
756  }
757  }
758 
759  // constant element size case
760  else if (ConstantElementSize) {
761 
762  for (int j=0; j<NumExportIDs; j++) {
763  jj = MaxElementSize*ExportLIDs[j];
764  for (int i=0; i<numVectors; i++)
765  for (k=0; k<MaxElementSize; k++)
766  *ptr++ = From[i][jj+k];
767  }
768  }
769 
770  // variable element size case
771  else {
772 
773  int thisSizeOfPacket = numVectors*MaxElementSize;
774  for (int j=0; j<NumExportIDs; j++) {
775  ptr = (int *) Exports + j*thisSizeOfPacket;
776  jj = FromFirstPointInElementList[ExportLIDs[j]];
777  int ElementSize = FromElementSizeList[ExportLIDs[j]];
778  for (int i=0; i<numVectors; i++)
779  for (k=0; k<ElementSize; k++)
780  *ptr++ = From[i][jj+k];
781  }
782  }
783  }
784 
785  return(0);
786 }
787 
788 //=========================================================================
790  int NumImportIDs,
791  int * ImportLIDs,
792  int LenImports,
793  char * Imports,
794  int & SizeOfPacket,
795  Epetra_Distributor & Distor,
796  Epetra_CombineMode CombineMode,
797  const Epetra_OffsetIndex * Indexor )
798 {
799  (void)Source;
800  (void)LenImports;
801  (void)SizeOfPacket;
802  (void)Distor;
803  (void)Indexor;
804  int jj, k;
805 
806  if( CombineMode != Add
807  && CombineMode != Zero
808  && CombineMode != Insert
809  && CombineMode != InsertAdd
810  && CombineMode != Average
811  && CombineMode != Epetra_Max
812  && CombineMode != Epetra_Min
813  && CombineMode != AbsMin
814  && CombineMode != AbsMax )
815  EPETRA_CHK_ERR(-1); //Unsupported CombinedMode, will default to Zero
816 
817  if (NumImportIDs<=0) return(0);
818 
819  int ** To = Pointers_;
820  int numVectors = NumVectors_;
821  int MaxElementSize = Map().MaxElementSize();
822  bool ConstantElementSize = Map().ConstantElementSize();
823 
824  int * ToFirstPointInElementList = 0;
825  int * ToElementSizeList = 0;
826 
827  if (!ConstantElementSize) {
828  ToFirstPointInElementList = Map().FirstPointInElementList();
829  ToElementSizeList = Map().ElementSizeList();
830  }
831 
832  int * ptr;
833  // Unpack it...
834 
835  ptr = (int *) Imports;
836 
837  // Point entry case
838  if (MaxElementSize==1) {
839 
840  if (numVectors==1) {
841  if (CombineMode==InsertAdd) for (int j=0; j<NumImportIDs; j++) To[0][ImportLIDs[j]] = 0.0; // Zero out first
842  if (CombineMode==Add || CombineMode==InsertAdd) for (int j=0; j<NumImportIDs; j++) To[0][ImportLIDs[j]] += *ptr++; // Add to existing value
843  else if(CombineMode==Insert) for (int j=0; j<NumImportIDs; j++) To[0][ImportLIDs[j]] = *ptr++;
844  else if(CombineMode==AbsMax) for (int j=0; j<NumImportIDs; j++) {To[0][ImportLIDs[j]] = EPETRA_MAX( To[0][ImportLIDs[j]],std::abs(*ptr)); ptr++;}
845  else if(CombineMode==AbsMin) for (int j=0; j<NumImportIDs; j++) {To[0][ImportLIDs[j]] = EPETRA_MIN( To[0][ImportLIDs[j]],std::abs(*ptr)); ptr++;}
846  else if(CombineMode==Epetra_Max) for (int j=0; j<NumImportIDs; j++) {To[0][ImportLIDs[j]] = EPETRA_MAX( To[0][ImportLIDs[j]],*ptr); ptr++;}
847  else if(CombineMode==Epetra_Min) for (int j=0; j<NumImportIDs; j++) {To[0][ImportLIDs[j]] = EPETRA_MIN( To[0][ImportLIDs[j]],*ptr); ptr++;}
848  else if(CombineMode==Average) for (int j=0; j<NumImportIDs; j++) {To[0][ImportLIDs[j]] += *ptr++; To[0][ImportLIDs[j]] *= 0.5;} // Not a true avg if >2 occurance of an ID
849  // Note: The following form of averaging is not a true average if more that one value is combined.
850  // This might be an issue in the future, but we leave this way for now.
851 /*
852  if (CombineMode==Add)
853  for (int j=0; j<NumImportIDs; j++) To[0][ImportLIDs[j]] += *ptr++; // Add to existing value
854  else if(CombineMode==Insert)
855  for (int j=0; j<NumImportIDs; j++) To[0][ImportLIDs[j]] = *ptr++;
856  else if(CombineMode==InsertAdd) {
857  for (int j=0; j<NumImportIDs; j++) To[0][ImportLIDs[j]] = 0.0;
858  for (int j=0; j<NumImportIDs; j++) To[0][ImportLIDs[j]] += *ptr++;
859  }
860  else if(CombineMode==AbsMax)
861  for (int j=0; j<NumImportIDs; j++) {
862  To[0][ImportLIDs[j]] = EPETRA_MAX( To[0][ImportLIDs[j]],std::abs(*ptr));
863  ptr++;
864  }
865  // Note: The following form of averaging is not a true average if more that one value is combined.
866  // This might be an issue in the future, but we leave this way for now.
867  else if(CombineMode==Average)
868  for (int j=0; j<NumImportIDs; j++) {To[0][ImportLIDs[j]] += *ptr++; To[0][ImportLIDs[j]] *= 0.5;}
869 */
870  }
871 
872  else { // numVectors>1
873 
874  for (int j=0; j<NumImportIDs; j++) {
875  jj = ImportLIDs[j];
876  for (int i=0; i<numVectors; i++) {
877  if (CombineMode==InsertAdd) To[i][jj] = 0.0; // Zero out if needed
878  if (CombineMode==Add || CombineMode==InsertAdd) To[i][jj] += *ptr++; // Add to existing value
879  else if (CombineMode==Insert) To[i][jj] = *ptr++; // Insert values
880  else if (CombineMode==AbsMax) {To[i][jj] = EPETRA_MAX( To[i][jj], std::abs(*ptr)); ptr++; } // max of absolutes
881  else if (CombineMode==AbsMin) {To[i][jj] = EPETRA_MIN( To[i][jj], std::abs(*ptr)); ptr++; } // max of absolutes
882  else if (CombineMode==Epetra_Max) {To[i][jj] = EPETRA_MAX( To[i][jj], *ptr); ptr++; } // simple max
883  else if (CombineMode==Epetra_Min) {To[i][jj] = EPETRA_MIN( To[i][jj], *ptr); ptr++; } // simple min
884  else if (CombineMode==Average){To[i][jj] += *ptr++; To[i][jj] *= 0.5;}} // Not a true avg if >2 occurance of an ID
885 
886  }
887 /*
888  if (CombineMode==Add) {
889  for (int j=0; j<NumImportIDs; j++) {
890  jj = ImportLIDs[j];
891  for (int i=0; i<numVectors; i++)
892  To[i][jj] += *ptr++; // Add to existing value
893  }
894  }
895  else if(CombineMode==Insert) {
896  for (int j=0; j<NumImportIDs; j++) {
897  jj = ImportLIDs[j];
898  for (int i=0; i<numVectors; i++)
899  To[i][jj] = *ptr++;
900  }
901  }
902  else if(CombineMode==InsertAdd) {
903  for (int j=0; j<NumImportIDs; j++) {
904  jj = ImportLIDs[j];
905  for (int i=0; i<numVectors; i++)
906  To[i][jj] = 0.0;
907  }
908  for (int j=0; j<NumImportIDs; j++) {
909  jj = ImportLIDs[j];
910  for (int i=0; i<numVectors; i++)
911  To[i][jj] += *ptr++;
912  }
913  }
914  else if(CombineMode==AbsMax) {
915  for (int j=0; j<NumImportIDs; j++) {
916  jj = ImportLIDs[j];
917  for (int i=0; i<numVectors; i++) {
918  To[i][jj] = EPETRA_MAX( To[i][jj], std::abs(*ptr) );
919  ptr++;
920  }
921  }
922  }
923  // Note: The following form of averaging is not a true average if more that one value is combined.
924  // This might be an issue in the future, but we leave this way for now.
925  else if(CombineMode==Average) {
926  for (int j=0; j<NumImportIDs; j++) {
927  jj = ImportLIDs[j];
928  for (int i=0; i<numVectors; i++)
929  { To[i][jj] += *ptr++; To[i][jj] *= 0.5;}
930  }
931  }
932 */
933  }
934  }
935 
936  // constant element size case
937 
938  else if (ConstantElementSize) {
939 
940  for (int j=0; j<NumImportIDs; j++) {
941  jj = MaxElementSize*ImportLIDs[j];
942  for (int i=0; i<numVectors; i++) {
943  if (CombineMode==InsertAdd) for (k=0; k<MaxElementSize; k++) To[i][jj+k] = 0.0; // Zero out if needed
944  if (CombineMode==Add || CombineMode==InsertAdd) for (k=0; k<MaxElementSize; k++) To[i][jj+k] += *ptr++; // Add to existing value
945  else if (CombineMode==Insert) for (k=0; k<MaxElementSize; k++) To[i][jj+k] = *ptr++; // Insert values
946  else if (CombineMode==AbsMax) {for (k=0; k<MaxElementSize; k++) { To[i][jj+k] = EPETRA_MAX( To[i][jj+k], std::abs(*ptr)); ptr++; }} // max of absolutes
947  else if (CombineMode==AbsMin) {for (k=0; k<MaxElementSize; k++) { To[i][jj+k] = EPETRA_MIN( To[i][jj+k], std::abs(*ptr)); ptr++; }} // max of absolutes
948  else if (CombineMode==Epetra_Max) {for (k=0; k<MaxElementSize; k++) { To[i][jj+k] = EPETRA_MAX( To[i][jj+k], *ptr); ptr++; }} // simple max
949  else if (CombineMode==Epetra_Min) {for (k=0; k<MaxElementSize; k++) { To[i][jj+k] = EPETRA_MIN( To[i][jj+k], *ptr); ptr++; }} // simple min
950  else if (CombineMode==Average) {for (k=0; k<MaxElementSize; k++) { To[i][jj+k] += *ptr++; To[i][jj+k] *= 0.5;}} // Not a true avg if >2 occurance of an ID
951  }
952  }
953 /*
954  if (CombineMode==Add) {
955  for (int j=0; j<NumImportIDs; j++) {
956  jj = MaxElementSize*ImportLIDs[j];
957  for (int i=0; i<numVectors; i++)
958  for (k=0; k<MaxElementSize; k++)
959  To[i][jj+k] += *ptr++; // Add to existing value
960  }
961  }
962  else if(CombineMode==Insert) {
963  for (int j=0; j<NumImportIDs; j++) {
964  jj = MaxElementSize*ImportLIDs[j];
965  for (int i=0; i<numVectors; i++)
966  for (k=0; k<MaxElementSize; k++)
967  To[i][jj+k] = *ptr++;
968  }
969  }
970  else if(CombineMode==InsertAdd) {
971  for (int j=0; j<NumImportIDs; j++) {
972  jj = MaxElementSize*ImportLIDs[j];
973  for (int i=0; i<numVectors; i++)
974  for (k=0; k<MaxElementSize; k++)
975  To[i][jj+k] = 0.0;
976  }
977  for (int j=0; j<NumImportIDs; j++) {
978  jj = MaxElementSize*ImportLIDs[j];
979  for (int i=0; i<numVectors; i++)
980  for (k=0; k<MaxElementSize; k++)
981  To[i][jj+k] += *ptr++;
982  }
983  }
984  else if(CombineMode==AbsMax) {
985  for (int j=0; j<NumImportIDs; j++) {
986  jj = MaxElementSize*ImportLIDs[j];
987  for (int i=0; i<numVectors; i++)
988  for (k=0; k<MaxElementSize; k++) {
989  To[i][jj+k] = EPETRA_MAX( To[i][jj+k], std::abs(*ptr) );
990  ptr++;
991  }
992  }
993  }
994  // Note: The following form of averaging is not a true average if more that one value is combined.
995  // This might be an issue in the future, but we leave this way for now.
996  else if(CombineMode==Average) {
997  for (int j=0; j<NumImportIDs; j++) {
998  jj = MaxElementSize*ImportLIDs[j];
999  for (int i=0; i<numVectors; i++)
1000  for (k=0; k<MaxElementSize; k++)
1001  { To[i][jj+k] += *ptr++; To[i][jj+k] *= 0.5;}
1002  }
1003  }
1004 */
1005  }
1006 
1007  // variable element size case
1008 
1009  else {
1010  int thisSizeOfPacket = numVectors*MaxElementSize;
1011 
1012  for (int j=0; j<NumImportIDs; j++) {
1013  ptr = (int *) Imports + j*thisSizeOfPacket;
1014  jj = ToFirstPointInElementList[ImportLIDs[j]];
1015  int ElementSize = ToElementSizeList[ImportLIDs[j]];
1016  for (int i=0; i<numVectors; i++) {
1017  if (CombineMode==InsertAdd) for (k=0; k<ElementSize; k++) To[i][jj+k] = 0.0; // Zero out if needed
1018  if (CombineMode==Add || CombineMode==InsertAdd) for (k=0; k<ElementSize; k++) To[i][jj+k] += *ptr++; // Add to existing value
1019  else if (CombineMode==Insert) for (k=0; k<ElementSize; k++) To[i][jj+k] = *ptr++; // Insert values
1020  else if (CombineMode==AbsMax) {for (k=0; k<ElementSize; k++) { To[i][jj+k] = EPETRA_MAX( To[i][jj+k], std::abs(*ptr)); ptr++; }} // max of absolutes
1021  else if (CombineMode==Epetra_Max) {for (k=0; k<ElementSize; k++) { To[i][jj+k] = EPETRA_MAX( To[i][jj+k], *ptr); ptr++; }} // simple max
1022  else if (CombineMode==Epetra_Min) {for (k=0; k<ElementSize; k++) { To[i][jj+k] = EPETRA_MIN( To[i][jj+k], *ptr); ptr++; }} // simple min
1023  else if (CombineMode==Average) {for (k=0; k<ElementSize; k++) { To[i][jj+k] += *ptr++; To[i][jj+k] *= 0.5;}} // Not a true avg if >2 occurance of an ID
1024  }
1025  }
1026 /*
1027  if (CombineMode==Add) {
1028  for (int j=0; j<NumImportIDs; j++) {
1029  ptr = (double *) Imports + j*thisSizeOfPacket;
1030  jj = ToFirstPointInElementList[ImportLIDs[j]];
1031  int ElementSize = ToElementSizeList[ImportLIDs[j]];
1032  for (int i=0; i<numVectors; i++)
1033  for (k=0; k<ElementSize; k++)
1034  To[i][jj+k] += *ptr++; // Add to existing value
1035  }
1036  }
1037  else if(CombineMode==Insert){
1038  for (int j=0; j<NumImportIDs; j++) {
1039  ptr = (double *) Imports + j*thisSizeOfPacket;
1040  jj = ToFirstPointInElementList[ImportLIDs[j]];
1041  int ElementSize = ToElementSizeList[ImportLIDs[j]];
1042  for (int i=0; i<numVectors; i++)
1043  for (k=0; k<ElementSize; k++)
1044  To[i][jj+k] = *ptr++;
1045  }
1046  }
1047  else if(CombineMode==InsertAdd){
1048  for (int j=0; j<NumImportIDs; j++) {
1049  ptr = (double *) Imports + j*thisSizeOfPacket;
1050  jj = ToFirstPointInElementList[ImportLIDs[j]];
1051  int ElementSize = ToElementSizeList[ImportLIDs[j]];
1052  for (int i=0; i<numVectors; i++)
1053  for (k=0; k<ElementSize; k++)
1054  To[i][jj+k] = 0.0;
1055  }
1056  for (int j=0; j<NumImportIDs; j++) {
1057  ptr = (double *) Imports + j*thisSizeOfPacket;
1058  jj = ToFirstPointInElementList[ImportLIDs[j]];
1059  int ElementSize = ToElementSizeList[ImportLIDs[j]];
1060  for (int i=0; i<numVectors; i++)
1061  for (k=0; k<ElementSize; k++)
1062  To[i][jj+k] += *ptr++;
1063  }
1064  }
1065  else if(CombineMode==AbsMax){
1066  for (int j=0; j<NumImportIDs; j++) {
1067  ptr = (double *) Imports + j*thisSizeOfPacket;
1068  jj = ToFirstPointInElementList[ImportLIDs[j]];
1069  int ElementSize = ToElementSizeList[ImportLIDs[j]];
1070  for (int i=0; i<numVectors; i++)
1071  for (k=0; k<ElementSize; k++) {
1072  To[i][jj+k] = EPETRA_MAX( To[i][jj+k], std::abs(*ptr));
1073  ptr++;
1074  }
1075  }
1076  }
1077  // Note: The following form of averaging is not a true average if more that one value is combined.
1078  // This might be an issue in the future, but we leave this way for now.
1079  else if(CombineMode==Average) {
1080  for (int j=0; j<NumImportIDs; j++) {
1081  ptr = (double *) Imports + j*thisSizeOfPacket;
1082  jj = ToFirstPointInElementList[ImportLIDs[j]];
1083  int ElementSize = ToElementSizeList[ImportLIDs[j]];
1084  for (int i=0; i<numVectors; i++)
1085  for (k=0; k<ElementSize; k++)
1086  { To[i][jj+k] += *ptr++; To[i][jj+k] *= 0.5;}
1087  }
1088  }
1089 */
1090  }
1091 
1092  return(0);
1093 }
1094 
1095 //=========================================================================
1096 int Epetra_IntMultiVector::MinValue (int* Result) const {
1097 
1098  // Minimum value of each vector in MultiVector
1099 
1100  int ierr = 0;
1101 
1102  const int myLength = MyLength_;
1104 
1105  for (int i=0; i < NumVectors_; i++)
1106  {
1107  const int * const from = Pointers_[i];
1108  int MinVal = 2000000000; // 2 billion is close to largest 32 bit int
1109  if (myLength>0) MinVal = from[0];
1110 #ifdef EPETRA_HAVE_OMP
1111 #pragma omp parallel default(none) shared(MinVal)
1112 {
1113  int localMinVal = MinVal;
1114 #pragma omp for
1115  for (int j=0; j< myLength; j++) localMinVal = EPETRA_MIN(localMinVal,from[j]);
1116 #pragma omp critical
1117  {
1118  MinVal = EPETRA_MIN(MinVal,localMinVal);
1119  }
1120 }
1121  OrdinalTemp_[i] = MinVal;
1122 #else
1123  for (int j=0; j< myLength; j++) MinVal = EPETRA_MIN(MinVal,from[j]);
1124  OrdinalTemp_[i] = MinVal;
1125 #endif
1126  }
1127 
1128  if (myLength > 0) {
1129  for(int i=0; i<NumVectors_; ++i) {
1130  Result[i] = OrdinalTemp_[i];
1131  }
1132  }
1133 
1134  //If myLength == 0 and Comm_->NumProc() == 1, then Result has
1135  //not been referenced. Also, if vector contents are uninitialized
1136  //then Result contents are not well defined...
1137 
1138  if (Comm_->NumProc() == 1 || !DistributedGlobal()) return(ierr);
1139 
1140  //We're going to use MPI_Allgather to gather every proc's local-
1141  //min values onto every other proc. We'll use the last position
1142  //of the OrdinalTemp_ array to indicate whether this proc has
1143  //valid data that should be considered by other procs when forming
1144  //the global-min results.
1145 
1146  if (myLength == 0) OrdinalTemp_[NumVectors_] = 0;
1147  else OrdinalTemp_[NumVectors_] = 1;
1148 
1149  //Now proceed to handle the parallel case. We'll gather local-min
1150  //values from other procs and form a global-min. If any processor
1151  //has myLength>0, we'll end up with a valid result.
1152 
1153 #ifdef EPETRA_MPI
1154  const Epetra_MpiComm* epetrampicomm =
1155  dynamic_cast<const Epetra_MpiComm*>(Comm_);
1156  if (!epetrampicomm) {
1157  return(-2);
1158  }
1159 
1160  MPI_Comm mpicomm = epetrampicomm->GetMpiComm();
1161  int numProcs = epetrampicomm->NumProc();
1162  int* owork = new int[numProcs*(NumVectors_+1)];
1163 
1164  MPI_Allgather(OrdinalTemp_, NumVectors_+1, MPI_INT,
1165  owork, NumVectors_+1, MPI_INT, mpicomm);
1166 
1167  //if myLength==0, then our Result array currently contains
1168  //Epetra_MaxOrdinal from the local-min calculations above. In this
1169  //case we'll overwrite our Result array with values from the first
1170  //processor that sent valid data.
1171  bool overwrite = myLength == 0 ? true : false;
1172 
1173  int myPID = epetrampicomm->MyPID();
1174  int* owork_ptr = owork;
1175 
1176  for(int j=0; j<numProcs; ++j) {
1177 
1178  //skip data from self, and skip data from
1179  //procs with OrdinalTemp_[NumVectors_] == 0.
1180  if (j == myPID || owork_ptr[NumVectors_] == 0) {
1181  owork_ptr += NumVectors_+1;
1182  continue;
1183  }
1184 
1185  for(int i=0; i<NumVectors_; ++i) {
1186  int val = owork_ptr[i];
1187 
1188  //Set val into our Result array if overwrite is true (see above),
1189  //or if val is less than our current Result[i].
1190  if (overwrite || (Result[i] > val)) Result[i] = val;
1191  }
1192 
1193  //Now set overwrite to false so that we'll do the right thing
1194  //when processing data from subsequent processors.
1195  if (overwrite) overwrite = false;
1196 
1197  owork_ptr += NumVectors_+1;
1198  }
1199 
1200  delete [] owork;
1201 #endif
1202 
1203  // UpdateFlops(0); Strictly speaking there are not FLOPS in this routine
1204 
1205  return(ierr);
1206 }
1207 
1208 //=========================================================================
1209 int Epetra_IntMultiVector::MaxValue (int* Result) const {
1210 
1211  // Maximum value of each vector in MultiVector
1212 
1213  int ierr = 0;
1214 
1215  const int myLength = MyLength_;
1217 
1218  for (int i=0; i < NumVectors_; i++)
1219  {
1220  const int * const from = Pointers_[i];
1221  int MaxVal = -2000000000; // Negative 2 billion is close to smallest 32 bit int
1222  if (myLength>0) MaxVal = from[0];
1223 #ifdef EPETRA_HAVE_OMP
1224 #pragma omp parallel default(none) shared(MaxVal)
1225 {
1226  int localMaxVal = MaxVal;
1227 #pragma omp for
1228  for (int j=0; j< myLength; j++) localMaxVal = EPETRA_MAX(localMaxVal,from[j]);
1229 #pragma omp critical
1230  {
1231  MaxVal = EPETRA_MAX(MaxVal,localMaxVal);
1232  }
1233 }
1234  OrdinalTemp_[i] = MaxVal;
1235 #else
1236  for (int j=0; j< myLength; j++) MaxVal = EPETRA_MAX(MaxVal,from[j]);
1237  OrdinalTemp_[i] = MaxVal;
1238 #endif
1239  }
1240 
1241  if (myLength > 0) {
1242  for(int i=0; i<NumVectors_; ++i) {
1243  Result[i] = OrdinalTemp_[i];
1244  }
1245  }
1246 
1247  //If myLength == 0 and Comm_->NumProc() == 1, then Result has
1248  //not been referenced. Also, if vector contents are uninitialized
1249  //then Result contents are not well defined...
1250 
1251  if (Comm_->NumProc() == 1 || !DistributedGlobal()) return(ierr);
1252 
1253  //We're going to use MPI_Allgather to gather every proc's local-
1254  //max values onto every other proc. We'll use the last position
1255  //of the OrdinalTemp_ array to indicate whether this proc has
1256  //valid data that should be considered by other procs when forming
1257  //the global-max results.
1258 
1259  if (myLength == 0) OrdinalTemp_[NumVectors_] = 0.0;
1260  else OrdinalTemp_[NumVectors_] = 1.0;
1261 
1262  //Now proceed to handle the parallel case. We'll gather local-max
1263  //values from other procs and form a global-max. If any processor
1264  //has myLength>0, we'll end up with a valid result.
1265 
1266 #ifdef EPETRA_MPI
1267  const Epetra_MpiComm* epetrampicomm =
1268  dynamic_cast<const Epetra_MpiComm*>(Comm_);
1269  if (!epetrampicomm) {
1270  return(-2);
1271  }
1272 
1273  MPI_Comm mpicomm = epetrampicomm->GetMpiComm();
1274  int numProcs = epetrampicomm->NumProc();
1275  int* owork = new int[numProcs*(NumVectors_+1)];
1276 
1277  MPI_Allgather(OrdinalTemp_, NumVectors_+1, MPI_INT,
1278  owork, NumVectors_+1, MPI_INT, mpicomm);
1279 
1280  //if myLength==0, then our Result array currently contains
1281  //-Epetra_MaxOrdinal from the local-max calculations above. In this
1282  //case we'll overwrite our Result array with values from the first
1283  //processor that sent valid data.
1284  bool overwrite = myLength == 0 ? true : false;
1285 
1286  int myPID = epetrampicomm->MyPID();
1287  int* owork_ptr = owork;
1288 
1289  for(int j=0; j<numProcs; ++j) {
1290 
1291  //skip data from self, and skip data from
1292  //procs with OrdinalTemp_[NumVectors_] == 0.
1293  if (j == myPID || owork_ptr[NumVectors_] == 0) {
1294  owork_ptr += NumVectors_+1;
1295  continue;
1296  }
1297 
1298  for(int i=0; i<NumVectors_; ++i) {
1299  int val = owork_ptr[i];
1300 
1301  //Set val into our Result array if overwrite is true (see above),
1302  //or if val is larger than our current Result[i].
1303  if (overwrite || (Result[i] < val)) Result[i] = val;
1304  }
1305 
1306  //Now set overwrite to false so that we'll do the right thing
1307  //when processing data from subsequent processors.
1308  if (overwrite) overwrite = false;
1309 
1310  owork_ptr += NumVectors_+1;
1311  }
1312 
1313  delete [] owork;
1314 #endif
1315 
1316  // UpdateFlops(0); Strictly speaking there are not FLOPS in this routine
1317 
1318  return(ierr);
1319 }
1320 
1321 
1322 //=======================================================================
1324 
1325  // Epetra_IntMultiVector::operator () --- return non-const reference
1326 
1327  if (index < 0 || index >=NumVectors_)
1328  throw ReportError("Vector index = " + toString(index) + "is out of range. Number of Vectors = " + toString(NumVectors_), -1);
1329 
1330  UpdateIntVectors();
1331 
1332  // Create a new Epetra_IntVector that is a view of ith vector, if not already present
1333  if (IntVectors_[index]==0)
1334  IntVectors_[index] = new Epetra_IntVector(View, Map(), Pointers_[index]);
1335  return(IntVectors_[index]);
1336 }
1337 
1338 //=======================================================================
1340 
1341  // Epetra_IntMultiVector::operator () --- return non-const reference
1342 
1343  if (index < 0 || index >=NumVectors_)
1344  throw ReportError("Vector index = " + toString(index) + "is out of range. Number of Vectors = " + toString(NumVectors_), -1);
1345 
1346  UpdateIntVectors();
1347 
1348  if (IntVectors_[index]==0)
1349  IntVectors_[index] = new Epetra_IntVector(View, Map(), Pointers_[index]);
1350 
1351  const Epetra_IntVector * & temp = (const Epetra_IntVector * &) (IntVectors_[index]);
1352  return(temp);
1353 }
1354 
1355 //========================================================================
1357 
1358  // Check for special case of this=Source
1359  if (this != &Source) Assign(Source);
1360 
1361  return(*this);
1362 }
1363 
1364 //=========================================================================
1366 
1367  const int myLength = MyLength_;
1368  if (NumVectors_ != A.NumVectors())
1369  throw ReportError("Number of vectors incompatible in Assign. The this MultiVector has NumVectors = " + toString(NumVectors_)
1370  + ". The A MultiVector has NumVectors = " + toString(A.NumVectors()), -3);
1371  if (myLength != A.MyLength())
1372  throw ReportError("Length of MultiVectors incompatible in Assign. The this MultiVector has MyLength = " + toString(myLength)
1373  + ". The A MultiVector has MyLength = " + toString(A.MyLength()), -4);
1374 
1375  int ** A_Pointers = A.Pointers();
1376  for (int i = 0; i< NumVectors_; i++) {
1377  int * const to = Pointers_[i];
1378  const int * const from = A_Pointers[i];
1379 #ifdef EPETRA_HAVE_OMP
1380 #pragma omp parallel for default(none)
1381 #endif
1382  for (int j=0; j<myLength; j++) to[j] = from[j];
1383  }
1384  return;
1385  }
1386 
1387 //=========================================================================
1389 
1390  // Global reduction on each entry of a Replicated Local MultiVector
1391  const int myLength = MyLength_;
1392  int * source = 0;
1393  if (myLength>0) source = new int[myLength*NumVectors_];
1394  int * target = 0;
1395  bool packed = (ConstantStride_ && (Stride_==myLength));
1396  if (packed) {
1397  for (int i=0; i<myLength*NumVectors_; i++) source[i] = Values_[i];
1398  target = Values_;
1399  }
1400  else {
1401  int * tmp1 = source;
1402  for (int i = 0; i < NumVectors_; i++) {
1403  int * tmp2 = Pointers_[i];
1404  for (int j=0; j< myLength; j++) *tmp1++ = *tmp2++;
1405  }
1406  if (myLength>0) target = new int[myLength*NumVectors_];
1407  }
1408 
1409  Comm_->SumAll(source, target, myLength*NumVectors_);
1410  if (myLength>0) delete [] source;
1411  if (!packed) {
1412  int * tmp2 = target;
1413  for (int i = 0; i < NumVectors_; i++) {
1414  int * tmp1 = Pointers_[i];
1415  for (int j=0; j< myLength; j++) *tmp1++ = *tmp2++;
1416  }
1417  if (myLength>0) delete [] target;
1418  }
1419  // UpdateFlops(0); No serial Flops in this function
1420  return(0);
1421 }
1422 
1423 //=======================================================================
1424 int Epetra_IntMultiVector::ResetView(int ** ArrayOfPointers) {
1425 
1426  if (!UserAllocated_) {
1427  EPETRA_CHK_ERR(-1); // Can't reset view if multivector was not allocated as a view
1428  }
1429 
1430  for (int i = 0; i< NumVectors_; i++) Pointers_[i] = ArrayOfPointers[i];
1431  DoView();
1432 
1433  return(0);
1434 }
1435 
1436 //=======================================================================
1437 void Epetra_IntMultiVector::Print(std::ostream& os) const {
1438  int MyPID = Map().Comm().MyPID();
1439  int NumProc = Map().Comm().NumProc();
1440 
1441  for (int iproc=0; iproc < NumProc; iproc++) {
1442  if (MyPID==iproc) {
1443  int NumVectors1 = NumVectors();
1444  int NumMyElements1 =Map(). NumMyElements();
1445  int MaxElementSize1 = Map().MaxElementSize();
1446  int * FirstPointInElementList1 = NULL;
1447  if (MaxElementSize1!=1) FirstPointInElementList1 = Map().FirstPointInElementList();
1448  int ** A_Pointers = Pointers();
1449 
1450  if (MyPID==0) {
1451  os.width(8);
1452  os << " MyPID"; os << " ";
1453  os.width(12);
1454  if (MaxElementSize1==1)
1455  os << "GID ";
1456  else
1457  os << " GID/Point";
1458  for (int j = 0; j < NumVectors1 ; j++)
1459  {
1460  os.width(20);
1461  os << "Value ";
1462  }
1463  os << std::endl;
1464  }
1465  for (int i=0; i < NumMyElements1; i++) {
1466  for (int ii=0; ii< Map().ElementSize(i); ii++) {
1467  int iii;
1468  os.width(10);
1469  os << MyPID; os << " ";
1470  os.width(10);
1471  if (MaxElementSize1==1) {
1472  if(Map().GlobalIndicesInt())
1473  {
1474 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
1475  int * MyGlobalElements1 = Map().MyGlobalElements();
1476  os << MyGlobalElements1[i] << " ";
1477 #else
1478  throw ReportError("Epetra_IntMultiVector::Print: ERROR, GlobalIndicesInt but no API for it.",-1);
1479 #endif
1480  }
1481  else if(Map().GlobalIndicesLongLong())
1482  {
1483 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1484  long long * MyGlobalElements1 = Map().MyGlobalElements64();
1485  os << MyGlobalElements1[i] << " ";
1486 #else
1487  throw ReportError("Epetra_IntMultiVector::Print: ERROR, GlobalIndicesLongLong but no API for it.",-1);
1488 #endif
1489  }
1490  else
1491  throw ReportError("Epetra_IntMultiVector::Print ERROR, Don't know map global index type.",-1);
1492 
1493  iii = i;
1494  }
1495  else {
1496  if(Map().GlobalIndicesInt())
1497  {
1498 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
1499  int * MyGlobalElements1 = Map().MyGlobalElements();
1500  os << MyGlobalElements1[i]<< "/" << ii << " ";
1501 #else
1502  throw ReportError("Epetra_IntMultiVector::Print: ERROR, GlobalIndicesInt but no API for it.",-1);
1503 #endif
1504  }
1505  else if(Map().GlobalIndicesLongLong())
1506  {
1507 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1508  long long * MyGlobalElements1 = Map().MyGlobalElements64();
1509  os << MyGlobalElements1[i]<< "/" << ii << " ";
1510 #else
1511  throw ReportError("Epetra_IntMultiVector::Print: ERROR, GlobalIndicesLongLong but no API for it.",-1);
1512 #endif
1513  }
1514  else
1515  throw ReportError("Epetra_IntMultiVector::Print ERROR, Don't know map global index type.",-1);
1516 
1517  iii = FirstPointInElementList1[i]+ii;
1518  }
1519  for (int j = 0; j < NumVectors1 ; j++)
1520  {
1521  os.width(20);
1522  os << A_Pointers[j][iii];
1523  }
1524  os << std::endl;
1525  }
1526  }
1527  os << std::flush;
1528  }
1529 
1530  // Do a few global ops to give I/O a chance to complete
1531  Map().Comm().Barrier();
1532  Map().Comm().Barrier();
1533  Map().Comm().Barrier();
1534  }
1535  return;
1536 }
int ** Pointers() const
Get pointer to individual vector pointers.
Epetra_IntVector *& operator()(int i)
Vector access function.
int * FirstPointInElementList() const
Pointer to internal array containing a mapping between the local elements and the first local point n...
bool DistributedGlobal() const
Returns true if this multi-vector is distributed global, i.e., not local replicated.
int PutScalar(int OrdinalConstant)
Initialize all values in a multi-vector with constant value.
int CopyAndPermute(const Epetra_SrcDistObject &Source, int NumSameIDs, int NumPermuteIDs, int *PermuteToLIDs, int *PermuteFromLIDs, const Epetra_OffsetIndex *Indexor, Epetra_CombineMode CombineMode=Zero)
Perform ID copies and permutations that are on processor.
int ElementSize() const
Returns the size of elements in the map; only valid if map has constant element size.
MPI_Comm GetMpiComm() const
Get the MPI Communicator (identical to Comm() method; used when we know we are MPI.
Epetra_Distributor: The Epetra Gather/Scatter Setup Base Class.
Epetra_OffsetIndex: This class builds index for efficient mapping of data from one Epetra_CrsGraph ba...
Epetra_IntVector: A class for constructing and using dense integer vectors on a parallel computer...
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
int MyLength() const
Returns the local vector length on the calling processor of vectors in the multi-vector.
int MinValue(int *Result) const
Compute minimum value of each vector in multi-vector.
int NumProc() const
Returns total number of processes.
bool ConstantElementSize() const
Returns true if map has constant element size.
#define EPETRA_CHK_ERR(a)
int ExtractCopy(int *A, int MyLDA) const
Put multi-vector values into user-provided two-dimensional array.
Epetra_IntMultiVector(const Epetra_BlockMap &Map, int NumVectors, bool zeroOut=true)
Basic Epetra_IntMultiVector constuctor.
int * ElementSizeList() const
List of the element sizes corresponding to the array MyGlobalElements().
Epetra_IntMultiVector & operator=(const Epetra_IntMultiVector &Source)
= Operator.
int ExtractView(int **A, int *MyLDA) const
Set user-provided addresses of A and MyLDA.
int NumVectors() const
Returns the number of vectors in the multi-vector.
#define EPETRA_MIN(x, y)
Epetra_IntMultiVector: A class for constructing and using dense multi-vectors, vectors and matrices i...
Epetra_MpiComm: The Epetra MPI Communication Class.
virtual void Barrier() const =0
Epetra_Comm Barrier function.
int ReplaceGlobalValue(int GlobalRow, int VectorIndex, int OrdinalValue)
Replace current value at the specified (GlobalRow, VectorIndex) location with OrdinalValue.
virtual int MyPID() const =0
Return my process ID.
virtual ~Epetra_IntMultiVector()
Epetra_MultiVector destructor.
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const =0
Epetra_Comm Global Max function.
void Assign(const Epetra_IntMultiVector &rhs)
Epetra_IntVector ** IntVectors_
int ResetView(int **ArrayOfPointers)
Reset the view of an existing multivector to point to new user data.
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
Epetra_Comm Global Sum function.
int NumMyElements() const
Number of elements on the calling processor.
int SetSeed(unsigned int Seed_in)
Set seed for Random function.
std::string toString(const int &x) const
int FirstPointInElement(int LID) const
Returns the requested entry in the FirstPointInElementList; see FirstPointInElementList() for details...
int MaxValue(int *Result) const
Compute maximum value of each vector in multi-vector.
const Epetra_Comm * Comm_
Epetra_CompObject: Functionality and data that is common to all computational classes.
int ChangeGlobalValue(int_type GlobalBlockRow, int BlockRowOffset, int VectorIndex, int OrdinalValue, bool SumInto)
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
virtual void Print(std::ostream &os) const
Print method.
int PackAndPrepare(const Epetra_SrcDistObject &Source, int NumExportIDs, int *ExportLIDs, int &LenExports, char *&Exports, int &SizeOfPacket, int *Sizes, bool &VarSizes, Epetra_Distributor &Distor)
Perform any packing or preparation required for call to DoTransfer().
virtual int NumProc() const =0
Returns total number of processes.
Epetra_BlockMap Map_
int ChangeMyValue(int MyBlockRow, int BlockRowOffset, int VectorIndex, int OrdinalValue, bool SumInto)
int SumIntoGlobalValue(int GlobalRow, int VectorIndex, int OrdinalValue)
Adds OrdinalValue to existing value at the specified (GlobalRow, VectorIndex) location.
int MaxElementSize() const
Maximum element size across all processors.
Epetra_CombineMode
int CheckSizes(const Epetra_SrcDistObject &A)
Allows the source and target (this) objects to be compared for compatibility, return nonzero if not...
int SumIntoMyValue(int MyRow, int VectorIndex, int OrdinalValue)
Adds OrdinalValue to existing value at the specified (MyRow, VectorIndex) location.
int ReplaceMyValue(int MyRow, int VectorIndex, int OrdinalValue)
Replace current value at the specified (MyRow, VectorIndex) location with OrdinalValue.
virtual int ReportError(const std::string Message, int ErrorCode) const
Error reporting method.
int ReplaceMap(const Epetra_BlockMap &map)
Replace map, only if new map has same point-structure as current map.
int NumMyPoints() const
Number of local points for this map; equals the sum of all element sizes on the calling processor...
Epetra_SrcDistObject: A class for supporting flexible source distributed objects for import/export op...
long long * MyGlobalElements64() const
Epetra_DataAccess
int UnpackAndCombine(const Epetra_SrcDistObject &Source, int NumImportIDs, int *ImportLIDs, int LenImports, char *Imports, int &SizeOfPacket, Epetra_Distributor &Distor, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor)
Perform any unpacking and combining after call to DoTransfer().
#define EPETRA_MAX(x, y)
int MyPID() const
Return my process ID.
Epetra_DistObject: A class for constructing and using dense multi-vectors, vectors and matrices in pa...
const Epetra_BlockMap & Map() const
Returns the address of the Epetra_BlockMap for this multi-vector.