Epetra Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
60 Epetra_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
363 int 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
372 int 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
381 int 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
390 int 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
399 int 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
408 int 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
417 int 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
426 int 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 //=========================================================================
434 int 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 //=========================================================================
441 int 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 //=========================================================================
448 int 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 //=========================================================================
454 int 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 //=========================================================================
461 template<typename int_type>
462 int 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 //=========================================================================
473 int 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 
521 int 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.
541 
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 
556 int 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 
574 int 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 
587 int Epetra_MultiVector::ExtractView(double ***ArrayOfPointers) const {
588  *ArrayOfPointers = Pointers_;
589 
590  return(0);
591 }
592 
593 
594 //=========================================================================
595 int 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
978  }
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 //=========================================================================
1123 int 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);
1130  UpdateDoubleTemp();
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())
1149  Comm_->SumAll(DoubleTemp_, Result, NumVectors_);
1150  else
1151  for (int i=0; i< NumVectors_; ++i) Result[i] = DoubleTemp_[i];
1152 
1153  UpdateFlops(2*GlobalLength_*NumVectors_);
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
1245  UpdateFlops(GlobalLength_*NumVectors_);
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  }
1270  UpdateFlops(GlobalLength_*NumVectors_);
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  }
1297  UpdateFlops(GlobalLength_*NumVectors_);
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  }
1309  UpdateFlops(2*GlobalLength_*NumVectors_);
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  }
1321  UpdateFlops(2*GlobalLength_*NumVectors_);
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  }
1334  UpdateFlops(3*GlobalLength_*NumVectors_);
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_;
1358  if (NumVectors_ != A.NumVectors() || NumVectors_ != B.NumVectors()) EPETRA_CHK_ERR(-1);
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  }
1378  UpdateFlops(2*GlobalLength_*NumVectors_);
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  }
1392  UpdateFlops(2*GlobalLength_*NumVectors_);
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  }
1406  UpdateFlops(3*GlobalLength_*NumVectors_);
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  }
1423  UpdateFlops(3*GlobalLength_*NumVectors_);
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  }
1437  UpdateFlops(3*GlobalLength_*NumVectors_);
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  }
1451  UpdateFlops(4*GlobalLength_*NumVectors_);
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  }
1469  UpdateFlops(4*GlobalLength_*NumVectors_);
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  }
1484  UpdateFlops(4*GlobalLength_*NumVectors_);
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  }
1499  UpdateFlops(5*GlobalLength_*NumVectors_);
1500  }
1501  }
1502 
1503 
1504  return(0);
1505  }
1506 
1507 //=========================================================================
1508 int 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_;
1515  UpdateDoubleTemp();
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())
1537  Comm_->SumAll(DoubleTemp_, Result, NumVectors_);
1538  else
1539  for (int i=0; i< NumVectors_; ++i) Result[i] = DoubleTemp_[i];
1540 
1541  UpdateFlops(2*GlobalLength_*NumVectors_);
1542 
1543  return(0);
1544 }
1545 
1546 //=========================================================================
1547 int 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_;
1555  UpdateDoubleTemp();
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())
1568  Comm_->SumAll(DoubleTemp_, Result, NumVectors_);
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 
1574  UpdateFlops(2*GlobalLength_*NumVectors_);
1575 
1576  return(0);
1577 }
1578 
1579 //=========================================================================
1580 int Epetra_MultiVector::NormInf (double* Result) const {
1581 
1582  // Inf-norm of each vector in MultiVector
1583 
1584 
1585  int myLength = MyLength_;
1586  UpdateDoubleTemp();
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())
1615  Comm_->MaxAll(DoubleTemp_, Result, NumVectors_);
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 //=========================================================================
1624 int 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 
1638  UpdateDoubleTemp();
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()) {
1659  Comm_->SumAll(DoubleTemp_, Result, NumVectors_);
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 
1668  UpdateFlops(3*GlobalLength_*NumVectors_);
1669 
1670  return(0);
1671  }
1672 
1673 //=========================================================================
1674 int 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_;
1681  UpdateDoubleTemp();
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 //=========================================================================
1787 int 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_;
1794  UpdateDoubleTemp();
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 //=========================================================================
1900 int 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 
1911  UpdateDoubleTemp();
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())
1924  Comm_->SumAll(DoubleTemp_, Result, NumVectors_);
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 
1930  UpdateFlops(GlobalLength_*NumVectors_);
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  {
2056  UpdateFlops(2*GlobalLength_*n*k);
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  }
2124  UpdateFlops(GlobalLength_*NumVectors_);
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  }
2139  UpdateFlops(2*GlobalLength_*NumVectors_);
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  }
2156  UpdateFlops(2*GlobalLength_*NumVectors_);
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  }
2170  UpdateFlops(3*GlobalLength_*NumVectors_);
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  }
2188  UpdateFlops(3*GlobalLength_*NumVectors_);
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  }
2204  UpdateFlops(4*GlobalLength_*NumVectors_);
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  }
2247  UpdateFlops(GlobalLength_*NumVectors_);
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  }
2262  UpdateFlops(2*GlobalLength_*NumVectors_);
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  }
2279  UpdateFlops(2*GlobalLength_*NumVectors_);
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  }
2294  UpdateFlops(3*GlobalLength_*NumVectors_);
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  }
2312  UpdateFlops(3*GlobalLength_*NumVectors_);
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  }
2327  UpdateFlops(4*GlobalLength_*NumVectors_);
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 //=======================================================================
2435 int 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 //=======================================================================
2447 void 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 }
int ChangeMyValue(int MyBlockRow, int BlockRowOffset, int VectorIndex, double ScalarValue, bool SumInto)
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
int * FirstPointInElementList() const
Pointer to internal array containing a mapping between the local elements and the first local point n...
bool DistributedGlobal() const
Returns true if this multi-vector is distributed global, i.e., not local replicated.
int Abs(const Epetra_MultiVector &A)
Puts element-wise absolute values of input Multi-vector in target.
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)
int Random()
Set multi-vector values to random numbers.
int ReplaceMyValue(int MyRow, int VectorIndex, double ScalarValue)
Replace current value at the specified (MyRow, VectorIndex) location with ScalarValue.
int ElementSize() const
Returns the size of elements in the map; only valid if map has constant element size.
MPI_Comm GetMpiComm() const
Get the MPI Communicator (identical to Comm() method; used when we know we are MPI.
int MaxValue(double *Result) const
Compute maximum value of each vector in multi-vector.
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.
Epetra_Distributor: The Epetra Gather/Scatter Setup Base Class.
Epetra_OffsetIndex: This class builds index for efficient mapping of data from one Epetra_CrsGraph ba...
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
int NormWeighted(const Epetra_MultiVector &Weights, double *Result) const
Compute Weighted 2-norm (RMS Norm) of each vector in multi-vector.
bool UniqueGIDs() const
Returns true if map GIDs are 1-to-1.
int PutScalar(double ScalarConstant)
Initialize all values in a multi-vector with constant value.
int MyLength() const
Returns the local vector length on the calling processor of vectors in the multi-vector.
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 NumProc() const
Returns total number of processes.
bool ConstantElementSize() const
Returns true if map has constant element size.
virtual void Print(std::ostream &os) const
Print method.
double ** Pointers() const
Get pointer to individual vector pointers.
#define EPETRA_CHK_ERR(a)
int * ElementSizeList() const
List of the element sizes corresponding to the array MyGlobalElements().
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
#define EPETRA_MIN(x, y)
Epetra_MultiVector(const Epetra_BlockMap &Map, int NumVectors, bool zeroOut=true)
Basic Epetra_MultiVector constuctor.
Epetra_MpiComm: The Epetra MPI Communication Class.
virtual void Barrier() const =0
Epetra_Comm Barrier function.
int ExtractCopy(double *A, int MyLDA) const
Put multi-vector values into user-provided two-dimensional array.
int NumVectors() const
Returns the number of vectors in the multi-vector.
virtual int MyPID() const =0
Return my process ID.
int Dot(const Epetra_MultiVector &A, double *Result) const
Computes dot product of each corresponding pair of vectors.
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const =0
Epetra_Comm Global Max function.
void UpdateDoubleTemp() const
Epetra_Vector ** Vectors_
int SumIntoGlobalValue(int GlobalRow, int VectorIndex, double ScalarValue)
Adds ScalarValue to existing value at the specified (GlobalRow, VectorIndex) location.
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
Epetra_Comm Global Sum function.
int NumMyElements() const
Number of elements on the calling processor.
Epetra_MultiVector & operator=(const Epetra_MultiVector &Source)
= Operator.
int Scale(double ScalarValue)
Scale the current values of a multi-vector, this = ScalarValue*this.
virtual ~Epetra_MultiVector()
Epetra_MultiVector destructor.
int SetSeed(unsigned int Seed_in)
Set seed for Random function.
std::string toString(const int &x) const
long long GlobalLength64() const
Returns the 64-bit global vector length of vectors in the multi-vector.
int FirstPointInElement(int LID) const
Returns the requested entry in the FirstPointInElementList; see FirstPointInElementList() for details...
float ASUM(const int N, const float *X, const int INCX=1) const
Epetra_BLAS one norm function (SASUM).
Definition: Epetra_BLAS.cpp:65
const Epetra_Comm * Comm_
int NormInf(double *Result) const
Compute Inf-norm of each vector in multi-vector.
#define EPETRA_SGN(x)
int ResetView(double **ArrayOfPointers)
Reset the view of an existing multivector to point to new user data.
const double Epetra_MaxDouble
void UpdateVectors() const
int ReplaceGlobalValue(int GlobalRow, int VectorIndex, double ScalarValue)
Replace current value at the specified (GlobalRow, VectorIndex) location with ScalarValue.
Epetra_CompObject: Functionality and data that is common to all computational classes.
int CheckSizes(const Epetra_SrcDistObject &A)
Allows the source and target (this) objects to be compared for compatibility, return nonzero if not...
int SumIntoMyValue(int MyRow, int VectorIndex, double ScalarValue)
Adds ScalarValue to existing value at the specified (MyRow, VectorIndex) location.
int ExtractView(double **A, int *MyLDA) const
Set user-provided addresses of A and MyLDA.
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
int MeanValue(double *Result) const
Compute mean (average) value of each vector in multi-vector.
void SCAL(const int N, const float ALPHA, float *X, const int INCX=1) const
Epetra_BLAS vector scale function (SSCAL)
Definition: Epetra_BLAS.cpp:89
int MinValue(double *Result) const
Compute minimum value of each vector in multi-vector.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
int Stride() const
Returns the stride between vectors in the multi-vector (only meaningful if ConstantStride() is true)...
bool ConstantStride() const
Returns true if this multi-vector has constant stride between vectors.
const double Epetra_MinDouble
int MyPID() const
Return my process ID.
double * Values() const
Get pointer to MultiVector values.
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 Norm1(double *Result) const
Compute 1-norm of each vector in multi-vector.
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.
virtual int NumProc() const =0
Returns total number of processes.
int IAMAX(const int N, const float *X, const int INCX=1) const
Epetra_BLAS arg maximum of absolute value function (ISAMAX)
double RandomDouble()
Returns a random double on the interval (-1.0,1.0)
Definition: Epetra_Util.cpp:90
Epetra_BlockMap Map_
int ChangeGlobalValue(int_type GlobalBlockRow, int BlockRowOffset, int VectorIndex, double ScalarValue, bool SumInto)
int Norm2(double *Result) const
Compute 2-norm of each vector in multi-vector.
void Assign(const Epetra_MultiVector &rhs)
int MaxElementSize() const
Maximum element size across all processors.
Epetra_CombineMode
int Reciprocal(const Epetra_MultiVector &A)
Puts element-wise reciprocal values of input Multi-vector in target.
virtual int ReportError(const std::string Message, int ErrorCode) const
Error reporting method.
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).
Definition: Epetra_BLAS.cpp:73
int ReplaceMap(const Epetra_BlockMap &map)
Replace map, only if new map has same point-structure as current map.
int NumProc() const
Returns total number of processes.
int NumMyPoints() const
Number of local points for this map; equals the sum of all element sizes on the calling processor...
Epetra_SrcDistObject: A class for supporting flexible source distributed objects for import/export op...
long long * MyGlobalElements64() const
Epetra_DataAccess
int 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 n
#define EPETRA_MAX(x, y)
void UpdateFlops(int Flops_in) const
Increment Flop count for this object.
int MyPID() const
Return my process ID.
Epetra_DistObject: A class for constructing and using dense multi-vectors, vectors and matrices in pa...
Epetra_Vector *& operator()(int i)
Vector access function.
const Epetra_BlockMap & Map() const
Returns the address of the Epetra_BlockMap for this multi-vector.
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.