EpetraExt 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   (void)Indexor;
00393 
00394   const EpetraExt_BlockDiagMatrix & A = dynamic_cast<const EpetraExt_BlockDiagMatrix &>(Source);
00395 
00396   double *From=A.Values();
00397   double *To = Values_;
00398 
00399   int * ToFirstPointInElementList = 0;
00400   int * FromFirstPointInElementList = 0;
00401   int * FromElementSizeList = 0;
00402   int MaxElementSize = DataMap().MaxElementSize();
00403   bool ConstantElementSize = DataMap().ConstantElementSize();
00404 
00405   if (!ConstantElementSize) {
00406     ToFirstPointInElementList =   DataMap().FirstPointInElementList();
00407     FromFirstPointInElementList = A.DataMap().FirstPointInElementList();
00408     FromElementSizeList = A.DataMap().ElementSizeList();
00409   }
00410   int j, jj, jjj, k;
00411   
00412   int NumSameEntries;
00413 
00414   bool Case1 = false;
00415   bool Case2 = false;
00416   // bool Case3 = false;
00417 
00418   if (MaxElementSize==1) {
00419     Case1 = true;
00420     NumSameEntries = NumSameIDs;
00421   }
00422   else if (ConstantElementSize) {
00423     Case2 = true;
00424     NumSameEntries = NumSameIDs * MaxElementSize;
00425   }
00426   else {
00427     // Case3 = true;
00428     NumSameEntries = FromFirstPointInElementList[NumSameIDs];
00429   }
00430 
00431   // Short circuit for the case where the source and target vector is the same.
00432   if (To==From) NumSameEntries = 0;
00433   
00434   // Do copy first
00435   if (NumSameIDs>0)
00436     if (To!=From) {
00437   for (j=0; j<NumSameEntries; j++)
00438     To[j] = From[j];
00439     }
00440   // Do local permutation next
00441   if (NumPermuteIDs>0) {
00442   
00443     // Point entry case
00444     if (Case1) {
00445       
00446       for (j=0; j<NumPermuteIDs; j++) 
00447   To[PermuteToLIDs[j]] = From[PermuteFromLIDs[j]];
00448     }
00449     // constant element size case
00450     else if (Case2) {
00451       
00452       for (j=0; j<NumPermuteIDs; j++) {
00453   jj = MaxElementSize*PermuteToLIDs[j];
00454   jjj = MaxElementSize*PermuteFromLIDs[j];
00455     for (k=0; k<MaxElementSize; k++)
00456       To[jj+k] = From[jjj+k];
00457       }
00458     }
00459     
00460     // variable element size case
00461     else {
00462       
00463       for (j=0; j<NumPermuteIDs; j++) {
00464   jj = ToFirstPointInElementList[PermuteToLIDs[j]];
00465   jjj = FromFirstPointInElementList[PermuteFromLIDs[j]];
00466   int ElementSize = FromElementSizeList[PermuteFromLIDs[j]];
00467     for (k=0; k<ElementSize; k++)
00468       To[jj+k] = From[jjj+k];
00469       }
00470     }
00471   }
00472   return(0);
00473 }
00474 
00475 //=========================================================================
00476 // Perform any packing or preparation required for call to DoTransfer().
00477 int EpetraExt_BlockDiagMatrix::PackAndPrepare(const Epetra_SrcDistObject& Source,
00478                                            int NumExportIDs,
00479                                            int* ExportLIDs,
00480                                            int& LenExports,
00481                                            char*& Exports,
00482                                            int& SizeOfPacket,
00483                                            int* Sizes,
00484                                            bool & VarSizes,
00485                                            Epetra_Distributor& Distor){
00486   (void)Sizes;
00487   (void)VarSizes;
00488   (void)Distor;
00489   const EpetraExt_BlockDiagMatrix & A = dynamic_cast<const EpetraExt_BlockDiagMatrix &>(Source);
00490 
00491   int j, jj, k;
00492 
00493   double *From=A.Values();
00494   int MaxElementSize = DataMap().MaxElementSize();
00495   bool ConstantElementSize = DataMap().ConstantElementSize();
00496 
00497   int * FromFirstPointInElementList = 0;
00498   int * FromElementSizeList = 0;
00499 
00500   if (!ConstantElementSize) {
00501     FromFirstPointInElementList = A.DataMap().FirstPointInElementList();
00502     FromElementSizeList = A.DataMap().ElementSizeList();
00503   }
00504 
00505   SizeOfPacket = MaxElementSize * (int)sizeof(double); 
00506 
00507   if(NumExportIDs*SizeOfPacket>LenExports) {
00508     if (LenExports>0) delete [] Exports;
00509     LenExports = NumExportIDs*SizeOfPacket;
00510     Exports = new char[LenExports];
00511   }
00512 
00513   double * ptr;
00514 
00515   if (NumExportIDs>0) {
00516     ptr = (double *) Exports;
00517     
00518     // Point entry case
00519     if (MaxElementSize==1) for (j=0; j<NumExportIDs; j++) *ptr++ = From[ExportLIDs[j]];
00520 
00521     // constant element size case
00522     else if (ConstantElementSize) {      
00523       for (j=0; j<NumExportIDs; j++) {
00524   jj = MaxElementSize*ExportLIDs[j];
00525     for (k=0; k<MaxElementSize; k++)
00526       *ptr++ = From[jj+k];
00527       }
00528     }
00529     
00530     // variable element size case
00531     else {     
00532       SizeOfPacket = MaxElementSize;
00533       for (j=0; j<NumExportIDs; j++) {
00534   ptr = (double *) Exports + j*SizeOfPacket;
00535   jj = FromFirstPointInElementList[ExportLIDs[j]];
00536   int ElementSize = FromElementSizeList[ExportLIDs[j]];
00537     for (k=0; k<ElementSize; k++)
00538       *ptr++ = From[jj+k];
00539       }
00540     }
00541   }
00542 
00543   return(0);
00544 }
00545 
00546 
00547 //=========================================================================
00548 // Perform any unpacking and combining after call to DoTransfer().
00549 int EpetraExt_BlockDiagMatrix::UnpackAndCombine(const Epetra_SrcDistObject& Source, 
00550                                              int NumImportIDs,
00551                                              int* ImportLIDs, 
00552                                              int LenImports,
00553                                              char* Imports,
00554                                              int& SizeOfPacket, 
00555                                              Epetra_Distributor& Distor,
00556                                              Epetra_CombineMode CombineMode,
00557                                              const Epetra_OffsetIndex * Indexor){
00558   (void)Source;
00559   (void)LenImports;
00560   (void)Distor;
00561   (void)Indexor;
00562   int j, jj, k;
00563   
00564   if(    CombineMode != Add
00565       && CombineMode != Zero
00566       && CombineMode != Insert
00567       && CombineMode != Average
00568       && CombineMode != AbsMax )
00569     EPETRA_CHK_ERR(-1); //Unsupported CombinedMode, will default to Zero
00570 
00571   if (NumImportIDs<=0) return(0);
00572 
00573   double * To = Values_;
00574   int MaxElementSize = DataMap().MaxElementSize();
00575   bool ConstantElementSize = DataMap().ConstantElementSize();
00576 
00577   int * ToFirstPointInElementList = 0;
00578   int * ToElementSizeList = 0;
00579 
00580   if (!ConstantElementSize) {
00581     ToFirstPointInElementList = DataMap().FirstPointInElementList();
00582     ToElementSizeList = DataMap().ElementSizeList();
00583   }
00584   
00585   double * ptr;
00586   // Unpack it...
00587 
00588   ptr = (double *) Imports;
00589     
00590   // Point entry case
00591   if (MaxElementSize==1) {
00592       
00593       if (CombineMode==Add)
00594   for (j=0; j<NumImportIDs; j++) To[ImportLIDs[j]] += *ptr++; // Add to existing value
00595       else if(CombineMode==Insert)
00596   for (j=0; j<NumImportIDs; j++) To[ImportLIDs[j]] = *ptr++;
00597       else if(CombineMode==AbsMax)
00598         for (j=0; j<NumImportIDs; j++) {
00599     To[ImportLIDs[j]] = EPETRA_MAX( To[ImportLIDs[j]],std::abs(*ptr));
00600     ptr++;
00601   }
00602       // Note:  The following form of averaging is not a true average if more that one value is combined.
00603       //        This might be an issue in the future, but we leave this way for now.
00604       else if(CombineMode==Average)
00605   for (j=0; j<NumImportIDs; j++) {To[ImportLIDs[j]] += *ptr++; To[ImportLIDs[j]] /= 2;}
00606   }
00607 
00608   // constant element size case
00609 
00610   else if (ConstantElementSize) {
00611    
00612     if (CombineMode==Add) {
00613       for (j=0; j<NumImportIDs; j++) {
00614   jj = MaxElementSize*ImportLIDs[j];
00615     for (k=0; k<MaxElementSize; k++)
00616       To[jj+k] += *ptr++; // Add to existing value
00617       }
00618     }
00619     else if(CombineMode==Insert) {
00620       for (j=0; j<NumImportIDs; j++) {
00621   jj = MaxElementSize*ImportLIDs[j];
00622     for (k=0; k<MaxElementSize; k++)
00623       To[jj+k] = *ptr++;
00624       }
00625     }
00626     else if(CombineMode==AbsMax) {
00627       for (j=0; j<NumImportIDs; j++) {
00628   jj = MaxElementSize*ImportLIDs[j];
00629   for (k=0; k<MaxElementSize; k++) {
00630       To[jj+k] = EPETRA_MAX( To[jj+k], std::abs(*ptr));
00631       ptr++;
00632   }
00633       }
00634     }
00635     // Note:  The following form of averaging is not a true average if more that one value is combined.
00636     //        This might be an issue in the future, but we leave this way for now.
00637     else if(CombineMode==Average) {
00638       for (j=0; j<NumImportIDs; j++) {
00639   jj = MaxElementSize*ImportLIDs[j];
00640     for (k=0; k<MaxElementSize; k++)
00641       { To[jj+k] += *ptr++; To[jj+k] /= 2;}
00642       }
00643     }
00644   }
00645     
00646   // variable element size case
00647 
00648   else {
00649       
00650     SizeOfPacket = MaxElementSize;
00651 
00652     if (CombineMode==Add) {
00653       for (j=0; j<NumImportIDs; j++) {
00654   ptr = (double *) Imports + j*SizeOfPacket;
00655   jj = ToFirstPointInElementList[ImportLIDs[j]];
00656   int ElementSize = ToElementSizeList[ImportLIDs[j]];
00657     for (k=0; k<ElementSize; k++)
00658       To[jj+k] += *ptr++; // Add to existing value
00659       }
00660     }
00661     else  if(CombineMode==Insert){
00662       for (j=0; j<NumImportIDs; j++) {
00663   ptr = (double *) Imports + j*SizeOfPacket;
00664   jj = ToFirstPointInElementList[ImportLIDs[j]];
00665   int ElementSize = ToElementSizeList[ImportLIDs[j]];
00666     for (k=0; k<ElementSize; k++)
00667       To[jj+k] = *ptr++;
00668       }
00669     }
00670     else  if(CombineMode==AbsMax){
00671       for (j=0; j<NumImportIDs; j++) {
00672   ptr = (double *) Imports + j*SizeOfPacket;
00673   jj = ToFirstPointInElementList[ImportLIDs[j]];
00674   int ElementSize = ToElementSizeList[ImportLIDs[j]];
00675   for (k=0; k<ElementSize; k++) {
00676       To[jj+k] = EPETRA_MAX( To[jj+k], std::abs(*ptr));
00677       ptr++;
00678   }
00679       }
00680     }
00681     // Note:  The following form of averaging is not a true average if more that one value is combined.
00682     //        This might be an issue in the future, but we leave this way for now.
00683     else if(CombineMode==Average) {
00684       for (j=0; j<NumImportIDs; j++) {
00685   ptr = (double *) Imports + j*SizeOfPacket;
00686   jj = ToFirstPointInElementList[ImportLIDs[j]];
00687   int ElementSize = ToElementSizeList[ImportLIDs[j]];
00688     for (k=0; k<ElementSize; k++)
00689       { To[jj+k] += *ptr++; To[jj+k] /= 2;}
00690       }
00691     }
00692   }
00693   
00694   return(0);
00695 }
00696   
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines