Epetra Package Browser (Single Doxygen Collection) Development
Epetra_MultiVector.cpp
Go to the documentation of this file.
00001 
00002 //@HEADER
00003 // ************************************************************************
00004 //
00005 //               Epetra: Linear Algebra Services Package
00006 //                 Copyright 2011 Sandia Corporation
00007 //
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00039 //
00040 // ************************************************************************
00041 //@HEADER
00042 
00043 #include "Epetra_ConfigDefs.h"
00044 #include "Epetra_MultiVector.h"
00045 #include "Epetra_Vector.h"
00046 #include "Epetra_Comm.h"
00047 #ifdef EPETRA_MPI
00048 #include "Epetra_MpiComm.h"
00049 #endif
00050 #include "Epetra_BlockMap.h"
00051 #include "Epetra_Map.h"
00052 #include "Epetra_Import.h"
00053 #include "Epetra_Export.h"
00054 #include "Epetra_Distributor.h"
00055 
00056 //=============================================================================
00057 
00058 // Epetra_BlockMap Constructor
00059 
00060 Epetra_MultiVector::Epetra_MultiVector(const Epetra_BlockMap& map, int numVectors, bool zeroOut)
00061   : Epetra_DistObject(map, "Epetra::MultiVector"),
00062     Epetra_CompObject(),
00063     Values_(0),
00064     Pointers_(0),
00065     MyLength_(map.NumMyPoints()),
00066     GlobalLength_(map.NumGlobalPoints64()),
00067     NumVectors_(numVectors),
00068     UserAllocated_(false),
00069     ConstantStride_(true),
00070     Stride_(map.NumMyPoints()),
00071     Allocated_(false)
00072 {
00073   Util_.SetSeed(1);
00074 
00075     AllocateForCopy();
00076 
00077     for (int i = 0; i< NumVectors_; i++) Pointers_[i] = Values_+i*Stride_;
00078 
00079   if(zeroOut) PutScalar(0.0); // Fill all vectors with zero.
00080 }
00081 //==========================================================================
00082 
00083 // Copy Constructor
00084 
00085 Epetra_MultiVector::Epetra_MultiVector(const Epetra_MultiVector& Source)
00086   : Epetra_DistObject(Source),
00087     Epetra_CompObject(Source),
00088     Values_(0),
00089     Pointers_(0),
00090     MyLength_(Source.MyLength_),
00091     GlobalLength_(Source.GlobalLength_),
00092     NumVectors_(Source.NumVectors_),
00093     UserAllocated_(false),
00094     ConstantStride_(true),
00095     Stride_(0),
00096     Allocated_(false),
00097     Util_(Source.Util_)
00098 {
00099   AllocateForCopy();
00100 
00101   double ** Source_Pointers = Source.Pointers();
00102   for (int i = 0; i< NumVectors_; i++) Pointers_[i] = Source_Pointers[i];
00103 
00104   DoCopy();
00105 
00106 }
00107 //==========================================================================
00108 
00109 // This constructor copies in or makes view of a standard Fortran array
00110 
00111 Epetra_MultiVector::Epetra_MultiVector(Epetra_DataAccess CV, const Epetra_BlockMap& map,
00112                double *A, int MyLDA, int numVectors)
00113   : Epetra_DistObject(map, "Epetra::MultiVector"),
00114     Epetra_CompObject(),
00115     Values_(0),
00116     Pointers_(0),
00117     MyLength_(map.NumMyPoints()),
00118     GlobalLength_(map.NumGlobalPoints64()),
00119     NumVectors_(numVectors),
00120     UserAllocated_(false),
00121     ConstantStride_(true),
00122     Stride_(map.NumMyPoints()),
00123     Allocated_(false)
00124 {
00125   Util_.SetSeed(1);
00126 
00127   if (CV==Copy) AllocateForCopy();
00128   else AllocateForView();
00129 
00130   for (int i = 0; i< NumVectors_; i++) Pointers_[i] = A + i*MyLDA;
00131 
00132    if (CV==Copy) DoCopy();
00133    else DoView();
00134 
00135 }
00136 
00137 //==========================================================================
00138 
00139 // This constructor copies in or makes view of a C/C++ array of pointer
00140 
00141 Epetra_MultiVector::Epetra_MultiVector(Epetra_DataAccess CV, const Epetra_BlockMap& map,
00142                 double **ArrayOfPointers, int numVectors)
00143   : Epetra_DistObject(map, "Epetra::MultiVector"),
00144     Epetra_CompObject(),
00145     Values_(0),
00146     Pointers_(0),
00147     MyLength_(map.NumMyPoints()),
00148     GlobalLength_(map.NumGlobalPoints64()),
00149     NumVectors_(numVectors),
00150     UserAllocated_(false),
00151     ConstantStride_(true),
00152     Stride_(map.NumMyPoints()),
00153     Allocated_(false)
00154 {
00155   Util_.SetSeed(1);
00156 
00157   if (CV==Copy) AllocateForCopy();
00158   else AllocateForView();
00159 
00160   for (int i = 0; i< NumVectors_; i++) Pointers_[i] = ArrayOfPointers[i];
00161 
00162    if (CV==Copy) DoCopy();
00163    else DoView();
00164 
00165 }
00166 
00167 //==========================================================================
00168 
00169 // This constructor copies or makes view of selected vectors, specified in Indices,
00170 // from an existing MultiVector
00171 
00172 Epetra_MultiVector::Epetra_MultiVector(Epetra_DataAccess CV, const Epetra_MultiVector& Source,
00173                int *Indices, int numVectors)
00174   : Epetra_DistObject(Source.Map(), "Epetra::MultiVector"),
00175     Epetra_CompObject(),
00176     Values_(0),
00177     Pointers_(0),
00178     MyLength_(Source.MyLength_),
00179     GlobalLength_(Source.GlobalLength_),
00180     NumVectors_(numVectors),
00181     UserAllocated_(false),
00182     ConstantStride_(true),
00183     Stride_(0),
00184     Allocated_(false)
00185 {
00186   Util_.SetSeed(1);
00187 
00188   if (CV==Copy) AllocateForCopy();
00189   else AllocateForView();
00190 
00191   double ** Source_Pointers = Source.Pointers();
00192   for (int i = 0; i< NumVectors_; i++) Pointers_[i] = Source_Pointers[Indices[i]];
00193 
00194    if (CV==Copy) DoCopy();
00195    else DoView();
00196 
00197 }
00198 
00199 //==========================================================================
00200 
00201 // This interface copies or makes view of a range of vectors from an existing MultiVector
00202 
00203 Epetra_MultiVector::Epetra_MultiVector(Epetra_DataAccess CV, const Epetra_MultiVector& Source,
00204                int StartIndex, int numVectors)
00205   : Epetra_DistObject(Source.Map(), "Epetra::MultiVector"),
00206     Epetra_CompObject(),
00207     Values_(0),
00208     Pointers_(0),
00209     MyLength_(Source.MyLength_),
00210     GlobalLength_(Source.GlobalLength_),
00211     NumVectors_(numVectors),
00212     UserAllocated_(false),
00213     ConstantStride_(true),
00214     Stride_(0),
00215     Allocated_(false)
00216 {
00217   Util_.SetSeed(1);
00218 
00219   if (CV==Copy) AllocateForCopy();
00220   else AllocateForView();
00221 
00222   double ** Source_Pointers = Source.Pointers();
00223   for (int i = 0; i< NumVectors_; i++) Pointers_[i] = Source_Pointers[StartIndex+i];
00224 
00225    if (CV==Copy) DoCopy();
00226    else DoView();
00227 }
00228 
00229 //=========================================================================
00230 Epetra_MultiVector::~Epetra_MultiVector(){
00231 
00232   if (!Allocated_) return;
00233 
00234   delete [] Pointers_;
00235   if (!UserAllocated_ && Values_!=0) delete [] Values_;
00236 
00237   if (Vectors_!=0) {
00238     for (int i=0; i<NumVectors_; i++) if (Vectors_[i]!=0) delete Vectors_[i];
00239     delete [] Vectors_;
00240   }
00241 
00242 
00243   if (DoubleTemp_!=0) delete [] DoubleTemp_;
00244 
00245 }
00246 
00247 //=========================================================================
00248 int Epetra_MultiVector::AllocateForCopy(void)
00249 {
00250 
00251   if (Allocated_) return(0);
00252 
00253   if (NumVectors_<=0)
00254     throw ReportError("Number of Vectors = " + toString(NumVectors_) + ", but must be greater than zero", -1);
00255 
00256   Stride_ = Map_.NumMyPoints();
00257   if (Stride_>0) Values_ = new double[Stride_ * NumVectors_];
00258   Pointers_ = new double *[NumVectors_];
00259 
00260   DoubleTemp_ = 0;
00261   Vectors_ = 0;
00262 
00263   int randval = rand(); // Use POSIX standard random function
00264   if (DistributedGlobal_)
00265     Util_.SetSeed(2*Comm_->MyPID() + randval);
00266   else {
00267     int locrandval = randval;
00268     Comm_->MaxAll(&locrandval, &randval, 1);
00269     Util_.SetSeed(randval); // Replicated Local MultiVectors must have same seeds
00270   }
00271 
00272   Allocated_ = true;
00273   UserAllocated_ = false;
00274   return(0);
00275 }
00276 
00277 //=========================================================================
00278 int Epetra_MultiVector::DoCopy(void)
00279 {
00280   // On entry Pointers_ contains pointers to the incoming vectors.  These
00281   // pointers are the only unique piece of information for each of the
00282   // constructors.
00283 
00284   // \internal { Optimization of this function can impact performance since it
00285   //           involves a fair amount of memory traffic.}
00286 
00287   // On exit, Pointers_ is redefined to point to its own MultiVector vectors.
00288 
00289   for (int i = 0; i< NumVectors_; i++)
00290     {
00291       double * from = Pointers_[i];
00292       double * to = Values_+i*Stride_;
00293       Pointers_[i] = to;
00294       const int myLength = MyLength_;
00295 #ifdef EPETRA_HAVE_OMP
00296 #pragma omp parallel for default(none) shared(to,from)
00297       for (int j=0; j<myLength; j++) to[j] = from[j];
00298 #else
00299       memcpy(to, from, myLength*sizeof(double));
00300 #endif
00301     }
00302 
00303   return(0);
00304 }
00305 //=========================================================================
00306 int Epetra_MultiVector::AllocateForView(void)
00307 {
00308 
00309   if (NumVectors_<=0)
00310     throw ReportError("Number of Vectors = " + toString(NumVectors_) + ", but must be greater than zero", -1);
00311 
00312   Pointers_ = new double *[NumVectors_];
00313 
00314   DoubleTemp_ = 0;
00315   Vectors_ = 0;
00316 
00317   int randval = rand(); // Use POSIX standard random function
00318   if (DistributedGlobal_)
00319     Util_.SetSeed(2*Comm_->MyPID() + randval);
00320   else {
00321     int locrandval = randval;
00322     Comm_->MaxAll(&locrandval, &randval, 1);
00323     Util_.SetSeed(randval); // Replicated Local MultiVectors must have same seeds
00324   }
00325 
00326   Allocated_ = true;
00327   UserAllocated_ = true;
00328 
00329   return(0);
00330 }
00331 
00332 //=========================================================================
00333 int Epetra_MultiVector::DoView(void)
00334 {
00335   // On entry Pointers_ contains pointers to the incoming vectors.  These
00336   // pointers are the only unique piece of information for each of the
00337   // constructors.
00338 
00339 
00340   Values_ = Pointers_[0];
00341 
00342   if (NumVectors_==1) {
00343     Stride_ = Map_.NumMyPoints();
00344     ConstantStride_ = true;
00345     return(0);
00346   }
00347 
00348   // Remainder of code checks if MultiVector has regular stride
00349 
00350   Stride_ = (int)(Pointers_[1] - Pointers_[0]);
00351   ConstantStride_ = false;
00352 
00353   for (int i = 1; i < NumVectors_-1; i++) {
00354     if (Pointers_[i+1] - Pointers_[i] != Stride_) return(0);
00355   }
00356 
00357   ConstantStride_ = true;
00358 
00359   return(0);
00360 }
00361 //=========================================================================
00362 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00363 int Epetra_MultiVector::ReplaceGlobalValue(int GlobalRow, int VectorIndex, double ScalarValue) {
00364 
00365  // Use the more general method below
00366   EPETRA_CHK_ERR(ChangeGlobalValue<int>(GlobalRow, 0, VectorIndex, ScalarValue, false));
00367   return(0);
00368 }
00369 #endif
00370 //=========================================================================
00371 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00372 int Epetra_MultiVector::ReplaceGlobalValue(long long GlobalRow, int VectorIndex, double ScalarValue) {
00373 
00374  // Use the more general method below
00375   EPETRA_CHK_ERR(ChangeGlobalValue<long long>(GlobalRow, 0, VectorIndex, ScalarValue, false));
00376   return(0);
00377 }
00378 #endif
00379 //=========================================================================
00380 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00381 int Epetra_MultiVector::ReplaceGlobalValue(int GlobalBlockRow, int BlockRowOffset,
00382              int VectorIndex, double ScalarValue) {
00383   // Use the more general method below
00384   EPETRA_CHK_ERR(ChangeGlobalValue<int>(GlobalBlockRow, BlockRowOffset, VectorIndex, ScalarValue, false));
00385   return(0);
00386 }
00387 #endif
00388 //=========================================================================
00389 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00390 int Epetra_MultiVector::ReplaceGlobalValue(long long GlobalBlockRow, int BlockRowOffset,
00391              int VectorIndex, double ScalarValue) {
00392   // Use the more general method below
00393   EPETRA_CHK_ERR(ChangeGlobalValue<long long>(GlobalBlockRow, BlockRowOffset, VectorIndex, ScalarValue, false));
00394   return(0);
00395 }
00396 #endif
00397 //=========================================================================
00398 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00399 int Epetra_MultiVector::SumIntoGlobalValue(int GlobalRow, int VectorIndex, double ScalarValue) {
00400 
00401   // Use the more general method below
00402   EPETRA_CHK_ERR(ChangeGlobalValue<int>(GlobalRow, 0, VectorIndex, ScalarValue, true));
00403   return(0);
00404 }
00405 #endif
00406 //=========================================================================
00407 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00408 int Epetra_MultiVector::SumIntoGlobalValue(long long GlobalRow, int VectorIndex, double ScalarValue) {
00409 
00410   // Use the more general method below
00411   EPETRA_CHK_ERR(ChangeGlobalValue<long long>(GlobalRow, 0, VectorIndex, ScalarValue, true));
00412   return(0);
00413 }
00414 #endif
00415 //=========================================================================
00416 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00417 int Epetra_MultiVector::SumIntoGlobalValue(int GlobalBlockRow, int BlockRowOffset,
00418              int VectorIndex, double ScalarValue) {
00419   // Use the more general method below
00420   EPETRA_CHK_ERR(ChangeGlobalValue<int>(GlobalBlockRow, BlockRowOffset, VectorIndex, ScalarValue, true));
00421   return(0);
00422 }
00423 #endif
00424 //=========================================================================
00425 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00426 int Epetra_MultiVector::SumIntoGlobalValue(long long GlobalBlockRow, int BlockRowOffset,
00427              int VectorIndex, double ScalarValue) {
00428   // Use the more general method below
00429   EPETRA_CHK_ERR(ChangeGlobalValue<long long>(GlobalBlockRow, BlockRowOffset, VectorIndex, ScalarValue, true));
00430   return(0);
00431 }
00432 #endif
00433 //=========================================================================
00434 int Epetra_MultiVector::ReplaceMyValue(int MyRow, int VectorIndex, double ScalarValue) {
00435 
00436   // Use the more general method below
00437   EPETRA_CHK_ERR(ChangeMyValue(MyRow, 0, VectorIndex, ScalarValue, false));
00438   return(0);
00439 }
00440 //=========================================================================
00441 int Epetra_MultiVector::ReplaceMyValue(int MyBlockRow, int BlockRowOffset,
00442              int VectorIndex, double ScalarValue) {
00443   // Use the more general method below
00444   EPETRA_CHK_ERR(ChangeMyValue(MyBlockRow, BlockRowOffset, VectorIndex, ScalarValue, false));
00445   return(0);
00446 }
00447 //=========================================================================
00448 int Epetra_MultiVector::SumIntoMyValue(int MyRow, int VectorIndex, double ScalarValue) {
00449   // Use the more general method below
00450   EPETRA_CHK_ERR(ChangeMyValue(MyRow, 0, VectorIndex, ScalarValue, true));
00451   return(0);
00452 }
00453 //=========================================================================
00454 int Epetra_MultiVector::SumIntoMyValue(int MyBlockRow, int BlockRowOffset,
00455              int VectorIndex, double ScalarValue) {
00456   // Use the more general method below
00457   EPETRA_CHK_ERR(ChangeMyValue(MyBlockRow, BlockRowOffset, VectorIndex, ScalarValue, true));
00458   return(0);
00459 }
00460 //=========================================================================
00461 template<typename int_type>
00462 int Epetra_MultiVector::ChangeGlobalValue(int_type GlobalBlockRow, int BlockRowOffset,
00463              int VectorIndex, double ScalarValue, bool SumInto) {
00464 
00465   if(!Map().template GlobalIndicesIsType<int_type>())
00466   throw ReportError("Epetra_MultiVector::ChangeGlobalValues mismatch between argument types (int/long long) and map type.", -1);
00467 
00468   // Convert GID to LID and call LID version
00469   EPETRA_CHK_ERR(ChangeMyValue(Map().LID(GlobalBlockRow), BlockRowOffset, VectorIndex, ScalarValue, SumInto));
00470   return(0);
00471 }
00472 //=========================================================================
00473 int Epetra_MultiVector::ChangeMyValue(int MyBlockRow, int BlockRowOffset,
00474              int VectorIndex, double ScalarValue, bool SumInto) {
00475 
00476   if (!Map().MyLID(MyBlockRow)) EPETRA_CHK_ERR(1); // I don't own this one, return a warning flag
00477   if (VectorIndex>= NumVectors_) EPETRA_CHK_ERR(-1); // Consider this a real error
00478   if (BlockRowOffset<0 || BlockRowOffset>=Map().ElementSize(MyBlockRow)) EPETRA_CHK_ERR(-2); // Offset is out-of-range
00479 
00480   int entry = Map().FirstPointInElement(MyBlockRow);
00481 
00482   if (SumInto)
00483     Pointers_[VectorIndex][entry+BlockRowOffset] += ScalarValue;
00484   else
00485     Pointers_[VectorIndex][entry+BlockRowOffset] = ScalarValue;
00486 
00487   return(0);
00488 }
00489 //=========================================================================
00490 int Epetra_MultiVector::Random() {
00491   // Generate random numbers drawn from a uniform distribution on
00492   // the interval (-1,1) using a multiplicative congruential generator
00493   // with modulus 2^31 - 1.
00494   /*
00495   const double a = 16807.0, BigInt=2147483647.0, DbleOne=1.0, DbleTwo=2.0;
00496 
00497   for(int i=0; i < NumVectors_; i++)
00498     for (int j=0; j<MyLength_; j++){
00499       Seed_ = fmod( a*Seed_, BigInt );
00500       Pointers_[i][j] = DbleTwo*(Seed_/BigInt)-DbleOne;
00501     }
00502   */
00503 
00504    const int myLength = MyLength_;
00505   for(int i = 0; i < NumVectors_; i++) {
00506           double * const to = Pointers_[i];
00507 #ifdef EPETRA_HAVE_OMP_NONASSOCIATIVE
00508 #pragma omp parallel for default(none)
00509 #endif
00510     for(int j = 0; j < myLength; j++)
00511       to[j] = Util_.RandomDouble();
00512         }
00513 
00514   return(0);
00515 }
00516 
00517 //=========================================================================
00518 
00519 // Extract a copy of a Epetra_MultiVector.  Put in a user's Fortran-style array
00520 
00521 int Epetra_MultiVector::ExtractCopy(double *A, int MyLDA) const {
00522   if (NumVectors_>1 && Stride_ > MyLDA) EPETRA_CHK_ERR(-1); // LDA not big enough
00523 
00524   const int myLength = MyLength_;
00525   for (int i=0; i< NumVectors_; i++)
00526     {
00527       double * from = Pointers_[i];
00528       double * to = A + i*MyLDA;
00529       for (int j=0; j<myLength; j++) *to++ = *from++;
00530     }
00531 
00532   return(0);
00533 }
00534 
00535 //=========================================================================
00536 int Epetra_MultiVector::ReplaceMap(const Epetra_BlockMap& map)
00537 {
00538   // mfh 28 Mar 2013: We can't check for compatibility across the
00539   // whole communicator, unless we know that the current and new
00540   // Maps are nonnull on _all_ participating processes.
00541 
00542   // So, we'll check to make sure that the maps are the same size on this processor and then
00543   // just go with it.
00544   if(Map().NumMyElements() == map.NumMyElements() && Map().NumMyPoints() == map.NumMyPoints()) {
00545     Epetra_DistObject::Map_ = map;
00546     return(0);
00547   }
00548 
00549   return(-1);
00550 }
00551 
00552 //=========================================================================
00553 
00554 // Extract a copy of a Epetra_MultiVector.  Put in a user's array of pointers
00555 
00556 int Epetra_MultiVector::ExtractCopy(double **ArrayOfPointers) const {
00557   const int myLength = MyLength_;
00558   for (int i=0; i< NumVectors_; i++)
00559     {
00560       double * from = Pointers_[i];
00561       double * to = ArrayOfPointers[i];
00562       memcpy(to, from, myLength*sizeof(double));
00563     }
00564 
00565   return(0);
00566 }
00567 
00568 
00569 
00570 //=========================================================================
00571 
00572 // Extract a view of a Epetra_MultiVector.  Set up a user's Fortran-style array
00573 
00574 int Epetra_MultiVector::ExtractView(double **A, int *MyLDA) const {
00575   if (!ConstantStride_) EPETRA_CHK_ERR(-1);  // Can't make a Fortran-style view if not constant stride
00576   *MyLDA = Stride_; // Set user's LDA
00577   *A = Values_; // Set user's value pointer
00578   return(0);
00579 }
00580 
00581 
00582 
00583 //=========================================================================
00584 
00585 // Extract a view of a Epetra_MultiVector.  Put in a user's array of pointers
00586 
00587 int Epetra_MultiVector::ExtractView(double ***ArrayOfPointers) const {
00588   *ArrayOfPointers = Pointers_;
00589 
00590   return(0);
00591 }
00592 
00593 
00594 //=========================================================================
00595 int Epetra_MultiVector::PutScalar(double ScalarConstant) {
00596 
00597   // Fills MultiVector with the value ScalarConstant **/
00598 
00599   const int myLength = MyLength_;
00600   for (int i = 0; i < NumVectors_; i++) {
00601     double * const to = Pointers_[i];
00602 #ifdef EPETRA_HAVE_OMP
00603 #pragma omp parallel for default(none) shared(ScalarConstant)
00604 #endif
00605     for (int j=0; j<myLength; j++) to[j] = ScalarConstant;
00606   }
00607   return(0);
00608 }
00609 //=========================================================================
00610 int Epetra_MultiVector::CheckSizes(const Epetra_SrcDistObject& Source) {
00611 
00612   const Epetra_MultiVector & A = dynamic_cast<const Epetra_MultiVector &>(Source);
00613 
00614   if (NumVectors()!=A.NumVectors()) {EPETRA_CHK_ERR(-3)};
00615   return(0);
00616 }
00617 
00618 //=========================================================================
00619 int Epetra_MultiVector::CopyAndPermute(const Epetra_SrcDistObject& Source,
00620                                        int NumSameIDs,
00621                                        int NumPermuteIDs,
00622                                        int * PermuteToLIDs,
00623                                        int *PermuteFromLIDs,
00624                                        const Epetra_OffsetIndex * Indexor)
00625 {
00626   (void)Indexor;
00627 
00628   const Epetra_MultiVector & A = dynamic_cast<const Epetra_MultiVector &>(Source);
00629 
00630   double **From = A.Pointers();
00631   double **To = Pointers_;
00632   int numVectors = NumVectors_;
00633 
00634   int * ToFirstPointInElementList = 0;
00635   int * FromFirstPointInElementList = 0;
00636   int * FromElementSizeList = 0;
00637   int MaxElementSize = Map().MaxElementSize();
00638   bool ConstantElementSize = Map().ConstantElementSize();
00639 
00640   if (!ConstantElementSize) {
00641     ToFirstPointInElementList =   Map().FirstPointInElementList();
00642     FromFirstPointInElementList = A.Map().FirstPointInElementList();
00643     FromElementSizeList = A.Map().ElementSizeList();
00644   }
00645   int jj, jjj, k;
00646 
00647   int NumSameEntries;
00648 
00649   bool Case1 = false;
00650   bool Case2 = false;
00651   // bool Case3 = false;
00652 
00653   if (MaxElementSize==1) {
00654     Case1 = true;
00655     NumSameEntries = NumSameIDs;
00656   }
00657   else if (ConstantElementSize) {
00658     Case2 = true;
00659     NumSameEntries = NumSameIDs * MaxElementSize;
00660   }
00661   else {
00662     // Case3 = true;
00663     NumSameEntries = FromFirstPointInElementList[NumSameIDs];
00664   }
00665 
00666   // Short circuit for the case where the source and target vector is the same.
00667   if (To==From) NumSameEntries = 0;
00668 
00669   // Do copy first
00670   if (NumSameIDs>0)
00671     if (To!=From) {
00672       for (int i=0; i < numVectors; i++)
00673   for (int j=0; j<NumSameEntries; j++)
00674     To[i][j] = From[i][j];
00675     }
00676   // Do local permutation next
00677   if (NumPermuteIDs>0) {
00678 
00679     // Point entry case
00680     if (Case1) {
00681 
00682       if (numVectors==1)
00683   for (int j=0; j<NumPermuteIDs; j++)
00684     To[0][PermuteToLIDs[j]] = From[0][PermuteFromLIDs[j]];
00685 
00686       else {
00687   for (int j=0; j<NumPermuteIDs; j++) {
00688     jj = PermuteToLIDs[j];
00689     jjj = PermuteFromLIDs[j];
00690     for (int i=0; i<numVectors; i++)
00691       To[i][jj] = From[i][jjj];
00692   }
00693       }
00694     }
00695     // constant element size case
00696     else if (Case2) {
00697 
00698       for (int j=0; j<NumPermuteIDs; j++) {
00699   jj = MaxElementSize*PermuteToLIDs[j];
00700   jjj = MaxElementSize*PermuteFromLIDs[j];
00701   for (int i=0; i<numVectors; i++)
00702     for (k=0; k<MaxElementSize; k++)
00703       To[i][jj+k] = From[i][jjj+k];
00704       }
00705     }
00706 
00707     // variable element size case
00708     else {
00709 
00710       for (int j=0; j<NumPermuteIDs; j++) {
00711   jj = ToFirstPointInElementList[PermuteToLIDs[j]];
00712   jjj = FromFirstPointInElementList[PermuteFromLIDs[j]];
00713   int ElementSize = FromElementSizeList[PermuteFromLIDs[j]];
00714   for (int i=0; i<numVectors; i++)
00715     for (k=0; k<ElementSize; k++)
00716       To[i][jj+k] = From[i][jjj+k];
00717       }
00718     }
00719   }
00720   return(0);
00721 }
00722 
00723 //=========================================================================
00724 int Epetra_MultiVector::PackAndPrepare(const Epetra_SrcDistObject & Source,
00725                                        int NumExportIDs,
00726                                        int * ExportLIDs,
00727                                        int & LenExports,
00728                                        char * & Exports,
00729                                        int & SizeOfPacket,
00730                                        int * Sizes,
00731                                        bool & VarSizes,
00732                                        Epetra_Distributor & Distor)
00733 {
00734   (void)Sizes;
00735   (void)VarSizes;
00736   (void)Distor;
00737 
00738   const Epetra_MultiVector & A = dynamic_cast<const Epetra_MultiVector &>(Source);
00739   int jj, k;
00740 
00741   double **From = A.Pointers();
00742   int MaxElementSize = Map().MaxElementSize();
00743   int numVectors = NumVectors_;
00744   bool ConstantElementSize = Map().ConstantElementSize();
00745 
00746   int * FromFirstPointInElementList = 0;
00747   int * FromElementSizeList = 0;
00748 
00749   if (!ConstantElementSize) {
00750     FromFirstPointInElementList = A.Map().FirstPointInElementList();
00751     FromElementSizeList = A.Map().ElementSizeList();
00752   }
00753 
00754   double * DoubleExports = 0;
00755 
00756   SizeOfPacket = numVectors*MaxElementSize*(int)sizeof(double);
00757 
00758   if(SizeOfPacket*NumExportIDs>LenExports) {
00759     if (LenExports>0) delete [] Exports;
00760     LenExports = SizeOfPacket*NumExportIDs;
00761     DoubleExports = new double[numVectors*MaxElementSize*NumExportIDs];
00762     Exports = (char *) DoubleExports;
00763   }
00764 
00765   double * ptr;
00766 
00767   if (NumExportIDs>0) {
00768     ptr = (double *) Exports;
00769 
00770     // Point entry case
00771     if (MaxElementSize==1) {
00772 
00773       if (numVectors==1)
00774         for (int j=0; j<NumExportIDs; j++)
00775           *ptr++ = From[0][ExportLIDs[j]];
00776 
00777       else {
00778   for (int j=0; j<NumExportIDs; j++) {
00779     jj = ExportLIDs[j];
00780     for (int i=0; i<numVectors; i++)
00781       *ptr++ = From[i][jj];
00782   }
00783       }
00784     }
00785 
00786     // constant element size case
00787     else if (ConstantElementSize) {
00788 
00789       for (int j=0; j<NumExportIDs; j++) {
00790   jj = MaxElementSize*ExportLIDs[j];
00791   for (int i=0; i<numVectors; i++)
00792     for (k=0; k<MaxElementSize; k++)
00793       *ptr++ = From[i][jj+k];
00794       }
00795     }
00796 
00797     // variable element size case
00798     else {
00799 
00800       int thisSizeOfPacket = numVectors*MaxElementSize;
00801       for (int j=0; j<NumExportIDs; j++) {
00802   ptr = (double *) Exports + j*thisSizeOfPacket;
00803   jj = FromFirstPointInElementList[ExportLIDs[j]];
00804   int ElementSize = FromElementSizeList[ExportLIDs[j]];
00805   for (int i=0; i<numVectors; i++)
00806     for (k=0; k<ElementSize; k++)
00807       *ptr++ = From[i][jj+k];
00808       }
00809     }
00810   }
00811 
00812   return(0);
00813 }
00814 
00815 //=========================================================================
00816 int Epetra_MultiVector::UnpackAndCombine(const Epetra_SrcDistObject & Source,
00817                                          int NumImportIDs,
00818                                          int * ImportLIDs,
00819                                          int LenImports,
00820                                          char * Imports,
00821                                          int & SizeOfPacket,
00822                                          Epetra_Distributor & Distor,
00823                                          Epetra_CombineMode CombineMode,
00824                                          const Epetra_OffsetIndex * Indexor )
00825 {
00826   (void)Source;
00827   (void)LenImports;
00828   (void)SizeOfPacket;
00829   (void)Distor;
00830   (void)Indexor;
00831   int jj, k;
00832 
00833   if(    CombineMode != Add
00834       && CombineMode != Zero
00835       && CombineMode != Insert
00836       && CombineMode != InsertAdd
00837       && CombineMode != Average
00838       && CombineMode != Epetra_Max
00839       && CombineMode != Epetra_Min
00840       && CombineMode != AbsMin
00841       && CombineMode != AbsMax )
00842       EPETRA_CHK_ERR(-1); //Unsupported CombinedMode, will default to Zero
00843 
00844   if (NumImportIDs<=0) return(0);
00845 
00846   double ** To = Pointers_;
00847   int numVectors = NumVectors_;
00848   int MaxElementSize = Map().MaxElementSize();
00849   bool ConstantElementSize = Map().ConstantElementSize();
00850 
00851   int * ToFirstPointInElementList = 0;
00852   int * ToElementSizeList = 0;
00853 
00854   if (!ConstantElementSize) {
00855     ToFirstPointInElementList = Map().FirstPointInElementList();
00856     ToElementSizeList = Map().ElementSizeList();
00857   }
00858 
00859   double * ptr;
00860   // Unpack it...
00861 
00862   ptr = (double *) Imports;
00863 
00864   // Point entry case
00865   if (MaxElementSize==1) {
00866 
00867     if (numVectors==1) {
00868       if (CombineMode==InsertAdd) for (int j=0; j<NumImportIDs; j++) To[0][ImportLIDs[j]] = 0.0; // Zero out first
00869       if (CombineMode==Add || CombineMode==InsertAdd) for (int j=0; j<NumImportIDs; j++) To[0][ImportLIDs[j]] += *ptr++; // Add to existing value
00870       else if(CombineMode==Insert) for (int j=0; j<NumImportIDs; j++) To[0][ImportLIDs[j]] = *ptr++;
00871       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++;}
00872       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++;}
00873       else if(CombineMode==Epetra_Max) for (int j=0; j<NumImportIDs; j++) {To[0][ImportLIDs[j]] = EPETRA_MAX( To[0][ImportLIDs[j]],*ptr); ptr++;}
00874       else if(CombineMode==Epetra_Min) for (int j=0; j<NumImportIDs; j++) {To[0][ImportLIDs[j]] = EPETRA_MIN( To[0][ImportLIDs[j]],*ptr); ptr++;}
00875       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
00876       // Note:  The following form of averaging is not a true average if more that one value is combined.
00877       //        This might be an issue in the future, but we leave this way for now.
00878 /*
00879       if (CombineMode==Add)
00880   for (int j=0; j<NumImportIDs; j++) To[0][ImportLIDs[j]] += *ptr++; // Add to existing value
00881       else if(CombineMode==Insert)
00882   for (int j=0; j<NumImportIDs; j++) To[0][ImportLIDs[j]] = *ptr++;
00883       else if(CombineMode==InsertAdd) {
00884   for (int j=0; j<NumImportIDs; j++) To[0][ImportLIDs[j]] = 0.0;
00885   for (int j=0; j<NumImportIDs; j++) To[0][ImportLIDs[j]] += *ptr++;
00886       }
00887       else if(CombineMode==AbsMax)
00888         for (int j=0; j<NumImportIDs; j++) {
00889     To[0][ImportLIDs[j]] = EPETRA_MAX( To[0][ImportLIDs[j]],std::abs(*ptr));
00890     ptr++;
00891   }
00892       // Note:  The following form of averaging is not a true average if more that one value is combined.
00893       //        This might be an issue in the future, but we leave this way for now.
00894       else if(CombineMode==Average)
00895   for (int j=0; j<NumImportIDs; j++) {To[0][ImportLIDs[j]] += *ptr++; To[0][ImportLIDs[j]] *= 0.5;}
00896 */
00897     }
00898 
00899     else {  // numVectors>1
00900 
00901       for (int j=0; j<NumImportIDs; j++) {
00902         jj = ImportLIDs[j];
00903         for (int i=0; i<numVectors; i++) {
00904         if (CombineMode==InsertAdd) To[i][jj] = 0.0; // Zero out if needed
00905         if (CombineMode==Add || CombineMode==InsertAdd) To[i][jj] += *ptr++; // Add to existing value
00906         else if (CombineMode==Insert) To[i][jj] = *ptr++; // Insert values
00907         else if (CombineMode==AbsMax) {To[i][jj] = EPETRA_MAX( To[i][jj], std::abs(*ptr)); ptr++; } // max of absolutes
00908         else if (CombineMode==AbsMin) {To[i][jj] = EPETRA_MIN( To[i][jj], std::abs(*ptr)); ptr++; } // max of absolutes
00909         else if (CombineMode==Epetra_Max)    {To[i][jj] = EPETRA_MAX( To[i][jj], *ptr); ptr++; } // simple max
00910         else if (CombineMode==Epetra_Min)    {To[i][jj] = EPETRA_MIN( To[i][jj], *ptr); ptr++; } // simple min
00911         else if (CombineMode==Average){To[i][jj] += *ptr++; To[i][jj] *= 0.5;}} // Not a true avg if >2 occurance of an ID
00912 
00913       }
00914 /*
00915       if (CombineMode==Add) {
00916   for (int j=0; j<NumImportIDs; j++) {
00917     jj = ImportLIDs[j];
00918     for (int i=0; i<numVectors; i++)
00919       To[i][jj] += *ptr++; // Add to existing value
00920   }
00921       }
00922       else if(CombineMode==Insert) {
00923   for (int j=0; j<NumImportIDs; j++) {
00924     jj = ImportLIDs[j];
00925     for (int i=0; i<numVectors; i++)
00926       To[i][jj] = *ptr++;
00927   }
00928       }
00929       else if(CombineMode==InsertAdd) {
00930   for (int j=0; j<NumImportIDs; j++) {
00931     jj = ImportLIDs[j];
00932     for (int i=0; i<numVectors; i++)
00933       To[i][jj] = 0.0;
00934   }
00935   for (int j=0; j<NumImportIDs; j++) {
00936     jj = ImportLIDs[j];
00937     for (int i=0; i<numVectors; i++)
00938       To[i][jj] += *ptr++;
00939   }
00940       }
00941       else if(CombineMode==AbsMax) {
00942         for (int j=0; j<NumImportIDs; j++) {
00943           jj = ImportLIDs[j];
00944           for (int i=0; i<numVectors; i++) {
00945             To[i][jj] = EPETRA_MAX( To[i][jj], std::abs(*ptr) );
00946       ptr++;
00947     }
00948         }
00949       }
00950       // Note:  The following form of averaging is not a true average if more that one value is combined.
00951       //        This might be an issue in the future, but we leave this way for now.
00952       else if(CombineMode==Average) {
00953   for (int j=0; j<NumImportIDs; j++) {
00954     jj = ImportLIDs[j];
00955     for (int i=0; i<numVectors; i++)
00956       { To[i][jj] += *ptr++;  To[i][jj] *= 0.5;}
00957   }
00958       }
00959 */
00960     }
00961   }
00962 
00963   // constant element size case
00964 
00965   else if (ConstantElementSize) {
00966 
00967     for (int j=0; j<NumImportIDs; j++) {
00968       jj = MaxElementSize*ImportLIDs[j];
00969       for (int i=0; i<numVectors; i++) {
00970         if (CombineMode==InsertAdd) for (k=0; k<MaxElementSize; k++) To[i][jj+k] = 0.0; // Zero out if needed
00971         if (CombineMode==Add || CombineMode==InsertAdd) for (k=0; k<MaxElementSize; k++) To[i][jj+k] += *ptr++; // Add to existing value
00972         else if (CombineMode==Insert) for (k=0; k<MaxElementSize; k++) To[i][jj+k] = *ptr++; // Insert values
00973         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
00974         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
00975         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
00976         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
00977         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
00978      }
00979    }
00980 /*
00981     if (CombineMode==Add) {
00982       for (int j=0; j<NumImportIDs; j++) {
00983   jj = MaxElementSize*ImportLIDs[j];
00984   for (int i=0; i<numVectors; i++)
00985     for (k=0; k<MaxElementSize; k++)
00986       To[i][jj+k] += *ptr++; // Add to existing value
00987       }
00988     }
00989     else if(CombineMode==Insert) {
00990       for (int j=0; j<NumImportIDs; j++) {
00991   jj = MaxElementSize*ImportLIDs[j];
00992   for (int i=0; i<numVectors; i++)
00993     for (k=0; k<MaxElementSize; k++)
00994       To[i][jj+k] = *ptr++;
00995       }
00996     }
00997     else if(CombineMode==InsertAdd) {
00998       for (int j=0; j<NumImportIDs; j++) {
00999   jj = MaxElementSize*ImportLIDs[j];
01000   for (int i=0; i<numVectors; i++)
01001     for (k=0; k<MaxElementSize; k++)
01002       To[i][jj+k] = 0.0;
01003       }
01004       for (int j=0; j<NumImportIDs; j++) {
01005   jj = MaxElementSize*ImportLIDs[j];
01006   for (int i=0; i<numVectors; i++)
01007     for (k=0; k<MaxElementSize; k++)
01008       To[i][jj+k] += *ptr++;
01009       }
01010     }
01011     else if(CombineMode==AbsMax) {
01012       for (int j=0; j<NumImportIDs; j++) {
01013   jj = MaxElementSize*ImportLIDs[j];
01014   for (int i=0; i<numVectors; i++)
01015     for (k=0; k<MaxElementSize; k++) {
01016       To[i][jj+k] = EPETRA_MAX( To[i][jj+k], std::abs(*ptr) );
01017       ptr++;
01018     }
01019       }
01020     }
01021     // Note:  The following form of averaging is not a true average if more that one value is combined.
01022     //        This might be an issue in the future, but we leave this way for now.
01023     else if(CombineMode==Average) {
01024       for (int j=0; j<NumImportIDs; j++) {
01025   jj = MaxElementSize*ImportLIDs[j];
01026   for (int i=0; i<numVectors; i++)
01027     for (k=0; k<MaxElementSize; k++)
01028       { To[i][jj+k] += *ptr++; To[i][jj+k] *= 0.5;}
01029       }
01030     }
01031 */
01032   }
01033 
01034   // variable element size case
01035 
01036   else {
01037     int thisSizeOfPacket = numVectors*MaxElementSize;
01038 
01039     for (int j=0; j<NumImportIDs; j++) {
01040       ptr = (double *) Imports + j*thisSizeOfPacket;
01041       jj = ToFirstPointInElementList[ImportLIDs[j]];
01042       int ElementSize = ToElementSizeList[ImportLIDs[j]];
01043       for (int i=0; i<numVectors; i++) {
01044         if (CombineMode==InsertAdd) for (k=0; k<ElementSize; k++) To[i][jj+k] = 0.0; // Zero out if needed
01045         if (CombineMode==Add || CombineMode==InsertAdd) for (k=0; k<ElementSize; k++) To[i][jj+k] += *ptr++; // Add to existing value
01046         else if (CombineMode==Insert) for (k=0; k<ElementSize; k++) To[i][jj+k] = *ptr++; // Insert values
01047         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
01048         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
01049         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
01050         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
01051      }
01052    }
01053 /*
01054     if (CombineMode==Add) {
01055       for (int j=0; j<NumImportIDs; j++) {
01056   ptr = (double *) Imports + j*thisSizeOfPacket;
01057   jj = ToFirstPointInElementList[ImportLIDs[j]];
01058   int ElementSize = ToElementSizeList[ImportLIDs[j]];
01059   for (int i=0; i<numVectors; i++)
01060     for (k=0; k<ElementSize; k++)
01061       To[i][jj+k] += *ptr++; // Add to existing value
01062       }
01063     }
01064     else  if(CombineMode==Insert){
01065       for (int j=0; j<NumImportIDs; j++) {
01066   ptr = (double *) Imports + j*thisSizeOfPacket;
01067   jj = ToFirstPointInElementList[ImportLIDs[j]];
01068   int ElementSize = ToElementSizeList[ImportLIDs[j]];
01069   for (int i=0; i<numVectors; i++)
01070     for (k=0; k<ElementSize; k++)
01071       To[i][jj+k] = *ptr++;
01072       }
01073     }
01074     else  if(CombineMode==InsertAdd){
01075       for (int j=0; j<NumImportIDs; j++) {
01076   ptr = (double *) Imports + j*thisSizeOfPacket;
01077   jj = ToFirstPointInElementList[ImportLIDs[j]];
01078   int ElementSize = ToElementSizeList[ImportLIDs[j]];
01079   for (int i=0; i<numVectors; i++)
01080     for (k=0; k<ElementSize; k++)
01081       To[i][jj+k] = 0.0;
01082       }
01083       for (int j=0; j<NumImportIDs; j++) {
01084   ptr = (double *) Imports + j*thisSizeOfPacket;
01085   jj = ToFirstPointInElementList[ImportLIDs[j]];
01086   int ElementSize = ToElementSizeList[ImportLIDs[j]];
01087   for (int i=0; i<numVectors; i++)
01088     for (k=0; k<ElementSize; k++)
01089       To[i][jj+k] += *ptr++;
01090       }
01091     }
01092     else  if(CombineMode==AbsMax){
01093       for (int j=0; j<NumImportIDs; j++) {
01094   ptr = (double *) Imports + j*thisSizeOfPacket;
01095   jj = ToFirstPointInElementList[ImportLIDs[j]];
01096   int ElementSize = ToElementSizeList[ImportLIDs[j]];
01097   for (int i=0; i<numVectors; i++)
01098     for (k=0; k<ElementSize; k++) {
01099       To[i][jj+k] = EPETRA_MAX( To[i][jj+k], std::abs(*ptr));
01100       ptr++;
01101     }
01102       }
01103     }
01104     // Note:  The following form of averaging is not a true average if more that one value is combined.
01105     //        This might be an issue in the future, but we leave this way for now.
01106     else if(CombineMode==Average) {
01107       for (int j=0; j<NumImportIDs; j++) {
01108   ptr = (double *) Imports + j*thisSizeOfPacket;
01109   jj = ToFirstPointInElementList[ImportLIDs[j]];
01110   int ElementSize = ToElementSizeList[ImportLIDs[j]];
01111   for (int i=0; i<numVectors; i++)
01112     for (k=0; k<ElementSize; k++)
01113       { To[i][jj+k] += *ptr++; To[i][jj+k] *= 0.5;}
01114       }
01115     }
01116 */
01117   }
01118 
01119   return(0);
01120 }
01121 
01122 //=========================================================================
01123 int Epetra_MultiVector::Dot(const Epetra_MultiVector& A, double *Result) const {
01124 
01125   // Dot product of two MultiVectors
01126 
01127   if (NumVectors_ != A.NumVectors()) EPETRA_CHK_ERR(-1);
01128   const int myLength = MyLength_;
01129   if (myLength != A.MyLength()) EPETRA_CHK_ERR(-2);
01130   UpdateDoubleTemp();
01131 
01132   double **A_Pointers = A.Pointers();
01133 
01134 #ifdef EPETRA_HAVE_OMP_NONASSOCIATIVE
01135   for (int i=0; i < NumVectors_; i++)
01136     {
01137       const double * const from = Pointers_[i];
01138       const double * const fromA = A_Pointers[i];
01139       double sum = 0.0;
01140 #pragma omp parallel for reduction (+:sum) default(none)
01141       for (int j=0; j < myLength; j++) sum += from[j] * fromA[j];
01142       DoubleTemp_[i] = sum;
01143     }
01144 #else
01145   for (int i=0; i < NumVectors_; i++) DoubleTemp_[i] = DOT(myLength, Pointers_[i], A_Pointers[i]);
01146 #endif
01147 
01148   if (DistributedGlobal())
01149     Comm_->SumAll(DoubleTemp_, Result, NumVectors_);
01150   else
01151     for (int i=0; i< NumVectors_; ++i) Result[i] = DoubleTemp_[i];
01152 
01153   UpdateFlops(2*GlobalLength_*NumVectors_);
01154 
01155   return(0);
01156 }
01157 //=========================================================================
01158 int Epetra_MultiVector::Abs(const Epetra_MultiVector& A) {
01159 
01160   // this[i][j] = std::abs(A[i][j])
01161 
01162   const int myLength = MyLength_;
01163   if (NumVectors_ != A.NumVectors()) EPETRA_CHK_ERR(-1);
01164   if (myLength != A.MyLength()) EPETRA_CHK_ERR(-2);
01165 
01166   double **A_Pointers = A.Pointers();
01167 
01168   for (int i=0; i < NumVectors_; i++) {
01169     double * const to = Pointers_[i];
01170     const double * const from = A_Pointers[i];
01171 #ifdef EPETRA_HAVE_OMP
01172 #pragma omp parallel for default(none)
01173 #endif
01174     for (int j=0; j < myLength; j++) to[j] = std::abs(from[j]);
01175   }
01176 
01177   return(0);
01178 }
01179 //=========================================================================
01180 int Epetra_MultiVector::Reciprocal(const Epetra_MultiVector& A) {
01181 
01182   // this[i][j] = 1.0/(A[i][j])
01183 
01184   int ierr = 0;
01185 #ifndef EPETRA_HAVE_OMP
01186   int localierr = 0;
01187 #endif
01188   const int myLength = MyLength_;
01189   if (NumVectors_ != A.NumVectors()) EPETRA_CHK_ERR(-1);
01190   if (myLength != A.MyLength()) EPETRA_CHK_ERR(-2);
01191 
01192   double **A_Pointers = A.Pointers();
01193 
01194   for (int i=0; i < NumVectors_; i++) {
01195     double * const to = Pointers_[i];
01196     const double * const from = A_Pointers[i];
01197 #ifdef EPETRA_HAVE_OMP
01198 #pragma omp parallel default(none) shared(ierr)
01199 {
01200     int localierr = 0;
01201 #pragma omp for
01202 #endif
01203     for (int j=0; j < myLength; j++) {
01204       double value = from[j];
01205       // Set error to 1 to signal that zero rowsum found (supercedes ierr = 2)
01206       if (std::abs(value)<Epetra_MinDouble) {
01207         if (value==0.0) localierr = 1;
01208         else if (localierr!=1) localierr = 2;
01209         to[j] = EPETRA_SGN(value) * Epetra_MaxDouble;
01210       }
01211       else
01212         to[j] = 1.0/value;
01213     }
01214 #ifdef EPETRA_HAVE_OMP
01215 #pragma omp critical
01216 #endif
01217 {
01218    if (localierr==1) ierr = 1;
01219    else if (localierr==2 && ierr!=1) ierr = 2;
01220 }
01221 #ifdef EPETRA_HAVE_OMP
01222 }
01223 #endif
01224   }
01225   EPETRA_CHK_ERR(ierr);
01226   return(0);
01227 }
01228   //=========================================================================
01229   int Epetra_MultiVector::Scale (double ScalarValue) {
01230 
01231     // scales a MultiVector in place by a scalar
01232 
01233 
01234   const int myLength = MyLength_;
01235 #ifdef EPETRA_HAVE_OMP
01236     for (int i = 0; i < NumVectors_; i++) {
01237       double * const to = Pointers_[i];
01238 #pragma omp parallel for default(none) shared(ScalarValue)
01239       for (int j = 0; j < myLength; j++) to[j] = ScalarValue * to[j];
01240     }
01241 #else
01242     for (int i = 0; i < NumVectors_; i++)
01243       SCAL(myLength, ScalarValue, Pointers_[i]);
01244 #endif
01245     UpdateFlops(GlobalLength_*NumVectors_);
01246 
01247     return(0);
01248   }
01249 
01250   //=========================================================================
01251   int Epetra_MultiVector::Scale (double ScalarA, const Epetra_MultiVector& A) {
01252 
01253     // scales a MultiVector by a scalar and put in the this:
01254     // this = ScalarA * A
01255 
01256   const int myLength = MyLength_;
01257   if (NumVectors_ != A.NumVectors()) EPETRA_CHK_ERR(-1);
01258   if (myLength != A.MyLength()) EPETRA_CHK_ERR(-2);
01259 
01260     double **A_Pointers = (double**)A.Pointers();
01261 
01262     for (int i = 0; i < NumVectors_; i++) {
01263       double * const to = Pointers_[i];
01264       const double * const from = A_Pointers[i];
01265 #ifdef EPETRA_HAVE_OMP
01266 #pragma omp parallel for default(none) shared(ScalarA)
01267 #endif
01268       for (int j = 0; j < myLength; j++) to[j] = ScalarA * from[j];
01269     }
01270     UpdateFlops(GlobalLength_*NumVectors_);
01271 
01272     return(0);
01273   }
01274 
01275   //=========================================================================
01276   int Epetra_MultiVector::Update(double ScalarA, const Epetra_MultiVector& A, double ScalarThis) {
01277 
01278 
01279     // linear combination of two MultiVectors: this = ScalarThis * this + ScalarA * A
01280 
01281   const int myLength = MyLength_;
01282   if (NumVectors_ != A.NumVectors()) EPETRA_CHK_ERR(-1);
01283   if (myLength != A.MyLength()) EPETRA_CHK_ERR(-2);
01284 
01285     double **A_Pointers = (double**)A.Pointers();
01286 
01287     if (ScalarThis==0.0)
01288       {
01289   for (int i = 0; i < NumVectors_; i++) {
01290           double * const to = Pointers_[i];
01291           const double * const from = A_Pointers[i];
01292 #ifdef EPETRA_HAVE_OMP
01293 #pragma omp parallel for default(none) shared(ScalarA)
01294 #endif
01295     for (int j = 0; j < myLength; j++) to[j] = ScalarA * from[j];
01296         }
01297   UpdateFlops(GlobalLength_*NumVectors_);
01298       }
01299     else if (ScalarThis==1.0)
01300       {
01301   for (int i = 0; i < NumVectors_; i++) {
01302           double * const to = Pointers_[i];
01303           const double * const from = A_Pointers[i];
01304 #ifdef EPETRA_HAVE_OMP
01305 #pragma omp parallel for default(none) shared(ScalarA)
01306 #endif
01307     for (int j = 0; j < myLength; j++) to[j] = to[j] + ScalarA * from[j];
01308         }
01309   UpdateFlops(2*GlobalLength_*NumVectors_);
01310       }
01311     else if (ScalarA==1.0)
01312       {
01313   for (int i = 0; i < NumVectors_; i++) {
01314           double * const to = Pointers_[i];
01315           const double * const from = A_Pointers[i];
01316 #ifdef EPETRA_HAVE_OMP
01317 #pragma omp parallel for default(none) shared(ScalarThis)
01318 #endif
01319     for (int j = 0; j < myLength; j++) to[j] = ScalarThis * to[j] + from[j];
01320         }
01321   UpdateFlops(2*GlobalLength_*NumVectors_);
01322       }
01323     else
01324       {
01325   for (int i = 0; i < NumVectors_; i++) {
01326           double * const to = Pointers_[i];
01327           const double * const from = A_Pointers[i];
01328 #ifdef EPETRA_HAVE_OMP
01329 #pragma omp parallel for default(none) shared(ScalarA,ScalarThis)
01330 #endif
01331     for (int j = 0; j < myLength; j++) to[j] = ScalarThis * to[j] +
01332                               ScalarA *  from[j];
01333         }
01334   UpdateFlops(3*GlobalLength_*NumVectors_);
01335       }
01336 
01337     return(0);
01338   }
01339 
01340 //=========================================================================
01341 int Epetra_MultiVector::Update(double ScalarA, const Epetra_MultiVector& A,
01342           double ScalarB, const Epetra_MultiVector& B, double ScalarThis) {
01343 
01344 
01345   // linear combination of three MultiVectors:
01346   // this = ScalarThis * this + ScalarA * A + ScalarB * B
01347 
01348   if (ScalarA==0.0) {
01349     EPETRA_CHK_ERR(Update(ScalarB, B, ScalarThis));
01350     return(0);
01351   }
01352   if (ScalarB==0.0) {
01353     EPETRA_CHK_ERR(Update(ScalarA, A, ScalarThis));
01354     return(0);
01355   }
01356 
01357   const int myLength = MyLength_;
01358   if (NumVectors_ != A.NumVectors() || NumVectors_ != B.NumVectors()) EPETRA_CHK_ERR(-1);
01359   if (myLength != A.MyLength() || myLength != B.MyLength()) EPETRA_CHK_ERR(-2);
01360 
01361     double **A_Pointers = (double**)A.Pointers();
01362     double **B_Pointers = (double**)B.Pointers();
01363 
01364     if (ScalarThis==0.0)
01365       {
01366   if (ScalarA==1.0)
01367     {
01368       for (int i = 0; i < NumVectors_; i++) {
01369               double * const to = Pointers_[i];
01370               const double * const fromA = A_Pointers[i];
01371               const double * const fromB = B_Pointers[i];
01372 #ifdef EPETRA_HAVE_OMP
01373 #pragma omp parallel for default(none) shared(ScalarB)
01374 #endif
01375         for (int j = 0; j < myLength; j++) to[j] =           fromA[j] +
01376                             ScalarB * fromB[j];
01377             }
01378       UpdateFlops(2*GlobalLength_*NumVectors_);
01379     }
01380   else if (ScalarB==1.0)
01381     {
01382       for (int i = 0; i < NumVectors_; i++) {
01383               double * const to = Pointers_[i];
01384               const double * const fromA = A_Pointers[i];
01385               const double * const fromB = B_Pointers[i];
01386 #ifdef EPETRA_HAVE_OMP
01387 #pragma omp parallel for default(none) shared(ScalarA)
01388 #endif
01389         for (int j = 0; j < myLength; j++) to[j] = ScalarA * fromA[j] +
01390                                       fromB[j];
01391             }
01392       UpdateFlops(2*GlobalLength_*NumVectors_);
01393     }
01394   else
01395     {
01396       for (int i = 0; i < NumVectors_; i++) {
01397               double * const to = Pointers_[i];
01398               const double * const fromA = A_Pointers[i];
01399               const double * const fromB = B_Pointers[i];
01400 #ifdef EPETRA_HAVE_OMP
01401 #pragma omp parallel for default(none) shared(ScalarA,ScalarB)
01402 #endif
01403         for (int j = 0; j < myLength; j++) to[j] = ScalarA * fromA[j] +
01404                             ScalarB * fromB[j];
01405             }
01406       UpdateFlops(3*GlobalLength_*NumVectors_);
01407     }
01408       }
01409     else if (ScalarThis==1.0)
01410       {
01411   if (ScalarA==1.0)
01412     {
01413       for (int i = 0; i < NumVectors_; i++) {
01414               double * const to = Pointers_[i];
01415               const double * const fromA = A_Pointers[i];
01416               const double * const fromB = B_Pointers[i];
01417 #ifdef EPETRA_HAVE_OMP
01418 #pragma omp parallel for default(none) shared(ScalarB)
01419 #endif
01420         for (int j = 0; j < myLength; j++) to[j] +=           fromA[j] +
01421                              ScalarB * fromB[j];
01422             }
01423       UpdateFlops(3*GlobalLength_*NumVectors_);
01424     }
01425   else if (ScalarB==1.0)
01426     {
01427       for (int i = 0; i < NumVectors_; i++) {
01428               double * const to = Pointers_[i];
01429               const double * const fromA = A_Pointers[i];
01430               const double * const fromB = B_Pointers[i];
01431 #ifdef EPETRA_HAVE_OMP
01432 #pragma omp parallel for default(none) shared(ScalarA)
01433 #endif
01434         for (int j = 0; j < myLength; j++) to[j] += ScalarA * fromA[j] +
01435                                        fromB[j];
01436             }
01437       UpdateFlops(3*GlobalLength_*NumVectors_);
01438     }
01439   else
01440     {
01441       for (int i = 0; i < NumVectors_; i++) {
01442               double * const to = Pointers_[i];
01443               const double * const fromA = A_Pointers[i];
01444               const double * const fromB = B_Pointers[i];
01445 #ifdef EPETRA_HAVE_OMP
01446 #pragma omp parallel for default(none) shared(ScalarA,ScalarB)
01447 #endif
01448         for (int j = 0; j < myLength; j++) to[j] += ScalarA * fromA[j] +
01449                              ScalarB * fromB[j];
01450             }
01451       UpdateFlops(4*GlobalLength_*NumVectors_);
01452     }
01453       }
01454     else
01455       {
01456   if (ScalarA==1.0)
01457     {
01458       for (int i = 0; i < NumVectors_; i++) {
01459               double * const to = Pointers_[i];
01460               const double * const fromA = A_Pointers[i];
01461               const double * const fromB = B_Pointers[i];
01462 #ifdef EPETRA_HAVE_OMP
01463 #pragma omp parallel for default(none) shared(ScalarA,ScalarB,ScalarThis)
01464 #endif
01465         for (int j = 0; j < myLength; j++) to[j] =  ScalarThis *    to[j] +
01466                                        fromA[j] +
01467                              ScalarB * fromB[j];
01468             }
01469       UpdateFlops(4*GlobalLength_*NumVectors_);
01470     }
01471   else if (ScalarB==1.0)
01472     {
01473       for (int i = 0; i < NumVectors_; i++) {
01474               double * const to = Pointers_[i];
01475               const double * const fromA = A_Pointers[i];
01476               const double * const fromB = B_Pointers[i];
01477 #ifdef EPETRA_HAVE_OMP
01478 #pragma omp parallel for default(none) shared(ScalarA,ScalarThis)
01479 #endif
01480         for (int j = 0; j < myLength; j++) to[j] =  ScalarThis *    to[j] +
01481                              ScalarA * fromA[j] +
01482                                        fromB[j];
01483             }
01484       UpdateFlops(4*GlobalLength_*NumVectors_);
01485     }
01486   else
01487     {
01488       for (int i = 0; i < NumVectors_; i++) {
01489               double * const to = Pointers_[i];
01490               const double * const fromA = A_Pointers[i];
01491               const double * const fromB = B_Pointers[i];
01492 #ifdef EPETRA_HAVE_OMP
01493 #pragma omp parallel for default(none) shared(ScalarA,ScalarB,ScalarThis)
01494 #endif
01495         for (int j = 0; j < myLength; j++) to[j] =  ScalarThis *    to[j] +
01496                              ScalarA * fromA[j] +
01497                              ScalarB * fromB[j];
01498             }
01499       UpdateFlops(5*GlobalLength_*NumVectors_);
01500     }
01501       }
01502 
01503 
01504     return(0);
01505   }
01506 
01507 //=========================================================================
01508 int  Epetra_MultiVector::Norm1 (double* Result) const {
01509 
01510   // 1-norm of each vector in MultiVector
01511 
01512   if (!Map().UniqueGIDs()) {EPETRA_CHK_ERR(-1);}
01513 
01514   const int myLength = MyLength_;
01515   UpdateDoubleTemp();
01516 #ifdef EPETRA_HAVE_OMP_NONASSOCIATIVE
01517   for (int i=0; i < NumVectors_; i++)
01518     {
01519       const double * const from = Pointers_[i];
01520       double asum = 0.0;
01521 #pragma omp parallel default(none) shared(asum)
01522 {
01523       double localasum = 0.0;
01524 #pragma omp for
01525       for (int j=0; j< myLength; j++) localasum += std::abs(from[j]);
01526 #pragma omp critical
01527       asum += localasum;
01528 }
01529       DoubleTemp_[i] = asum;
01530   }
01531 #else
01532 
01533   for (int i=0; i < NumVectors_; i++) DoubleTemp_[i] = ASUM(myLength, Pointers_[i]);
01534 #endif
01535 
01536   if (DistributedGlobal())
01537     Comm_->SumAll(DoubleTemp_, Result, NumVectors_);
01538   else
01539     for (int i=0; i< NumVectors_; ++i) Result[i] = DoubleTemp_[i];
01540 
01541   UpdateFlops(2*GlobalLength_*NumVectors_);
01542 
01543   return(0);
01544 }
01545 
01546 //=========================================================================
01547 int  Epetra_MultiVector::Norm2 (double* Result) const {
01548 
01549   // 2-norm of each vector in MultiVector
01550 
01551 
01552   if (!Map().UniqueGIDs()) {EPETRA_CHK_ERR(-1);}
01553 
01554   const int myLength = MyLength_;
01555   UpdateDoubleTemp();
01556 
01557   for (int i=0; i < NumVectors_; i++)
01558     {
01559       const double * const from = Pointers_[i];
01560       double sum = 0.0;
01561 #ifdef EPETRA_HAVE_OMP_NONASSOCIATIVE
01562 #pragma omp parallel for reduction (+:sum) default(none)
01563 #endif
01564       for (int j=0; j < myLength; j++) sum += from[j] * from[j];
01565       DoubleTemp_[i] = sum;
01566     }
01567   if (DistributedGlobal())
01568     Comm_->SumAll(DoubleTemp_, Result, NumVectors_);
01569   else
01570     for (int i=0; i< NumVectors_; ++i) Result[i] = DoubleTemp_[i];
01571 
01572   for (int i=0; i < NumVectors_; i++) Result[i] = std::sqrt(Result[i]);
01573 
01574   UpdateFlops(2*GlobalLength_*NumVectors_);
01575 
01576   return(0);
01577 }
01578 
01579 //=========================================================================
01580 int  Epetra_MultiVector::NormInf (double* Result) const {
01581 
01582   // Inf-norm of each vector in MultiVector
01583 
01584 
01585   const int myLength = MyLength_;
01586   UpdateDoubleTemp();
01587 
01588   for (int i=0; i < NumVectors_; i++)
01589     {
01590       DoubleTemp_[i] = 0.0;
01591       double normval = 0.0;
01592       const double * const from = Pointers_[i];
01593       if (myLength>0) normval = from[0];
01594 #ifdef EPETRA_HAVE_OMP
01595 #pragma omp parallel default(none) shared(normval)
01596 {
01597       double localnormval = 0.0;
01598 #pragma omp for
01599       for (int j=0; j< myLength; j++) {
01600          localnormval = EPETRA_MAX(localnormval,std::abs(from[j]));
01601       }
01602 #pragma omp critical
01603       {
01604       normval = EPETRA_MAX(normval,localnormval);
01605       }
01606 }
01607       DoubleTemp_[i] = normval;
01608 #else
01609       (void) normval; // silence unused value warning in non-OpenMP build
01610       int jj = IAMAX(myLength, Pointers_[i]);
01611       if (jj>-1) DoubleTemp_[i] = std::abs(Pointers_[i][jj]);
01612 #endif
01613     }
01614   if (DistributedGlobal())
01615     Comm_->MaxAll(DoubleTemp_, Result, NumVectors_);
01616   else
01617     for (int i=0; i< NumVectors_; ++i) Result[i] = DoubleTemp_[i];
01618 
01619   // UpdateFlops(0);  Strictly speaking there are not FLOPS in this routine
01620   return(0);
01621 }
01622 
01623 //=========================================================================
01624 int  Epetra_MultiVector::NormWeighted (const Epetra_MultiVector& Weights, double* Result) const {
01625 
01626   // Weighted 2-norm of each vector in MultiVector
01627 
01628   // If only one vector in Weights, we assume it will be used as the weights for all vectors
01629 
01630   const int myLength = MyLength_;
01631   bool OneW = false;
01632   if (Weights.NumVectors()==1) OneW = true;
01633   else if (NumVectors_ != Weights.NumVectors()) EPETRA_CHK_ERR(-1);
01634 
01635   if (myLength != Weights.MyLength()) EPETRA_CHK_ERR(-2);
01636 
01637 
01638   UpdateDoubleTemp();
01639 
01640   double *W = Weights.Values();
01641   double **W_Pointers = Weights.Pointers();
01642 
01643   for (int i=0; i < NumVectors_; i++)
01644     {
01645       if (!OneW) W = W_Pointers[i]; // If Weights has the same number of vectors as this, use each weight vector
01646       double sum = 0.0;
01647       const double * const from = Pointers_[i];
01648 #ifdef EPETRA_HAVE_OMP_NONASSOCIATIVE
01649 #pragma omp parallel for reduction (+:sum) default(none) shared(W)
01650 #endif
01651       for (int j=0; j < myLength; j++) {
01652         double tmp = from[j]/W[j];
01653         sum += tmp * tmp;
01654       }
01655       DoubleTemp_[i] = sum;
01656     }
01657   double OneOverN;
01658   if (DistributedGlobal()) {
01659     Comm_->SumAll(DoubleTemp_, Result, NumVectors_);
01660     OneOverN = 1.0 / (double) GlobalLength_;
01661   }
01662   else {
01663     for (int i=0; i< NumVectors_; ++i) Result[i] = DoubleTemp_[i];
01664     OneOverN = 1.0 / (double) myLength;
01665   }
01666   for (int i=0; i < NumVectors_; i++) Result[i] = std::sqrt(Result[i]*OneOverN);
01667 
01668   UpdateFlops(3*GlobalLength_*NumVectors_);
01669 
01670   return(0);
01671   }
01672 
01673 //=========================================================================
01674 int  Epetra_MultiVector::MinValue (double* Result) const {
01675 
01676   // Minimum value of each vector in MultiVector
01677 
01678   int ierr = 0;
01679 
01680   const int myLength = MyLength_;
01681   UpdateDoubleTemp();
01682 
01683   for (int i=0; i < NumVectors_; i++)
01684     {
01685       const double * const from = Pointers_[i];
01686       double MinVal = Epetra_MaxDouble;
01687       if (myLength>0) MinVal = from[0];
01688 #ifdef EPETRA_HAVE_OMP
01689 #pragma omp parallel default(none) shared(MinVal)
01690 {
01691       double localMinVal = MinVal;
01692 #pragma omp for
01693       for (int j=0; j< myLength; j++) localMinVal = EPETRA_MIN(localMinVal,from[j]);
01694 #pragma omp critical
01695       {
01696       MinVal = EPETRA_MIN(MinVal,localMinVal);
01697       }
01698 }
01699       DoubleTemp_[i] = MinVal;
01700 #else
01701       for (int j=0; j< myLength; j++) MinVal = EPETRA_MIN(MinVal,from[j]);
01702       DoubleTemp_[i] = MinVal;
01703 #endif
01704     }
01705 
01706   if (myLength > 0) {
01707     for(int i=0; i<NumVectors_; ++i) {
01708       Result[i] = DoubleTemp_[i];
01709     }
01710   }
01711 
01712   //If myLength == 0 and Comm_->NumProc() == 1, then Result has
01713   //not been referenced. Also, if vector contents are uninitialized
01714   //then Result contents are not well defined...
01715 
01716   if (Comm_->NumProc() == 1 || !DistributedGlobal()) return(ierr);
01717 
01718   //We're going to use MPI_Allgather to gather every proc's local-
01719   //min values onto every other proc. We'll use the last position
01720   //of the DoubleTemp_ array to indicate whether this proc has
01721   //valid data that should be considered by other procs when forming
01722   //the global-min results.
01723 
01724   if (myLength == 0) DoubleTemp_[NumVectors_] = 0.0;
01725   else DoubleTemp_[NumVectors_] = 1.0;
01726 
01727   //Now proceed to handle the parallel case. We'll gather local-min
01728   //values from other procs and form a global-min. If any processor
01729   //has myLength>0, we'll end up with a valid result.
01730 
01731 #ifdef EPETRA_MPI
01732   const Epetra_MpiComm* epetrampicomm =
01733     dynamic_cast<const Epetra_MpiComm*>(Comm_);
01734   if (!epetrampicomm) {
01735     return(-2);
01736   }
01737 
01738   MPI_Comm mpicomm = epetrampicomm->GetMpiComm();
01739   int numProcs = epetrampicomm->NumProc();
01740   double* dwork = new double[numProcs*(NumVectors_+1)];
01741 
01742   MPI_Allgather(DoubleTemp_, NumVectors_+1, MPI_DOUBLE,
01743                 dwork, NumVectors_+1, MPI_DOUBLE, mpicomm);
01744 
01745   //if myLength==0, then our Result array currently contains
01746   //Epetra_MaxDouble from the local-min calculations above. In this
01747   //case we'll overwrite our Result array with values from the first
01748   //processor that sent valid data.
01749   bool overwrite = myLength == 0 ? true : false;
01750 
01751   int myPID = epetrampicomm->MyPID();
01752   double* dwork_ptr = dwork;
01753 
01754   for(int j=0; j<numProcs; ++j) {
01755 
01756     //skip data from self, and skip data from
01757     //procs with DoubleTemp_[NumVectors_] == 0.0.
01758     if (j == myPID || dwork_ptr[NumVectors_] < 0.5) {
01759       dwork_ptr += NumVectors_+1;
01760       continue;
01761     }
01762 
01763     for(int i=0; i<NumVectors_; ++i) {
01764       double val = dwork_ptr[i];
01765 
01766       //Set val into our Result array if overwrite is true (see above),
01767       //or if val is less than our current Result[i].
01768       if (overwrite || (Result[i] > val)) Result[i] = val;
01769     }
01770 
01771     //Now set overwrite to false so that we'll do the right thing
01772     //when processing data from subsequent processors.
01773     if (overwrite) overwrite = false;
01774 
01775     dwork_ptr += NumVectors_+1;
01776   }
01777 
01778   delete [] dwork;
01779 #endif
01780 
01781   // UpdateFlops(0);  Strictly speaking there are not FLOPS in this routine
01782 
01783   return(ierr);
01784 }
01785 
01786 //=========================================================================
01787 int  Epetra_MultiVector::MaxValue (double* Result) const {
01788 
01789   // Maximum value of each vector in MultiVector
01790 
01791   int ierr = 0;
01792 
01793   const int myLength = MyLength_;
01794   UpdateDoubleTemp();
01795 
01796   for (int i=0; i < NumVectors_; i++)
01797     {
01798       const double * const from = Pointers_[i];
01799       double MaxVal = -Epetra_MaxDouble;
01800       if (myLength>0) MaxVal = from[0];
01801 #ifdef EPETRA_HAVE_OMP
01802 #pragma omp parallel default(none) shared(MaxVal)
01803 {
01804       double localMaxVal = MaxVal;
01805 #pragma omp for
01806       for (int j=0; j< myLength; j++) localMaxVal = EPETRA_MAX(localMaxVal,from[j]);
01807 #pragma omp critical
01808       {
01809       MaxVal = EPETRA_MAX(MaxVal,localMaxVal);
01810       }
01811 }
01812       DoubleTemp_[i] = MaxVal;
01813 #else
01814       for (int j=0; j< myLength; j++) MaxVal = EPETRA_MAX(MaxVal,from[j]);
01815       DoubleTemp_[i] = MaxVal;
01816 #endif
01817     }
01818 
01819   if (myLength > 0) {
01820     for(int i=0; i<NumVectors_; ++i) {
01821       Result[i] = DoubleTemp_[i];
01822     }
01823   }
01824 
01825   //If myLength == 0 and Comm_->NumProc() == 1, then Result has
01826   //not been referenced. Also, if vector contents are uninitialized
01827   //then Result contents are not well defined...
01828 
01829   if (Comm_->NumProc() == 1  || !DistributedGlobal()) return(ierr);
01830 
01831   //We're going to use MPI_Allgather to gather every proc's local-
01832   //max values onto every other proc. We'll use the last position
01833   //of the DoubleTemp_ array to indicate whether this proc has
01834   //valid data that should be considered by other procs when forming
01835   //the global-max results.
01836 
01837   if (myLength == 0) DoubleTemp_[NumVectors_] = 0.0;
01838   else DoubleTemp_[NumVectors_] = 1.0;
01839 
01840   //Now proceed to handle the parallel case. We'll gather local-max
01841   //values from other procs and form a global-max. If any processor
01842   //has myLength>0, we'll end up with a valid result.
01843 
01844 #ifdef EPETRA_MPI
01845   const Epetra_MpiComm* epetrampicomm =
01846     dynamic_cast<const Epetra_MpiComm*>(Comm_);
01847   if (!epetrampicomm) {
01848     return(-2);
01849   }
01850 
01851   MPI_Comm mpicomm = epetrampicomm->GetMpiComm();
01852   int numProcs = epetrampicomm->NumProc();
01853   double* dwork = new double[numProcs*(NumVectors_+1)];
01854 
01855   MPI_Allgather(DoubleTemp_, NumVectors_+1, MPI_DOUBLE,
01856                 dwork, NumVectors_+1, MPI_DOUBLE, mpicomm);
01857 
01858   //if myLength==0, then our Result array currently contains
01859   //-Epetra_MaxDouble from the local-max calculations above. In this
01860   //case we'll overwrite our Result array with values from the first
01861   //processor that sent valid data.
01862   bool overwrite = myLength == 0 ? true : false;
01863 
01864   int myPID = epetrampicomm->MyPID();
01865   double* dwork_ptr = dwork;
01866 
01867   for(int j=0; j<numProcs; ++j) {
01868 
01869     //skip data from self, and skip data from
01870     //procs with DoubleTemp_[NumVectors_] == 0.0.
01871     if (j == myPID || dwork_ptr[NumVectors_] < 0.5) {
01872       dwork_ptr += NumVectors_+1;
01873       continue;
01874     }
01875 
01876     for(int i=0; i<NumVectors_; ++i) {
01877       double val = dwork_ptr[i];
01878 
01879       //Set val into our Result array if overwrite is true (see above),
01880       //or if val is larger than our current Result[i].
01881       if (overwrite || (Result[i] < val)) Result[i] = val;
01882     }
01883 
01884     //Now set overwrite to false so that we'll do the right thing
01885     //when processing data from subsequent processors.
01886     if (overwrite) overwrite = false;
01887 
01888     dwork_ptr += NumVectors_+1;
01889   }
01890 
01891   delete [] dwork;
01892 #endif
01893 
01894   // UpdateFlops(0);  Strictly speaking there are not FLOPS in this routine
01895 
01896   return(ierr);
01897 }
01898 
01899 //=========================================================================
01900 int  Epetra_MultiVector::MeanValue (double* Result) const {
01901 
01902   // Mean value of each vector in MultiVector
01903 
01904   const int myLength = MyLength_;
01905 
01906   if (!Map().UniqueGIDs()) {EPETRA_CHK_ERR(-1);}
01907 
01908   double fGlobalLength = 1.0/EPETRA_MAX((double) GlobalLength_, 1.0);
01909 
01910 
01911   UpdateDoubleTemp();
01912 
01913   for (int i=0; i < NumVectors_; i++)
01914     {
01915       double sum = 0.0;
01916       const double * const from = Pointers_[i];
01917 #ifdef EPETRA_HAVE_OMP_NONASSOCIATIVE
01918 #pragma omp parallel for reduction (+:sum) default(none)
01919 #endif
01920       for (int j=0; j < myLength; j++) sum += from[j];
01921       DoubleTemp_[i] = sum;
01922     }
01923   if (DistributedGlobal())
01924     Comm_->SumAll(DoubleTemp_, Result, NumVectors_);
01925   else
01926     for (int i=0; i< NumVectors_; ++i) Result[i] = DoubleTemp_[i];
01927 
01928   for (int i=0; i < NumVectors_; i++) Result[i] = Result[i]*fGlobalLength;
01929 
01930   UpdateFlops(GlobalLength_*NumVectors_);
01931 
01932   return(0);
01933 }
01934 
01935   //=========================================================================
01936   int  Epetra_MultiVector::Multiply (char TransA, char TransB, double ScalarAB,
01937           const Epetra_MultiVector& A,
01938           const Epetra_MultiVector& B,
01939           double ScalarThis ) {
01940 
01941     // This routine performs a variety of matrix-matrix multiply operations, interpreting
01942     // the Epetra_MultiVector (this-aka C , A and B) as 2D matrices.  Variations are due to
01943     // the fact that A, B and C can be local replicated or global distributed
01944     // Epetra_MultiVectors and that we may or may not operate with the transpose of
01945     // A and B.  Possible cases are:
01946 
01947     //                                       Num
01948     //     OPERATIONS                        case  Notes
01949     // 1) C(local) = A^X(local) * B^X(local)  4   (X=Trans or Not, No comm needed)
01950     // 2) C(local) = A^T(distr) * B  (distr)  1   (2D dot product, replicate C)
01951     // 3) C(distr) = A  (distr) * B^X(local)  2   (2D vector update, no comm needed)
01952 
01953     // Note that the following operations are not meaningful for
01954     // 1D distributions:
01955 
01956     // 1) C(local) = A^T(distr) * B^T(distr)  1
01957     // 2) C(local) = A  (distr) * B^X(distr)  2
01958     // 3) C(distr) = A^X(local) * B^X(local)  4
01959     // 4) C(distr) = A^X(local) * B^X(distr)  4
01960     // 5) C(distr) = A^T(distr) * B^X(local)  2
01961     // 6) C(local) = A^X(distr) * B^X(local)  4
01962     // 7) C(distr) = A^X(distr) * B^X(local)  4
01963     // 8) C(local) = A^X(local) * B^X(distr)  4
01964 
01965     // Total of 32 case (2^5).
01966 
01967 
01968     //if (!ConstantStride_    ||
01971 
01972     // Check for compatible dimensions
01973 
01974     int A_nrows = (TransA=='T') ? A.NumVectors() : A.MyLength();
01975     int A_ncols = (TransA=='T') ? A.MyLength() : A.NumVectors();
01976     int B_nrows = (TransB=='T') ? B.NumVectors() : B.MyLength();
01977     int B_ncols = (TransB=='T') ? B.MyLength() : B.NumVectors();
01978 
01979     double Scalar_local = ScalarThis; // local copy of Scalar
01980   const int myLength = MyLength_;
01981 
01982     if( myLength      != A_nrows     ||   // RAB: 2002/01/25: Minor reformat to allow
01983     A_ncols        != B_nrows     ||   //   setting breakpoint at error return.
01984     NumVectors_    != B_ncols  )
01985     EPETRA_CHK_ERR(-2); // Return error
01986 
01987     bool A_is_local = (!A.DistributedGlobal());
01988     bool B_is_local = (!B.DistributedGlobal());
01989     bool C_is_local = (!DistributedGlobal_);
01990     bool Case1 = ( A_is_local &&  B_is_local &&  C_is_local);  // Case 1 above
01991     bool Case2 = (!A_is_local && !B_is_local &&  C_is_local && TransA=='T' );// Case 2
01992     bool Case3 = (!A_is_local &&  B_is_local && !C_is_local && TransA=='N');// Case 3
01993 
01994     if (Case2 && (!A.Map().UniqueGIDs() || !B.Map().UniqueGIDs())) {EPETRA_CHK_ERR(-4);}
01995     if (Case3 && (!A.Map().UniqueGIDs() || !Map().UniqueGIDs())) {EPETRA_CHK_ERR(-5);}
01996 
01997     // Test for meaningful cases
01998 
01999     if (Case1 || Case2 || Case3)
02000       {
02001   if (ScalarThis!=0.0 && Case2)
02002     {
02003       const int MyPID = Comm_->MyPID();
02004       if (MyPID!=0) Scalar_local = 0.0;
02005     }
02006 
02007         // Check if A, B, C have constant stride, if not then make temp copy (strided)
02008 
02009         Epetra_MultiVector * A_tmp, * B_tmp, *C_tmp;
02010         if (!ConstantStride_) C_tmp = new Epetra_MultiVector(*this);
02011         else C_tmp = this;
02012 
02013         if (!A.ConstantStride()) A_tmp = new Epetra_MultiVector(A);
02014         else A_tmp = (Epetra_MultiVector *) &A;
02015 
02016         if (!B.ConstantStride()) B_tmp = new Epetra_MultiVector(B);
02017         else B_tmp = (Epetra_MultiVector *) &B;
02018 
02019 
02020   int m = myLength;
02021   int n = NumVectors_;
02022   int k = A_ncols;
02023   int lda = EPETRA_MAX(A_tmp->Stride(),1); // The reference BLAS implementation requires lda, ldb, ldc > 0 even if m, n, or k = 0
02024   int ldb = EPETRA_MAX(B_tmp->Stride(),1);
02025   int ldc = EPETRA_MAX(C_tmp->Stride(),1);
02026   double *Ap = A_tmp->Values();
02027   double *Bp = B_tmp->Values();
02028   double *Cp = C_tmp->Values();
02029 
02030   GEMM(TransA, TransB,  m, n, k, ScalarAB,
02031        Ap, lda, Bp, ldb, Scalar_local, Cp, ldc);
02032 
02033   // FLOP Counts
02034   //                                       Num
02035   //     OPERATIONS                        case  Notes
02036   // 1) C(local) = A^X(local) * B^X(local)  4   (X=Trans or Not, No comm needed)
02037   // 2) C(local) = A^T(distr) * B  (distr)  1   (2D dot product, replicate C)
02038   // 3) C(distr) = A  (distr) * B^X(local)  2   (2D vector update, no comm needed)
02039 
02040   // For Case 1 we only count the local operations, since we are interested in serial
02041   // cost.  Computation on other processors is redundant.
02042   if (Case1)
02043     {
02044       UpdateFlops(2*m*n*k);
02045       if (ScalarAB!=1.0) UpdateFlops(m*n);
02046       if (ScalarThis==1.0) UpdateFlops(m*n); else if (ScalarThis!=0.0) UpdateFlops(2*m*n);
02047     }
02048   else if (Case2)
02049     {
02050       UpdateFlops(2*m*n*A.GlobalLength64());
02051       if (ScalarAB!=1.0) UpdateFlops(m*n);
02052       if (ScalarThis==1.0) UpdateFlops(m*n); else if (ScalarThis!=0.0) UpdateFlops(2*m*n);
02053     }
02054   else
02055     {
02056       UpdateFlops(2*GlobalLength_*n*k);
02057       if (ScalarAB!=1.0) UpdateFlops(GlobalLength_*n);
02058       if (ScalarThis==1.0) UpdateFlops(GlobalLength_*n);
02059       else if (ScalarThis!=0.0) UpdateFlops(2*GlobalLength_*n);
02060     }
02061 
02062   // If A or B were not strided, dispose of extra copies.
02063   if (!A.ConstantStride()) delete A_tmp;
02064   if (!B.ConstantStride()) delete B_tmp;
02065 
02066   // If C was not strided, copy from strided version and delete
02067   if (!ConstantStride_)
02068     {
02069       C_tmp->ExtractCopy(Pointers_);
02070       delete C_tmp;
02071     }
02072 
02073   // If Case 2 then sum up C and distribute it to all processors.
02074 
02075   if (Case2) {EPETRA_CHK_ERR(Reduce());}
02076 
02077   return(0);
02078 
02079       }
02080     else {EPETRA_CHK_ERR(-3)}; // Return error: not a supported operation
02081 
02082   return(0);
02083   }
02084 
02085 
02086 //=========================================================================
02087 int Epetra_MultiVector::Multiply(double ScalarAB, const Epetra_MultiVector& A, const Epetra_MultiVector& B,
02088            double ScalarThis) {
02089 
02090 
02091   // Hadamard product of two MultiVectors:
02092   // this = ScalarThis * this + ScalarAB * A * B (element-wise)
02093 
02094   if (ScalarAB==0.0) {
02095     EPETRA_CHK_ERR(Scale(ScalarThis));
02096     return(0);
02097   }
02098   const int myLength = MyLength_;
02099 
02100   if (A.NumVectors() != 1 && A.NumVectors() != B.NumVectors()) EPETRA_CHK_ERR(-1); // A must have one column or be the same as B.
02101   if (NumVectors_ != B.NumVectors()) EPETRA_CHK_ERR(-2);
02102   if (myLength != A.MyLength() || myLength != B.MyLength()) EPETRA_CHK_ERR(-3);
02103 
02104   double **A_Pointers = (double**)A.Pointers();
02105   double **B_Pointers = (double**)B.Pointers();
02106 
02107   int IncA = 1;
02108   if (A.NumVectors() == 1 ) IncA = 0;
02109 
02110     if (ScalarThis==0.0) {
02111       if (ScalarAB==1.0)
02112   {
02113     for (int i = 0; i < NumVectors_; i++) {
02114       const double * const Aptr = A_Pointers[i*IncA];
02115             const double * const Bptr = B_Pointers[i];
02116             double * const to = Pointers_[i];
02117 #ifdef EPETRA_HAVE_OMP
02118 #pragma omp parallel for default(none)
02119 #endif
02120       for (int j = 0; j < myLength; j++) {
02121               to[j] =  Aptr[j] * Bptr[j];
02122             }
02123     }
02124     UpdateFlops(GlobalLength_*NumVectors_);
02125   }
02126       else
02127   {
02128     for (int i = 0; i < NumVectors_; i++) {
02129       const double * const Aptr = A_Pointers[i*IncA];
02130             const double * const Bptr = B_Pointers[i];
02131             double * const to = Pointers_[i];
02132 #ifdef EPETRA_HAVE_OMP
02133 #pragma omp parallel for default(none) shared(ScalarAB)
02134 #endif
02135       for (int j = 0; j < myLength; j++) {
02136               to[j] = ScalarAB * Aptr[j] * Bptr[j];
02137             }
02138     }
02139     UpdateFlops(2*GlobalLength_*NumVectors_);
02140   }
02141     }
02142     else if (ScalarThis==1.0) {
02143       if (ScalarAB==1.0)
02144   {
02145     for (int i = 0; i < NumVectors_; i++) {
02146       const double * const Aptr = A_Pointers[i*IncA];
02147             const double * const Bptr = B_Pointers[i];
02148             double * const to = Pointers_[i];
02149 #ifdef EPETRA_HAVE_OMP
02150 #pragma omp parallel for default(none)
02151 #endif
02152       for (int j = 0; j < myLength; j++) {
02153               to[j] +=  Aptr[j] * Bptr[j];
02154             }
02155     }
02156     UpdateFlops(2*GlobalLength_*NumVectors_);
02157   }
02158       else {
02159     for (int i = 0; i < NumVectors_; i++) {
02160       const double * const Aptr = A_Pointers[i*IncA];
02161             const double * const Bptr = B_Pointers[i];
02162             double * const to = Pointers_[i];
02163 #ifdef EPETRA_HAVE_OMP
02164 #pragma omp parallel for default(none) shared(ScalarAB)
02165 #endif
02166       for (int j = 0; j < myLength; j++) {
02167               to[j] += ScalarAB * Aptr[j] * Bptr[j];
02168             }
02169           }
02170           UpdateFlops(3*GlobalLength_*NumVectors_);
02171       }
02172     }
02173     else { // if (ScalarThis!=1.0 && ScalarThis !=0 )
02174       if (ScalarAB==1.0)
02175   {
02176     for (int i = 0; i < NumVectors_; i++) {
02177       const double * const Aptr = A_Pointers[i*IncA];
02178             const double * const Bptr = B_Pointers[i];
02179             double * const to = Pointers_[i];
02180 #ifdef EPETRA_HAVE_OMP
02181 #pragma omp parallel for default(none) shared(ScalarThis)
02182 #endif
02183       for (int j = 0; j < myLength; j++) {
02184               to[j] =  ScalarThis * to[j] +
02185                   Aptr[j] * Bptr[j];
02186             }
02187     }
02188     UpdateFlops(3*GlobalLength_*NumVectors_);
02189   }
02190       else
02191   {
02192     for (int i = 0; i < NumVectors_; i++) {
02193       const double * const Aptr = A_Pointers[i*IncA];
02194             const double * const Bptr = B_Pointers[i];
02195             double * const to = Pointers_[i];
02196 #ifdef EPETRA_HAVE_OMP
02197 #pragma omp parallel for default(none) shared(ScalarThis,ScalarAB)
02198 #endif
02199       for (int j = 0; j < myLength; j++) {
02200               to[j] = ScalarThis * to[j] +
02201                      ScalarAB * Aptr[j] * Bptr[j];
02202             }
02203     }
02204     UpdateFlops(4*GlobalLength_*NumVectors_);
02205   }
02206     }
02207   return(0);
02208 }
02209 //=========================================================================
02210 int Epetra_MultiVector::ReciprocalMultiply(double ScalarAB, const Epetra_MultiVector& A, const Epetra_MultiVector& B,
02211            double ScalarThis) {
02212 
02213 
02214   // Hadamard product of two MultiVectors:
02215   // this = ScalarThis * this + ScalarAB * B / A (element-wise)
02216 
02217   if (ScalarAB==0.0) {
02218     EPETRA_CHK_ERR(Scale(ScalarThis));
02219     return(0);
02220   }
02221   const int myLength = MyLength_;
02222 
02223   if (A.NumVectors() != 1 && A.NumVectors() != B.NumVectors()) EPETRA_CHK_ERR(-1); // A must have one column or be the same as B.
02224   if (NumVectors_ != B.NumVectors()) EPETRA_CHK_ERR(-2);
02225   if (myLength != A.MyLength() || myLength != B.MyLength()) EPETRA_CHK_ERR(-3);
02226 
02227   double **A_Pointers = (double**)A.Pointers();
02228   double **B_Pointers = (double**)B.Pointers();
02229 
02230   int IncA = 1;
02231   if (A.NumVectors() == 1 ) IncA = 0;
02232 
02233     if (ScalarThis==0.0) {
02234       if (ScalarAB==1.0)
02235   {
02236     for (int i = 0; i < NumVectors_; i++) {
02237       const double * const Aptr = A_Pointers[i*IncA];
02238             const double * const Bptr = B_Pointers[i];
02239             double * const to = Pointers_[i];
02240 #ifdef EPETRA_HAVE_OMP
02241 #pragma omp parallel for default(none)
02242 #endif
02243       for (int j = 0; j < myLength; j++) {
02244               to[j] = Bptr[j] / Aptr[j];
02245             }
02246     }
02247     UpdateFlops(GlobalLength_*NumVectors_);
02248   }
02249       else
02250   {
02251     for (int i = 0; i < NumVectors_; i++) {
02252       const double * const Aptr = A_Pointers[i*IncA];
02253             const double * const Bptr = B_Pointers[i];
02254             double * const to = Pointers_[i];
02255 #ifdef EPETRA_HAVE_OMP
02256 #pragma omp parallel for default(none) shared(ScalarAB)
02257 #endif
02258       for (int j = 0; j < myLength; j++) {
02259               to[j] = ScalarAB * Bptr[j] / Aptr[j];
02260             }
02261     }
02262     UpdateFlops(2*GlobalLength_*NumVectors_);
02263   }
02264     }
02265     else if (ScalarThis==1.0) {
02266       if (ScalarAB==1.0)
02267   {
02268     for (int i = 0; i < NumVectors_; i++) {
02269       const double * const Aptr = A_Pointers[i*IncA];
02270             const double * const Bptr = B_Pointers[i];
02271             double * const to = Pointers_[i];
02272 #ifdef EPETRA_HAVE_OMP
02273 #pragma omp parallel for default(none)
02274 #endif
02275       for (int j = 0; j < myLength; j++) {
02276               to[j] +=  Bptr[j] / Aptr[j];
02277             }
02278     }
02279     UpdateFlops(2*GlobalLength_*NumVectors_);
02280   }
02281       else
02282   {
02283     for (int i = 0; i < NumVectors_; i++) {
02284       const double * const Aptr = A_Pointers[i*IncA];
02285             const double * const Bptr = B_Pointers[i];
02286             double * const to = Pointers_[i];
02287 #ifdef EPETRA_HAVE_OMP
02288 #pragma omp parallel for default(none) shared(ScalarAB)
02289 #endif
02290       for (int j = 0; j < myLength; j++) {
02291               to[j] += ScalarAB * Bptr[j] / Aptr[j];
02292             }
02293     }
02294     UpdateFlops(3*GlobalLength_*NumVectors_);
02295         }
02296     }
02297     else { // if (ScalarThis!=1.0 && ScalarThis !=0 )
02298       if (ScalarAB==1.0)
02299   {
02300     for (int i = 0; i < NumVectors_; i++) {
02301       const double * const Aptr = A_Pointers[i*IncA];
02302             const double * const Bptr = B_Pointers[i];
02303             double * const to = Pointers_[i];
02304 #ifdef EPETRA_HAVE_OMP
02305 #pragma omp parallel for default(none) shared(ScalarThis)
02306 #endif
02307       for (int j = 0; j < myLength; j++) {
02308               to[j] =  ScalarThis * to[j] +
02309                    Bptr[j] / Aptr[j];
02310             }
02311     }
02312     UpdateFlops(3*GlobalLength_*NumVectors_);
02313   }
02314       else {
02315     for (int i = 0; i < NumVectors_; i++) {
02316       const double * const Aptr = A_Pointers[i*IncA];
02317             const double * const Bptr = B_Pointers[i];
02318             double * const to = Pointers_[i];
02319 #ifdef EPETRA_HAVE_OMP
02320 #pragma omp parallel for default(none) shared(ScalarAB,ScalarThis)
02321 #endif
02322       for (int j = 0; j < myLength; j++) {
02323               to[j] = ScalarThis * to[j] + ScalarAB *
02324                 Bptr[j] / Aptr[j];
02325             }
02326     }
02327           UpdateFlops(4*GlobalLength_*NumVectors_);
02328       }
02329     }
02330   return(0);
02331 }
02332 
02333 
02334 //=======================================================================
02335 Epetra_Vector *& Epetra_MultiVector::operator () (int index)  {
02336 
02337   //  Epetra_MultiVector::operator () --- return non-const reference
02338 
02339   if (index < 0 || index >=NumVectors_)
02340     throw ReportError("Vector index = " + toString(index) + "is out of range. Number of Vectors = " + toString(NumVectors_), -1);
02341 
02342   UpdateVectors();
02343 
02344   // Create a new Epetra_Vector that is a view of ith vector, if not already present
02345   if (Vectors_[index]==0)
02346     Vectors_[index] = new Epetra_Vector(View, Map(), Pointers_[index]);
02347   return(Vectors_[index]);
02348 }
02349 
02350 //=======================================================================
02351 const Epetra_Vector *& Epetra_MultiVector::operator () (int index) const {
02352 
02353   //  Epetra_MultiVector::operator () --- return non-const reference
02354 
02355   if (index < 0 || index >=NumVectors_)
02356     throw ReportError("Vector index = " + toString(index) + "is out of range. Number of Vectors = " + toString(NumVectors_), -1);
02357 
02358   UpdateVectors();
02359 
02360   if (Vectors_[index]==0)
02361     Vectors_[index] = new Epetra_Vector(View, Map(), Pointers_[index]);
02362 
02363   const Epetra_Vector * & temp = (const Epetra_Vector * &) (Vectors_[index]);
02364   return(temp);
02365 }
02366 
02367 //========================================================================
02368 Epetra_MultiVector& Epetra_MultiVector::operator = (const Epetra_MultiVector& Source) {
02369 
02370   // Check for special case of this=Source
02371   if (this != &Source) Assign(Source);
02372 
02373   return(*this);
02374 }
02375 
02376 //=========================================================================
02377 void Epetra_MultiVector::Assign(const Epetra_MultiVector& A) {
02378 
02379   const int myLength = MyLength_;
02380   if (NumVectors_ != A.NumVectors())
02381     throw ReportError("Number of vectors incompatible in Assign.  The this MultiVector has NumVectors = " + toString(NumVectors_)
02382           + ".  The A MultiVector has NumVectors = " + toString(A.NumVectors()), -3);
02383   if (myLength != A.MyLength())
02384     throw ReportError("Length of MultiVectors incompatible in Assign.  The this MultiVector has MyLength = " + toString(myLength)
02385           + ".  The A MultiVector has MyLength = " + toString(A.MyLength()), -4);
02386 
02387   double ** A_Pointers = A.Pointers();
02388   for (int i = 0; i< NumVectors_; i++) {
02389       double * const to = Pointers_[i];
02390       const double * const from = A_Pointers[i];
02391 #ifdef EPETRA_HAVE_OMP
02392 #pragma omp parallel for default(none)
02393 #endif
02394       for (int j=0; j<myLength; j++) to[j] = from[j];
02395     }
02396     return;
02397   }
02398 
02399   //=========================================================================
02400   int  Epetra_MultiVector::Reduce() {
02401 
02402     // Global reduction on each entry of a Replicated Local MultiVector
02403   const int myLength = MyLength_;
02404     double * source = 0;
02405       if (myLength>0) source = new double[myLength*NumVectors_];
02406     double * target = 0;
02407     bool packed = (ConstantStride_ && (Stride_==myLength));
02408     if (packed) {
02409        for (int i=0; i<myLength*NumVectors_; i++) source[i] = Values_[i];
02410        target = Values_;
02411     }
02412     else {
02413        double * tmp1 = source;
02414        for (int i = 0; i < NumVectors_; i++) {
02415          double * tmp2 = Pointers_[i];
02416          for (int j=0; j< myLength; j++) *tmp1++ = *tmp2++;
02417        }
02418        if (myLength>0) target = new double[myLength*NumVectors_];
02419     }
02420 
02421     Comm_->SumAll(source, target, myLength*NumVectors_);
02422     if (myLength>0) delete [] source;
02423     if (!packed) {
02424        double * tmp2 = target;
02425        for (int i = 0; i < NumVectors_; i++) {
02426          double * tmp1 = Pointers_[i];
02427          for (int j=0; j< myLength; j++) *tmp1++ = *tmp2++;
02428        }
02429        if (myLength>0) delete [] target;
02430     }
02431     // UpdateFlops(0);  No serial Flops in this function
02432     return(0);
02433   }
02434 //=======================================================================
02435 int Epetra_MultiVector::ResetView(double ** ArrayOfPointers) {
02436 
02437   if (!UserAllocated_) {
02438     EPETRA_CHK_ERR(-1); // Can't reset view if multivector was not allocated as a view
02439   }
02440 
02441   for (int i = 0; i< NumVectors_; i++) Pointers_[i] = ArrayOfPointers[i];
02442   DoView();
02443 
02444   return(0);
02445   }
02446 //=======================================================================
02447 void Epetra_MultiVector::Print(std::ostream& os) const {
02448   int MyPID = Map().Comm().MyPID();
02449   int NumProc = Map().Comm().NumProc();
02450 
02451   for (int iproc=0; iproc < NumProc; iproc++) {
02452     if (MyPID==iproc) {
02453       int NumVectors1 = NumVectors();
02454       int NumMyElements1 =Map(). NumMyElements();
02455       int MaxElementSize1 = Map().MaxElementSize();
02456       int * FirstPointInElementList1 = NULL;
02457       if (MaxElementSize1!=1) FirstPointInElementList1 = Map().FirstPointInElementList();
02458       double ** A_Pointers = Pointers();
02459 
02460       if (MyPID==0) {
02461   os.width(8);
02462   os <<  "     MyPID"; os << "    ";
02463   os.width(12);
02464   if (MaxElementSize1==1)
02465     os <<  "GID  ";
02466   else
02467     os <<  "     GID/Point";
02468   for (int j = 0; j < NumVectors1 ; j++)
02469     {
02470       os.width(20);
02471       os <<  "Value  ";
02472     }
02473   os << std::endl;
02474       }
02475       for (int i=0; i < NumMyElements1; i++) {
02476   for (int ii=0; ii< Map().ElementSize(i); ii++) {
02477        int iii;
02478     os.width(10);
02479     os <<  MyPID; os << "    ";
02480     os.width(10);
02481     if (MaxElementSize1==1) {
02482       if(Map().GlobalIndicesInt())
02483       {
02484 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
02485       int * MyGlobalElements1 = Map().MyGlobalElements();
02486       os << MyGlobalElements1[i] << "    ";
02487 #else
02488             throw ReportError("Epetra_MultiVector::Print: ERROR, GlobalIndicesInt but no API for it.",-1);
02489 #endif
02490       }
02491       else if(Map().GlobalIndicesLongLong())
02492       {
02493 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
02494       long long * MyGlobalElements1 = Map().MyGlobalElements64();
02495       os << MyGlobalElements1[i] << "    ";
02496 #else
02497             throw ReportError("Epetra_MultiVector::Print: ERROR, GlobalIndicesLongLong but no API for it.",-1);
02498 #endif
02499       }
02500       else
02501         throw ReportError("Epetra_MultiVector::Print ERROR, Don't know map global index type.",-1);
02502 
02503        iii = i;
02504        }
02505     else {
02506       if(Map().GlobalIndicesInt())
02507       {
02508 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
02509       int * MyGlobalElements1 = Map().MyGlobalElements();
02510       os <<  MyGlobalElements1[i]<< "/" << ii << "    ";
02511 #else
02512             throw ReportError("Epetra_MultiVector::Print: ERROR, GlobalIndicesInt but no API for it.",-1);
02513 #endif
02514       }
02515       else if(Map().GlobalIndicesLongLong())
02516       {
02517 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
02518       long long * MyGlobalElements1 = Map().MyGlobalElements64();
02519       os <<  MyGlobalElements1[i]<< "/" << ii << "    ";
02520 #else
02521             throw ReportError("Epetra_MultiVector::Print: ERROR, GlobalIndicesLongLong but no API for it.",-1);
02522 #endif
02523       }
02524       else
02525         throw ReportError("Epetra_MultiVector::Print ERROR, Don't know map global index type.",-1);
02526 
02527          iii = FirstPointInElementList1[i]+ii;
02528        }
02529     for (int j = 0; j < NumVectors1 ; j++)
02530       {
02531         os.width(20);
02532         os <<  A_Pointers[j][iii];
02533       }
02534     os << std::endl;
02535   }
02536       }
02537       os << std::flush;
02538     }
02539 
02540     // Do a few global ops to give I/O a chance to complete
02541     Map().Comm().Barrier();
02542     Map().Comm().Barrier();
02543     Map().Comm().Barrier();
02544   }
02545   return;
02546 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines