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