Epetra Package Browser (Single Doxygen Collection)  Development
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 
1154 
1155  return(0);
1156 }
1157 //=========================================================================
1159 
1160  // this[i][j] = std::abs(A[i][j])
1161 
1162  int myLength = MyLength_;
1163  if (NumVectors_ != A.NumVectors()) EPETRA_CHK_ERR(-1);
1164  if (myLength != A.MyLength()) EPETRA_CHK_ERR(-2);
1165 
1166  double **A_Pointers = A.Pointers();
1167 
1168  for (int i=0; i < NumVectors_; i++) {
1169  double * to = Pointers_[i];
1170  const double * from = A_Pointers[i];
1171 #ifdef EPETRA_HAVE_OMP
1172 #pragma omp parallel for default(none) shared(myLength,to,from)
1173 #endif
1174  for (int j=0; j < myLength; j++) to[j] = std::abs(from[j]);
1175  }
1176 
1177  return(0);
1178 }
1179 //=========================================================================
1181 
1182  // this[i][j] = 1.0/(A[i][j])
1183 
1184  int ierr = 0;
1185 #ifndef EPETRA_HAVE_OMP
1186  int localierr = 0;
1187 #endif
1188  int myLength = MyLength_;
1189  if (NumVectors_ != A.NumVectors()) EPETRA_CHK_ERR(-1);
1190  if (myLength != A.MyLength()) EPETRA_CHK_ERR(-2);
1191 
1192  double **A_Pointers = A.Pointers();
1193 
1194  for (int i=0; i < NumVectors_; i++) {
1195  double * to = Pointers_[i];
1196  const double * from = A_Pointers[i];
1197 #ifdef EPETRA_HAVE_OMP
1198 #pragma omp parallel default(none) shared(ierr,myLength,to,from)
1199 {
1200  int localierr = 0;
1201 #pragma omp for
1202 #endif
1203  for (int j=0; j < myLength; j++) {
1204  double value = from[j];
1205  // Set error to 1 to signal that zero rowsum found (supercedes ierr = 2)
1206  if (std::abs(value)<Epetra_MinDouble) {
1207  if (value==0.0) localierr = 1;
1208  else if (localierr!=1) localierr = 2;
1209  to[j] = EPETRA_SGN(value) * Epetra_MaxDouble;
1210  }
1211  else
1212  to[j] = 1.0/value;
1213  }
1214 #ifdef EPETRA_HAVE_OMP
1215 #pragma omp critical
1216 #endif
1217 {
1218  if (localierr==1) ierr = 1;
1219  else if (localierr==2 && ierr!=1) ierr = 2;
1220 }
1221 #ifdef EPETRA_HAVE_OMP
1222 }
1223 #endif
1224  }
1225  EPETRA_CHK_ERR(ierr);
1226  return(0);
1227 }
1228  //=========================================================================
1229  int Epetra_MultiVector::Scale (double ScalarValue) {
1230 
1231  // scales a MultiVector in place by a scalar
1232 
1233 
1234  int myLength = MyLength_;
1235 #ifdef EPETRA_HAVE_OMP
1236  for (int i = 0; i < NumVectors_; i++) {
1237  double * to = Pointers_[i];
1238 #pragma omp parallel for default(none) shared(ScalarValue,myLength,to)
1239  for (int j = 0; j < myLength; j++) to[j] = ScalarValue * to[j];
1240  }
1241 #else
1242  for (int i = 0; i < NumVectors_; i++)
1243  SCAL(myLength, ScalarValue, Pointers_[i]);
1244 #endif
1246 
1247  return(0);
1248  }
1249 
1250  //=========================================================================
1251  int Epetra_MultiVector::Scale (double ScalarA, const Epetra_MultiVector& A) {
1252 
1253  // scales a MultiVector by a scalar and put in the this:
1254  // this = ScalarA * A
1255 
1256  int myLength = MyLength_;
1257  if (NumVectors_ != A.NumVectors()) EPETRA_CHK_ERR(-1);
1258  if (myLength != A.MyLength()) EPETRA_CHK_ERR(-2);
1259 
1260  double **A_Pointers = (double**)A.Pointers();
1261 
1262  for (int i = 0; i < NumVectors_; i++) {
1263  double * to = Pointers_[i];
1264  const double * from = A_Pointers[i];
1265 #ifdef EPETRA_HAVE_OMP
1266 #pragma omp parallel for default(none) shared(ScalarA,myLength,to,from)
1267 #endif
1268  for (int j = 0; j < myLength; j++) to[j] = ScalarA * from[j];
1269  }
1271 
1272  return(0);
1273  }
1274 
1275  //=========================================================================
1276  int Epetra_MultiVector::Update(double ScalarA, const Epetra_MultiVector& A, double ScalarThis) {
1277 
1278 
1279  // linear combination of two MultiVectors: this = ScalarThis * this + ScalarA * A
1280 
1281  int myLength = MyLength_;
1282  if (NumVectors_ != A.NumVectors()) EPETRA_CHK_ERR(-1);
1283  if (myLength != A.MyLength()) EPETRA_CHK_ERR(-2);
1284 
1285  double **A_Pointers = (double**)A.Pointers();
1286 
1287  if (ScalarThis==0.0)
1288  {
1289  for (int i = 0; i < NumVectors_; i++) {
1290  double * to = Pointers_[i];
1291  const double * from = A_Pointers[i];
1292 #ifdef EPETRA_HAVE_OMP
1293 #pragma omp parallel for default(none) shared(ScalarA,myLength,to,from)
1294 #endif
1295  for (int j = 0; j < myLength; j++) to[j] = ScalarA * from[j];
1296  }
1298  }
1299  else if (ScalarThis==1.0)
1300  {
1301  for (int i = 0; i < NumVectors_; i++) {
1302  double * to = Pointers_[i];
1303  const double * from = A_Pointers[i];
1304 #ifdef EPETRA_HAVE_OMP
1305 #pragma omp parallel for default(none) shared(ScalarA,to,from,myLength)
1306 #endif
1307  for (int j = 0; j < myLength; j++) to[j] = to[j] + ScalarA * from[j];
1308  }
1310  }
1311  else if (ScalarA==1.0)
1312  {
1313  for (int i = 0; i < NumVectors_; i++) {
1314  double * to = Pointers_[i];
1315  const double * from = A_Pointers[i];
1316 #ifdef EPETRA_HAVE_OMP
1317 #pragma omp parallel for default(none) shared(ScalarThis,myLength,to,from)
1318 #endif
1319  for (int j = 0; j < myLength; j++) to[j] = ScalarThis * to[j] + from[j];
1320  }
1322  }
1323  else
1324  {
1325  for (int i = 0; i < NumVectors_; i++) {
1326  double * to = Pointers_[i];
1327  const double * from = A_Pointers[i];
1328 #ifdef EPETRA_HAVE_OMP
1329 #pragma omp parallel for default(none) shared(ScalarA,ScalarThis,to,from,myLength)
1330 #endif
1331  for (int j = 0; j < myLength; j++) to[j] = ScalarThis * to[j] +
1332  ScalarA * from[j];
1333  }
1335  }
1336 
1337  return(0);
1338  }
1339 
1340 //=========================================================================
1341 int Epetra_MultiVector::Update(double ScalarA, const Epetra_MultiVector& A,
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  }
1379  }
1380  else if (ScalarB==1.0)
1381  {
1382  for (int i = 0; i < NumVectors_; i++) {
1383  double * to = Pointers_[i];
1384  const double * fromA = A_Pointers[i];
1385  const double * fromB = B_Pointers[i];
1386 #ifdef EPETRA_HAVE_OMP
1387 #pragma omp parallel for default(none) shared(ScalarA,myLength,to,fromA,fromB)
1388 #endif
1389  for (int j = 0; j < myLength; j++) to[j] = ScalarA * fromA[j] +
1390  fromB[j];
1391  }
1393  }
1394  else
1395  {
1396  for (int i = 0; i < NumVectors_; i++) {
1397  double * to = Pointers_[i];
1398  const double * fromA = A_Pointers[i];
1399  const double * fromB = B_Pointers[i];
1400 #ifdef EPETRA_HAVE_OMP
1401 #pragma omp parallel for default(none) shared(ScalarA,ScalarB,myLength,to,fromA,fromB)
1402 #endif
1403  for (int j = 0; j < myLength; j++) to[j] = ScalarA * fromA[j] +
1404  ScalarB * fromB[j];
1405  }
1407  }
1408  }
1409  else if (ScalarThis==1.0)
1410  {
1411  if (ScalarA==1.0)
1412  {
1413  for (int i = 0; i < NumVectors_; i++) {
1414  double * to = Pointers_[i];
1415  const double * fromA = A_Pointers[i];
1416  const double * fromB = B_Pointers[i];
1417 #ifdef EPETRA_HAVE_OMP
1418 #pragma omp parallel for default(none) shared(ScalarB,myLength,to,fromA,fromB)
1419 #endif
1420  for (int j = 0; j < myLength; j++) to[j] += fromA[j] +
1421  ScalarB * fromB[j];
1422  }
1424  }
1425  else if (ScalarB==1.0)
1426  {
1427  for (int i = 0; i < NumVectors_; i++) {
1428  double * to = Pointers_[i];
1429  const double * fromA = A_Pointers[i];
1430  const double * fromB = B_Pointers[i];
1431 #ifdef EPETRA_HAVE_OMP
1432 #pragma omp parallel for default(none) shared(ScalarA,myLength,to,fromA,fromB)
1433 #endif
1434  for (int j = 0; j < myLength; j++) to[j] += ScalarA * fromA[j] +
1435  fromB[j];
1436  }
1438  }
1439  else
1440  {
1441  for (int i = 0; i < NumVectors_; i++) {
1442  double * to = Pointers_[i];
1443  const double * fromA = A_Pointers[i];
1444  const double * fromB = B_Pointers[i];
1445 #ifdef EPETRA_HAVE_OMP
1446 #pragma omp parallel for default(none) shared(ScalarA,ScalarB,myLength,to,fromA,fromB)
1447 #endif
1448  for (int j = 0; j < myLength; j++) to[j] += ScalarA * fromA[j] +
1449  ScalarB * fromB[j];
1450  }
1452  }
1453  }
1454  else
1455  {
1456  if (ScalarA==1.0)
1457  {
1458  for (int i = 0; i < NumVectors_; i++) {
1459  double * to = Pointers_[i];
1460  const double * fromA = A_Pointers[i];
1461  const double * fromB = B_Pointers[i];
1462 #ifdef EPETRA_HAVE_OMP
1463 #pragma omp parallel for default(none) shared(ScalarA,ScalarB,ScalarThis,myLength,to,fromA,fromB)
1464 #endif
1465  for (int j = 0; j < myLength; j++) to[j] = ScalarThis * to[j] +
1466  fromA[j] +
1467  ScalarB * fromB[j];
1468  }
1470  }
1471  else if (ScalarB==1.0)
1472  {
1473  for (int i = 0; i < NumVectors_; i++) {
1474  double * to = Pointers_[i];
1475  const double * fromA = A_Pointers[i];
1476  const double * fromB = B_Pointers[i];
1477 #ifdef EPETRA_HAVE_OMP
1478 #pragma omp parallel for default(none) shared(ScalarA,ScalarThis,myLength,to,fromA,fromB)
1479 #endif
1480  for (int j = 0; j < myLength; j++) to[j] = ScalarThis * to[j] +
1481  ScalarA * fromA[j] +
1482  fromB[j];
1483  }
1485  }
1486  else
1487  {
1488  for (int i = 0; i < NumVectors_; i++) {
1489  double * to = Pointers_[i];
1490  const double * fromA = A_Pointers[i];
1491  const double * fromB = B_Pointers[i];
1492 #ifdef EPETRA_HAVE_OMP
1493 #pragma omp parallel for default(none) shared(ScalarA,ScalarB,ScalarThis,myLength,to,fromA,fromB)
1494 #endif
1495  for (int j = 0; j < myLength; j++) to[j] = ScalarThis * to[j] +
1496  ScalarA * fromA[j] +
1497  ScalarB * fromB[j];
1498  }
1500  }
1501  }
1502 
1503 
1504  return(0);
1505  }
1506 
1507 //=========================================================================
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 
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 
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 
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 
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  }
2125  }
2126  else
2127  {
2128  for (int i = 0; i < NumVectors_; i++) {
2129  const double * Aptr = A_Pointers[i*IncA];
2130  const double * Bptr = B_Pointers[i];
2131  double * to = Pointers_[i];
2132 #ifdef EPETRA_HAVE_OMP
2133 #pragma omp parallel for default(none) shared(ScalarAB,myLength,to,Aptr,Bptr)
2134 #endif
2135  for (int j = 0; j < myLength; j++) {
2136  to[j] = ScalarAB * Aptr[j] * Bptr[j];
2137  }
2138  }
2140  }
2141  }
2142  else if (ScalarThis==1.0) {
2143  if (ScalarAB==1.0)
2144  {
2145  for (int i = 0; i < NumVectors_; i++) {
2146  const double * Aptr = A_Pointers[i*IncA];
2147  const double * Bptr = B_Pointers[i];
2148  double * to = Pointers_[i];
2149 #ifdef EPETRA_HAVE_OMP
2150 #pragma omp parallel for default(none) shared(myLength,to,Aptr,Bptr)
2151 #endif
2152  for (int j = 0; j < myLength; j++) {
2153  to[j] += Aptr[j] * Bptr[j];
2154  }
2155  }
2157  }
2158  else {
2159  for (int i = 0; i < NumVectors_; i++) {
2160  const double * Aptr = A_Pointers[i*IncA];
2161  const double * Bptr = B_Pointers[i];
2162  double * to = Pointers_[i];
2163 #ifdef EPETRA_HAVE_OMP
2164 #pragma omp parallel for default(none) shared(ScalarAB,myLength,to,Aptr,Bptr)
2165 #endif
2166  for (int j = 0; j < myLength; j++) {
2167  to[j] += ScalarAB * Aptr[j] * Bptr[j];
2168  }
2169  }
2171  }
2172  }
2173  else { // if (ScalarThis!=1.0 && ScalarThis !=0 )
2174  if (ScalarAB==1.0)
2175  {
2176  for (int i = 0; i < NumVectors_; i++) {
2177  const double * Aptr = A_Pointers[i*IncA];
2178  const double * Bptr = B_Pointers[i];
2179  double * to = Pointers_[i];
2180 #ifdef EPETRA_HAVE_OMP
2181 #pragma omp parallel for default(none) shared(ScalarThis,myLength,to,Aptr,Bptr)
2182 #endif
2183  for (int j = 0; j < myLength; j++) {
2184  to[j] = ScalarThis * to[j] +
2185  Aptr[j] * Bptr[j];
2186  }
2187  }
2189  }
2190  else
2191  {
2192  for (int i = 0; i < NumVectors_; i++) {
2193  const double * Aptr = A_Pointers[i*IncA];
2194  const double * Bptr = B_Pointers[i];
2195  double * to = Pointers_[i];
2196 #ifdef EPETRA_HAVE_OMP
2197 #pragma omp parallel for default(none) shared(ScalarThis,ScalarAB,myLength,to,Aptr,Bptr)
2198 #endif
2199  for (int j = 0; j < myLength; j++) {
2200  to[j] = ScalarThis * to[j] +
2201  ScalarAB * Aptr[j] * Bptr[j];
2202  }
2203  }
2205  }
2206  }
2207  return(0);
2208 }
2209 //=========================================================================
2211  double ScalarThis) {
2212 
2213 
2214  // Hadamard product of two MultiVectors:
2215  // this = ScalarThis * this + ScalarAB * B / A (element-wise)
2216 
2217  if (ScalarAB==0.0) {
2218  EPETRA_CHK_ERR(Scale(ScalarThis));
2219  return(0);
2220  }
2221  int myLength = MyLength_;
2222 
2223  if (A.NumVectors() != 1 && A.NumVectors() != B.NumVectors()) EPETRA_CHK_ERR(-1); // A must have one column or be the same as B.
2224  if (NumVectors_ != B.NumVectors()) EPETRA_CHK_ERR(-2);
2225  if (myLength != A.MyLength() || myLength != B.MyLength()) EPETRA_CHK_ERR(-3);
2226 
2227  double **A_Pointers = (double**)A.Pointers();
2228  double **B_Pointers = (double**)B.Pointers();
2229 
2230  int IncA = 1;
2231  if (A.NumVectors() == 1 ) IncA = 0;
2232 
2233  if (ScalarThis==0.0) {
2234  if (ScalarAB==1.0)
2235  {
2236  for (int i = 0; i < NumVectors_; i++) {
2237  const double * Aptr = A_Pointers[i*IncA];
2238  const double * Bptr = B_Pointers[i];
2239  double * to = Pointers_[i];
2240 #ifdef EPETRA_HAVE_OMP
2241 #pragma omp parallel for default(none) shared(myLength,to,Aptr,Bptr)
2242 #endif
2243  for (int j = 0; j < myLength; j++) {
2244  to[j] = Bptr[j] / Aptr[j];
2245  }
2246  }
2248  }
2249  else
2250  {
2251  for (int i = 0; i < NumVectors_; i++) {
2252  const double * Aptr = A_Pointers[i*IncA];
2253  const double * Bptr = B_Pointers[i];
2254  double * to = Pointers_[i];
2255 #ifdef EPETRA_HAVE_OMP
2256 #pragma omp parallel for default(none) shared(ScalarAB,myLength,to,Aptr,Bptr)
2257 #endif
2258  for (int j = 0; j < myLength; j++) {
2259  to[j] = ScalarAB * Bptr[j] / Aptr[j];
2260  }
2261  }
2263  }
2264  }
2265  else if (ScalarThis==1.0) {
2266  if (ScalarAB==1.0)
2267  {
2268  for (int i = 0; i < NumVectors_; i++) {
2269  const double * Aptr = A_Pointers[i*IncA];
2270  const double * Bptr = B_Pointers[i];
2271  double * to = Pointers_[i];
2272 #ifdef EPETRA_HAVE_OMP
2273 #pragma omp parallel for default(none) shared(myLength,to,Aptr,Bptr)
2274 #endif
2275  for (int j = 0; j < myLength; j++) {
2276  to[j] += Bptr[j] / Aptr[j];
2277  }
2278  }
2280  }
2281  else
2282  {
2283  for (int i = 0; i < NumVectors_; i++) {
2284  const double * Aptr = A_Pointers[i*IncA];
2285  const double * Bptr = B_Pointers[i];
2286  double * to = Pointers_[i];
2287 #ifdef EPETRA_HAVE_OMP
2288 #pragma omp parallel for default(none) shared(ScalarAB,myLength,to,Aptr,Bptr)
2289 #endif
2290  for (int j = 0; j < myLength; j++) {
2291  to[j] += ScalarAB * Bptr[j] / Aptr[j];
2292  }
2293  }
2295  }
2296  }
2297  else { // if (ScalarThis!=1.0 && ScalarThis !=0 )
2298  if (ScalarAB==1.0)
2299  {
2300  for (int i = 0; i < NumVectors_; i++) {
2301  const double * Aptr = A_Pointers[i*IncA];
2302  const double * Bptr = B_Pointers[i];
2303  double * to = Pointers_[i];
2304 #ifdef EPETRA_HAVE_OMP
2305 #pragma omp parallel for default(none) shared(ScalarThis,myLength,to,Aptr,Bptr)
2306 #endif
2307  for (int j = 0; j < myLength; j++) {
2308  to[j] = ScalarThis * to[j] +
2309  Bptr[j] / Aptr[j];
2310  }
2311  }
2313  }
2314  else {
2315  for (int i = 0; i < NumVectors_; i++) {
2316  const double * Aptr = A_Pointers[i*IncA];
2317  const double * Bptr = B_Pointers[i];
2318  double * to = Pointers_[i];
2319 #ifdef EPETRA_HAVE_OMP
2320 #pragma omp parallel for default(none) shared(ScalarAB,ScalarThis,myLength,to,Aptr,Bptr)
2321 #endif
2322  for (int j = 0; j < myLength; j++) {
2323  to[j] = ScalarThis * to[j] + ScalarAB *
2324  Bptr[j] / Aptr[j];
2325  }
2326  }
2328  }
2329  }
2330  return(0);
2331 }
2332 
2333 
2334 //=======================================================================
2336 
2337  // Epetra_MultiVector::operator () --- return non-const reference
2338 
2339  if (index < 0 || index >=NumVectors_)
2340  throw ReportError("Vector index = " + toString(index) + "is out of range. Number of Vectors = " + toString(NumVectors_), -1);
2341 
2342  UpdateVectors();
2343 
2344  // Create a new Epetra_Vector that is a view of ith vector, if not already present
2345  if (Vectors_[index]==0)
2346  Vectors_[index] = new Epetra_Vector(View, Map(), Pointers_[index]);
2347  return(Vectors_[index]);
2348 }
2349 
2350 //=======================================================================
2352 
2353  // Epetra_MultiVector::operator () --- return non-const reference
2354 
2355  if (index < 0 || index >=NumVectors_)
2356  throw ReportError("Vector index = " + toString(index) + "is out of range. Number of Vectors = " + toString(NumVectors_), -1);
2357 
2358  UpdateVectors();
2359 
2360  if (Vectors_[index]==0)
2361  Vectors_[index] = new Epetra_Vector(View, Map(), Pointers_[index]);
2362 
2363  const Epetra_Vector * & temp = (const Epetra_Vector * &) (Vectors_[index]);
2364  return(temp);
2365 }
2366 
2367 //========================================================================
2369 
2370  // Check for special case of this=Source
2371  if (this != &Source) Assign(Source);
2372 
2373  return(*this);
2374 }
2375 
2376 //=========================================================================
2378 
2379  int myLength = MyLength_;
2380  if (NumVectors_ != A.NumVectors())
2381  throw ReportError("Number of vectors incompatible in Assign. The this MultiVector has NumVectors = " + toString(NumVectors_)
2382  + ". The A MultiVector has NumVectors = " + toString(A.NumVectors()), -3);
2383  if (myLength != A.MyLength())
2384  throw ReportError("Length of MultiVectors incompatible in Assign. The this MultiVector has MyLength = " + toString(myLength)
2385  + ". The A MultiVector has MyLength = " + toString(A.MyLength()), -4);
2386 
2387  double ** A_Pointers = A.Pointers();
2388  for (int i = 0; i< NumVectors_; i++) {
2389  double * to = Pointers_[i];
2390  const double * from = A_Pointers[i];
2391 #ifdef EPETRA_HAVE_OMP
2392 #pragma omp parallel for default(none) shared(myLength,to,from)
2393 #endif
2394  for (int j=0; j<myLength; j++) to[j] = from[j];
2395  }
2396  return;
2397  }
2398 
2399  //=========================================================================
2401 
2402  // Global reduction on each entry of a Replicated Local MultiVector
2403  const int myLength = MyLength_;
2404  double * source = 0;
2405  if (myLength>0) source = new double[myLength*NumVectors_];
2406  double * target = 0;
2407  bool packed = (ConstantStride_ && (Stride_==myLength));
2408  if (packed) {
2409  for (int i=0; i<myLength*NumVectors_; i++) source[i] = Values_[i];
2410  target = Values_;
2411  }
2412  else {
2413  double * tmp1 = source;
2414  for (int i = 0; i < NumVectors_; i++) {
2415  double * tmp2 = Pointers_[i];
2416  for (int j=0; j< myLength; j++) *tmp1++ = *tmp2++;
2417  }
2418  if (myLength>0) target = new double[myLength*NumVectors_];
2419  }
2420 
2421  Comm_->SumAll(source, target, myLength*NumVectors_);
2422  if (myLength>0) delete [] source;
2423  if (!packed) {
2424  double * tmp2 = target;
2425  for (int i = 0; i < NumVectors_; i++) {
2426  double * tmp1 = Pointers_[i];
2427  for (int j=0; j< myLength; j++) *tmp1++ = *tmp2++;
2428  }
2429  if (myLength>0) delete [] target;
2430  }
2431  // UpdateFlops(0); No serial Flops in this function
2432  return(0);
2433  }
2434 //=======================================================================
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)
int Norm2(double *Result) const
Compute 2-norm of each vector in multi-vector.
int Norm1(double *Result) const
Compute 1-norm of each vector in multi-vector.
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
double * Values() const
Get pointer to MultiVector values.
int NumProc() const
Returns total number of processes.
const Epetra_BlockMap & Map() const
Returns the address of the Epetra_BlockMap for this multi-vector.
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
int Abs(const Epetra_MultiVector &A)
Puts element-wise absolute values of input Multi-vector in target.
int Dot(const Epetra_MultiVector &A, double *Result) const
Computes dot product of each corresponding pair of vectors.
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.
void UpdateFlops(int Flops_in) const
Increment Flop count for this object.
bool UniqueGIDs() const
Returns true if map GIDs are 1-to-1.
bool ConstantStride() const
Returns true if this multi-vector has constant stride between vectors.
int CopyAndPermute(const Epetra_SrcDistObject &Source, int NumSameIDs, int NumPermuteIDs, int *PermuteToLIDs, int *PermuteFromLIDs, const Epetra_OffsetIndex *Indexor, Epetra_CombineMode CombineMode=Zero)
Perform ID copies and permutations that are on processor.
int MyLength() const
Returns the local vector length on the calling processor of vectors in the multi-vector.
Epetra_Distributor: The Epetra Gather/Scatter Setup Base Class.
int MeanValue(double *Result) const
Compute mean (average) value of each vector in multi-vector.
virtual void Print(std::ostream &os) const
Print method.
bool ConstantElementSize() const
Returns true if map has constant element size.
int NormWeighted(const Epetra_MultiVector &Weights, double *Result) const
Compute Weighted 2-norm (RMS Norm) of each vector in multi-vector.
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
Epetra_OffsetIndex: This class builds index for efficient mapping of data from one Epetra_CrsGraph ba...
int IAMAX(const int N, const float *X, const int INCX=1) const
Epetra_BLAS arg maximum of absolute value function (ISAMAX)
float ASUM(const int N, const float *X, const int INCX=1) const
Epetra_BLAS one norm function (SASUM).
Definition: Epetra_BLAS.cpp:65
int PutScalar(double ScalarConstant)
Initialize all values in a multi-vector with constant value.
int UnpackAndCombine(const Epetra_SrcDistObject &Source, int NumImportIDs, int *ImportLIDs, int LenImports, char *Imports, int &SizeOfPacket, Epetra_Distributor &Distor, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor)
Perform any unpacking and combining after call to DoTransfer().
#define EPETRA_CHK_ERR(a)
int NormInf(double *Result) const
Compute Inf-norm of each vector in multi-vector.
int Stride() const
Returns the stride between vectors in the multi-vector (only meaningful if ConstantStride() is true)...
void UpdateDoubleTemp() const
int MaxValue(double *Result) const
Compute maximum value of each vector in multi-vector.
int * FirstPointInElementList() const
Pointer to internal array containing a mapping between the local elements and the first local point n...
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
int ExtractCopy(double *A, int MyLDA) const
Put multi-vector values into user-provided two-dimensional array.
#define EPETRA_MIN(x, y)
virtual int ReportError(const std::string Message, int ErrorCode) const
Error reporting method.
Epetra_MultiVector(const Epetra_BlockMap &Map, int NumVectors, bool zeroOut=true)
Basic Epetra_MultiVector constuctor.
Epetra_MpiComm: The Epetra MPI Communication Class.
int FirstPointInElement(int LID) const
Returns the requested entry in the FirstPointInElementList; see FirstPointInElementList() for details...
virtual void Barrier() const =0
Epetra_Comm Barrier function.
int MinValue(double *Result) const
Compute minimum value of each vector in multi-vector.
long long * MyGlobalElements64() const
virtual int MyPID() const =0
Return my process ID.
int * ElementSizeList() const
List of the element sizes corresponding to the array MyGlobalElements().
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const =0
Epetra_Comm Global Max function.
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.
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.
const Epetra_Comm * Comm_
#define EPETRA_SGN(x)
int ResetView(double **ArrayOfPointers)
Reset the view of an existing multivector to point to new user data.
int NumVectors() const
Returns the number of vectors in the multi-vector.
int NumMyElements() const
Number of elements on the calling processor.
const double Epetra_MaxDouble
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.
MPI_Comm GetMpiComm() const
Get the MPI Communicator (identical to Comm() method; used when we know we are MPI.
int CheckSizes(const Epetra_SrcDistObject &A)
Allows the source and target (this) objects to be compared for compatibility, return nonzero if not...
int ElementSize() const
Returns the size of elements in the map; only valid if map has constant element size.
std::string toString(const int &x) const
int SumIntoMyValue(int MyRow, int VectorIndex, double ScalarValue)
Adds ScalarValue to existing value at the specified (MyRow, VectorIndex) location.
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
long long GlobalLength64() const
Returns the 64-bit global vector length of vectors in the multi-vector.
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 MyPID() const
Return my process ID.
const double Epetra_MinDouble
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 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.
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)
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.
int ExtractView(double **A, int *MyLDA) const
Set user-provided addresses of A and MyLDA.
int ReplaceMap(const Epetra_BlockMap &map)
Replace map, only if new map has same point-structure as current map.
Epetra_SrcDistObject: A class for supporting flexible source distributed objects for import/export op...
int NumMyPoints() const
Number of local points for this map; equals the sum of all element sizes on the calling processor...
Epetra_DataAccess
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 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().
double ** Pointers() const
Get pointer to individual vector pointers.
#define EPETRA_MAX(x, y)
bool DistributedGlobal() const
Returns true if this multi-vector is distributed global, i.e., not local replicated.
Epetra_DistObject: A class for constructing and using dense multi-vectors, vectors and matrices in pa...
Epetra_Vector *& operator()(int i)
Vector access function.
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.