Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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"
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
60Epetra_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 int myLength = MyLength_;
295#ifdef EPETRA_HAVE_OMP
296#pragma omp parallel for default(none) shared(to,from,myLength)
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
363int 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
372int 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
381int 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
390int 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
399int 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
408int 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
417int 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
426int 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//=========================================================================
434int 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//=========================================================================
441int 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//=========================================================================
448int 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//=========================================================================
454int 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//=========================================================================
461template<typename int_type>
462int 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//=========================================================================
473int 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
494int 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
529int 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
547int 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
560int Epetra_IntMultiVector::ExtractView(int ***ArrayOfPointers) const {
561 *ArrayOfPointers = Pointers_;
562
563 return(0);
564}
565
566
567//=========================================================================
568int Epetra_IntMultiVector::PutScalar(int ScalarConstant) {
569
570 // Fills MultiVector with the value ScalarConstant **/
571
572 int myLength = MyLength_;
573 for (int i = 0; i < NumVectors_; i++) {
574 int * to = Pointers_[i];
575#ifdef EPETRA_HAVE_OMP
576#pragma omp parallel for default(none) shared(ScalarConstant,myLength,to)
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//=========================================================================
1096int Epetra_IntMultiVector::MinValue (int* Result) const {
1097
1098 // Minimum value of each vector in MultiVector
1099
1100 int ierr = 0;
1101
1102 int myLength = MyLength_;
1104
1105 for (int i=0; i < NumVectors_; i++)
1106 {
1107 const int * 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,myLength,from)
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//=========================================================================
1209int Epetra_IntMultiVector::MaxValue (int* Result) const {
1210
1211 // Maximum value of each vector in MultiVector
1212
1213 int ierr = 0;
1214
1215 int myLength = MyLength_;
1217
1218 for (int i=0; i < NumVectors_; i++)
1219 {
1220 const int * 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,myLength,from)
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
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
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 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 * to = Pointers_[i];
1378 const int * from = A_Pointers[i];
1379#ifdef EPETRA_HAVE_OMP
1380#pragma omp parallel for default(none) shared(myLength,to,from)
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//=======================================================================
1424int 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//=======================================================================
1437void 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}
Epetra_CombineMode
@ Insert
@ AbsMin
@ AbsMax
@ InsertAdd
@ Epetra_Max
@ Epetra_Min
@ Average
@ Zero
@ Epetra_AddLocalAlso
#define EPETRA_MIN(x, y)
#define EPETRA_MAX(x, y)
#define EPETRA_CHK_ERR(a)
Epetra_DataAccess
@ View
@ Copy
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
int MaxElementSize() const
Maximum element size across all processors.
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
long long * MyGlobalElements64() const
int * ElementSizeList() const
List of the element sizes corresponding to the array MyGlobalElements().
int NumMyPoints() const
Number of local points for this map; equals the sum of all element sizes on the calling processor.
int ElementSize() const
Returns the size of elements in the map; only valid if map has constant element size.
int * FirstPointInElementList() const
Pointer to internal array containing a mapping between the local elements and the first local point n...
bool GlobalIndicesInt() const
Returns true if map create with int NumGlobalElements.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
bool ConstantElementSize() const
Returns true if map has constant element size.
int NumMyElements() const
Number of elements on the calling processor.
int FirstPointInElement(int LID) const
Returns the requested entry in the FirstPointInElementList; see FirstPointInElementList() for details...
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const =0
Epetra_Comm Global Max function.
virtual int NumProc() const =0
Returns total number of processes.
virtual int MyPID() const =0
Return my process ID.
virtual void Barrier() const =0
Epetra_Comm Barrier function.
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
Epetra_Comm Global Sum function.
Epetra_CompObject: Functionality and data that is common to all computational classes.
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.
bool DistributedGlobal() const
Returns true if this multi-vector is distributed global, i.e., not local replicated.
const Epetra_Comm * Comm_
Epetra_BlockMap Map_
Epetra_Distributor: The Epetra Gather/Scatter Setup Base Class.
Epetra_IntMultiVector: A class for constructing and using dense multi-vectors, vectors and matrices i...
int MaxValue(int *Result) const
Compute maximum value of each vector in multi-vector.
int ChangeGlobalValue(int_type GlobalBlockRow, int BlockRowOffset, int VectorIndex, int OrdinalValue, bool SumInto)
int ReplaceMap(const Epetra_BlockMap &map)
Replace map, only if new map has same point-structure as current map.
int MinValue(int *Result) const
Compute minimum value of each vector in multi-vector.
int ReplaceMyValue(int MyRow, int VectorIndex, int OrdinalValue)
Replace current value at the specified (MyRow, VectorIndex) location with OrdinalValue.
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().
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().
int ExtractCopy(int *A, int MyLDA) const
Put multi-vector values into user-provided two-dimensional array.
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 SumIntoGlobalValue(int GlobalRow, int VectorIndex, int OrdinalValue)
Adds OrdinalValue to existing value at the specified (GlobalRow, VectorIndex) location.
int MyLength() const
Returns the local vector length on the calling processor of vectors in the multi-vector.
Epetra_IntMultiVector & operator=(const Epetra_IntMultiVector &Source)
= Operator.
void Assign(const Epetra_IntMultiVector &rhs)
virtual ~Epetra_IntMultiVector()
Epetra_MultiVector destructor.
int ChangeMyValue(int MyBlockRow, int BlockRowOffset, int VectorIndex, int OrdinalValue, bool SumInto)
int PutScalar(int OrdinalConstant)
Initialize all values in a multi-vector with constant value.
int ExtractView(int **A, int *MyLDA) const
Set user-provided addresses of A and MyLDA.
int SumIntoMyValue(int MyRow, int VectorIndex, int OrdinalValue)
Adds OrdinalValue to existing value at the specified (MyRow, VectorIndex) location.
int NumVectors() const
Returns the number of vectors in the multi-vector.
int ResetView(int **ArrayOfPointers)
Reset the view of an existing multivector to point to new user data.
int CheckSizes(const Epetra_SrcDistObject &A)
Allows the source and target (this) objects to be compared for compatibility, return nonzero if not.
Epetra_IntMultiVector(const Epetra_BlockMap &Map, int NumVectors, bool zeroOut=true)
Basic Epetra_IntMultiVector constuctor.
virtual void Print(std::ostream &os) const
Print method.
Epetra_IntVector *& operator()(int i)
Vector access function.
int ReplaceGlobalValue(int GlobalRow, int VectorIndex, int OrdinalValue)
Replace current value at the specified (GlobalRow, VectorIndex) location with OrdinalValue.
Epetra_IntVector ** IntVectors_
int ** Pointers() const
Get pointer to individual vector pointers.
Epetra_IntVector: A class for constructing and using dense integer vectors on a parallel computer.
Epetra_MpiComm: The Epetra MPI Communication Class.
MPI_Comm GetMpiComm() const
Get the MPI Communicator (identical to Comm() method; used when we know we are MPI.
int NumProc() const
Returns total number of processes.
int MyPID() const
Return my process ID.
virtual int ReportError(const std::string Message, int ErrorCode) const
Error reporting method.
std::string toString(const int &x) const
Epetra_OffsetIndex: This class builds index for efficient mapping of data from one Epetra_CrsGraph ba...
Epetra_SrcDistObject: A class for supporting flexible source distributed objects for import/export op...
int SetSeed(unsigned int Seed_in)
Set seed for Random function.