EpetraExt Package Browser (Single Doxygen Collection) Development
EpetraExt_BlockDiagMatrix.cpp
Go to the documentation of this file.
00001 /*
00002 //@HEADER
00003 // ***********************************************************************
00004 //
00005 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00006 //                 Copyright (2011) Sandia Corporation
00007 //
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00039 //
00040 // ***********************************************************************
00041 //@HEADER
00042 */
00043 
00044 #include "EpetraExt_BlockDiagMatrix.h"
00045 #include "Epetra_MultiVector.h"
00046 #include "Epetra_Comm.h"
00047 #include "Epetra_LAPACK.h"
00048 #include "Epetra_Distributor.h"
00049 
00050 #define AM_MULTIPLY 0
00051 #define AM_INVERT   1
00052 #define AM_FACTOR   2 
00053 
00054 //=========================================================================
00055 // Constructor
00056 EpetraExt_BlockDiagMatrix::EpetraExt_BlockDiagMatrix(const Epetra_BlockMap& Map,bool zero_out)
00057   : Epetra_DistObject(Map, "EpetraExt::BlockDiagMatrix"),
00058     HasComputed_(false),
00059     ApplyMode_(AM_MULTIPLY),
00060     DataMap_(0),
00061     Values_(0),
00062     Pivots_(0)
00063 {
00064   Allocate();
00065   if(zero_out) PutScalar(0.0);
00066 }
00067 
00068 
00069 //=========================================================================
00070 // Destructor
00071 EpetraExt_BlockDiagMatrix::~EpetraExt_BlockDiagMatrix(){
00072   if(DataMap_) delete DataMap_;
00073   if(Pivots_) delete [] Pivots_;
00074   if(Values_) delete [] Values_;
00075 }
00076 
00077 
00078 //=========================================================================
00079 // Copy constructor.  
00080 EpetraExt_BlockDiagMatrix::EpetraExt_BlockDiagMatrix(const EpetraExt_BlockDiagMatrix& Source)
00081   : Epetra_DistObject(Source),
00082     HasComputed_(false),
00083     ApplyMode_(AM_MULTIPLY),
00084     Values_(0),
00085     Pivots_(0)
00086 {
00087   DataMap_=new Epetra_BlockMap(*Source.DataMap_);
00088   Pivots_=new int[NumMyUnknowns()];
00089   Values_=new double[DataMap_->NumMyPoints()];
00090   DoCopy(Source);
00091 }
00092 
00093 //=========================================================================
00094 // Allocate
00095 void EpetraExt_BlockDiagMatrix::Allocate(){
00096 
00097   int DataSize=0, NumBlocks=NumMyBlocks();
00098   Pivots_=new int[NumMyUnknowns()];
00099   int *ElementSize=new int[NumBlocks];
00100   
00101   for(int i=0;i<NumBlocks;i++) {
00102     ElementSize[i]=BlockSize(i)*BlockSize(i);
00103     DataSize+=ElementSize[i];
00104   }
00105   
00106 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00107   if(Map().GlobalIndicesInt()) {
00108     DataMap_=new Epetra_BlockMap(-1,Map().NumMyElements(),Map().MyGlobalElements(),ElementSize,0,Map().Comm());
00109   }
00110   else
00111 #endif
00112 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00113   if(Map().GlobalIndicesLongLong()) {
00114     DataMap_=new Epetra_BlockMap((long long) -1,Map().NumMyElements(),Map().MyGlobalElements64(),ElementSize,0,Map().Comm());
00115   }
00116   else
00117 #endif
00118     throw "EpetraExt_BlockDiagMatrix::Allocate: GlobalIndices type unknown";
00119 
00120   Values_=new double[DataSize];  
00121   delete [] ElementSize;
00122 }
00123 
00124 
00125 //=========================================================================
00127 int EpetraExt_BlockDiagMatrix::SetParameters(Teuchos::ParameterList & List){
00128   List_=List;
00129 
00130   // Inverse storage mode
00131   std::string dummy("multiply");
00132   std::string ApplyMode=List_.get("apply mode",dummy);
00133   if(ApplyMode==std::string("multiply"))    ApplyMode_=AM_MULTIPLY;
00134   else if(ApplyMode==std::string("invert")) ApplyMode_=AM_INVERT;
00135   else if(ApplyMode==std::string("factor")) ApplyMode_=AM_FACTOR;
00136   else EPETRA_CHK_ERR(-1);
00137 
00138   return 0;
00139 }
00140 
00141 
00142 //=========================================================================
00143 void EpetraExt_BlockDiagMatrix::PutScalar(double value){
00144   int MaxData=NumData();
00145   for (int i=0;i<MaxData;i++) Values_[i]=value;
00146 }
00147 
00148 //=========================================================================
00149 // Assignment operator: Needs the same maps
00150 EpetraExt_BlockDiagMatrix& EpetraExt_BlockDiagMatrix::operator=(const EpetraExt_BlockDiagMatrix& Source){
00151   DoCopy(Source);
00152   return(*this);
00153 }
00154 //=========================================================================
00155 int EpetraExt_BlockDiagMatrix::DoCopy(const EpetraExt_BlockDiagMatrix& Source){
00156   // Need the same map
00157   if(!Map().SameAs(Source.Map()) || !DataMap_->SameAs(*Source.DataMap_))
00158     throw ReportError("Maps of DistBlockMatrices incompatible in assignment",-1);
00159 
00160   int MaxData=Source.NumData();
00161 
00162   for(int i=0;i<MaxData;i++)                 Values_[i]=Source.Values_[i];
00163   for(int i=0;i<Source.NumMyUnknowns();i++)  Pivots_[i]=Source.Pivots_[i];
00164 
00165   List_=Source.List_;
00166   ApplyMode_=Source.ApplyMode_;
00167   HasComputed_=Source.HasComputed_;
00168 
00169   return 0;
00170 }
00171 
00172 
00173 //=========================================================================
00175 int EpetraExt_BlockDiagMatrix::Compute(){
00176   int info;
00177 
00178   if(ApplyMode_==AM_MULTIPLY)
00179     // Multiply mode - noop    
00180     return 0;
00181   else {
00182     // Factorization - Needed for both 'factor' and 'invert' modes
00183     int NumBlocks=NumMyBlocks();
00184     for(int i=0;i<NumBlocks;i++){
00185       int Nb=BlockSize(i);
00186       if(Nb==1) {
00187         // Optimize for Block Size 1        
00188         Values_[DataMap_->FirstPointInElement(i)]=1.0/Values_[DataMap_->FirstPointInElement(i)];
00189       }
00190       else if(Nb==2) {
00191         // Optimize for Block Size 2
00192         double * v=&Values_[DataMap_->FirstPointInElement(i)];          
00193         double d=1/(v[0]*v[3]-v[1]*v[2]);
00194         double v0old=v[0];
00195         v[0]=v[3]*d;
00196         v[1]=-v[1]*d;
00197         v[2]=-v[2]*d;
00198         v[3]=v0old*d;
00199       }
00200       else{
00201         // "Large" Block - Use LAPACK
00202         LAPACK.GETRF(Nb,Nb,&Values_[DataMap_->FirstPointInElement(i)],Nb,&Pivots_[Map().FirstPointInElement(i)],&info);
00203         if(info) EPETRA_CHK_ERR(-2);
00204       }
00205     }
00206     
00207     if(ApplyMode_==AM_INVERT){
00208       // Invert - Use the factorization and invert the blocks in-place
00209       int lwork=3*DataMap_->MaxMyElementSize();
00210       std::vector<double> work(lwork);
00211       for(int i=0;i<NumBlocks;i++){
00212         int Nb=BlockSize(i);
00213         if(Nb==1 || Nb==2){
00214           // Optimize for Block Size 1 and 2
00215           // No need to do anything - factorization has already inverted the value
00216         }
00217         else{
00218           // "Large" Block - Use LAPACK
00219           LAPACK.GETRI(Nb,&Values_[DataMap_->FirstPointInElement(i)],Nb,&Pivots_[Map().FirstPointInElement(i)],&work[0],&lwork,&info);
00220           if(info) EPETRA_CHK_ERR(-3);
00221         }
00222       }
00223     }      
00224   }
00225   HasComputed_=true;
00226   return 0;
00227 }
00228 
00229 
00230 //=========================================================================
00231 int EpetraExt_BlockDiagMatrix::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
00232   int info;
00233   // Sanity Checks
00234   int NumVectors=X.NumVectors();  
00235   if(NumVectors!=Y.NumVectors())
00236     EPETRA_CHK_ERR(-1);
00237   if(!HasComputed_ && (ApplyMode_==AM_INVERT || ApplyMode_==AM_FACTOR))
00238     EPETRA_CHK_ERR(-2);
00239   
00240   //NTS: MultiVector's MyLength and [] Operators are  "points" level operators
00241   //not a "block/element" level operators.
00242 
00243   const int *vlist=DataMap_->FirstPointInElementList();
00244   const int *xlist=Map().FirstPointInElementList();
00245   const int *blocksize=Map().ElementSizeList();
00246   
00247   if(ApplyMode_==AM_MULTIPLY || ApplyMode_==AM_INVERT){
00248     // Multiply & Invert mode have the same apply
00249     int NumBlocks=NumMyBlocks();
00250     for(int i=0;i<NumBlocks;i++){
00251       int Nb=blocksize[i];
00252       int vidx0=vlist[i];
00253       int xidx0=xlist[i];
00254       for(int j=0;j<NumVectors;j++){  
00255   if(Nb==1) {
00256     // Optimize for size = 1
00257     Y[j][xidx0]=Values_[vidx0]*X[j][xidx0];
00258   }
00259   else if(Nb==2){
00260     // Optimize for size = 2
00261     Y[j][xidx0  ]=Values_[vidx0  ]*X[j][xidx0] + Values_[vidx0+2]*X[j][xidx0+1];
00262     Y[j][xidx0+1]=Values_[vidx0+1]*X[j][xidx0] + Values_[vidx0+3]*X[j][xidx0+1];
00263   }
00264   else{
00265     // "Large" Block - Use BLAS
00266     //void  GEMV (const char TRANS, const int M, const int N, const double ALPHA, const double *A, const int LDA, const double *X, const double BETA, double *Y, const int INCX=1, const int INCY=1) const 
00267     GEMV('N',Nb,Nb,1.0,&Values_[vidx0],Nb,&X[j][xidx0],0.0,&Y[j][xidx0]);
00268   }
00269       }   
00270     }
00271   }
00272   else{
00273     // Factorization mode has a different apply
00274     int NumBlocks=NumMyBlocks();
00275     for(int i=0;i<NumBlocks;i++){
00276       int Nb=blocksize[i];
00277       int vidx0=vlist[i];
00278       int xidx0=xlist[i];      
00279       for(int j=0;j<NumVectors;j++){
00280   if(Nb==1) {
00281     // Optimize for size = 1 - use the inverse
00282     Y[j][xidx0]=Values_[vidx0]*X[j][xidx0];
00283   }
00284   else if(Nb==2){
00285     // Optimize for size = 2 - use the inverse
00286     Y[j][xidx0  ]=Values_[vidx0  ]*X[j][xidx0] + Values_[vidx0+2]*X[j][xidx0+1];
00287     Y[j][xidx0+1]=Values_[vidx0+1]*X[j][xidx0] + Values_[vidx0+3]*X[j][xidx0+1];
00288   }
00289   else{
00290     // "Large" Block - use LAPACK
00291     //    void  GETRS (const char TRANS, const int N, const int NRHS, const double *A, const int LDA, const int *IPIV, double *X, const int LDX, int *INFO) const 
00292     for(int k=0;k<Nb;k++) Y[j][xidx0+k]=X[j][xidx0+k];
00293     LAPACK.GETRS('N',Nb,1,&Values_[vidx0],Nb,&Pivots_[xidx0],&Y[j][xidx0],Nb,&info);
00294     if(info) EPETRA_CHK_ERR(info);
00295   }
00296       }
00297     }    
00298   }  
00299   return 0;
00300 }
00301 
00302 
00303 
00304 
00305 //=========================================================================
00306 // Print method
00307 void EpetraExt_BlockDiagMatrix::Print(std::ostream & os) const{
00308   int MyPID = DataMap_->Comm().MyPID();
00309   int NumProc = DataMap_->Comm().NumProc();
00310   
00311   for (int iproc=0; iproc < NumProc; iproc++) {
00312     if (MyPID==iproc) {
00313       int NumMyElements1 =DataMap_->NumMyElements();
00314       int MaxElementSize1 = DataMap_->MaxElementSize();
00315     int * MyGlobalElements1_int = 0;
00316     long long * MyGlobalElements1_LL = 0;
00317 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00318   if(DataMap_->GlobalIndicesInt()) {
00319       MyGlobalElements1_int = DataMap_->MyGlobalElements();
00320   }
00321   else
00322 #endif
00323 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00324   if(DataMap_->GlobalIndicesLongLong()) {
00325       MyGlobalElements1_LL = DataMap_->MyGlobalElements64();
00326   }
00327   else
00328 #endif
00329     throw "EpetraExt_BlockDiagMatrix::Print: GlobalIndices type unknown";
00330 
00331       int * FirstPointInElementList1;
00332       if (MaxElementSize1!=1) FirstPointInElementList1 = DataMap_->FirstPointInElementList();
00333 
00334       if (MyPID==0) {
00335   os.width(8);
00336   os <<  "     MyPID"; os << "    ";
00337   os.width(12);
00338   if (MaxElementSize1==1)
00339     os <<  "GID  ";
00340   else
00341     os <<  "     GID/Point";
00342   os.width(20);
00343   os <<  "Values ";
00344   os << std::endl;
00345       }
00346       for (int i=0; i < NumMyElements1; i++) {
00347   for (int ii=0; ii< DataMap_->ElementSize(i); ii++) {
00348     int iii;
00349     os.width(10);
00350     os <<  MyPID; os << "    ";
00351     os.width(10);
00352     if (MaxElementSize1==1) {
00353       os << (MyGlobalElements1_int ? MyGlobalElements1_int[i] : MyGlobalElements1_LL[i]) << "    ";
00354       iii = i;
00355     }
00356     else {
00357       os <<  (MyGlobalElements1_int ? MyGlobalElements1_int[i] : MyGlobalElements1_LL[i]) << "/" << ii << "    ";
00358       iii = FirstPointInElementList1[i]+ii;
00359     }
00360           os.width(20);
00361           os <<  Values_[iii];
00362     os << std::endl;
00363   }
00364       }
00365       os << std::flush; 
00366     }
00367 
00368     // Do a few global ops to give I/O a chance to complete
00369     Map().Comm().Barrier();
00370     Map().Comm().Barrier();
00371     Map().Comm().Barrier();
00372   }
00373   return;
00374 }
00375 
00376 
00377 //=========================================================================
00378 // Allows the source and target (\e this) objects to be compared for compatibility, return nonzero if not.
00379 int EpetraExt_BlockDiagMatrix::CheckSizes(const Epetra_SrcDistObject& Source){
00380   return &Map() == &Source.Map();
00381 }
00382 
00383 
00384  //=========================================================================
00385 // Perform ID copies and permutations that are on processor.
00386 int EpetraExt_BlockDiagMatrix::CopyAndPermute(const Epetra_SrcDistObject& Source,
00387                                            int NumSameIDs, 
00388                                            int NumPermuteIDs,
00389                                            int * PermuteToLIDs,
00390                                            int * PermuteFromLIDs,
00391                                            const Epetra_OffsetIndex * Indexor,
00392                                            Epetra_CombineMode CombineMode){
00393   (void)Indexor;
00394 
00395   const EpetraExt_BlockDiagMatrix & A = dynamic_cast<const EpetraExt_BlockDiagMatrix &>(Source);
00396 
00397   double *From=A.Values();
00398   double *To = Values_;
00399 
00400   int * ToFirstPointInElementList = 0;
00401   int * FromFirstPointInElementList = 0;
00402   int * FromElementSizeList = 0;
00403   int MaxElementSize = DataMap().MaxElementSize();
00404   bool ConstantElementSize = DataMap().ConstantElementSize();
00405 
00406   if (!ConstantElementSize) {
00407     ToFirstPointInElementList =   DataMap().FirstPointInElementList();
00408     FromFirstPointInElementList = A.DataMap().FirstPointInElementList();
00409     FromElementSizeList = A.DataMap().ElementSizeList();
00410   }
00411   int j, jj, jjj, k;
00412   
00413   int NumSameEntries;
00414 
00415   bool Case1 = false;
00416   bool Case2 = false;
00417   // bool Case3 = false;
00418 
00419   if (MaxElementSize==1) {
00420     Case1 = true;
00421     NumSameEntries = NumSameIDs;
00422   }
00423   else if (ConstantElementSize) {
00424     Case2 = true;
00425     NumSameEntries = NumSameIDs * MaxElementSize;
00426   }
00427   else {
00428     // Case3 = true;
00429     NumSameEntries = FromFirstPointInElementList[NumSameIDs];
00430   }
00431 
00432   // Short circuit for the case where the source and target vector is the same.
00433   if (To==From) NumSameEntries = 0;
00434   
00435   // Do copy first
00436   if (NumSameIDs>0)
00437     if (To!=From) {
00438       if (CombineMode==Epetra_AddLocalAlso) for (int j=0; j<NumSameEntries; j++) To[j] += From[j]; // Add to existing value
00439       else for (int j=0; j<NumSameEntries; j++) To[j] = From[j];
00440     }
00441   // Do local permutation next
00442   if (NumPermuteIDs>0) {
00443   
00444     // Point entry case
00445     if (Case1) {
00446       
00447       if (CombineMode==Epetra_AddLocalAlso) for (int j=0; j<NumPermuteIDs; j++) To[PermuteToLIDs[j]] += From[PermuteFromLIDs[j]]; // Add to existing value
00448       else for (int j=0; j<NumPermuteIDs; j++) To[PermuteToLIDs[j]] = From[PermuteFromLIDs[j]];
00449     }
00450     // constant element size case
00451     else if (Case2) {
00452       
00453       for (j=0; j<NumPermuteIDs; j++) {
00454   jj = MaxElementSize*PermuteToLIDs[j];
00455   jjj = MaxElementSize*PermuteFromLIDs[j];
00456       if (CombineMode==Epetra_AddLocalAlso) for (k=0; k<MaxElementSize; k++) To[jj+k] += From[jjj+k]; // Add to existing value
00457       else for (k=0; k<MaxElementSize; k++) To[jj+k] = From[jjj+k];
00458       }
00459     }
00460     
00461     // variable element size case
00462     else {
00463       
00464       for (j=0; j<NumPermuteIDs; j++) {
00465   jj = ToFirstPointInElementList[PermuteToLIDs[j]];
00466   jjj = FromFirstPointInElementList[PermuteFromLIDs[j]];
00467   int ElementSize = FromElementSizeList[PermuteFromLIDs[j]];
00468       if (CombineMode==Epetra_AddLocalAlso) for (k=0; k<ElementSize; k++) To[jj+k] += From[jjj+k]; // Add to existing value
00469       else for (k=0; k<ElementSize; k++) To[jj+k] = From[jjj+k];
00470       }
00471     }
00472   }
00473   return(0);
00474 }
00475 
00476 //=========================================================================
00477 // Perform any packing or preparation required for call to DoTransfer().
00478 int EpetraExt_BlockDiagMatrix::PackAndPrepare(const Epetra_SrcDistObject& Source,
00479                                            int NumExportIDs,
00480                                            int* ExportLIDs,
00481                                            int& LenExports,
00482                                            char*& Exports,
00483                                            int& SizeOfPacket,
00484                                            int* Sizes,
00485                                            bool & VarSizes,
00486                                            Epetra_Distributor& Distor){
00487   (void)Sizes;
00488   (void)VarSizes;
00489   (void)Distor;
00490   const EpetraExt_BlockDiagMatrix & A = dynamic_cast<const EpetraExt_BlockDiagMatrix &>(Source);
00491 
00492   int j, jj, k;
00493 
00494   double *From=A.Values();
00495   int MaxElementSize = DataMap().MaxElementSize();
00496   bool ConstantElementSize = DataMap().ConstantElementSize();
00497 
00498   int * FromFirstPointInElementList = 0;
00499   int * FromElementSizeList = 0;
00500 
00501   if (!ConstantElementSize) {
00502     FromFirstPointInElementList = A.DataMap().FirstPointInElementList();
00503     FromElementSizeList = A.DataMap().ElementSizeList();
00504   }
00505 
00506   SizeOfPacket = MaxElementSize * (int)sizeof(double); 
00507 
00508   if(NumExportIDs*SizeOfPacket>LenExports) {
00509     if (LenExports>0) delete [] Exports;
00510     LenExports = NumExportIDs*SizeOfPacket;
00511     Exports = new char[LenExports];
00512   }
00513 
00514   double * ptr;
00515 
00516   if (NumExportIDs>0) {
00517     ptr = (double *) Exports;
00518     
00519     // Point entry case
00520     if (MaxElementSize==1) for (j=0; j<NumExportIDs; j++) *ptr++ = From[ExportLIDs[j]];
00521 
00522     // constant element size case
00523     else if (ConstantElementSize) {      
00524       for (j=0; j<NumExportIDs; j++) {
00525   jj = MaxElementSize*ExportLIDs[j];
00526     for (k=0; k<MaxElementSize; k++)
00527       *ptr++ = From[jj+k];
00528       }
00529     }
00530     
00531     // variable element size case
00532     else {     
00533       SizeOfPacket = MaxElementSize;
00534       for (j=0; j<NumExportIDs; j++) {
00535   ptr = (double *) Exports + j*SizeOfPacket;
00536   jj = FromFirstPointInElementList[ExportLIDs[j]];
00537   int ElementSize = FromElementSizeList[ExportLIDs[j]];
00538     for (k=0; k<ElementSize; k++)
00539       *ptr++ = From[jj+k];
00540       }
00541     }
00542   }
00543 
00544   return(0);
00545 }
00546 
00547 
00548 //=========================================================================
00549 // Perform any unpacking and combining after call to DoTransfer().
00550 int EpetraExt_BlockDiagMatrix::UnpackAndCombine(const Epetra_SrcDistObject& Source, 
00551                                              int NumImportIDs,
00552                                              int* ImportLIDs, 
00553                                              int LenImports,
00554                                              char* Imports,
00555                                              int& SizeOfPacket, 
00556                                              Epetra_Distributor& Distor,
00557                                              Epetra_CombineMode CombineMode,
00558                                              const Epetra_OffsetIndex * Indexor){
00559   (void)Source;
00560   (void)LenImports;
00561   (void)Distor;
00562   (void)Indexor;
00563   int j, jj, k;
00564   
00565   if(    CombineMode != Add
00566       && CombineMode != Zero
00567       && CombineMode != Insert
00568       && CombineMode != Average
00569       && CombineMode != AbsMax )
00570     EPETRA_CHK_ERR(-1); //Unsupported CombinedMode, will default to Zero
00571 
00572   if (NumImportIDs<=0) return(0);
00573 
00574   double * To = Values_;
00575   int MaxElementSize = DataMap().MaxElementSize();
00576   bool ConstantElementSize = DataMap().ConstantElementSize();
00577 
00578   int * ToFirstPointInElementList = 0;
00579   int * ToElementSizeList = 0;
00580 
00581   if (!ConstantElementSize) {
00582     ToFirstPointInElementList = DataMap().FirstPointInElementList();
00583     ToElementSizeList = DataMap().ElementSizeList();
00584   }
00585   
00586   double * ptr;
00587   // Unpack it...
00588 
00589   ptr = (double *) Imports;
00590     
00591   // Point entry case
00592   if (MaxElementSize==1) {
00593       
00594       if (CombineMode==Add)
00595   for (j=0; j<NumImportIDs; j++) To[ImportLIDs[j]] += *ptr++; // Add to existing value
00596       else if(CombineMode==Insert)
00597   for (j=0; j<NumImportIDs; j++) To[ImportLIDs[j]] = *ptr++;
00598       else if(CombineMode==AbsMax)
00599         for (j=0; j<NumImportIDs; j++) {
00600     To[ImportLIDs[j]] = EPETRA_MAX( To[ImportLIDs[j]],std::abs(*ptr));
00601     ptr++;
00602   }
00603       // Note:  The following form of averaging is not a true average if more that one value is combined.
00604       //        This might be an issue in the future, but we leave this way for now.
00605       else if(CombineMode==Average)
00606   for (j=0; j<NumImportIDs; j++) {To[ImportLIDs[j]] += *ptr++; To[ImportLIDs[j]] /= 2;}
00607   }
00608 
00609   // constant element size case
00610 
00611   else if (ConstantElementSize) {
00612    
00613     if (CombineMode==Add) {
00614       for (j=0; j<NumImportIDs; j++) {
00615   jj = MaxElementSize*ImportLIDs[j];
00616     for (k=0; k<MaxElementSize; k++)
00617       To[jj+k] += *ptr++; // Add to existing value
00618       }
00619     }
00620     else if(CombineMode==Insert) {
00621       for (j=0; j<NumImportIDs; j++) {
00622   jj = MaxElementSize*ImportLIDs[j];
00623     for (k=0; k<MaxElementSize; k++)
00624       To[jj+k] = *ptr++;
00625       }
00626     }
00627     else if(CombineMode==AbsMax) {
00628       for (j=0; j<NumImportIDs; j++) {
00629   jj = MaxElementSize*ImportLIDs[j];
00630   for (k=0; k<MaxElementSize; k++) {
00631       To[jj+k] = EPETRA_MAX( To[jj+k], std::abs(*ptr));
00632       ptr++;
00633   }
00634       }
00635     }
00636     // Note:  The following form of averaging is not a true average if more that one value is combined.
00637     //        This might be an issue in the future, but we leave this way for now.
00638     else if(CombineMode==Average) {
00639       for (j=0; j<NumImportIDs; j++) {
00640   jj = MaxElementSize*ImportLIDs[j];
00641     for (k=0; k<MaxElementSize; k++)
00642       { To[jj+k] += *ptr++; To[jj+k] /= 2;}
00643       }
00644     }
00645   }
00646     
00647   // variable element size case
00648 
00649   else {
00650       
00651     SizeOfPacket = MaxElementSize;
00652 
00653     if (CombineMode==Add) {
00654       for (j=0; j<NumImportIDs; j++) {
00655   ptr = (double *) Imports + j*SizeOfPacket;
00656   jj = ToFirstPointInElementList[ImportLIDs[j]];
00657   int ElementSize = ToElementSizeList[ImportLIDs[j]];
00658     for (k=0; k<ElementSize; k++)
00659       To[jj+k] += *ptr++; // Add to existing value
00660       }
00661     }
00662     else  if(CombineMode==Insert){
00663       for (j=0; j<NumImportIDs; j++) {
00664   ptr = (double *) Imports + j*SizeOfPacket;
00665   jj = ToFirstPointInElementList[ImportLIDs[j]];
00666   int ElementSize = ToElementSizeList[ImportLIDs[j]];
00667     for (k=0; k<ElementSize; k++)
00668       To[jj+k] = *ptr++;
00669       }
00670     }
00671     else  if(CombineMode==AbsMax){
00672       for (j=0; j<NumImportIDs; j++) {
00673   ptr = (double *) Imports + j*SizeOfPacket;
00674   jj = ToFirstPointInElementList[ImportLIDs[j]];
00675   int ElementSize = ToElementSizeList[ImportLIDs[j]];
00676   for (k=0; k<ElementSize; k++) {
00677       To[jj+k] = EPETRA_MAX( To[jj+k], std::abs(*ptr));
00678       ptr++;
00679   }
00680       }
00681     }
00682     // Note:  The following form of averaging is not a true average if more that one value is combined.
00683     //        This might be an issue in the future, but we leave this way for now.
00684     else if(CombineMode==Average) {
00685       for (j=0; j<NumImportIDs; j++) {
00686   ptr = (double *) Imports + j*SizeOfPacket;
00687   jj = ToFirstPointInElementList[ImportLIDs[j]];
00688   int ElementSize = ToElementSizeList[ImportLIDs[j]];
00689     for (k=0; k<ElementSize; k++)
00690       { To[jj+k] += *ptr++; To[jj+k] /= 2;}
00691       }
00692     }
00693   }
00694   
00695   return(0);
00696 }
00697   
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines