Epetra_MultiVector.cpp

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

Generated on Thu Sep 18 12:37:58 2008 for Epetra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1