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++) {
00788     To[0][ImportLIDs[j]] = EPETRA_MAX( To[0][ImportLIDs[j]],std::abs(*ptr));
00789     ptr++;
00790   }
00791       // Note:  The following form of averaging is not a true average if more that one value is combined.
00792       //        This might be an issue in the future, but we leave this way for now.
00793       else if(CombineMode==Average)
00794   for (j=0; j<NumImportIDs; j++) {To[0][ImportLIDs[j]] += *ptr++; To[0][ImportLIDs[j]] *= 0.5;}
00795     }
00796 
00797     else {  // NumVectors>1
00798 
00799       if (CombineMode==Add) {
00800   for (j=0; j<NumImportIDs; j++) {
00801     jj = ImportLIDs[j];
00802     for (i=0; i<NumVectors; i++)
00803       To[i][jj] += *ptr++; // Add to existing value
00804   }
00805       }
00806       else if(CombineMode==Insert) {
00807   for (j=0; j<NumImportIDs; j++) {
00808     jj = ImportLIDs[j];
00809     for (i=0; i<NumVectors; i++)
00810       To[i][jj] = *ptr++;
00811   }
00812       }
00813       else if(CombineMode==InsertAdd) {
00814   for (j=0; j<NumImportIDs; j++) {
00815     jj = ImportLIDs[j];
00816     for (i=0; i<NumVectors; i++)
00817       To[i][jj] = 0.0;
00818   }
00819   for (j=0; j<NumImportIDs; j++) {
00820     jj = ImportLIDs[j];
00821     for (i=0; i<NumVectors; i++)
00822       To[i][jj] += *ptr++;
00823   }
00824       }
00825       else if(CombineMode==AbsMax) {
00826         for (j=0; j<NumImportIDs; j++) {
00827           jj = ImportLIDs[j];
00828           for (i=0; i<NumVectors; i++) {
00829             To[i][jj] = EPETRA_MAX( To[i][jj], std::abs(*ptr) );
00830       ptr++;
00831     }
00832         }
00833       }
00834       // Note:  The following form of averaging is not a true average if more that one value is combined.
00835       //        This might be an issue in the future, but we leave this way for now.
00836       else if(CombineMode==Average) {
00837   for (j=0; j<NumImportIDs; j++) {
00838     jj = ImportLIDs[j];
00839     for (i=0; i<NumVectors; i++)
00840       { To[i][jj] += *ptr++;  To[i][jj] *= 0.5;}
00841   }
00842       }
00843     }
00844   }
00845 
00846   // constant element size case
00847 
00848   else if (ConstantElementSize) {
00849    
00850     if (CombineMode==Add) {
00851       for (j=0; j<NumImportIDs; j++) {
00852   jj = MaxElementSize*ImportLIDs[j];
00853   for (i=0; i<NumVectors; i++)
00854     for (k=0; k<MaxElementSize; k++)
00855       To[i][jj+k] += *ptr++; // Add to existing value
00856       }
00857     }
00858     else if(CombineMode==Insert) {
00859       for (j=0; j<NumImportIDs; j++) {
00860   jj = MaxElementSize*ImportLIDs[j];
00861   for (i=0; i<NumVectors; i++)
00862     for (k=0; k<MaxElementSize; k++)
00863       To[i][jj+k] = *ptr++;
00864       }
00865     }
00866     else if(CombineMode==InsertAdd) {
00867       for (j=0; j<NumImportIDs; j++) {
00868   jj = MaxElementSize*ImportLIDs[j];
00869   for (i=0; i<NumVectors; i++)
00870     for (k=0; k<MaxElementSize; k++)
00871       To[i][jj+k] = 0.0;
00872       }
00873       for (j=0; j<NumImportIDs; j++) {
00874   jj = MaxElementSize*ImportLIDs[j];
00875   for (i=0; i<NumVectors; i++)
00876     for (k=0; k<MaxElementSize; k++)
00877       To[i][jj+k] += *ptr++;
00878       }
00879     }
00880     else if(CombineMode==AbsMax) {
00881       for (j=0; j<NumImportIDs; j++) {
00882   jj = MaxElementSize*ImportLIDs[j];
00883   for (i=0; i<NumVectors; i++)
00884     for (k=0; k<MaxElementSize; k++) {
00885       To[i][jj+k] = EPETRA_MAX( To[i][jj+k], std::abs(*ptr) );
00886       ptr++;
00887     }
00888       }
00889     }
00890     // Note:  The following form of averaging is not a true average if more that one value is combined.
00891     //        This might be an issue in the future, but we leave this way for now.
00892     else if(CombineMode==Average) {
00893       for (j=0; j<NumImportIDs; j++) {
00894   jj = MaxElementSize*ImportLIDs[j];
00895   for (i=0; i<NumVectors; i++)
00896     for (k=0; k<MaxElementSize; k++)
00897       { To[i][jj+k] += *ptr++; To[i][jj+k] *= 0.5;}
00898       }
00899     }
00900   }
00901     
00902   // variable element size case
00903 
00904   else {
00905     int thisSizeOfPacket = NumVectors*MaxElementSize;
00906 
00907     if (CombineMode==Add) {
00908       for (j=0; j<NumImportIDs; j++) {
00909   ptr = (double *) Imports + j*thisSizeOfPacket;
00910   jj = ToFirstPointInElementList[ImportLIDs[j]];
00911   int ElementSize = ToElementSizeList[ImportLIDs[j]];
00912   for (i=0; i<NumVectors; i++)
00913     for (k=0; k<ElementSize; k++)
00914       To[i][jj+k] += *ptr++; // Add to existing value
00915       }
00916     }
00917     else  if(CombineMode==Insert){
00918       for (j=0; j<NumImportIDs; j++) {
00919   ptr = (double *) Imports + j*thisSizeOfPacket;
00920   jj = ToFirstPointInElementList[ImportLIDs[j]];
00921   int ElementSize = ToElementSizeList[ImportLIDs[j]];
00922   for (i=0; i<NumVectors; i++)
00923     for (k=0; k<ElementSize; k++)
00924       To[i][jj+k] = *ptr++;
00925       }
00926     }
00927     else  if(CombineMode==InsertAdd){
00928       for (j=0; j<NumImportIDs; j++) {
00929   ptr = (double *) Imports + j*thisSizeOfPacket;
00930   jj = ToFirstPointInElementList[ImportLIDs[j]];
00931   int ElementSize = ToElementSizeList[ImportLIDs[j]];
00932   for (i=0; i<NumVectors; i++)
00933     for (k=0; k<ElementSize; k++)
00934       To[i][jj+k] = 0.0;
00935       }
00936       for (j=0; j<NumImportIDs; j++) {
00937   ptr = (double *) Imports + j*thisSizeOfPacket;
00938   jj = ToFirstPointInElementList[ImportLIDs[j]];
00939   int ElementSize = ToElementSizeList[ImportLIDs[j]];
00940   for (i=0; i<NumVectors; i++)
00941     for (k=0; k<ElementSize; k++)
00942       To[i][jj+k] += *ptr++;
00943       }
00944     }
00945     else  if(CombineMode==AbsMax){
00946       for (j=0; j<NumImportIDs; j++) {
00947   ptr = (double *) Imports + j*thisSizeOfPacket;
00948   jj = ToFirstPointInElementList[ImportLIDs[j]];
00949   int ElementSize = ToElementSizeList[ImportLIDs[j]];
00950   for (i=0; i<NumVectors; i++)
00951     for (k=0; k<ElementSize; k++) {
00952       To[i][jj+k] = EPETRA_MAX( To[i][jj+k], std::abs(*ptr));
00953       ptr++;
00954     }
00955       }
00956     }
00957     // Note:  The following form of averaging is not a true average if more that one value is combined.
00958     //        This might be an issue in the future, but we leave this way for now.
00959     else if(CombineMode==Average) {
00960       for (j=0; j<NumImportIDs; j++) {
00961   ptr = (double *) Imports + j*thisSizeOfPacket;
00962   jj = ToFirstPointInElementList[ImportLIDs[j]];
00963   int ElementSize = ToElementSizeList[ImportLIDs[j]];
00964   for (i=0; i<NumVectors; i++)
00965     for (k=0; k<ElementSize; k++)
00966       { To[i][jj+k] += *ptr++; To[i][jj+k] *= 0.5;}
00967       }
00968     }
00969   }
00970 
00971   return(0);
00972 }
00973 
00974 //=========================================================================
00975 int Epetra_MultiVector::Dot(const Epetra_MultiVector& A, double *Result) const {
00976 
00977   // Dot product of two MultiVectors 
00978 
00979   int i;
00980   if (NumVectors_ != A.NumVectors()) EPETRA_CHK_ERR(-1);
00981   if (MyLength_ != A.MyLength()) EPETRA_CHK_ERR(-2);
00982   UpdateDoubleTemp();
00983     
00984   double **A_Pointers = A.Pointers();
00985 
00986   for (i=0; i < NumVectors_; i++) DoubleTemp_[i] = DOT(MyLength_, Pointers_[i], A_Pointers[i]);
00987   
00988   Comm_->SumAll(DoubleTemp_, Result, NumVectors_);
00989   
00990   UpdateFlops(2*GlobalLength_*NumVectors_);
00991 
00992   return(0);
00993 }
00994 //=========================================================================
00995 int Epetra_MultiVector::Abs(const Epetra_MultiVector& A) {
00996 
00997   // this[i][j] = std::abs(A[i][j])
00998 
00999   int i, j;
01000   if (NumVectors_ != A.NumVectors()) EPETRA_CHK_ERR(-1);
01001   if (MyLength_ != A.MyLength()) EPETRA_CHK_ERR(-2);
01002 
01003   double **A_Pointers = A.Pointers();
01004 
01005   for (i=0; i < NumVectors_; i++) 
01006     for (j=0; j < MyLength_; j++)
01007       Pointers_[i][j] = std::abs(A_Pointers[i][j]);
01008 
01009   return(0);
01010 }
01011 //=========================================================================
01012 int Epetra_MultiVector::Reciprocal(const Epetra_MultiVector& A) {
01013 
01014   // this[i][j] = 1.0/(A[i][j])
01015 
01016   int ierr = 0;
01017   int i, j;
01018   if (NumVectors_ != A.NumVectors()) EPETRA_CHK_ERR(-1);
01019   if (MyLength_ != A.MyLength()) EPETRA_CHK_ERR(-2);
01020 
01021   double **A_Pointers = A.Pointers();
01022 
01023   for (i=0; i < NumVectors_; i++) 
01024     for (j=0; j < MyLength_; j++) {
01025       double value = A_Pointers[i][j];
01026       // Set error to 1 to signal that zero rowsum found (supercedes ierr = 2)     
01027       if (std::abs(value)<Epetra_MinDouble) {
01028   if (value==0.0) ierr = 1;
01029   else if (ierr!=1) ierr = 2;
01030   Pointers_[i][j] = EPETRA_SGN(value) * Epetra_MaxDouble;
01031       }
01032       else
01033   Pointers_[i][j] = 1.0/value;
01034     }
01035   EPETRA_CHK_ERR(ierr);
01036   return(0);
01037 }
01038   //=========================================================================
01039   int Epetra_MultiVector::Scale (double ScalarValue) {
01040 
01041     // scales a MultiVector in place by a scalar
01042   
01043     for (int i = 0; i < NumVectors_; i++)
01044       SCAL(MyLength_, ScalarValue, Pointers_[i]);
01045 
01046     UpdateFlops(GlobalLength_*NumVectors_);
01047 
01048     return(0);
01049   }
01050 
01051   //=========================================================================
01052   int Epetra_MultiVector::Scale (double ScalarA, const Epetra_MultiVector& A) {
01053 
01054     // scales a MultiVector by a scalar and put in the this:
01055     // this = ScalarA * A
01056 
01057   if (NumVectors_ != A.NumVectors()) EPETRA_CHK_ERR(-1);
01058   if (MyLength_ != A.MyLength()) EPETRA_CHK_ERR(-2);
01059 
01060     double **A_Pointers = (double**)A.Pointers();
01061     
01062     for (int i = 0; i < NumVectors_; i++)
01063       for (int j = 0; j < MyLength_; j++) Pointers_[i][j] = ScalarA * A_Pointers[i][j];
01064 
01065     UpdateFlops(GlobalLength_*NumVectors_);
01066 
01067     return(0);
01068   }
01069 
01070   //=========================================================================
01071   int Epetra_MultiVector::Update(double ScalarA, const Epetra_MultiVector& A, double ScalarThis) {
01072 
01073     int i, j;
01074 
01075     // linear combination of two MultiVectors: this = ScalarThis * this + ScalarA * A
01076 
01077   if (NumVectors_ != A.NumVectors()) EPETRA_CHK_ERR(-1);
01078   if (MyLength_ != A.MyLength()) EPETRA_CHK_ERR(-2);
01079 
01080     double **A_Pointers = (double**)A.Pointers();
01081 
01082     if (ScalarThis==0.0)
01083       {
01084   for (i = 0; i < NumVectors_; i++)
01085     for (j = 0; j < MyLength_; j++) Pointers_[i][j] = ScalarA * A_Pointers[i][j];
01086   UpdateFlops(GlobalLength_*NumVectors_);
01087       }
01088     else if (ScalarThis==1.0)
01089       {
01090   for (i = 0; i < NumVectors_; i++)
01091     for (j = 0; j < MyLength_; j++) Pointers_[i][j] = Pointers_[i][j] + ScalarA * A_Pointers[i][j];
01092   UpdateFlops(2*GlobalLength_*NumVectors_);
01093       }
01094     else if (ScalarA==1.0)
01095       {
01096   for (i = 0; i < NumVectors_; i++)
01097     for (j = 0; j < MyLength_; j++) Pointers_[i][j] = ScalarThis * Pointers_[i][j] + A_Pointers[i][j];
01098   UpdateFlops(2*GlobalLength_*NumVectors_);
01099       }
01100     else
01101       {
01102   for (i = 0; i < NumVectors_; i++)
01103     for (j = 0; j < MyLength_; j++) Pointers_[i][j] = ScalarThis * Pointers_[i][j] +
01104                               ScalarA *  A_Pointers[i][j];
01105   UpdateFlops(3*GlobalLength_*NumVectors_);
01106       }
01107 
01108     return(0);
01109   }
01110 
01111 //=========================================================================
01112 int Epetra_MultiVector::Update(double ScalarA, const Epetra_MultiVector& A, 
01113           double ScalarB, const Epetra_MultiVector& B, double ScalarThis) {
01114   
01115   int i, j;
01116   
01117   // linear combination of three MultiVectors: 
01118   // this = ScalarThis * this + ScalarA * A + ScalarB * B
01119   
01120   if (ScalarA==0.0) {
01121     EPETRA_CHK_ERR(Update(ScalarB, B, ScalarThis));
01122     return(0);
01123   }
01124   if (ScalarB==0.0) {
01125     EPETRA_CHK_ERR(Update(ScalarA, A, ScalarThis));
01126     return(0);
01127   }
01128          
01129   if (NumVectors_ != A.NumVectors() || NumVectors_ != B.NumVectors()) EPETRA_CHK_ERR(-1);
01130   if (MyLength_ != A.MyLength() || MyLength_ != B.MyLength()) EPETRA_CHK_ERR(-2);
01131   
01132     double **A_Pointers = (double**)A.Pointers();
01133     double **B_Pointers = (double**)B.Pointers();
01134 
01135     if (ScalarThis==0.0)
01136       {
01137   if (ScalarA==1.0)
01138     {
01139       for (i = 0; i < NumVectors_; i++)
01140         for (j = 0; j < MyLength_; j++) Pointers_[i][j] =           A_Pointers[i][j] + 
01141                             ScalarB * B_Pointers[i][j];
01142       UpdateFlops(2*GlobalLength_*NumVectors_);
01143     }
01144   else if (ScalarB==1.0)
01145     {
01146       for (i = 0; i < NumVectors_; i++)
01147         for (j = 0; j < MyLength_; j++) Pointers_[i][j] = ScalarA * A_Pointers[i][j] +
01148                                       B_Pointers[i][j];
01149       UpdateFlops(2*GlobalLength_*NumVectors_);
01150     }
01151   else
01152     {
01153       for (i = 0; i < NumVectors_; i++)
01154         for (j = 0; j < MyLength_; j++) Pointers_[i][j] = ScalarA * A_Pointers[i][j] + 
01155                             ScalarB * B_Pointers[i][j];
01156       UpdateFlops(3*GlobalLength_*NumVectors_);
01157     }
01158       }
01159     else if (ScalarThis==1.0)
01160       {
01161   if (ScalarA==1.0)
01162     {
01163       for (i = 0; i < NumVectors_; i++)
01164         for (j = 0; j < MyLength_; j++) Pointers_[i][j] +=           A_Pointers[i][j] + 
01165                              ScalarB * B_Pointers[i][j];
01166       UpdateFlops(3*GlobalLength_*NumVectors_);
01167     }
01168   else if (ScalarB==1.0)
01169     {
01170       for (i = 0; i < NumVectors_; i++)
01171         for (j = 0; j < MyLength_; j++) Pointers_[i][j] += ScalarA * A_Pointers[i][j] +
01172                                        B_Pointers[i][j];
01173       UpdateFlops(3*GlobalLength_*NumVectors_);
01174     }
01175   else
01176     {
01177       for (i = 0; i < NumVectors_; i++)
01178         for (j = 0; j < MyLength_; j++) Pointers_[i][j] += ScalarA * A_Pointers[i][j] + 
01179                              ScalarB * B_Pointers[i][j];
01180       UpdateFlops(4*GlobalLength_*NumVectors_);
01181     }
01182       }
01183     else
01184       {
01185   if (ScalarA==1.0)
01186     {
01187       for (i = 0; i < NumVectors_; i++)
01188         for (j = 0; j < MyLength_; j++) Pointers_[i][j] =  ScalarThis *    Pointers_[i][j] +
01189                                        A_Pointers[i][j] + 
01190                              ScalarB * B_Pointers[i][j];
01191       UpdateFlops(4*GlobalLength_*NumVectors_);
01192     }
01193   else if (ScalarB==1.0)
01194     {
01195       for (i = 0; i < NumVectors_; i++)
01196         for (j = 0; j < MyLength_; j++) Pointers_[i][j] =  ScalarThis *    Pointers_[i][j] +
01197                              ScalarA * A_Pointers[i][j] +
01198                                        B_Pointers[i][j];
01199       UpdateFlops(4*GlobalLength_*NumVectors_);
01200     }
01201   else
01202     {
01203       for (i = 0; i < NumVectors_; i++)
01204         for (j = 0; j < MyLength_; j++) Pointers_[i][j] =  ScalarThis *    Pointers_[i][j] +
01205                              ScalarA * A_Pointers[i][j] + 
01206                              ScalarB * B_Pointers[i][j];
01207       UpdateFlops(5*GlobalLength_*NumVectors_);
01208     }
01209       }
01210 
01211 
01212     return(0);
01213   }
01214 
01215 //=========================================================================
01216 int  Epetra_MultiVector::Norm1 (double* Result) const {
01217   
01218   // 1-norm of each vector in MultiVector 
01219     
01220   int i;
01221 
01222   UpdateDoubleTemp();
01223 
01224   for (i=0; i < NumVectors_; i++) DoubleTemp_[i] = ASUM(MyLength_, Pointers_[i]);
01225   
01226   Comm_->SumAll(DoubleTemp_, Result, NumVectors_);
01227   
01228   UpdateFlops(2*GlobalLength_*NumVectors_);
01229 
01230   return(0);
01231 }
01232 
01233 //=========================================================================
01234 int  Epetra_MultiVector::Norm2 (double* Result) const {
01235   
01236   // 2-norm of each vector in MultiVector 
01237   
01238   int i, j;
01239 
01240   UpdateDoubleTemp();
01241 
01242   for (i=0; i < NumVectors_; i++) 
01243     {
01244       double sum = 0.0;
01245       for (j=0; j < MyLength_; j++) sum += Pointers_[i][j] * Pointers_[i][j];
01246       DoubleTemp_[i] = sum;
01247     }
01248   Comm_->SumAll(DoubleTemp_, Result, NumVectors_);
01249   for (i=0; i < NumVectors_; i++) Result[i] = std::sqrt(Result[i]);
01250   
01251   UpdateFlops(2*GlobalLength_*NumVectors_);
01252   
01253   return(0);
01254 }
01255 
01256 //=========================================================================
01257 int  Epetra_MultiVector::NormInf (double* Result) const {
01258   
01259   // Inf-norm of each vector in MultiVector 
01260     
01261   int i, j;
01262 
01263   UpdateDoubleTemp();
01264 
01265   for (i=0; i < NumVectors_; i++) 
01266     {
01267       DoubleTemp_[i] = 0.0;
01268       j = IAMAX(MyLength_, Pointers_[i]);
01269       if (j>-1) DoubleTemp_[i] = std::abs(Pointers_[i][j]);
01270     }
01271   Comm_->MaxAll(DoubleTemp_, Result, NumVectors_);
01272   
01273   // UpdateFlops(0);  Strictly speaking there are not FLOPS in this routine  
01274   return(0);
01275 }
01276 
01277 //=========================================================================
01278 int  Epetra_MultiVector::NormWeighted (const Epetra_MultiVector& Weights, double* Result) const {
01279   
01280   // Weighted 2-norm of each vector in MultiVector 
01281 
01282   // If only one vector in Weights, we assume it will be used as the weights for all vectors
01283 
01284   int i, j;
01285   bool OneW = false;
01286   if (Weights.NumVectors()==1) OneW = true;
01287   else if (NumVectors_ != Weights.NumVectors()) EPETRA_CHK_ERR(-1);
01288 
01289   if (MyLength_ != Weights.MyLength()) EPETRA_CHK_ERR(-2);
01290 
01291 
01292   UpdateDoubleTemp();
01293 
01294   double *W = Weights.Values();
01295   double **W_Pointers = Weights.Pointers();
01296   
01297   for (i=0; i < NumVectors_; i++) 
01298     {
01299       if (!OneW) W = W_Pointers[i]; // If Weights has the same number of vectors as this, use each weight vector
01300       double sum = 0.0;
01301       for (j=0; j < MyLength_; j++) {
01302         double tmp = Pointers_[i][j]/W[j];
01303         sum += tmp * tmp;
01304       }
01305       DoubleTemp_[i] = sum;
01306     }
01307   Comm_->SumAll(DoubleTemp_, Result, NumVectors_);
01308   double OneOverN = 1.0 / (double) GlobalLength_;
01309   for (i=0; i < NumVectors_; i++) Result[i] = std::sqrt(Result[i]*OneOverN);
01310   
01311   UpdateFlops(3*GlobalLength_*NumVectors_);
01312   
01313   return(0);
01314   }
01315 
01316 //=========================================================================
01317 int  Epetra_MultiVector::MinValue (double* Result) const {
01318   
01319   // Minimum value of each vector in MultiVector 
01320   
01321   int i, j, ierr = 0;
01322 
01323   UpdateDoubleTemp();
01324 
01325   for (i=0; i < NumVectors_; i++) 
01326     {
01327       double MinVal = Epetra_MaxDouble;
01328       if (MyLength_>0) MinVal = Pointers_[i][0];
01329       for (j=0; j< MyLength_; j++) MinVal = EPETRA_MIN(MinVal,Pointers_[i][j]); 
01330       DoubleTemp_[i] = MinVal;
01331     }
01332 
01333   if (MyLength_ > 0) {
01334     for(i=0; i<NumVectors_; ++i) {
01335       Result[i] = DoubleTemp_[i];
01336     }
01337   }
01338 
01339   //If MyLength_ == 0 and Comm_->NumProc() == 1, then Result has
01340   //not been referenced. Also, if vector contents are uninitialized
01341   //then Result contents are not well defined...
01342 
01343   if (Comm_->NumProc() == 1) return(ierr);
01344 
01345   //We're going to use MPI_Allgather to gather every proc's local-
01346   //min values onto every other proc. We'll use the last position
01347   //of the DoubleTemp_ array to indicate whether this proc has
01348   //valid data that should be considered by other procs when forming
01349   //the global-min results.
01350 
01351   if (MyLength_ == 0) DoubleTemp_[NumVectors_] = 0.0;
01352   else DoubleTemp_[NumVectors_] = 1.0;
01353 
01354   //Now proceed to handle the parallel case. We'll gather local-min
01355   //values from other procs and form a global-min. If any processor
01356   //has MyLength_>0, we'll end up with a valid result.
01357 
01358 #ifdef EPETRA_MPI
01359   const Epetra_MpiComm* epetrampicomm =
01360     dynamic_cast<const Epetra_MpiComm*>(Comm_);
01361   if (!epetrampicomm) {
01362     return(-2);
01363   }
01364 
01365   MPI_Comm mpicomm = epetrampicomm->GetMpiComm();
01366   int numProcs = epetrampicomm->NumProc();
01367   double* dwork = new double[numProcs*(NumVectors_+1)];
01368 
01369   MPI_Allgather(DoubleTemp_, NumVectors_+1, MPI_DOUBLE,
01370                 dwork, NumVectors_+1, MPI_DOUBLE, mpicomm);
01371 
01372   //if MyLength_==0, then our Result array currently contains
01373   //Epetra_MaxDouble from the local-min calculations above. In this
01374   //case we'll overwrite our Result array with values from the first
01375   //processor that sent valid data.
01376   bool overwrite = MyLength_ == 0 ? true : false;
01377 
01378   int myPID = epetrampicomm->MyPID();
01379   double* dwork_ptr = dwork;
01380 
01381   for(j=0; j<numProcs; ++j) {
01382 
01383     //skip data from self, and skip data from
01384     //procs with DoubleTemp_[NumVectors_] == 0.0.
01385     if (j == myPID || dwork_ptr[NumVectors_] < 0.5) {
01386       dwork_ptr += NumVectors_+1;
01387       continue;
01388     }
01389 
01390     for(i=0; i<NumVectors_; ++i) {
01391       double val = dwork_ptr[i];
01392 
01393       //Set val into our Result array if overwrite is true (see above),
01394       //or if val is less than our current Result[i].
01395       if (overwrite || (Result[i] > val)) Result[i] = val;
01396     }
01397 
01398     //Now set overwrite to false so that we'll do the right thing
01399     //when processing data from subsequent processors.
01400     if (overwrite) overwrite = false;
01401 
01402     dwork_ptr += NumVectors_+1;
01403   }
01404 
01405   delete [] dwork;
01406 #endif
01407 
01408   // UpdateFlops(0);  Strictly speaking there are not FLOPS in this routine
01409   
01410   return(ierr);
01411 }
01412 
01413 //=========================================================================
01414 int  Epetra_MultiVector::MaxValue (double* Result) const {
01415   
01416   // Maximum value of each vector in MultiVector 
01417   
01418   int i, j, ierr = 0;
01419 
01420   UpdateDoubleTemp();
01421 
01422   for (i=0; i < NumVectors_; i++) 
01423     {
01424       double MaxVal = -Epetra_MaxDouble;
01425       if (MyLength_>0) MaxVal = Pointers_[i][0];
01426       for (j=0; j< MyLength_; j++) MaxVal = EPETRA_MAX(MaxVal,Pointers_[i][j]); 
01427       DoubleTemp_[i] = MaxVal;
01428     }
01429 
01430   if (MyLength_ > 0) {
01431     for(i=0; i<NumVectors_; ++i) {
01432       Result[i] = DoubleTemp_[i];
01433     }
01434   }
01435 
01436   //If MyLength_ == 0 and Comm_->NumProc() == 1, then Result has
01437   //not been referenced. Also, if vector contents are uninitialized
01438   //then Result contents are not well defined...
01439 
01440   if (Comm_->NumProc() == 1) return(ierr);
01441 
01442   //We're going to use MPI_Allgather to gather every proc's local-
01443   //max values onto every other proc. We'll use the last position
01444   //of the DoubleTemp_ array to indicate whether this proc has
01445   //valid data that should be considered by other procs when forming
01446   //the global-max results.
01447 
01448   if (MyLength_ == 0) DoubleTemp_[NumVectors_] = 0.0;
01449   else DoubleTemp_[NumVectors_] = 1.0;
01450 
01451   //Now proceed to handle the parallel case. We'll gather local-max
01452   //values from other procs and form a global-max. If any processor
01453   //has MyLength_>0, we'll end up with a valid result.
01454 
01455 #ifdef EPETRA_MPI
01456   const Epetra_MpiComm* epetrampicomm =
01457     dynamic_cast<const Epetra_MpiComm*>(Comm_);
01458   if (!epetrampicomm) {
01459     return(-2);
01460   }
01461 
01462   MPI_Comm mpicomm = epetrampicomm->GetMpiComm();
01463   int numProcs = epetrampicomm->NumProc();
01464   double* dwork = new double[numProcs*(NumVectors_+1)];
01465 
01466   MPI_Allgather(DoubleTemp_, NumVectors_+1, MPI_DOUBLE,
01467                 dwork, NumVectors_+1, MPI_DOUBLE, mpicomm);
01468 
01469   //if MyLength_==0, then our Result array currently contains
01470   //-Epetra_MaxDouble from the local-max calculations above. In this
01471   //case we'll overwrite our Result array with values from the first
01472   //processor that sent valid data.
01473   bool overwrite = MyLength_ == 0 ? true : false;
01474 
01475   int myPID = epetrampicomm->MyPID();
01476   double* dwork_ptr = dwork;
01477 
01478   for(j=0; j<numProcs; ++j) {
01479 
01480     //skip data from self, and skip data from
01481     //procs with DoubleTemp_[NumVectors_] == 0.0.
01482     if (j == myPID || dwork_ptr[NumVectors_] < 0.5) {
01483       dwork_ptr += NumVectors_+1;
01484       continue;
01485     }
01486 
01487     for(i=0; i<NumVectors_; ++i) {
01488       double val = dwork_ptr[i];
01489 
01490       //Set val into our Result array if overwrite is true (see above),
01491       //or if val is larger than our current Result[i].
01492       if (overwrite || (Result[i] < val)) Result[i] = val; 
01493     }
01494 
01495     //Now set overwrite to false so that we'll do the right thing
01496     //when processing data from subsequent processors.
01497     if (overwrite) overwrite = false;
01498 
01499     dwork_ptr += NumVectors_+1;
01500   }
01501 
01502   delete [] dwork;
01503 #endif
01504 
01505   // UpdateFlops(0);  Strictly speaking there are not FLOPS in this routine
01506   
01507   return(ierr);
01508 }
01509 
01510 //=========================================================================
01511 int  Epetra_MultiVector::MeanValue (double* Result) const {
01512   
01513   // Mean value of each vector in MultiVector 
01514   
01515   int i, j;
01516   double fGlobalLength = 1.0/EPETRA_MAX((double) GlobalLength_, 1.0);
01517   
01518 
01519   UpdateDoubleTemp();
01520 
01521   for (i=0; i < NumVectors_; i++) 
01522     {
01523       double sum = 0.0;
01524       for (j=0; j < MyLength_; j++) sum += Pointers_[i][j];
01525       DoubleTemp_[i] = sum;
01526     }
01527   Comm_->SumAll(DoubleTemp_, Result, NumVectors_);
01528   for (i=0; i < NumVectors_; i++) Result[i] = Result[i]*fGlobalLength;
01529   
01530   UpdateFlops(GlobalLength_*NumVectors_);
01531 
01532   return(0);
01533 }
01534 
01535   //=========================================================================
01536   int  Epetra_MultiVector::Multiply (char TransA, char TransB, double ScalarAB, 
01537           const Epetra_MultiVector& A, 
01538           const Epetra_MultiVector& B,
01539           double ScalarThis ) {
01540 
01541     // This routine performs a variety of matrix-matrix multiply operations, interpreting
01542     // the Epetra_MultiVector (this-aka C , A and B) as 2D matrices.  Variations are due to
01543     // the fact that A, B and C can be local replicated or global distributed
01544     // Epetra_MultiVectors and that we may or may not operate with the transpose of 
01545     // A and B.  Possible cases are:
01546 
01547     //                                       Num
01548     //     OPERATIONS                        case  Notes
01549     // 1) C(local) = A^X(local) * B^X(local)  4   (X=Trans or Not, No comm needed) 
01550     // 2) C(local) = A^T(distr) * B  (distr)  1   (2D dot product, replicate C)
01551     // 3) C(distr) = A  (distr) * B^X(local)  2   (2D vector update, no comm needed)
01552 
01553     // Note that the following operations are not meaningful for 
01554     // 1D distributions:
01555 
01556     // 1) C(local) = A^T(distr) * B^T(distr)  1
01557     // 2) C(local) = A  (distr) * B^X(distr)  2
01558     // 3) C(distr) = A^X(local) * B^X(local)  4
01559     // 4) C(distr) = A^X(local) * B^X(distr)  4
01560     // 5) C(distr) = A^T(distr) * B^X(local)  2
01561     // 6) C(local) = A^X(distr) * B^X(local)  4
01562     // 7) C(distr) = A^X(distr) * B^X(local)  4
01563     // 8) C(local) = A^X(local) * B^X(distr)  4
01564 
01565     // Total of 32 case (2^5).
01566 
01567     
01568     //if (!ConstantStride_    ||
01571 
01572     // Check for compatible dimensions
01573 
01574     int A_nrows = (TransA=='T') ? A.NumVectors() : A.MyLength();
01575     int A_ncols = (TransA=='T') ? A.MyLength() : A.NumVectors();
01576     int B_nrows = (TransB=='T') ? B.NumVectors() : B.MyLength();
01577     int B_ncols = (TransB=='T') ? B.MyLength() : B.NumVectors();
01578 
01579     double Scalar_local = ScalarThis; // local copy of Scalar
01580 
01581     if( MyLength_      != A_nrows     ||   // RAB: 2002/01/25: Minor reformat to allow
01582     A_ncols        != B_nrows     ||   //   setting breakpoint at error return.
01583     NumVectors_    != B_ncols  )
01584     EPETRA_CHK_ERR(-2); // Return error
01585 
01586     bool A_is_local = (!A.DistributedGlobal());
01587     bool B_is_local = (!B.DistributedGlobal());
01588     bool C_is_local = (!DistributedGlobal_);
01589     bool Case1 = ( A_is_local &&  B_is_local &&  C_is_local);  // Case 1 above
01590     bool Case2 = (!A_is_local && !B_is_local &&  C_is_local && TransA=='T' );// Case 2
01591     bool Case3 = (!A_is_local &&  B_is_local && !C_is_local && TransA=='N');// Case 3
01592   
01593     // Test for meaningful cases
01594 
01595     if (Case1 || Case2 || Case3)
01596       {
01597   if (ScalarThis!=0.0 && Case2)
01598     {
01599       const int MyPID = Comm_->MyPID();
01600       if (MyPID!=0) Scalar_local = 0.0;
01601     }
01602 
01603         // Check if A, B, C have constant stride, if not then make temp copy (strided)
01604 
01605         Epetra_MultiVector * A_tmp, * B_tmp, *C_tmp;
01606         if (!ConstantStride_) C_tmp = new Epetra_MultiVector(*this);
01607         else C_tmp = this;
01608           
01609         if (!A.ConstantStride()) A_tmp = new Epetra_MultiVector(A);
01610         else A_tmp = (Epetra_MultiVector *) &A;
01611     
01612         if (!B.ConstantStride()) B_tmp = new Epetra_MultiVector(B);
01613         else B_tmp = (Epetra_MultiVector *) &B;
01614       
01615     
01616   int m = MyLength_;
01617   int n = NumVectors_;
01618   int k = A_ncols;
01619   int lda = EPETRA_MAX(A_tmp->Stride(),1); // The reference BLAS implementation requires lda, ldb, ldc > 0 even if m, n, or k = 0
01620   int ldb = EPETRA_MAX(B_tmp->Stride(),1);
01621   int ldc = EPETRA_MAX(C_tmp->Stride(),1);
01622   double *Ap = A_tmp->Values();
01623   double *Bp = B_tmp->Values();
01624   double *Cp = C_tmp->Values();
01625    
01626   GEMM(TransA, TransB,  m, n, k, ScalarAB,
01627        Ap, lda, Bp, ldb, Scalar_local, Cp, ldc);
01628       
01629   // FLOP Counts
01630   //                                       Num
01631   //     OPERATIONS                        case  Notes 
01632   // 1) C(local) = A^X(local) * B^X(local)  4   (X=Trans or Not, No comm needed)
01633   // 2) C(local) = A^T(distr) * B  (distr)  1   (2D dot product, replicate C)      
01634   // 3) C(distr) = A  (distr) * B^X(local)  2   (2D vector update, no comm needed)
01635 
01636   // For Case 1 we only count the local operations, since we are interested in serial
01637   // cost.  Computation on other processors is redundant.
01638   if (Case1)
01639     {     
01640       UpdateFlops(2*m*n*k);
01641       if (ScalarAB!=1.0) UpdateFlops(m*n);
01642       if (ScalarThis==1.0) UpdateFlops(m*n); else if (ScalarThis!=0.0) UpdateFlops(2*m*n);
01643     }
01644   else if (Case2)
01645     {     
01646       UpdateFlops(2*m*n*A.GlobalLength());
01647       if (ScalarAB!=1.0) UpdateFlops(m*n);
01648       if (ScalarThis==1.0) UpdateFlops(m*n); else if (ScalarThis!=0.0) UpdateFlops(2*m*n);
01649     }
01650   else
01651     {
01652       UpdateFlops(2*GlobalLength_*n*k);
01653       if (ScalarAB!=1.0) UpdateFlops(GlobalLength_*n);
01654       if (ScalarThis==1.0) UpdateFlops(GlobalLength_*n);
01655       else if (ScalarThis!=0.0) UpdateFlops(2*GlobalLength_*n);
01656     }
01657 
01658   // If A or B were not strided, dispose of extra copies.
01659   if (!A.ConstantStride()) delete A_tmp;
01660   if (!B.ConstantStride()) delete B_tmp;
01661 
01662   // If C was not strided, copy from strided version and delete
01663   if (!ConstantStride_) 
01664     {
01665       C_tmp->ExtractCopy(Pointers_);
01666       delete C_tmp;
01667     }
01668 
01669   // If Case 2 then sum up C and distribute it to all processors.
01670 
01671   if (Case2) {EPETRA_CHK_ERR(Reduce());}
01672 
01673   return(0);
01674 
01675       }
01676     else {EPETRA_CHK_ERR(-3)}; // Return error: not a supported operation
01677 
01678   return(0);
01679   }
01680 
01681 
01682 //=========================================================================
01683 int Epetra_MultiVector::Multiply(double ScalarAB, const Epetra_MultiVector& A, const Epetra_MultiVector& B,
01684            double ScalarThis) {
01685   
01686   int i, j;
01687   
01688   // Hadamard product of two MultiVectors: 
01689   // this = ScalarThis * this + ScalarAB * A * B (element-wise)
01690   
01691   if (ScalarAB==0.0) {
01692     EPETRA_CHK_ERR(Scale(ScalarThis));
01693     return(0);
01694   }
01695          
01696   if (A.NumVectors() != 1 && A.NumVectors() != B.NumVectors()) EPETRA_CHK_ERR(-1); // A must have one column or be the same as B.
01697   if (NumVectors_ != B.NumVectors()) EPETRA_CHK_ERR(-2);
01698   if (MyLength_ != A.MyLength() || MyLength_ != B.MyLength()) EPETRA_CHK_ERR(-3);
01699   
01700   double **A_Pointers = (double**)A.Pointers();
01701   double **B_Pointers = (double**)B.Pointers();
01702 
01703   int IncA = 1;
01704   if (A.NumVectors() == 1 ) IncA = 0;
01705 
01706     if (ScalarThis==0.0) {
01707       if (ScalarAB==1.0)
01708   {
01709     for (i = 0; i < NumVectors_; i++) {
01710       double * Aptr = A_Pointers[i*IncA];
01711       for (j = 0; j < MyLength_; j++) {
01712               Pointers_[i][j] =  Aptr[j] * B_Pointers[i][j];
01713             }
01714     }
01715     UpdateFlops(GlobalLength_*NumVectors_);
01716   }
01717       else
01718   {
01719     for (i = 0; i < NumVectors_; i++) {
01720       double * Aptr = A_Pointers[i*IncA];
01721       for (j = 0; j < MyLength_; j++) {
01722               Pointers_[i][j] = ScalarAB * Aptr[j] * B_Pointers[i][j];
01723             }
01724     }
01725     UpdateFlops(2*GlobalLength_*NumVectors_);
01726   }
01727     }
01728     else if (ScalarThis==1.0) {
01729       if (ScalarAB==1.0)
01730   {
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] +=  Aptr[j] * B_Pointers[i][j];
01735             }
01736     }
01737     UpdateFlops(2*GlobalLength_*NumVectors_);
01738   }
01739       else {
01740     for (i = 0; i < NumVectors_; i++) {
01741       double * Aptr = A_Pointers[i*IncA];
01742       for (j = 0; j < MyLength_; j++) {
01743               Pointers_[i][j] += ScalarAB * Aptr[j] * B_Pointers[i][j];
01744             }
01745           }
01746           UpdateFlops(3*GlobalLength_*NumVectors_);
01747       }
01748     }
01749     else { // if (ScalarThis!=1.0 && ScalarThis !=0 ) 
01750       if (ScalarAB==1.0)
01751   {
01752     for (i = 0; i < NumVectors_; i++) {
01753       double * Aptr = A_Pointers[i*IncA];
01754       for (j = 0; j < MyLength_; j++) {
01755               Pointers_[i][j] =  ScalarThis * Pointers_[i][j] +
01756                   Aptr[j] * B_Pointers[i][j];
01757             }
01758     }
01759     UpdateFlops(3*GlobalLength_*NumVectors_);
01760   }
01761       else
01762   {
01763     for (i = 0; i < NumVectors_; i++) {
01764       double * Aptr = A_Pointers[i*IncA];
01765       for (j = 0; j < MyLength_; j++) {
01766               Pointers_[i][j] = ScalarThis * Pointers_[i][j] +
01767                      ScalarAB * Aptr[j] * B_Pointers[i][j];
01768             }
01769     }
01770     UpdateFlops(4*GlobalLength_*NumVectors_);
01771   }
01772     }
01773   return(0);
01774 }
01775 //=========================================================================
01776 int Epetra_MultiVector::ReciprocalMultiply(double ScalarAB, const Epetra_MultiVector& A, const Epetra_MultiVector& B,
01777            double ScalarThis) {
01778   
01779   int i, j;
01780   
01781   // Hadamard product of two MultiVectors: 
01782   // this = ScalarThis * this + ScalarAB * B / A (element-wise)
01783   
01784   if (ScalarAB==0.0) {
01785     EPETRA_CHK_ERR(Scale(ScalarThis));
01786     return(0);
01787   }
01788          
01789   if (A.NumVectors() != 1 && A.NumVectors() != B.NumVectors()) EPETRA_CHK_ERR(-1); // A must have one column or be the same as B.
01790   if (NumVectors_ != B.NumVectors()) EPETRA_CHK_ERR(-2);
01791   if (MyLength_ != A.MyLength() || MyLength_ != B.MyLength()) EPETRA_CHK_ERR(-3);
01792   
01793   double **A_Pointers = (double**)A.Pointers();
01794   double **B_Pointers = (double**)B.Pointers();
01795 
01796   int IncA = 1;
01797   if (A.NumVectors() == 1 ) IncA = 0;
01798 
01799     if (ScalarThis==0.0) {
01800       if (ScalarAB==1.0)
01801   {
01802     for (i = 0; i < NumVectors_; i++) {
01803       double * Aptr = A_Pointers[i*IncA];
01804       for (j = 0; j < MyLength_; j++) {
01805               Pointers_[i][j] = B_Pointers[i][j] / Aptr[j];
01806             }
01807     }
01808     UpdateFlops(GlobalLength_*NumVectors_);
01809   }
01810       else
01811   {
01812     for (i = 0; i < NumVectors_; i++) {
01813       double * Aptr = A_Pointers[i*IncA];
01814       for (j = 0; j < MyLength_; j++) {
01815               Pointers_[i][j] = ScalarAB * B_Pointers[i][j] / Aptr[j];
01816             }
01817     }
01818     UpdateFlops(2*GlobalLength_*NumVectors_);
01819   }
01820     }
01821     else if (ScalarThis==1.0) {
01822       if (ScalarAB==1.0)
01823   {
01824     for (i = 0; i < NumVectors_; i++) {
01825       double * Aptr = A_Pointers[i*IncA];
01826       for (j = 0; j < MyLength_; j++) {
01827               Pointers_[i][j] +=  B_Pointers[i][j] / Aptr[j];
01828             }
01829     }
01830     UpdateFlops(2*GlobalLength_*NumVectors_);
01831   }
01832       else
01833   {
01834     for (i = 0; i < NumVectors_; i++) {
01835       double * Aptr = A_Pointers[i*IncA];
01836       for (j = 0; j < MyLength_; j++) {
01837               Pointers_[i][j] += ScalarAB * B_Pointers[i][j] / Aptr[j];
01838             }
01839     }
01840     UpdateFlops(3*GlobalLength_*NumVectors_);
01841         }
01842     }
01843     else { // if (ScalarThis!=1.0 && ScalarThis !=0 ) 
01844       if (ScalarAB==1.0)
01845   {
01846     for (i = 0; i < NumVectors_; i++) {
01847       double * Aptr = A_Pointers[i*IncA];
01848       for (j = 0; j < MyLength_; j++) {
01849               Pointers_[i][j] =  ScalarThis * Pointers_[i][j] +
01850                    B_Pointers[i][j] / Aptr[j];
01851             }
01852     }
01853     UpdateFlops(3*GlobalLength_*NumVectors_);
01854   }
01855       else {
01856     for (i = 0; i < NumVectors_; i++) {
01857       double * Aptr = A_Pointers[i*IncA];
01858       for (j = 0; j < MyLength_; j++) {
01859               Pointers_[i][j] = ScalarThis * Pointers_[i][j] + ScalarAB * 
01860                 B_Pointers[i][j] / Aptr[j];
01861             }
01862     }
01863           UpdateFlops(4*GlobalLength_*NumVectors_);
01864       }
01865     }
01866   return(0);
01867 }
01868 
01869 
01870 //=======================================================================
01871 Epetra_Vector *& Epetra_MultiVector::operator () (int index)  {
01872   
01873   //  Epetra_MultiVector::operator () --- return non-const reference 
01874   
01875   if (index < 0 || index >=NumVectors_) 
01876     throw ReportError("Vector index = " + toString(index) + "is out of range. Number of Vectors = " + toString(NumVectors_), -1);
01877 
01878   UpdateVectors();
01879 
01880   // Create a new Epetra_Vector that is a view of ith vector, if not already present
01881   if (Vectors_[index]==0)
01882     Vectors_[index] = new Epetra_Vector(View, Map(), Pointers_[index]);
01883   return(Vectors_[index]);
01884 }
01885 
01886 //=======================================================================
01887 const Epetra_Vector *& Epetra_MultiVector::operator () (int index) const {
01888   
01889   //  Epetra_MultiVector::operator () --- return non-const reference 
01890 
01891   if (index < 0 || index >=NumVectors_) 
01892     throw ReportError("Vector index = " + toString(index) + "is out of range. Number of Vectors = " + toString(NumVectors_), -1);
01893 
01894   UpdateVectors();
01895 
01896   if (Vectors_[index]==0)
01897     Vectors_[index] = new Epetra_Vector(View, Map(), Pointers_[index]);
01898 
01899   const Epetra_Vector * & temp = (const Epetra_Vector * &) (Vectors_[index]);
01900   return(temp);
01901 }
01902 
01903 //========================================================================
01904 Epetra_MultiVector& Epetra_MultiVector::operator = (const Epetra_MultiVector& Source) {
01905   
01906   // Check for special case of this=Source
01907   if (this != &Source) Assign(Source);
01908   
01909   return(*this);
01910 }
01911 
01912 //=========================================================================
01913 void Epetra_MultiVector::Assign(const Epetra_MultiVector& A) {
01914   
01915   if (NumVectors_ != A.NumVectors())
01916     throw ReportError("Number of vectors incompatible in Assign.  The this MultiVector has NumVectors = " + toString(NumVectors_)
01917           + ".  The A MultiVector has NumVectors = " + toString(A.NumVectors()), -3);
01918   if (MyLength_ != A.MyLength())
01919     throw ReportError("Length of MultiVectors incompatible in Assign.  The this MultiVector has MyLength = " + toString(MyLength_)
01920           + ".  The A MultiVector has MyLength = " + toString(A.MyLength()), -4);
01921   
01922   double ** A_Pointers = A.Pointers();
01923   for (int i = 0; i< NumVectors_; i++)
01924       for (int j=0; j<MyLength_; j++) Pointers_[i][j] = A_Pointers[i][j];
01925     return;    
01926   }
01927 
01928   //=========================================================================
01929   int  Epetra_MultiVector::Reduce() {
01930 
01931     // Global reduction on each entry of a Replicated Local MultiVector
01932     int i, j;
01933     double * source = 0;
01934       if (MyLength_>0) source = new double[MyLength_*NumVectors_];
01935     double * target = 0;
01936     bool packed = (ConstantStride_ && (Stride_==MyLength_));
01937     if (packed) {
01938        for (i=0; i<MyLength_*NumVectors_; i++) source[i] = Values_[i];
01939        target = Values_;
01940     }
01941     else {
01942        double * tmp1 = source;
01943        for (i = 0; i < NumVectors_; i++) {
01944          double * tmp2 = Pointers_[i];
01945          for (j=0; j< MyLength_; j++) *tmp1++ = *tmp2++;
01946        }
01947        if (MyLength_>0) target = new double[MyLength_*NumVectors_];
01948     }
01949 
01950     Comm_->SumAll(source, target, MyLength_*NumVectors_);
01951     if (MyLength_>0) delete [] source;
01952     if (!packed) {
01953        double * tmp2 = target;
01954        for (i = 0; i < NumVectors_; i++) {
01955          double * tmp1 = Pointers_[i];
01956          for (j=0; j< MyLength_; j++) *tmp1++ = *tmp2++;
01957        }
01958        if (MyLength_>0) delete [] target;
01959     }
01960     // UpdateFlops(0);  No serial Flops in this function
01961     return(0);
01962   }
01963 //=======================================================================
01964 int Epetra_MultiVector::ResetView(double ** ArrayOfPointers) {
01965   
01966   if (!UserAllocated_) {
01967     EPETRA_CHK_ERR(-1); // Can't reset view if multivector was not allocated as a view
01968   }
01969 
01970   for (int i = 0; i< NumVectors_; i++) Pointers_[i] = ArrayOfPointers[i];
01971   DoView();
01972   
01973   return(0);
01974   }
01975 //=======================================================================
01976 void Epetra_MultiVector::Print(ostream& os) const {
01977   int MyPID = Map().Comm().MyPID();
01978   int NumProc = Map().Comm().NumProc();
01979   
01980   for (int iproc=0; iproc < NumProc; iproc++) {
01981     if (MyPID==iproc) {
01982       int NumVectors1 = NumVectors();
01983       int NumMyElements1 =Map(). NumMyElements();
01984       int MaxElementSize1 = Map().MaxElementSize();
01985       int * MyGlobalElements1 = Map().MyGlobalElements();
01986       int * FirstPointInElementList1;
01987       if (MaxElementSize1!=1) FirstPointInElementList1 = Map().FirstPointInElementList();
01988       double ** A_Pointers = Pointers();
01989 
01990       if (MyPID==0) {
01991   os.width(8);
01992   os <<  "     MyPID"; os << "    ";
01993   os.width(12);
01994   if (MaxElementSize1==1)
01995     os <<  "GID  ";
01996   else
01997     os <<  "     GID/Point";
01998   for (int j = 0; j < NumVectors1 ; j++)
01999     {   
02000       os.width(20);
02001       os <<  "Value  ";
02002     }
02003   os << endl;
02004       }
02005       for (int i=0; i < NumMyElements1; i++) {
02006   for (int ii=0; ii< Map().ElementSize(i); ii++) {
02007        int iii;
02008     os.width(10);
02009     os <<  MyPID; os << "    ";
02010     os.width(10);
02011     if (MaxElementSize1==1) {
02012       os << MyGlobalElements1[i] << "    ";
02013        iii = i;
02014        }
02015     else {
02016       os <<  MyGlobalElements1[i]<< "/" << ii << "    ";
02017          iii = FirstPointInElementList1[i]+ii;
02018        }
02019     for (int j = 0; j < NumVectors1 ; j++)
02020       {   
02021         os.width(20);
02022         os <<  A_Pointers[j][iii];
02023       }
02024     os << endl;
02025   }
02026       }
02027       os << flush; 
02028     }
02029 
02030     // Do a few global ops to give I/O a chance to complete
02031     Map().Comm().Barrier();
02032     Map().Comm().Barrier();
02033     Map().Comm().Barrier();
02034   }
02035   return;
02036 }

Generated on Wed May 12 21:41:05 2010 for Epetra Package Browser (Single Doxygen Collection) by  doxygen 1.4.7