Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Epetra_MultiVector.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_MultiVector.h"
45#include "Epetra_Vector.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_MultiVector::Epetra_MultiVector(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.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 double ** 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 double *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 double **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 double ** 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 MultiVector
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 double ** 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 (Vectors_!=0) {
238 for (int i=0; i<NumVectors_; i++) if (Vectors_[i]!=0) delete Vectors_[i];
239 delete [] Vectors_;
240 }
241
242
243 if (DoubleTemp_!=0) delete [] DoubleTemp_;
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 double[Stride_ * NumVectors_];
258 Pointers_ = new double *[NumVectors_];
259
260 DoubleTemp_ = 0;
261 Vectors_ = 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 double * from = Pointers_[i];
292 double * 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(double));
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 double *[NumVectors_];
313
314 DoubleTemp_ = 0;
315 Vectors_ = 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_MultiVector::ReplaceGlobalValue(int GlobalRow, int VectorIndex, double ScalarValue) {
364
365 // Use the more general method below
366 EPETRA_CHK_ERR(ChangeGlobalValue<int>(GlobalRow, 0, VectorIndex, ScalarValue, false));
367 return(0);
368}
369#endif
370//=========================================================================
371#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
372int Epetra_MultiVector::ReplaceGlobalValue(long long GlobalRow, int VectorIndex, double ScalarValue) {
373
374 // Use the more general method below
375 EPETRA_CHK_ERR(ChangeGlobalValue<long long>(GlobalRow, 0, VectorIndex, ScalarValue, false));
376 return(0);
377}
378#endif
379//=========================================================================
380#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
381int Epetra_MultiVector::ReplaceGlobalValue(int GlobalBlockRow, int BlockRowOffset,
382 int VectorIndex, double ScalarValue) {
383 // Use the more general method below
384 EPETRA_CHK_ERR(ChangeGlobalValue<int>(GlobalBlockRow, BlockRowOffset, VectorIndex, ScalarValue, false));
385 return(0);
386}
387#endif
388//=========================================================================
389#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
390int Epetra_MultiVector::ReplaceGlobalValue(long long GlobalBlockRow, int BlockRowOffset,
391 int VectorIndex, double ScalarValue) {
392 // Use the more general method below
393 EPETRA_CHK_ERR(ChangeGlobalValue<long long>(GlobalBlockRow, BlockRowOffset, VectorIndex, ScalarValue, false));
394 return(0);
395}
396#endif
397//=========================================================================
398#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
399int Epetra_MultiVector::SumIntoGlobalValue(int GlobalRow, int VectorIndex, double ScalarValue) {
400
401 // Use the more general method below
402 EPETRA_CHK_ERR(ChangeGlobalValue<int>(GlobalRow, 0, VectorIndex, ScalarValue, true));
403 return(0);
404}
405#endif
406//=========================================================================
407#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
408int Epetra_MultiVector::SumIntoGlobalValue(long long GlobalRow, int VectorIndex, double ScalarValue) {
409
410 // Use the more general method below
411 EPETRA_CHK_ERR(ChangeGlobalValue<long long>(GlobalRow, 0, VectorIndex, ScalarValue, true));
412 return(0);
413}
414#endif
415//=========================================================================
416#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
417int Epetra_MultiVector::SumIntoGlobalValue(int GlobalBlockRow, int BlockRowOffset,
418 int VectorIndex, double ScalarValue) {
419 // Use the more general method below
420 EPETRA_CHK_ERR(ChangeGlobalValue<int>(GlobalBlockRow, BlockRowOffset, VectorIndex, ScalarValue, true));
421 return(0);
422}
423#endif
424//=========================================================================
425#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
426int Epetra_MultiVector::SumIntoGlobalValue(long long GlobalBlockRow, int BlockRowOffset,
427 int VectorIndex, double ScalarValue) {
428 // Use the more general method below
429 EPETRA_CHK_ERR(ChangeGlobalValue<long long>(GlobalBlockRow, BlockRowOffset, VectorIndex, ScalarValue, true));
430 return(0);
431}
432#endif
433//=========================================================================
434int Epetra_MultiVector::ReplaceMyValue(int MyRow, int VectorIndex, double ScalarValue) {
435
436 // Use the more general method below
437 EPETRA_CHK_ERR(ChangeMyValue(MyRow, 0, VectorIndex, ScalarValue, false));
438 return(0);
439}
440//=========================================================================
441int Epetra_MultiVector::ReplaceMyValue(int MyBlockRow, int BlockRowOffset,
442 int VectorIndex, double ScalarValue) {
443 // Use the more general method below
444 EPETRA_CHK_ERR(ChangeMyValue(MyBlockRow, BlockRowOffset, VectorIndex, ScalarValue, false));
445 return(0);
446}
447//=========================================================================
448int Epetra_MultiVector::SumIntoMyValue(int MyRow, int VectorIndex, double ScalarValue) {
449 // Use the more general method below
450 EPETRA_CHK_ERR(ChangeMyValue(MyRow, 0, VectorIndex, ScalarValue, true));
451 return(0);
452}
453//=========================================================================
454int Epetra_MultiVector::SumIntoMyValue(int MyBlockRow, int BlockRowOffset,
455 int VectorIndex, double ScalarValue) {
456 // Use the more general method below
457 EPETRA_CHK_ERR(ChangeMyValue(MyBlockRow, BlockRowOffset, VectorIndex, ScalarValue, true));
458 return(0);
459}
460//=========================================================================
461template<typename int_type>
462int Epetra_MultiVector::ChangeGlobalValue(int_type GlobalBlockRow, int BlockRowOffset,
463 int VectorIndex, double ScalarValue, bool SumInto) {
464
465 if(!Map().template GlobalIndicesIsType<int_type>())
466 throw ReportError("Epetra_MultiVector::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, ScalarValue, SumInto));
470 return(0);
471}
472//=========================================================================
473int Epetra_MultiVector::ChangeMyValue(int MyBlockRow, int BlockRowOffset,
474 int VectorIndex, double ScalarValue, 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] += ScalarValue;
484 else
485 Pointers_[VectorIndex][entry+BlockRowOffset] = ScalarValue;
486
487 return(0);
488}
489//=========================================================================
491 // Generate random numbers drawn from a uniform distribution on
492 // the interval (-1,1) using a multiplicative congruential generator
493 // with modulus 2^31 - 1.
494 /*
495 const double a = 16807.0, BigInt=2147483647.0, DbleOne=1.0, DbleTwo=2.0;
496
497 for(int i=0; i < NumVectors_; i++)
498 for (int j=0; j<MyLength_; j++){
499 Seed_ = fmod( a*Seed_, BigInt );
500 Pointers_[i][j] = DbleTwo*(Seed_/BigInt)-DbleOne;
501 }
502 */
503
504 const int myLength = MyLength_;
505 for(int i = 0; i < NumVectors_; i++) {
506 double * const to = Pointers_[i];
507#ifdef EPETRA_HAVE_OMP_NONASSOCIATIVE
508#pragma omp parallel for default(none)
509#endif
510 for(int j = 0; j < myLength; j++)
511 to[j] = Util_.RandomDouble();
512 }
513
514 return(0);
515}
516
517//=========================================================================
518
519// Extract a copy of a Epetra_MultiVector. Put in a user's Fortran-style array
520
521int Epetra_MultiVector::ExtractCopy(double *A, int MyLDA) const {
522 if (NumVectors_>1 && Stride_ > MyLDA) EPETRA_CHK_ERR(-1); // LDA not big enough
523
524 const int myLength = MyLength_;
525 for (int i=0; i< NumVectors_; i++)
526 {
527 double * from = Pointers_[i];
528 double * to = A + i*MyLDA;
529 for (int j=0; j<myLength; j++) *to++ = *from++;
530 }
531
532 return(0);
533}
534
535//=========================================================================
537{
538 // mfh 28 Mar 2013: We can't check for compatibility across the
539 // whole communicator, unless we know that the current and new
540 // Maps are nonnull on _all_ participating processes.
542 // So, we'll check to make sure that the maps are the same size on this processor and then
543 // just go with it.
544 if(Map().NumMyElements() == map.NumMyElements() && Map().NumMyPoints() == map.NumMyPoints()) {
546 return(0);
547 }
548
549 return(-1);
550}
551
552//=========================================================================
553
554// Extract a copy of a Epetra_MultiVector. Put in a user's array of pointers
555
556int Epetra_MultiVector::ExtractCopy(double **ArrayOfPointers) const {
557 const int myLength = MyLength_;
558 for (int i=0; i< NumVectors_; i++)
559 {
560 double * from = Pointers_[i];
561 double * to = ArrayOfPointers[i];
562 memcpy(to, from, myLength*sizeof(double));
563 }
564
565 return(0);
566}
567
568
569
570//=========================================================================
571
572// Extract a view of a Epetra_MultiVector. Set up a user's Fortran-style array
573
574int Epetra_MultiVector::ExtractView(double **A, int *MyLDA) const {
575 if (!ConstantStride_) EPETRA_CHK_ERR(-1); // Can't make a Fortran-style view if not constant stride
576 *MyLDA = Stride_; // Set user's LDA
577 *A = Values_; // Set user's value pointer
578 return(0);
579}
580
581
582
583//=========================================================================
584
585// Extract a view of a Epetra_MultiVector. Put in a user's array of pointers
586
587int Epetra_MultiVector::ExtractView(double ***ArrayOfPointers) const {
588 *ArrayOfPointers = Pointers_;
589
590 return(0);
591}
592
593
594//=========================================================================
595int Epetra_MultiVector::PutScalar(double ScalarConstant) {
596
597 // Fills MultiVector with the value ScalarConstant **/
598
599 int myLength = MyLength_;
600 for (int i = 0; i < NumVectors_; i++) {
601 double * to = Pointers_[i];
602#ifdef EPETRA_HAVE_OMP
603#pragma omp parallel for default(none) shared(ScalarConstant,myLength,to)
604#endif
605 for (int j=0; j<myLength; j++) to[j] = ScalarConstant;
606 }
607 return(0);
608}
609//=========================================================================
611
612 const Epetra_MultiVector & A = dynamic_cast<const Epetra_MultiVector &>(Source);
613
614 if (NumVectors()!=A.NumVectors()) {EPETRA_CHK_ERR(-3)};
615 return(0);
616}
617
618//=========================================================================
620 int NumSameIDs,
621 int NumPermuteIDs,
622 int * PermuteToLIDs,
623 int *PermuteFromLIDs,
624 const Epetra_OffsetIndex * Indexor,
625 Epetra_CombineMode CombineMode)
626{
627 (void)Indexor;
628
629 const Epetra_MultiVector & A = dynamic_cast<const Epetra_MultiVector &>(Source);
630
631 double **From = A.Pointers();
632 double **To = Pointers_;
633 int numVectors = NumVectors_;
634
635 int * ToFirstPointInElementList = 0;
636 int * FromFirstPointInElementList = 0;
637 int * FromElementSizeList = 0;
638 int MaxElementSize = Map().MaxElementSize();
639 bool ConstantElementSize = Map().ConstantElementSize();
640
641 if (!ConstantElementSize) {
642 ToFirstPointInElementList = Map().FirstPointInElementList();
643 FromFirstPointInElementList = A.Map().FirstPointInElementList();
644 FromElementSizeList = A.Map().ElementSizeList();
645 }
646 int jj, jjj, k;
647
648 int NumSameEntries;
649
650 bool Case1 = false;
651 bool Case2 = false;
652 // bool Case3 = false;
653
654 if (MaxElementSize==1) {
655 Case1 = true;
656 NumSameEntries = NumSameIDs;
657 }
658 else if (ConstantElementSize) {
659 Case2 = true;
660 NumSameEntries = NumSameIDs * MaxElementSize;
661 }
662 else {
663 // Case3 = true;
664 NumSameEntries = FromFirstPointInElementList[NumSameIDs];
665 }
666
667 // Short circuit for the case where the source and target vector is the same.
668 if (To==From) NumSameEntries = 0;
669
670 // Do copy first
671 if (NumSameIDs>0)
672 if (To!=From) {
673 for (int i=0; i < numVectors; i++) {
674 if (CombineMode==Epetra_AddLocalAlso) for (int j=0; j<NumSameEntries; j++) To[i][j] += From[i][j]; // Add to existing value
675 else for (int j=0; j<NumSameEntries; j++) To[i][j] = From[i][j];
676 }
677 }
678 // Do local permutation next
679 if (NumPermuteIDs>0) {
680
681 // Point entry case
682 if (Case1) {
683
684 if (numVectors==1) {
685 if (CombineMode==Epetra_AddLocalAlso) for (int j=0; j<NumPermuteIDs; j++) To[0][PermuteToLIDs[j]] += From[0][PermuteFromLIDs[j]]; // Add to existing value
686 else for (int j=0; j<NumPermuteIDs; j++) To[0][PermuteToLIDs[j]] = From[0][PermuteFromLIDs[j]];
687 }
688 else {
689 for (int j=0; j<NumPermuteIDs; j++) {
690 jj = PermuteToLIDs[j];
691 jjj = PermuteFromLIDs[j];
692 if (CombineMode==Epetra_AddLocalAlso) for (int i=0; i<numVectors; i++) To[i][jj] += From[i][jjj]; // Add to existing value
693 else for (int i=0; i<numVectors; i++) To[i][jj] = From[i][jjj];
694 }
695 }
696 }
697 // constant element size case
698 else if (Case2) {
699
700 for (int j=0; j<NumPermuteIDs; j++) {
701 jj = MaxElementSize*PermuteToLIDs[j];
702 jjj = MaxElementSize*PermuteFromLIDs[j];
703 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
704 else for(int i=0; i<numVectors; i++) for (k=0; k<MaxElementSize; k++) To[i][jj+k] = From[i][jjj+k];
705 }
706 }
707
708 // variable element size case
709 else {
710
711 for (int j=0; j<NumPermuteIDs; j++) {
712 jj = ToFirstPointInElementList[PermuteToLIDs[j]];
713 jjj = FromFirstPointInElementList[PermuteFromLIDs[j]];
714 int ElementSize = FromElementSizeList[PermuteFromLIDs[j]];
715 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
716 else for (int i=0; i<numVectors; i++) for (k=0; k<ElementSize; k++) To[i][jj+k] = From[i][jjj+k];
717 }
718 }
719 }
720 return(0);
721}
722
723//=========================================================================
725 int NumExportIDs,
726 int * ExportLIDs,
727 int & LenExports,
728 char * & Exports,
729 int & SizeOfPacket,
730 int * Sizes,
731 bool & VarSizes,
732 Epetra_Distributor & Distor)
733{
734 (void)Sizes;
735 (void)VarSizes;
736 (void)Distor;
737
738 const Epetra_MultiVector & A = dynamic_cast<const Epetra_MultiVector &>(Source);
739 int jj, k;
740
741 double **From = A.Pointers();
742 int MaxElementSize = Map().MaxElementSize();
743 int numVectors = NumVectors_;
744 bool ConstantElementSize = Map().ConstantElementSize();
745
746 int * FromFirstPointInElementList = 0;
747 int * FromElementSizeList = 0;
748
749 if (!ConstantElementSize) {
750 FromFirstPointInElementList = A.Map().FirstPointInElementList();
751 FromElementSizeList = A.Map().ElementSizeList();
752 }
753
754 double * DoubleExports = 0;
755
756 SizeOfPacket = numVectors*MaxElementSize*(int)sizeof(double);
757
758 if(SizeOfPacket*NumExportIDs>LenExports) {
759 if (LenExports>0) delete [] Exports;
760 LenExports = SizeOfPacket*NumExportIDs;
761 DoubleExports = new double[numVectors*MaxElementSize*NumExportIDs];
762 Exports = (char *) DoubleExports;
763 }
764
765 double * ptr;
766
767 if (NumExportIDs>0) {
768 ptr = (double *) Exports;
769
770 // Point entry case
771 if (MaxElementSize==1) {
772
773 if (numVectors==1)
774 for (int j=0; j<NumExportIDs; j++)
775 *ptr++ = From[0][ExportLIDs[j]];
776
777 else {
778 for (int j=0; j<NumExportIDs; j++) {
779 jj = ExportLIDs[j];
780 for (int i=0; i<numVectors; i++)
781 *ptr++ = From[i][jj];
782 }
783 }
784 }
785
786 // constant element size case
787 else if (ConstantElementSize) {
788
789 for (int j=0; j<NumExportIDs; j++) {
790 jj = MaxElementSize*ExportLIDs[j];
791 for (int i=0; i<numVectors; i++)
792 for (k=0; k<MaxElementSize; k++)
793 *ptr++ = From[i][jj+k];
794 }
795 }
796
797 // variable element size case
798 else {
799
800 int thisSizeOfPacket = numVectors*MaxElementSize;
801 for (int j=0; j<NumExportIDs; j++) {
802 ptr = (double *) Exports + j*thisSizeOfPacket;
803 jj = FromFirstPointInElementList[ExportLIDs[j]];
804 int ElementSize = FromElementSizeList[ExportLIDs[j]];
805 for (int i=0; i<numVectors; i++)
806 for (k=0; k<ElementSize; k++)
807 *ptr++ = From[i][jj+k];
808 }
809 }
810 }
811
812 return(0);
813}
814
815//=========================================================================
817 int NumImportIDs,
818 int * ImportLIDs,
819 int LenImports,
820 char * Imports,
821 int & SizeOfPacket,
822 Epetra_Distributor & Distor,
823 Epetra_CombineMode CombineMode,
824 const Epetra_OffsetIndex * Indexor )
825{
826 (void)Source;
827 (void)LenImports;
828 (void)SizeOfPacket;
829 (void)Distor;
830 (void)Indexor;
831 int jj, k;
832
833 if( CombineMode != Add
834 && CombineMode != Zero
835 && CombineMode != Insert
836 && CombineMode != InsertAdd
837 && CombineMode != Average
838 && CombineMode != Epetra_Max
839 && CombineMode != Epetra_Min
840 && CombineMode != AbsMin
841 && CombineMode != AbsMax )
842 EPETRA_CHK_ERR(-1); //Unsupported CombinedMode, will default to Zero
843
844 if (NumImportIDs<=0) return(0);
845
846 double ** To = Pointers_;
847 int numVectors = NumVectors_;
848 int MaxElementSize = Map().MaxElementSize();
849 bool ConstantElementSize = Map().ConstantElementSize();
850
851 int * ToFirstPointInElementList = 0;
852 int * ToElementSizeList = 0;
853
854 if (!ConstantElementSize) {
855 ToFirstPointInElementList = Map().FirstPointInElementList();
856 ToElementSizeList = Map().ElementSizeList();
857 }
858
859 double * ptr;
860 // Unpack it...
861
862 ptr = (double *) Imports;
863
864 // Point entry case
865 if (MaxElementSize==1) {
866
867 if (numVectors==1) {
868 if (CombineMode==InsertAdd) for (int j=0; j<NumImportIDs; j++) To[0][ImportLIDs[j]] = 0.0; // Zero out first
869 if (CombineMode==Add || CombineMode==InsertAdd) for (int j=0; j<NumImportIDs; j++) To[0][ImportLIDs[j]] += *ptr++; // Add to existing value
870 else if(CombineMode==Insert) for (int j=0; j<NumImportIDs; j++) To[0][ImportLIDs[j]] = *ptr++;
871 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++;}
872 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++;}
873 else if(CombineMode==Epetra_Max) for (int j=0; j<NumImportIDs; j++) {To[0][ImportLIDs[j]] = EPETRA_MAX( To[0][ImportLIDs[j]],*ptr); ptr++;}
874 else if(CombineMode==Epetra_Min) for (int j=0; j<NumImportIDs; j++) {To[0][ImportLIDs[j]] = EPETRA_MIN( To[0][ImportLIDs[j]],*ptr); ptr++;}
875 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
876 // Note: The following form of averaging is not a true average if more that one value is combined.
877 // This might be an issue in the future, but we leave this way for now.
878/*
879 if (CombineMode==Add)
880 for (int j=0; j<NumImportIDs; j++) To[0][ImportLIDs[j]] += *ptr++; // Add to existing value
881 else if(CombineMode==Insert)
882 for (int j=0; j<NumImportIDs; j++) To[0][ImportLIDs[j]] = *ptr++;
883 else if(CombineMode==InsertAdd) {
884 for (int j=0; j<NumImportIDs; j++) To[0][ImportLIDs[j]] = 0.0;
885 for (int j=0; j<NumImportIDs; j++) To[0][ImportLIDs[j]] += *ptr++;
886 }
887 else if(CombineMode==AbsMax)
888 for (int j=0; j<NumImportIDs; j++) {
889 To[0][ImportLIDs[j]] = EPETRA_MAX( To[0][ImportLIDs[j]],std::abs(*ptr));
890 ptr++;
891 }
892 // Note: The following form of averaging is not a true average if more that one value is combined.
893 // This might be an issue in the future, but we leave this way for now.
894 else if(CombineMode==Average)
895 for (int j=0; j<NumImportIDs; j++) {To[0][ImportLIDs[j]] += *ptr++; To[0][ImportLIDs[j]] *= 0.5;}
896*/
897 }
898
899 else { // numVectors>1
900
901 for (int j=0; j<NumImportIDs; j++) {
902 jj = ImportLIDs[j];
903 for (int i=0; i<numVectors; i++) {
904 if (CombineMode==InsertAdd) To[i][jj] = 0.0; // Zero out if needed
905 if (CombineMode==Add || CombineMode==InsertAdd) To[i][jj] += *ptr++; // Add to existing value
906 else if (CombineMode==Insert) To[i][jj] = *ptr++; // Insert values
907 else if (CombineMode==AbsMax) {To[i][jj] = EPETRA_MAX( To[i][jj], std::abs(*ptr)); ptr++; } // max of absolutes
908 else if (CombineMode==AbsMin) {To[i][jj] = EPETRA_MIN( To[i][jj], std::abs(*ptr)); ptr++; } // max of absolutes
909 else if (CombineMode==Epetra_Max) {To[i][jj] = EPETRA_MAX( To[i][jj], *ptr); ptr++; } // simple max
910 else if (CombineMode==Epetra_Min) {To[i][jj] = EPETRA_MIN( To[i][jj], *ptr); ptr++; } // simple min
911 else if (CombineMode==Average){To[i][jj] += *ptr++; To[i][jj] *= 0.5;}} // Not a true avg if >2 occurance of an ID
912
913 }
914/*
915 if (CombineMode==Add) {
916 for (int j=0; j<NumImportIDs; j++) {
917 jj = ImportLIDs[j];
918 for (int i=0; i<numVectors; i++)
919 To[i][jj] += *ptr++; // Add to existing value
920 }
921 }
922 else if(CombineMode==Insert) {
923 for (int j=0; j<NumImportIDs; j++) {
924 jj = ImportLIDs[j];
925 for (int i=0; i<numVectors; i++)
926 To[i][jj] = *ptr++;
927 }
928 }
929 else if(CombineMode==InsertAdd) {
930 for (int j=0; j<NumImportIDs; j++) {
931 jj = ImportLIDs[j];
932 for (int i=0; i<numVectors; i++)
933 To[i][jj] = 0.0;
934 }
935 for (int j=0; j<NumImportIDs; j++) {
936 jj = ImportLIDs[j];
937 for (int i=0; i<numVectors; i++)
938 To[i][jj] += *ptr++;
939 }
940 }
941 else if(CombineMode==AbsMax) {
942 for (int j=0; j<NumImportIDs; j++) {
943 jj = ImportLIDs[j];
944 for (int i=0; i<numVectors; i++) {
945 To[i][jj] = EPETRA_MAX( To[i][jj], std::abs(*ptr) );
946 ptr++;
947 }
948 }
949 }
950 // Note: The following form of averaging is not a true average if more that one value is combined.
951 // This might be an issue in the future, but we leave this way for now.
952 else if(CombineMode==Average) {
953 for (int j=0; j<NumImportIDs; j++) {
954 jj = ImportLIDs[j];
955 for (int i=0; i<numVectors; i++)
956 { To[i][jj] += *ptr++; To[i][jj] *= 0.5;}
957 }
958 }
959*/
960 }
961 }
962
963 // constant element size case
964
965 else if (ConstantElementSize) {
966
967 for (int j=0; j<NumImportIDs; j++) {
968 jj = MaxElementSize*ImportLIDs[j];
969 for (int i=0; i<numVectors; i++) {
970 if (CombineMode==InsertAdd) for (k=0; k<MaxElementSize; k++) To[i][jj+k] = 0.0; // Zero out if needed
971 if (CombineMode==Add || CombineMode==InsertAdd) for (k=0; k<MaxElementSize; k++) To[i][jj+k] += *ptr++; // Add to existing value
972 else if (CombineMode==Insert) for (k=0; k<MaxElementSize; k++) To[i][jj+k] = *ptr++; // Insert values
973 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
974 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
975 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
976 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
977 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
979 }
980/*
981 if (CombineMode==Add) {
982 for (int j=0; j<NumImportIDs; j++) {
983 jj = MaxElementSize*ImportLIDs[j];
984 for (int i=0; i<numVectors; i++)
985 for (k=0; k<MaxElementSize; k++)
986 To[i][jj+k] += *ptr++; // Add to existing value
987 }
988 }
989 else if(CombineMode==Insert) {
990 for (int j=0; j<NumImportIDs; j++) {
991 jj = MaxElementSize*ImportLIDs[j];
992 for (int i=0; i<numVectors; i++)
993 for (k=0; k<MaxElementSize; k++)
994 To[i][jj+k] = *ptr++;
995 }
996 }
997 else if(CombineMode==InsertAdd) {
998 for (int j=0; j<NumImportIDs; j++) {
999 jj = MaxElementSize*ImportLIDs[j];
1000 for (int i=0; i<numVectors; i++)
1001 for (k=0; k<MaxElementSize; k++)
1002 To[i][jj+k] = 0.0;
1003 }
1004 for (int j=0; j<NumImportIDs; j++) {
1005 jj = MaxElementSize*ImportLIDs[j];
1006 for (int i=0; i<numVectors; i++)
1007 for (k=0; k<MaxElementSize; k++)
1008 To[i][jj+k] += *ptr++;
1009 }
1010 }
1011 else if(CombineMode==AbsMax) {
1012 for (int j=0; j<NumImportIDs; j++) {
1013 jj = MaxElementSize*ImportLIDs[j];
1014 for (int i=0; i<numVectors; i++)
1015 for (k=0; k<MaxElementSize; k++) {
1016 To[i][jj+k] = EPETRA_MAX( To[i][jj+k], std::abs(*ptr) );
1017 ptr++;
1018 }
1019 }
1020 }
1021 // Note: The following form of averaging is not a true average if more that one value is combined.
1022 // This might be an issue in the future, but we leave this way for now.
1023 else if(CombineMode==Average) {
1024 for (int j=0; j<NumImportIDs; j++) {
1025 jj = MaxElementSize*ImportLIDs[j];
1026 for (int i=0; i<numVectors; i++)
1027 for (k=0; k<MaxElementSize; k++)
1028 { To[i][jj+k] += *ptr++; To[i][jj+k] *= 0.5;}
1029 }
1030 }
1031*/
1032 }
1033
1034 // variable element size case
1035
1036 else {
1037 int thisSizeOfPacket = numVectors*MaxElementSize;
1038
1039 for (int j=0; j<NumImportIDs; j++) {
1040 ptr = (double *) Imports + j*thisSizeOfPacket;
1041 jj = ToFirstPointInElementList[ImportLIDs[j]];
1042 int ElementSize = ToElementSizeList[ImportLIDs[j]];
1043 for (int i=0; i<numVectors; i++) {
1044 if (CombineMode==InsertAdd) for (k=0; k<ElementSize; k++) To[i][jj+k] = 0.0; // Zero out if needed
1045 if (CombineMode==Add || CombineMode==InsertAdd) for (k=0; k<ElementSize; k++) To[i][jj+k] += *ptr++; // Add to existing value
1046 else if (CombineMode==Insert) for (k=0; k<ElementSize; k++) To[i][jj+k] = *ptr++; // Insert values
1047 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
1048 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
1049 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
1050 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
1051 }
1052 }
1053/*
1054 if (CombineMode==Add) {
1055 for (int j=0; j<NumImportIDs; j++) {
1056 ptr = (double *) Imports + j*thisSizeOfPacket;
1057 jj = ToFirstPointInElementList[ImportLIDs[j]];
1058 int ElementSize = ToElementSizeList[ImportLIDs[j]];
1059 for (int i=0; i<numVectors; i++)
1060 for (k=0; k<ElementSize; k++)
1061 To[i][jj+k] += *ptr++; // Add to existing value
1062 }
1063 }
1064 else if(CombineMode==Insert){
1065 for (int j=0; j<NumImportIDs; j++) {
1066 ptr = (double *) Imports + j*thisSizeOfPacket;
1067 jj = ToFirstPointInElementList[ImportLIDs[j]];
1068 int ElementSize = ToElementSizeList[ImportLIDs[j]];
1069 for (int i=0; i<numVectors; i++)
1070 for (k=0; k<ElementSize; k++)
1071 To[i][jj+k] = *ptr++;
1072 }
1073 }
1074 else if(CombineMode==InsertAdd){
1075 for (int j=0; j<NumImportIDs; j++) {
1076 ptr = (double *) Imports + j*thisSizeOfPacket;
1077 jj = ToFirstPointInElementList[ImportLIDs[j]];
1078 int ElementSize = ToElementSizeList[ImportLIDs[j]];
1079 for (int i=0; i<numVectors; i++)
1080 for (k=0; k<ElementSize; k++)
1081 To[i][jj+k] = 0.0;
1082 }
1083 for (int j=0; j<NumImportIDs; j++) {
1084 ptr = (double *) Imports + j*thisSizeOfPacket;
1085 jj = ToFirstPointInElementList[ImportLIDs[j]];
1086 int ElementSize = ToElementSizeList[ImportLIDs[j]];
1087 for (int i=0; i<numVectors; i++)
1088 for (k=0; k<ElementSize; k++)
1089 To[i][jj+k] += *ptr++;
1090 }
1091 }
1092 else if(CombineMode==AbsMax){
1093 for (int j=0; j<NumImportIDs; j++) {
1094 ptr = (double *) Imports + j*thisSizeOfPacket;
1095 jj = ToFirstPointInElementList[ImportLIDs[j]];
1096 int ElementSize = ToElementSizeList[ImportLIDs[j]];
1097 for (int i=0; i<numVectors; i++)
1098 for (k=0; k<ElementSize; k++) {
1099 To[i][jj+k] = EPETRA_MAX( To[i][jj+k], std::abs(*ptr));
1100 ptr++;
1101 }
1102 }
1103 }
1104 // Note: The following form of averaging is not a true average if more that one value is combined.
1105 // This might be an issue in the future, but we leave this way for now.
1106 else if(CombineMode==Average) {
1107 for (int j=0; j<NumImportIDs; j++) {
1108 ptr = (double *) Imports + j*thisSizeOfPacket;
1109 jj = ToFirstPointInElementList[ImportLIDs[j]];
1110 int ElementSize = ToElementSizeList[ImportLIDs[j]];
1111 for (int i=0; i<numVectors; i++)
1112 for (k=0; k<ElementSize; k++)
1113 { To[i][jj+k] += *ptr++; To[i][jj+k] *= 0.5;}
1114 }
1115 }
1116*/
1117 }
1118
1119 return(0);
1120}
1121
1122//=========================================================================
1123int Epetra_MultiVector::Dot(const Epetra_MultiVector& A, double *Result) const {
1124
1125 // Dot product of two MultiVectors
1126
1127 if (NumVectors_ != A.NumVectors()) EPETRA_CHK_ERR(-1);
1128 const int myLength = MyLength_;
1129 if (myLength != A.MyLength()) EPETRA_CHK_ERR(-2);
1131
1132 double **A_Pointers = A.Pointers();
1133
1134#ifdef EPETRA_HAVE_OMP_NONASSOCIATIVE
1135 for (int i=0; i < NumVectors_; i++)
1136 {
1137 const double * const from = Pointers_[i];
1138 const double * const fromA = A_Pointers[i];
1139 double sum = 0.0;
1140#pragma omp parallel for reduction (+:sum) default(none)
1141 for (int j=0; j < myLength; j++) sum += from[j] * fromA[j];
1142 DoubleTemp_[i] = sum;
1143 }
1144#else
1145 for (int i=0; i < NumVectors_; i++) DoubleTemp_[i] = DOT(myLength, Pointers_[i], A_Pointers[i]);
1146#endif
1147
1148 if (DistributedGlobal())
1150 else
1151 for (int i=0; i< NumVectors_; ++i) Result[i] = DoubleTemp_[i];
1152
1154
1155 return(0);
1156}
1157//=========================================================================
1159
1160 // this[i][j] = std::abs(A[i][j])
1161
1162 int myLength = MyLength_;
1163 if (NumVectors_ != A.NumVectors()) EPETRA_CHK_ERR(-1);
1164 if (myLength != A.MyLength()) EPETRA_CHK_ERR(-2);
1165
1166 double **A_Pointers = A.Pointers();
1167
1168 for (int i=0; i < NumVectors_; i++) {
1169 double * to = Pointers_[i];
1170 const double * from = A_Pointers[i];
1171#ifdef EPETRA_HAVE_OMP
1172#pragma omp parallel for default(none) shared(myLength,to,from)
1173#endif
1174 for (int j=0; j < myLength; j++) to[j] = std::abs(from[j]);
1175 }
1176
1177 return(0);
1178}
1179//=========================================================================
1181
1182 // this[i][j] = 1.0/(A[i][j])
1183
1184 int ierr = 0;
1185#ifndef EPETRA_HAVE_OMP
1186 int localierr = 0;
1187#endif
1188 int myLength = MyLength_;
1189 if (NumVectors_ != A.NumVectors()) EPETRA_CHK_ERR(-1);
1190 if (myLength != A.MyLength()) EPETRA_CHK_ERR(-2);
1191
1192 double **A_Pointers = A.Pointers();
1193
1194 for (int i=0; i < NumVectors_; i++) {
1195 double * to = Pointers_[i];
1196 const double * from = A_Pointers[i];
1197#ifdef EPETRA_HAVE_OMP
1198#pragma omp parallel default(none) shared(ierr,myLength,to,from)
1199{
1200 int localierr = 0;
1201#pragma omp for
1202#endif
1203 for (int j=0; j < myLength; j++) {
1204 double value = from[j];
1205 // Set error to 1 to signal that zero rowsum found (supercedes ierr = 2)
1206 if (std::abs(value)<Epetra_MinDouble) {
1207 if (value==0.0) localierr = 1;
1208 else if (localierr!=1) localierr = 2;
1209 to[j] = EPETRA_SGN(value) * Epetra_MaxDouble;
1210 }
1211 else
1212 to[j] = 1.0/value;
1213 }
1214#ifdef EPETRA_HAVE_OMP
1215#pragma omp critical
1216#endif
1217{
1218 if (localierr==1) ierr = 1;
1219 else if (localierr==2 && ierr!=1) ierr = 2;
1220}
1221#ifdef EPETRA_HAVE_OMP
1222}
1223#endif
1224 }
1225 EPETRA_CHK_ERR(ierr);
1226 return(0);
1227}
1228 //=========================================================================
1229 int Epetra_MultiVector::Scale (double ScalarValue) {
1230
1231 // scales a MultiVector in place by a scalar
1232
1233
1234 int myLength = MyLength_;
1235#ifdef EPETRA_HAVE_OMP
1236 for (int i = 0; i < NumVectors_; i++) {
1237 double * to = Pointers_[i];
1238#pragma omp parallel for default(none) shared(ScalarValue,myLength,to)
1239 for (int j = 0; j < myLength; j++) to[j] = ScalarValue * to[j];
1240 }
1241#else
1242 for (int i = 0; i < NumVectors_; i++)
1243 SCAL(myLength, ScalarValue, Pointers_[i]);
1244#endif
1246
1247 return(0);
1248 }
1249
1250 //=========================================================================
1251 int Epetra_MultiVector::Scale (double ScalarA, const Epetra_MultiVector& A) {
1252
1253 // scales a MultiVector by a scalar and put in the this:
1254 // this = ScalarA * A
1255
1256 int myLength = MyLength_;
1257 if (NumVectors_ != A.NumVectors()) EPETRA_CHK_ERR(-1);
1258 if (myLength != A.MyLength()) EPETRA_CHK_ERR(-2);
1259
1260 double **A_Pointers = (double**)A.Pointers();
1261
1262 for (int i = 0; i < NumVectors_; i++) {
1263 double * to = Pointers_[i];
1264 const double * from = A_Pointers[i];
1265#ifdef EPETRA_HAVE_OMP
1266#pragma omp parallel for default(none) shared(ScalarA,myLength,to,from)
1267#endif
1268 for (int j = 0; j < myLength; j++) to[j] = ScalarA * from[j];
1269 }
1271
1272 return(0);
1273 }
1274
1275 //=========================================================================
1276 int Epetra_MultiVector::Update(double ScalarA, const Epetra_MultiVector& A, double ScalarThis) {
1277
1278
1279 // linear combination of two MultiVectors: this = ScalarThis * this + ScalarA * A
1280
1281 int myLength = MyLength_;
1282 if (NumVectors_ != A.NumVectors()) EPETRA_CHK_ERR(-1);
1283 if (myLength != A.MyLength()) EPETRA_CHK_ERR(-2);
1284
1285 double **A_Pointers = (double**)A.Pointers();
1286
1287 if (ScalarThis==0.0)
1288 {
1289 for (int i = 0; i < NumVectors_; i++) {
1290 double * to = Pointers_[i];
1291 const double * from = A_Pointers[i];
1292#ifdef EPETRA_HAVE_OMP
1293#pragma omp parallel for default(none) shared(ScalarA,myLength,to,from)
1294#endif
1295 for (int j = 0; j < myLength; j++) to[j] = ScalarA * from[j];
1296 }
1298 }
1299 else if (ScalarThis==1.0)
1300 {
1301 for (int i = 0; i < NumVectors_; i++) {
1302 double * to = Pointers_[i];
1303 const double * from = A_Pointers[i];
1304#ifdef EPETRA_HAVE_OMP
1305#pragma omp parallel for default(none) shared(ScalarA,to,from,myLength)
1306#endif
1307 for (int j = 0; j < myLength; j++) to[j] = to[j] + ScalarA * from[j];
1308 }
1310 }
1311 else if (ScalarA==1.0)
1312 {
1313 for (int i = 0; i < NumVectors_; i++) {
1314 double * to = Pointers_[i];
1315 const double * from = A_Pointers[i];
1316#ifdef EPETRA_HAVE_OMP
1317#pragma omp parallel for default(none) shared(ScalarThis,myLength,to,from)
1318#endif
1319 for (int j = 0; j < myLength; j++) to[j] = ScalarThis * to[j] + from[j];
1320 }
1322 }
1323 else
1324 {
1325 for (int i = 0; i < NumVectors_; i++) {
1326 double * to = Pointers_[i];
1327 const double * from = A_Pointers[i];
1328#ifdef EPETRA_HAVE_OMP
1329#pragma omp parallel for default(none) shared(ScalarA,ScalarThis,to,from,myLength)
1330#endif
1331 for (int j = 0; j < myLength; j++) to[j] = ScalarThis * to[j] +
1332 ScalarA * from[j];
1333 }
1335 }
1336
1337 return(0);
1338 }
1339
1340//=========================================================================
1342 double ScalarB, const Epetra_MultiVector& B, double ScalarThis) {
1343
1344
1345 // linear combination of three MultiVectors:
1346 // this = ScalarThis * this + ScalarA * A + ScalarB * B
1347
1348 if (ScalarA==0.0) {
1349 EPETRA_CHK_ERR(Update(ScalarB, B, ScalarThis));
1350 return(0);
1351 }
1352 if (ScalarB==0.0) {
1353 EPETRA_CHK_ERR(Update(ScalarA, A, ScalarThis));
1354 return(0);
1355 }
1356
1357 int myLength = MyLength_;
1359 if (myLength != A.MyLength() || myLength != B.MyLength()) EPETRA_CHK_ERR(-2);
1360
1361 double **A_Pointers = (double**)A.Pointers();
1362 double **B_Pointers = (double**)B.Pointers();
1363
1364 if (ScalarThis==0.0)
1365 {
1366 if (ScalarA==1.0)
1367 {
1368 for (int i = 0; i < NumVectors_; i++) {
1369 double * to = Pointers_[i];
1370 const double * fromA = A_Pointers[i];
1371 const double * fromB = B_Pointers[i];
1372#ifdef EPETRA_HAVE_OMP
1373#pragma omp parallel for default(none) shared(ScalarB,myLength,to,fromA,fromB)
1374#endif
1375 for (int j = 0; j < myLength; j++) to[j] = fromA[j] +
1376 ScalarB * fromB[j];
1377 }
1379 }
1380 else if (ScalarB==1.0)
1381 {
1382 for (int i = 0; i < NumVectors_; i++) {
1383 double * to = Pointers_[i];
1384 const double * fromA = A_Pointers[i];
1385 const double * fromB = B_Pointers[i];
1386#ifdef EPETRA_HAVE_OMP
1387#pragma omp parallel for default(none) shared(ScalarA,myLength,to,fromA,fromB)
1388#endif
1389 for (int j = 0; j < myLength; j++) to[j] = ScalarA * fromA[j] +
1390 fromB[j];
1391 }
1393 }
1394 else
1395 {
1396 for (int i = 0; i < NumVectors_; i++) {
1397 double * to = Pointers_[i];
1398 const double * fromA = A_Pointers[i];
1399 const double * fromB = B_Pointers[i];
1400#ifdef EPETRA_HAVE_OMP
1401#pragma omp parallel for default(none) shared(ScalarA,ScalarB,myLength,to,fromA,fromB)
1402#endif
1403 for (int j = 0; j < myLength; j++) to[j] = ScalarA * fromA[j] +
1404 ScalarB * fromB[j];
1405 }
1407 }
1408 }
1409 else if (ScalarThis==1.0)
1410 {
1411 if (ScalarA==1.0)
1412 {
1413 for (int i = 0; i < NumVectors_; i++) {
1414 double * to = Pointers_[i];
1415 const double * fromA = A_Pointers[i];
1416 const double * fromB = B_Pointers[i];
1417#ifdef EPETRA_HAVE_OMP
1418#pragma omp parallel for default(none) shared(ScalarB,myLength,to,fromA,fromB)
1419#endif
1420 for (int j = 0; j < myLength; j++) to[j] += fromA[j] +
1421 ScalarB * fromB[j];
1422 }
1424 }
1425 else if (ScalarB==1.0)
1426 {
1427 for (int i = 0; i < NumVectors_; i++) {
1428 double * to = Pointers_[i];
1429 const double * fromA = A_Pointers[i];
1430 const double * fromB = B_Pointers[i];
1431#ifdef EPETRA_HAVE_OMP
1432#pragma omp parallel for default(none) shared(ScalarA,myLength,to,fromA,fromB)
1433#endif
1434 for (int j = 0; j < myLength; j++) to[j] += ScalarA * fromA[j] +
1435 fromB[j];
1436 }
1438 }
1439 else
1440 {
1441 for (int i = 0; i < NumVectors_; i++) {
1442 double * to = Pointers_[i];
1443 const double * fromA = A_Pointers[i];
1444 const double * fromB = B_Pointers[i];
1445#ifdef EPETRA_HAVE_OMP
1446#pragma omp parallel for default(none) shared(ScalarA,ScalarB,myLength,to,fromA,fromB)
1447#endif
1448 for (int j = 0; j < myLength; j++) to[j] += ScalarA * fromA[j] +
1449 ScalarB * fromB[j];
1450 }
1452 }
1453 }
1454 else
1455 {
1456 if (ScalarA==1.0)
1457 {
1458 for (int i = 0; i < NumVectors_; i++) {
1459 double * to = Pointers_[i];
1460 const double * fromA = A_Pointers[i];
1461 const double * fromB = B_Pointers[i];
1462#ifdef EPETRA_HAVE_OMP
1463#pragma omp parallel for default(none) shared(ScalarA,ScalarB,ScalarThis,myLength,to,fromA,fromB)
1464#endif
1465 for (int j = 0; j < myLength; j++) to[j] = ScalarThis * to[j] +
1466 fromA[j] +
1467 ScalarB * fromB[j];
1468 }
1470 }
1471 else if (ScalarB==1.0)
1472 {
1473 for (int i = 0; i < NumVectors_; i++) {
1474 double * to = Pointers_[i];
1475 const double * fromA = A_Pointers[i];
1476 const double * fromB = B_Pointers[i];
1477#ifdef EPETRA_HAVE_OMP
1478#pragma omp parallel for default(none) shared(ScalarA,ScalarThis,myLength,to,fromA,fromB)
1479#endif
1480 for (int j = 0; j < myLength; j++) to[j] = ScalarThis * to[j] +
1481 ScalarA * fromA[j] +
1482 fromB[j];
1483 }
1485 }
1486 else
1487 {
1488 for (int i = 0; i < NumVectors_; i++) {
1489 double * to = Pointers_[i];
1490 const double * fromA = A_Pointers[i];
1491 const double * fromB = B_Pointers[i];
1492#ifdef EPETRA_HAVE_OMP
1493#pragma omp parallel for default(none) shared(ScalarA,ScalarB,ScalarThis,myLength,to,fromA,fromB)
1494#endif
1495 for (int j = 0; j < myLength; j++) to[j] = ScalarThis * to[j] +
1496 ScalarA * fromA[j] +
1497 ScalarB * fromB[j];
1498 }
1500 }
1501 }
1502
1503
1504 return(0);
1505 }
1506
1507//=========================================================================
1508int Epetra_MultiVector::Norm1 (double* Result) const {
1509
1510 // 1-norm of each vector in MultiVector
1511
1512 if (!Map().UniqueGIDs()) {EPETRA_CHK_ERR(-1);}
1513
1514 int myLength = MyLength_;
1516#ifdef EPETRA_HAVE_OMP_NONASSOCIATIVE
1517 for (int i=0; i < NumVectors_; i++)
1518 {
1519 const double * from = Pointers_[i];
1520 double asum = 0.0;
1521#pragma omp parallel default(none) shared(asum,myLength)
1522{
1523 double localasum = 0.0;
1524#pragma omp for
1525 for (int j=0; j< myLength; j++) localasum += std::abs(from[j]);
1526#pragma omp critical
1527 asum += localasum;
1528}
1529 DoubleTemp_[i] = asum;
1530 }
1531#else
1532
1533 for (int i=0; i < NumVectors_; i++) DoubleTemp_[i] = ASUM(myLength, Pointers_[i]);
1534#endif
1535
1536 if (DistributedGlobal())
1538 else
1539 for (int i=0; i< NumVectors_; ++i) Result[i] = DoubleTemp_[i];
1540
1542
1543 return(0);
1544}
1545
1546//=========================================================================
1547int Epetra_MultiVector::Norm2 (double* Result) const {
1548
1549 // 2-norm of each vector in MultiVector
1550
1551
1552 if (!Map().UniqueGIDs()) {EPETRA_CHK_ERR(-1);}
1553
1554 const int myLength = MyLength_;
1556
1557 for (int i=0; i < NumVectors_; i++)
1558 {
1559 const double * const from = Pointers_[i];
1560 double sum = 0.0;
1561#ifdef EPETRA_HAVE_OMP_NONASSOCIATIVE
1562#pragma omp parallel for reduction (+:sum) default(none)
1563#endif
1564 for (int j=0; j < myLength; j++) sum += from[j] * from[j];
1565 DoubleTemp_[i] = sum;
1566 }
1567 if (DistributedGlobal())
1569 else
1570 for (int i=0; i< NumVectors_; ++i) Result[i] = DoubleTemp_[i];
1571
1572 for (int i=0; i < NumVectors_; i++) Result[i] = std::sqrt(Result[i]);
1573
1575
1576 return(0);
1577}
1578
1579//=========================================================================
1580int Epetra_MultiVector::NormInf (double* Result) const {
1581
1582 // Inf-norm of each vector in MultiVector
1583
1584
1585 int myLength = MyLength_;
1587
1588 for (int i=0; i < NumVectors_; i++)
1589 {
1590 DoubleTemp_[i] = 0.0;
1591 double normval = 0.0;
1592 const double * from = Pointers_[i];
1593 if (myLength>0) normval = from[0];
1594#ifdef EPETRA_HAVE_OMP
1595#pragma omp parallel default(none) shared(normval,myLength,from)
1596{
1597 double localnormval = 0.0;
1598#pragma omp for
1599 for (int j=0; j< myLength; j++) {
1600 localnormval = EPETRA_MAX(localnormval,std::abs(from[j]));
1601 }
1602#pragma omp critical
1603 {
1604 normval = EPETRA_MAX(normval,localnormval);
1605 }
1606}
1607 DoubleTemp_[i] = normval;
1608#else
1609 (void) normval; // silence unused value warning in non-OpenMP build
1610 int jj = IAMAX(myLength, Pointers_[i]);
1611 if (jj>-1) DoubleTemp_[i] = std::abs(Pointers_[i][jj]);
1612#endif
1613 }
1614 if (DistributedGlobal())
1616 else
1617 for (int i=0; i< NumVectors_; ++i) Result[i] = DoubleTemp_[i];
1618
1619 // UpdateFlops(0); Strictly speaking there are not FLOPS in this routine
1620 return(0);
1621}
1622
1623//=========================================================================
1624int Epetra_MultiVector::NormWeighted (const Epetra_MultiVector& Weights, double* Result) const {
1625
1626 // Weighted 2-norm of each vector in MultiVector
1627
1628 // If only one vector in Weights, we assume it will be used as the weights for all vectors
1629
1630 int myLength = MyLength_;
1631 bool OneW = false;
1632 if (Weights.NumVectors()==1) OneW = true;
1633 else if (NumVectors_ != Weights.NumVectors()) EPETRA_CHK_ERR(-1);
1634
1635 if (myLength != Weights.MyLength()) EPETRA_CHK_ERR(-2);
1636
1637
1639
1640 double *W = Weights.Values();
1641 double **W_Pointers = Weights.Pointers();
1642
1643 for (int i=0; i < NumVectors_; i++)
1644 {
1645 if (!OneW) W = W_Pointers[i]; // If Weights has the same number of vectors as this, use each weight vector
1646 double sum = 0.0;
1647 const double * from = Pointers_[i];
1648#ifdef EPETRA_HAVE_OMP_NONASSOCIATIVE
1649#pragma omp parallel for reduction (+:sum) default(none) shared(W)
1650#endif
1651 for (int j=0; j < myLength; j++) {
1652 double tmp = from[j]/W[j];
1653 sum += tmp * tmp;
1654 }
1655 DoubleTemp_[i] = sum;
1656 }
1657 double OneOverN;
1658 if (DistributedGlobal()) {
1660 OneOverN = 1.0 / (double) GlobalLength_;
1661 }
1662 else {
1663 for (int i=0; i< NumVectors_; ++i) Result[i] = DoubleTemp_[i];
1664 OneOverN = 1.0 / (double) myLength;
1665 }
1666 for (int i=0; i < NumVectors_; i++) Result[i] = std::sqrt(Result[i]*OneOverN);
1667
1669
1670 return(0);
1671 }
1672
1673//=========================================================================
1674int Epetra_MultiVector::MinValue (double* Result) const {
1675
1676 // Minimum value of each vector in MultiVector
1677
1678 int ierr = 0;
1679
1680 int myLength = MyLength_;
1682
1683 for (int i=0; i < NumVectors_; i++)
1684 {
1685 const double * from = Pointers_[i];
1686 double MinVal = Epetra_MaxDouble;
1687 if (myLength>0) MinVal = from[0];
1688#ifdef EPETRA_HAVE_OMP
1689#pragma omp parallel default(none) shared(MinVal,from,myLength)
1690{
1691 double localMinVal = MinVal;
1692#pragma omp for
1693 for (int j=0; j< myLength; j++) localMinVal = EPETRA_MIN(localMinVal,from[j]);
1694#pragma omp critical
1695 {
1696 MinVal = EPETRA_MIN(MinVal,localMinVal);
1697 }
1698}
1699 DoubleTemp_[i] = MinVal;
1700#else
1701 for (int j=0; j< myLength; j++) MinVal = EPETRA_MIN(MinVal,from[j]);
1702 DoubleTemp_[i] = MinVal;
1703#endif
1704 }
1705
1706 if (myLength > 0) {
1707 for(int i=0; i<NumVectors_; ++i) {
1708 Result[i] = DoubleTemp_[i];
1709 }
1710 }
1711
1712 //If myLength == 0 and Comm_->NumProc() == 1, then Result has
1713 //not been referenced. Also, if vector contents are uninitialized
1714 //then Result contents are not well defined...
1715
1716 if (Comm_->NumProc() == 1 || !DistributedGlobal()) return(ierr);
1717
1718 //We're going to use MPI_Allgather to gather every proc's local-
1719 //min values onto every other proc. We'll use the last position
1720 //of the DoubleTemp_ array to indicate whether this proc has
1721 //valid data that should be considered by other procs when forming
1722 //the global-min results.
1723
1724 if (myLength == 0) DoubleTemp_[NumVectors_] = 0.0;
1725 else DoubleTemp_[NumVectors_] = 1.0;
1726
1727 //Now proceed to handle the parallel case. We'll gather local-min
1728 //values from other procs and form a global-min. If any processor
1729 //has myLength>0, we'll end up with a valid result.
1730
1731#ifdef EPETRA_MPI
1732 const Epetra_MpiComm* epetrampicomm =
1733 dynamic_cast<const Epetra_MpiComm*>(Comm_);
1734 if (!epetrampicomm) {
1735 return(-2);
1736 }
1737
1738 MPI_Comm mpicomm = epetrampicomm->GetMpiComm();
1739 int numProcs = epetrampicomm->NumProc();
1740 double* dwork = new double[numProcs*(NumVectors_+1)];
1741
1742 MPI_Allgather(DoubleTemp_, NumVectors_+1, MPI_DOUBLE,
1743 dwork, NumVectors_+1, MPI_DOUBLE, mpicomm);
1744
1745 //if myLength==0, then our Result array currently contains
1746 //Epetra_MaxDouble from the local-min calculations above. In this
1747 //case we'll overwrite our Result array with values from the first
1748 //processor that sent valid data.
1749 bool overwrite = myLength == 0 ? true : false;
1750
1751 int myPID = epetrampicomm->MyPID();
1752 double* dwork_ptr = dwork;
1753
1754 for(int j=0; j<numProcs; ++j) {
1755
1756 //skip data from self, and skip data from
1757 //procs with DoubleTemp_[NumVectors_] == 0.0.
1758 if (j == myPID || dwork_ptr[NumVectors_] < 0.5) {
1759 dwork_ptr += NumVectors_+1;
1760 continue;
1761 }
1762
1763 for(int i=0; i<NumVectors_; ++i) {
1764 double val = dwork_ptr[i];
1765
1766 //Set val into our Result array if overwrite is true (see above),
1767 //or if val is less than our current Result[i].
1768 if (overwrite || (Result[i] > val)) Result[i] = val;
1769 }
1770
1771 //Now set overwrite to false so that we'll do the right thing
1772 //when processing data from subsequent processors.
1773 if (overwrite) overwrite = false;
1774
1775 dwork_ptr += NumVectors_+1;
1776 }
1777
1778 delete [] dwork;
1779#endif
1780
1781 // UpdateFlops(0); Strictly speaking there are not FLOPS in this routine
1782
1783 return(ierr);
1784}
1785
1786//=========================================================================
1787int Epetra_MultiVector::MaxValue (double* Result) const {
1788
1789 // Maximum value of each vector in MultiVector
1790
1791 int ierr = 0;
1792
1793 int myLength = MyLength_;
1795
1796 for (int i=0; i < NumVectors_; i++)
1797 {
1798 const double * from = Pointers_[i];
1799 double MaxVal = -Epetra_MaxDouble;
1800 if (myLength>0) MaxVal = from[0];
1801#ifdef EPETRA_HAVE_OMP
1802#pragma omp parallel default(none) shared(MaxVal,myLength,from)
1803{
1804 double localMaxVal = MaxVal;
1805#pragma omp for
1806 for (int j=0; j< myLength; j++) localMaxVal = EPETRA_MAX(localMaxVal,from[j]);
1807#pragma omp critical
1808 {
1809 MaxVal = EPETRA_MAX(MaxVal,localMaxVal);
1810 }
1811}
1812 DoubleTemp_[i] = MaxVal;
1813#else
1814 for (int j=0; j< myLength; j++) MaxVal = EPETRA_MAX(MaxVal,from[j]);
1815 DoubleTemp_[i] = MaxVal;
1816#endif
1817 }
1818
1819 if (myLength > 0) {
1820 for(int i=0; i<NumVectors_; ++i) {
1821 Result[i] = DoubleTemp_[i];
1822 }
1823 }
1824
1825 //If myLength == 0 and Comm_->NumProc() == 1, then Result has
1826 //not been referenced. Also, if vector contents are uninitialized
1827 //then Result contents are not well defined...
1828
1829 if (Comm_->NumProc() == 1 || !DistributedGlobal()) return(ierr);
1830
1831 //We're going to use MPI_Allgather to gather every proc's local-
1832 //max values onto every other proc. We'll use the last position
1833 //of the DoubleTemp_ array to indicate whether this proc has
1834 //valid data that should be considered by other procs when forming
1835 //the global-max results.
1836
1837 if (myLength == 0) DoubleTemp_[NumVectors_] = 0.0;
1838 else DoubleTemp_[NumVectors_] = 1.0;
1839
1840 //Now proceed to handle the parallel case. We'll gather local-max
1841 //values from other procs and form a global-max. If any processor
1842 //has myLength>0, we'll end up with a valid result.
1843
1844#ifdef EPETRA_MPI
1845 const Epetra_MpiComm* epetrampicomm =
1846 dynamic_cast<const Epetra_MpiComm*>(Comm_);
1847 if (!epetrampicomm) {
1848 return(-2);
1849 }
1850
1851 MPI_Comm mpicomm = epetrampicomm->GetMpiComm();
1852 int numProcs = epetrampicomm->NumProc();
1853 double* dwork = new double[numProcs*(NumVectors_+1)];
1854
1855 MPI_Allgather(DoubleTemp_, NumVectors_+1, MPI_DOUBLE,
1856 dwork, NumVectors_+1, MPI_DOUBLE, mpicomm);
1857
1858 //if myLength==0, then our Result array currently contains
1859 //-Epetra_MaxDouble from the local-max calculations above. In this
1860 //case we'll overwrite our Result array with values from the first
1861 //processor that sent valid data.
1862 bool overwrite = myLength == 0 ? true : false;
1863
1864 int myPID = epetrampicomm->MyPID();
1865 double* dwork_ptr = dwork;
1866
1867 for(int j=0; j<numProcs; ++j) {
1868
1869 //skip data from self, and skip data from
1870 //procs with DoubleTemp_[NumVectors_] == 0.0.
1871 if (j == myPID || dwork_ptr[NumVectors_] < 0.5) {
1872 dwork_ptr += NumVectors_+1;
1873 continue;
1874 }
1875
1876 for(int i=0; i<NumVectors_; ++i) {
1877 double val = dwork_ptr[i];
1878
1879 //Set val into our Result array if overwrite is true (see above),
1880 //or if val is larger than our current Result[i].
1881 if (overwrite || (Result[i] < val)) Result[i] = val;
1882 }
1883
1884 //Now set overwrite to false so that we'll do the right thing
1885 //when processing data from subsequent processors.
1886 if (overwrite) overwrite = false;
1887
1888 dwork_ptr += NumVectors_+1;
1889 }
1890
1891 delete [] dwork;
1892#endif
1893
1894 // UpdateFlops(0); Strictly speaking there are not FLOPS in this routine
1895
1896 return(ierr);
1897}
1898
1899//=========================================================================
1900int Epetra_MultiVector::MeanValue (double* Result) const {
1901
1902 // Mean value of each vector in MultiVector
1903
1904 const int myLength = MyLength_;
1905
1906 if (!Map().UniqueGIDs()) {EPETRA_CHK_ERR(-1);}
1907
1908 double fGlobalLength = 1.0/EPETRA_MAX((double) GlobalLength_, 1.0);
1909
1910
1912
1913 for (int i=0; i < NumVectors_; i++)
1914 {
1915 double sum = 0.0;
1916 const double * const from = Pointers_[i];
1917#ifdef EPETRA_HAVE_OMP_NONASSOCIATIVE
1918#pragma omp parallel for reduction (+:sum) default(none)
1919#endif
1920 for (int j=0; j < myLength; j++) sum += from[j];
1921 DoubleTemp_[i] = sum;
1922 }
1923 if (DistributedGlobal())
1925 else
1926 for (int i=0; i< NumVectors_; ++i) Result[i] = DoubleTemp_[i];
1927
1928 for (int i=0; i < NumVectors_; i++) Result[i] = Result[i]*fGlobalLength;
1929
1931
1932 return(0);
1933}
1934
1935 //=========================================================================
1936 int Epetra_MultiVector::Multiply (char TransA, char TransB, double ScalarAB,
1937 const Epetra_MultiVector& A,
1938 const Epetra_MultiVector& B,
1939 double ScalarThis ) {
1940
1941 // This routine performs a variety of matrix-matrix multiply operations, interpreting
1942 // the Epetra_MultiVector (this-aka C , A and B) as 2D matrices. Variations are due to
1943 // the fact that A, B and C can be local replicated or global distributed
1944 // Epetra_MultiVectors and that we may or may not operate with the transpose of
1945 // A and B. Possible cases are:
1946
1947 // Num
1948 // OPERATIONS case Notes
1949 // 1) C(local) = A^X(local) * B^X(local) 4 (X=Trans or Not, No comm needed)
1950 // 2) C(local) = A^T(distr) * B (distr) 1 (2D dot product, replicate C)
1951 // 3) C(distr) = A (distr) * B^X(local) 2 (2D vector update, no comm needed)
1952
1953 // Note that the following operations are not meaningful for
1954 // 1D distributions:
1955
1956 // 1) C(local) = A^T(distr) * B^T(distr) 1
1957 // 2) C(local) = A (distr) * B^X(distr) 2
1958 // 3) C(distr) = A^X(local) * B^X(local) 4
1959 // 4) C(distr) = A^X(local) * B^X(distr) 4
1960 // 5) C(distr) = A^T(distr) * B^X(local) 2
1961 // 6) C(local) = A^X(distr) * B^X(local) 4
1962 // 7) C(distr) = A^X(distr) * B^X(local) 4
1963 // 8) C(local) = A^X(local) * B^X(distr) 4
1964
1965 // Total of 32 case (2^5).
1966
1967
1968 //if (!ConstantStride_ ||
1971
1972 // Check for compatible dimensions
1973
1974 int A_nrows = (TransA=='T') ? A.NumVectors() : A.MyLength();
1975 int A_ncols = (TransA=='T') ? A.MyLength() : A.NumVectors();
1976 int B_nrows = (TransB=='T') ? B.NumVectors() : B.MyLength();
1977 int B_ncols = (TransB=='T') ? B.MyLength() : B.NumVectors();
1978
1979 double Scalar_local = ScalarThis; // local copy of Scalar
1980 const int myLength = MyLength_;
1981
1982 if( myLength != A_nrows || // RAB: 2002/01/25: Minor reformat to allow
1983 A_ncols != B_nrows || // setting breakpoint at error return.
1984 NumVectors_ != B_ncols )
1985 EPETRA_CHK_ERR(-2); // Return error
1986
1987 bool A_is_local = (!A.DistributedGlobal());
1988 bool B_is_local = (!B.DistributedGlobal());
1989 bool C_is_local = (!DistributedGlobal());
1990 bool Case1 = ( A_is_local && B_is_local && C_is_local); // Case 1 above
1991 bool Case2 = (!A_is_local && !B_is_local && C_is_local && TransA=='T' );// Case 2
1992 bool Case3 = (!A_is_local && B_is_local && !C_is_local && TransA=='N');// Case 3
1993
1994 if (Case2 && (!A.Map().UniqueGIDs() || !B.Map().UniqueGIDs())) {EPETRA_CHK_ERR(-4);}
1995 if (Case3 && (!A.Map().UniqueGIDs() || !Map().UniqueGIDs())) {EPETRA_CHK_ERR(-5);}
1996
1997 // Test for meaningful cases
1998
1999 if (Case1 || Case2 || Case3)
2000 {
2001 if (ScalarThis!=0.0 && Case2)
2002 {
2003 const int MyPID = Comm_->MyPID();
2004 if (MyPID!=0) Scalar_local = 0.0;
2005 }
2006
2007 // Check if A, B, C have constant stride, if not then make temp copy (strided)
2008
2009 Epetra_MultiVector * A_tmp, * B_tmp, *C_tmp;
2010 if (!ConstantStride_) C_tmp = new Epetra_MultiVector(*this);
2011 else C_tmp = this;
2012
2013 if (!A.ConstantStride()) A_tmp = new Epetra_MultiVector(A);
2014 else A_tmp = (Epetra_MultiVector *) &A;
2015
2016 if (!B.ConstantStride()) B_tmp = new Epetra_MultiVector(B);
2017 else B_tmp = (Epetra_MultiVector *) &B;
2018
2019
2020 int m = myLength;
2021 int n = NumVectors_;
2022 int k = A_ncols;
2023 int lda = EPETRA_MAX(A_tmp->Stride(),1); // The reference BLAS implementation requires lda, ldb, ldc > 0 even if m, n, or k = 0
2024 int ldb = EPETRA_MAX(B_tmp->Stride(),1);
2025 int ldc = EPETRA_MAX(C_tmp->Stride(),1);
2026 double *Ap = A_tmp->Values();
2027 double *Bp = B_tmp->Values();
2028 double *Cp = C_tmp->Values();
2029
2030 GEMM(TransA, TransB, m, n, k, ScalarAB,
2031 Ap, lda, Bp, ldb, Scalar_local, Cp, ldc);
2032
2033 // FLOP Counts
2034 // Num
2035 // OPERATIONS case Notes
2036 // 1) C(local) = A^X(local) * B^X(local) 4 (X=Trans or Not, No comm needed)
2037 // 2) C(local) = A^T(distr) * B (distr) 1 (2D dot product, replicate C)
2038 // 3) C(distr) = A (distr) * B^X(local) 2 (2D vector update, no comm needed)
2039
2040 // For Case 1 we only count the local operations, since we are interested in serial
2041 // cost. Computation on other processors is redundant.
2042 if (Case1)
2043 {
2044 UpdateFlops(2*m*n*k);
2045 if (ScalarAB!=1.0) UpdateFlops(m*n);
2046 if (ScalarThis==1.0) UpdateFlops(m*n); else if (ScalarThis!=0.0) UpdateFlops(2*m*n);
2047 }
2048 else if (Case2)
2049 {
2050 UpdateFlops(2*m*n*A.GlobalLength64());
2051 if (ScalarAB!=1.0) UpdateFlops(m*n);
2052 if (ScalarThis==1.0) UpdateFlops(m*n); else if (ScalarThis!=0.0) UpdateFlops(2*m*n);
2053 }
2054 else
2055 {
2057 if (ScalarAB!=1.0) UpdateFlops(GlobalLength_*n);
2058 if (ScalarThis==1.0) UpdateFlops(GlobalLength_*n);
2059 else if (ScalarThis!=0.0) UpdateFlops(2*GlobalLength_*n);
2060 }
2061
2062 // If A or B were not strided, dispose of extra copies.
2063 if (!A.ConstantStride()) delete A_tmp;
2064 if (!B.ConstantStride()) delete B_tmp;
2065
2066 // If C was not strided, copy from strided version and delete
2067 if (!ConstantStride_)
2068 {
2069 C_tmp->ExtractCopy(Pointers_);
2070 delete C_tmp;
2071 }
2072
2073 // If Case 2 then sum up C and distribute it to all processors.
2074
2075 if (Case2) {EPETRA_CHK_ERR(Reduce());}
2076
2077 return(0);
2078
2079 }
2080 else {EPETRA_CHK_ERR(-3)}; // Return error: not a supported operation
2081
2082 return(0);
2083 }
2084
2085
2086//=========================================================================
2088 double ScalarThis) {
2089
2090
2091 // Hadamard product of two MultiVectors:
2092 // this = ScalarThis * this + ScalarAB * A * B (element-wise)
2093
2094 if (ScalarAB==0.0) {
2095 EPETRA_CHK_ERR(Scale(ScalarThis));
2096 return(0);
2097 }
2098 int myLength = MyLength_;
2099
2100 if (A.NumVectors() != 1 && A.NumVectors() != B.NumVectors()) EPETRA_CHK_ERR(-1); // A must have one column or be the same as B.
2101 if (NumVectors_ != B.NumVectors()) EPETRA_CHK_ERR(-2);
2102 if (myLength != A.MyLength() || myLength != B.MyLength()) EPETRA_CHK_ERR(-3);
2103
2104 double **A_Pointers = (double**)A.Pointers();
2105 double **B_Pointers = (double**)B.Pointers();
2106
2107 int IncA = 1;
2108 if (A.NumVectors() == 1 ) IncA = 0;
2109
2110 if (ScalarThis==0.0) {
2111 if (ScalarAB==1.0)
2112 {
2113 for (int i = 0; i < NumVectors_; i++) {
2114 const double * Aptr = A_Pointers[i*IncA];
2115 const double * Bptr = B_Pointers[i];
2116 double * to = Pointers_[i];
2117#ifdef EPETRA_HAVE_OMP
2118#pragma omp parallel for default(none) shared(myLength,to,Aptr,Bptr)
2119#endif
2120 for (int j = 0; j < myLength; j++) {
2121 to[j] = Aptr[j] * Bptr[j];
2122 }
2123 }
2125 }
2126 else
2127 {
2128 for (int i = 0; i < NumVectors_; i++) {
2129 const double * Aptr = A_Pointers[i*IncA];
2130 const double * Bptr = B_Pointers[i];
2131 double * to = Pointers_[i];
2132#ifdef EPETRA_HAVE_OMP
2133#pragma omp parallel for default(none) shared(ScalarAB,myLength,to,Aptr,Bptr)
2134#endif
2135 for (int j = 0; j < myLength; j++) {
2136 to[j] = ScalarAB * Aptr[j] * Bptr[j];
2137 }
2138 }
2140 }
2141 }
2142 else if (ScalarThis==1.0) {
2143 if (ScalarAB==1.0)
2144 {
2145 for (int i = 0; i < NumVectors_; i++) {
2146 const double * Aptr = A_Pointers[i*IncA];
2147 const double * Bptr = B_Pointers[i];
2148 double * to = Pointers_[i];
2149#ifdef EPETRA_HAVE_OMP
2150#pragma omp parallel for default(none) shared(myLength,to,Aptr,Bptr)
2151#endif
2152 for (int j = 0; j < myLength; j++) {
2153 to[j] += Aptr[j] * Bptr[j];
2154 }
2155 }
2157 }
2158 else {
2159 for (int i = 0; i < NumVectors_; i++) {
2160 const double * Aptr = A_Pointers[i*IncA];
2161 const double * Bptr = B_Pointers[i];
2162 double * to = Pointers_[i];
2163#ifdef EPETRA_HAVE_OMP
2164#pragma omp parallel for default(none) shared(ScalarAB,myLength,to,Aptr,Bptr)
2165#endif
2166 for (int j = 0; j < myLength; j++) {
2167 to[j] += ScalarAB * Aptr[j] * Bptr[j];
2168 }
2169 }
2171 }
2172 }
2173 else { // if (ScalarThis!=1.0 && ScalarThis !=0 )
2174 if (ScalarAB==1.0)
2175 {
2176 for (int i = 0; i < NumVectors_; i++) {
2177 const double * Aptr = A_Pointers[i*IncA];
2178 const double * Bptr = B_Pointers[i];
2179 double * to = Pointers_[i];
2180#ifdef EPETRA_HAVE_OMP
2181#pragma omp parallel for default(none) shared(ScalarThis,myLength,to,Aptr,Bptr)
2182#endif
2183 for (int j = 0; j < myLength; j++) {
2184 to[j] = ScalarThis * to[j] +
2185 Aptr[j] * Bptr[j];
2186 }
2187 }
2189 }
2190 else
2191 {
2192 for (int i = 0; i < NumVectors_; i++) {
2193 const double * Aptr = A_Pointers[i*IncA];
2194 const double * Bptr = B_Pointers[i];
2195 double * to = Pointers_[i];
2196#ifdef EPETRA_HAVE_OMP
2197#pragma omp parallel for default(none) shared(ScalarThis,ScalarAB,myLength,to,Aptr,Bptr)
2198#endif
2199 for (int j = 0; j < myLength; j++) {
2200 to[j] = ScalarThis * to[j] +
2201 ScalarAB * Aptr[j] * Bptr[j];
2202 }
2203 }
2205 }
2206 }
2207 return(0);
2208}
2209//=========================================================================
2211 double ScalarThis) {
2212
2213
2214 // Hadamard product of two MultiVectors:
2215 // this = ScalarThis * this + ScalarAB * B / A (element-wise)
2216
2217 if (ScalarAB==0.0) {
2218 EPETRA_CHK_ERR(Scale(ScalarThis));
2219 return(0);
2220 }
2221 int myLength = MyLength_;
2222
2223 if (A.NumVectors() != 1 && A.NumVectors() != B.NumVectors()) EPETRA_CHK_ERR(-1); // A must have one column or be the same as B.
2224 if (NumVectors_ != B.NumVectors()) EPETRA_CHK_ERR(-2);
2225 if (myLength != A.MyLength() || myLength != B.MyLength()) EPETRA_CHK_ERR(-3);
2226
2227 double **A_Pointers = (double**)A.Pointers();
2228 double **B_Pointers = (double**)B.Pointers();
2229
2230 int IncA = 1;
2231 if (A.NumVectors() == 1 ) IncA = 0;
2232
2233 if (ScalarThis==0.0) {
2234 if (ScalarAB==1.0)
2235 {
2236 for (int i = 0; i < NumVectors_; i++) {
2237 const double * Aptr = A_Pointers[i*IncA];
2238 const double * Bptr = B_Pointers[i];
2239 double * to = Pointers_[i];
2240#ifdef EPETRA_HAVE_OMP
2241#pragma omp parallel for default(none) shared(myLength,to,Aptr,Bptr)
2242#endif
2243 for (int j = 0; j < myLength; j++) {
2244 to[j] = Bptr[j] / Aptr[j];
2245 }
2246 }
2248 }
2249 else
2250 {
2251 for (int i = 0; i < NumVectors_; i++) {
2252 const double * Aptr = A_Pointers[i*IncA];
2253 const double * Bptr = B_Pointers[i];
2254 double * to = Pointers_[i];
2255#ifdef EPETRA_HAVE_OMP
2256#pragma omp parallel for default(none) shared(ScalarAB,myLength,to,Aptr,Bptr)
2257#endif
2258 for (int j = 0; j < myLength; j++) {
2259 to[j] = ScalarAB * Bptr[j] / Aptr[j];
2260 }
2261 }
2263 }
2264 }
2265 else if (ScalarThis==1.0) {
2266 if (ScalarAB==1.0)
2267 {
2268 for (int i = 0; i < NumVectors_; i++) {
2269 const double * Aptr = A_Pointers[i*IncA];
2270 const double * Bptr = B_Pointers[i];
2271 double * to = Pointers_[i];
2272#ifdef EPETRA_HAVE_OMP
2273#pragma omp parallel for default(none) shared(myLength,to,Aptr,Bptr)
2274#endif
2275 for (int j = 0; j < myLength; j++) {
2276 to[j] += Bptr[j] / Aptr[j];
2277 }
2278 }
2280 }
2281 else
2282 {
2283 for (int i = 0; i < NumVectors_; i++) {
2284 const double * Aptr = A_Pointers[i*IncA];
2285 const double * Bptr = B_Pointers[i];
2286 double * to = Pointers_[i];
2287#ifdef EPETRA_HAVE_OMP
2288#pragma omp parallel for default(none) shared(ScalarAB,myLength,to,Aptr,Bptr)
2289#endif
2290 for (int j = 0; j < myLength; j++) {
2291 to[j] += ScalarAB * Bptr[j] / Aptr[j];
2292 }
2293 }
2295 }
2296 }
2297 else { // if (ScalarThis!=1.0 && ScalarThis !=0 )
2298 if (ScalarAB==1.0)
2299 {
2300 for (int i = 0; i < NumVectors_; i++) {
2301 const double * Aptr = A_Pointers[i*IncA];
2302 const double * Bptr = B_Pointers[i];
2303 double * to = Pointers_[i];
2304#ifdef EPETRA_HAVE_OMP
2305#pragma omp parallel for default(none) shared(ScalarThis,myLength,to,Aptr,Bptr)
2306#endif
2307 for (int j = 0; j < myLength; j++) {
2308 to[j] = ScalarThis * to[j] +
2309 Bptr[j] / Aptr[j];
2310 }
2311 }
2313 }
2314 else {
2315 for (int i = 0; i < NumVectors_; i++) {
2316 const double * Aptr = A_Pointers[i*IncA];
2317 const double * Bptr = B_Pointers[i];
2318 double * to = Pointers_[i];
2319#ifdef EPETRA_HAVE_OMP
2320#pragma omp parallel for default(none) shared(ScalarAB,ScalarThis,myLength,to,Aptr,Bptr)
2321#endif
2322 for (int j = 0; j < myLength; j++) {
2323 to[j] = ScalarThis * to[j] + ScalarAB *
2324 Bptr[j] / Aptr[j];
2325 }
2326 }
2328 }
2329 }
2330 return(0);
2331}
2332
2333
2334//=======================================================================
2336
2337 // Epetra_MultiVector::operator () --- return non-const reference
2338
2339 if (index < 0 || index >=NumVectors_)
2340 throw ReportError("Vector index = " + toString(index) + "is out of range. Number of Vectors = " + toString(NumVectors_), -1);
2341
2342 UpdateVectors();
2343
2344 // Create a new Epetra_Vector that is a view of ith vector, if not already present
2345 if (Vectors_[index]==0)
2346 Vectors_[index] = new Epetra_Vector(View, Map(), Pointers_[index]);
2347 return(Vectors_[index]);
2348}
2349
2350//=======================================================================
2352
2353 // Epetra_MultiVector::operator () --- return non-const reference
2354
2355 if (index < 0 || index >=NumVectors_)
2356 throw ReportError("Vector index = " + toString(index) + "is out of range. Number of Vectors = " + toString(NumVectors_), -1);
2357
2358 UpdateVectors();
2359
2360 if (Vectors_[index]==0)
2361 Vectors_[index] = new Epetra_Vector(View, Map(), Pointers_[index]);
2362
2363 const Epetra_Vector * & temp = (const Epetra_Vector * &) (Vectors_[index]);
2364 return(temp);
2365}
2366
2367//========================================================================
2369
2370 // Check for special case of this=Source
2371 if (this != &Source) Assign(Source);
2372
2373 return(*this);
2374}
2375
2376//=========================================================================
2378
2379 int myLength = MyLength_;
2380 if (NumVectors_ != A.NumVectors())
2381 throw ReportError("Number of vectors incompatible in Assign. The this MultiVector has NumVectors = " + toString(NumVectors_)
2382 + ". The A MultiVector has NumVectors = " + toString(A.NumVectors()), -3);
2383 if (myLength != A.MyLength())
2384 throw ReportError("Length of MultiVectors incompatible in Assign. The this MultiVector has MyLength = " + toString(myLength)
2385 + ". The A MultiVector has MyLength = " + toString(A.MyLength()), -4);
2386
2387 double ** A_Pointers = A.Pointers();
2388 for (int i = 0; i< NumVectors_; i++) {
2389 double * to = Pointers_[i];
2390 const double * from = A_Pointers[i];
2391#ifdef EPETRA_HAVE_OMP
2392#pragma omp parallel for default(none) shared(myLength,to,from)
2393#endif
2394 for (int j=0; j<myLength; j++) to[j] = from[j];
2395 }
2396 return;
2397 }
2398
2399 //=========================================================================
2401
2402 // Global reduction on each entry of a Replicated Local MultiVector
2403 const int myLength = MyLength_;
2404 double * source = 0;
2405 if (myLength>0) source = new double[myLength*NumVectors_];
2406 double * target = 0;
2407 bool packed = (ConstantStride_ && (Stride_==myLength));
2408 if (packed) {
2409 for (int i=0; i<myLength*NumVectors_; i++) source[i] = Values_[i];
2410 target = Values_;
2411 }
2412 else {
2413 double * tmp1 = source;
2414 for (int i = 0; i < NumVectors_; i++) {
2415 double * tmp2 = Pointers_[i];
2416 for (int j=0; j< myLength; j++) *tmp1++ = *tmp2++;
2417 }
2418 if (myLength>0) target = new double[myLength*NumVectors_];
2419 }
2420
2421 Comm_->SumAll(source, target, myLength*NumVectors_);
2422 if (myLength>0) delete [] source;
2423 if (!packed) {
2424 double * tmp2 = target;
2425 for (int i = 0; i < NumVectors_; i++) {
2426 double * tmp1 = Pointers_[i];
2427 for (int j=0; j< myLength; j++) *tmp1++ = *tmp2++;
2428 }
2429 if (myLength>0) delete [] target;
2430 }
2431 // UpdateFlops(0); No serial Flops in this function
2432 return(0);
2433 }
2434//=======================================================================
2435int Epetra_MultiVector::ResetView(double ** ArrayOfPointers) {
2436
2437 if (!UserAllocated_) {
2438 EPETRA_CHK_ERR(-1); // Can't reset view if multivector was not allocated as a view
2439 }
2440
2441 for (int i = 0; i< NumVectors_; i++) Pointers_[i] = ArrayOfPointers[i];
2442 DoView();
2443
2444 return(0);
2445 }
2446//=======================================================================
2447void Epetra_MultiVector::Print(std::ostream& os) const {
2448 int MyPID = Map().Comm().MyPID();
2449 int NumProc = Map().Comm().NumProc();
2450
2451 for (int iproc=0; iproc < NumProc; iproc++) {
2452 if (MyPID==iproc) {
2453 int NumVectors1 = NumVectors();
2454 int NumMyElements1 =Map(). NumMyElements();
2455 int MaxElementSize1 = Map().MaxElementSize();
2456 int * FirstPointInElementList1 = NULL;
2457 if (MaxElementSize1!=1) FirstPointInElementList1 = Map().FirstPointInElementList();
2458 double ** A_Pointers = Pointers();
2459
2460 if (MyPID==0) {
2461 os.width(8);
2462 os << " MyPID"; os << " ";
2463 os.width(12);
2464 if (MaxElementSize1==1)
2465 os << "GID ";
2466 else
2467 os << " GID/Point";
2468 for (int j = 0; j < NumVectors1 ; j++)
2469 {
2470 os.width(20);
2471 os << "Value ";
2472 }
2473 os << std::endl;
2474 }
2475 for (int i=0; i < NumMyElements1; i++) {
2476 for (int ii=0; ii< Map().ElementSize(i); ii++) {
2477 int iii;
2478 os.width(10);
2479 os << MyPID; os << " ";
2480 os.width(10);
2481 if (MaxElementSize1==1) {
2482 if(Map().GlobalIndicesInt())
2483 {
2484#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
2485 int * MyGlobalElements1 = Map().MyGlobalElements();
2486 os << MyGlobalElements1[i] << " ";
2487#else
2488 throw ReportError("Epetra_MultiVector::Print: ERROR, GlobalIndicesInt but no API for it.",-1);
2489#endif
2490 }
2491 else if(Map().GlobalIndicesLongLong())
2492 {
2493#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2494 long long * MyGlobalElements1 = Map().MyGlobalElements64();
2495 os << MyGlobalElements1[i] << " ";
2496#else
2497 throw ReportError("Epetra_MultiVector::Print: ERROR, GlobalIndicesLongLong but no API for it.",-1);
2498#endif
2499 }
2500 else
2501 throw ReportError("Epetra_MultiVector::Print ERROR, Don't know map global index type.",-1);
2502
2503 iii = i;
2504 }
2505 else {
2506 if(Map().GlobalIndicesInt())
2507 {
2508#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
2509 int * MyGlobalElements1 = Map().MyGlobalElements();
2510 os << MyGlobalElements1[i]<< "/" << ii << " ";
2511#else
2512 throw ReportError("Epetra_MultiVector::Print: ERROR, GlobalIndicesInt but no API for it.",-1);
2513#endif
2514 }
2515 else if(Map().GlobalIndicesLongLong())
2516 {
2517#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2518 long long * MyGlobalElements1 = Map().MyGlobalElements64();
2519 os << MyGlobalElements1[i]<< "/" << ii << " ";
2520#else
2521 throw ReportError("Epetra_MultiVector::Print: ERROR, GlobalIndicesLongLong but no API for it.",-1);
2522#endif
2523 }
2524 else
2525 throw ReportError("Epetra_MultiVector::Print ERROR, Don't know map global index type.",-1);
2526
2527 iii = FirstPointInElementList1[i]+ii;
2528 }
2529 for (int j = 0; j < NumVectors1 ; j++)
2530 {
2531 os.width(20);
2532 os << A_Pointers[j][iii];
2533 }
2534 os << std::endl;
2535 }
2536 }
2537 os << std::flush;
2538 }
2539
2540 // Do a few global ops to give I/O a chance to complete
2541 Map().Comm().Barrier();
2542 Map().Comm().Barrier();
2543 Map().Comm().Barrier();
2544 }
2545 return;
2546}
Epetra_CombineMode
@ InsertAdd
@ Epetra_Max
@ Epetra_Min
@ Average
@ Epetra_AddLocalAlso
#define EPETRA_SGN(x)
#define EPETRA_MIN(x, y)
#define EPETRA_MAX(x, y)
#define EPETRA_CHK_ERR(a)
const double Epetra_MinDouble
const double Epetra_MaxDouble
Epetra_DataAccess
void SCAL(const int N, const float ALPHA, float *X, const int INCX=1) const
Epetra_BLAS vector scale function (SSCAL)
float DOT(const int N, const float *X, const float *Y, const int INCX=1, const int INCY=1) const
Epetra_BLAS dot product function (SDOT).
float ASUM(const int N, const float *X, const int INCX=1) const
Epetra_BLAS one norm function (SASUM).
int IAMAX(const int N, const float *X, const int INCX=1) const
Epetra_BLAS arg maximum of absolute value function (ISAMAX)
void GEMM(const char TRANSA, const char TRANSB, const int M, const int N, const int K, const float ALPHA, const float *A, const int LDA, const float *B, const int LDB, const float BETA, float *C, const int LDC) const
Epetra_BLAS matrix-matrix multiply function (SGEMM)
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 UniqueGIDs() const
Returns true if map GIDs are 1-to-1.
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.
void UpdateFlops(int Flops_in) const
Increment Flop count for this object.
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_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.
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
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 NormInf(double *Result) const
Compute Inf-norm of each vector in multi-vector.
int Scale(double ScalarValue)
Scale the current values of a multi-vector, this = ScalarValue*this.
int NumVectors() const
Returns the number of vectors in the multi-vector.
int Reciprocal(const Epetra_MultiVector &A)
Puts element-wise reciprocal values of input Multi-vector in target.
Epetra_Vector *& operator()(int i)
Vector access function.
virtual ~Epetra_MultiVector()
Epetra_MultiVector destructor.
int Abs(const Epetra_MultiVector &A)
Puts element-wise absolute values of input Multi-vector in target.
int SumIntoGlobalValue(int GlobalRow, int VectorIndex, double ScalarValue)
Adds ScalarValue to existing value at the specified (GlobalRow, VectorIndex) location.
int MinValue(double *Result) const
Compute minimum value of each vector in multi-vector.
int Dot(const Epetra_MultiVector &A, double *Result) const
Computes dot product of each corresponding pair of vectors.
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().
int MyLength() const
Returns the local vector length on the calling processor of vectors in the multi-vector.
int ExtractView(double **A, int *MyLDA) const
Set user-provided addresses of A and MyLDA.
Epetra_MultiVector & operator=(const Epetra_MultiVector &Source)
= Operator.
int NormWeighted(const Epetra_MultiVector &Weights, double *Result) const
Compute Weighted 2-norm (RMS Norm) of each vector in multi-vector.
int CheckSizes(const Epetra_SrcDistObject &A)
Allows the source and target (this) objects to be compared for compatibility, return nonzero if not.
int ResetView(double **ArrayOfPointers)
Reset the view of an existing multivector to point to new user data.
bool ConstantStride() const
Returns true if this multi-vector has constant stride between vectors.
int ReplaceGlobalValue(int GlobalRow, int VectorIndex, double ScalarValue)
Replace current value at the specified (GlobalRow, VectorIndex) location with ScalarValue.
int ChangeGlobalValue(int_type GlobalBlockRow, int BlockRowOffset, int VectorIndex, double ScalarValue, bool SumInto)
int Stride() const
Returns the stride between vectors in the multi-vector (only meaningful if ConstantStride() is true).
double * Values() const
Get pointer to MultiVector values.
void Assign(const Epetra_MultiVector &rhs)
int Multiply(char TransA, char TransB, double ScalarAB, const Epetra_MultiVector &A, const Epetra_MultiVector &B, double ScalarThis)
Matrix-Matrix multiplication, this = ScalarThis*this + ScalarAB*A*B.
int Norm1(double *Result) const
Compute 1-norm of each vector in multi-vector.
int ChangeMyValue(int MyBlockRow, int BlockRowOffset, int VectorIndex, double ScalarValue, bool SumInto)
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
Update multi-vector values with scaled values of A, this = ScalarThis*this + ScalarA*A.
int ReplaceMyValue(int MyRow, int VectorIndex, double ScalarValue)
Replace current value at the specified (MyRow, VectorIndex) location with ScalarValue.
double ** Pointers() const
Get pointer to individual vector pointers.
int MeanValue(double *Result) const
Compute mean (average) value of each vector in multi-vector.
int SumIntoMyValue(int MyRow, int VectorIndex, double ScalarValue)
Adds ScalarValue to existing value at the specified (MyRow, VectorIndex) location.
int Random()
Set multi-vector values to random numbers.
int Norm2(double *Result) const
Compute 2-norm of each vector in multi-vector.
int MaxValue(double *Result) const
Compute maximum value of each vector in multi-vector.
int ExtractCopy(double *A, int MyLDA) const
Put multi-vector values into user-provided two-dimensional array.
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().
long long GlobalLength64() const
Returns the 64-bit global vector length of vectors in the multi-vector.
int PutScalar(double ScalarConstant)
Initialize all values in a multi-vector with constant value.
int ReciprocalMultiply(double ScalarAB, const Epetra_MultiVector &A, const Epetra_MultiVector &B, double ScalarThis)
Multiply a Epetra_MultiVector by the reciprocal of another, element-by-element.
int ReplaceMap(const Epetra_BlockMap &map)
Replace map, only if new map has same point-structure as current map.
Epetra_Vector ** Vectors_
Epetra_MultiVector(const Epetra_BlockMap &Map, int NumVectors, bool zeroOut=true)
Basic Epetra_MultiVector constuctor.
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.
double RandomDouble()
Returns a random double on the interval (-1.0,1.0)
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.