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   DataMap_=new Epetra_BlockMap(-1,Map().NumMyElements(),Map().MyGlobalElements(),ElementSize,0,Map().Comm());
00107   Values_=new double[DataSize];  
00108   delete [] ElementSize;
00109 }
00110 
00111 
00112 //=========================================================================
00114 int EpetraExt_BlockDiagMatrix::SetParameters(Teuchos::ParameterList & List){
00115   List_=List;
00116 
00117   // Inverse storage mode
00118   string dummy("multiply");
00119   string ApplyMode=List_.get("apply mode",dummy);
00120   if(ApplyMode==string("multiply"))    ApplyMode_=AM_MULTIPLY;
00121   else if(ApplyMode==string("invert")) ApplyMode_=AM_INVERT;
00122   else if(ApplyMode==string("factor")) ApplyMode_=AM_FACTOR;
00123   else EPETRA_CHK_ERR(-1);
00124 
00125   return 0;
00126 }
00127 
00128 
00129 //=========================================================================
00130 void EpetraExt_BlockDiagMatrix::PutScalar(double value){
00131   int MaxData=NumData();
00132   for (int i=0;i<MaxData;i++) Values_[i]=value;
00133 }
00134 
00135 //=========================================================================
00136 // Assignment operator: Needs the same maps
00137 EpetraExt_BlockDiagMatrix& EpetraExt_BlockDiagMatrix::operator=(const EpetraExt_BlockDiagMatrix& Source){
00138   DoCopy(Source);
00139   return(*this);
00140 }
00141 //=========================================================================
00142 int EpetraExt_BlockDiagMatrix::DoCopy(const EpetraExt_BlockDiagMatrix& Source){
00143   // Need the same map
00144   if(!Map().SameAs(Source.Map()) || !DataMap_->SameAs(*Source.DataMap_))
00145     throw ReportError("Maps of DistBlockMatrices incompatible in assignment",-1);
00146 
00147   int MaxData=Source.NumData();
00148 
00149   for(int i=0;i<MaxData;i++)                 Values_[i]=Source.Values_[i];
00150   for(int i=0;i<Source.NumMyUnknowns();i++)  Pivots_[i]=Source.Pivots_[i];
00151 
00152   List_=Source.List_;
00153   ApplyMode_=Source.ApplyMode_;
00154   HasComputed_=Source.HasComputed_;
00155 
00156   return 0;
00157 }
00158 
00159 
00160 //=========================================================================
00162 int EpetraExt_BlockDiagMatrix::Compute(){
00163   int info;
00164 
00165   if(ApplyMode_==AM_MULTIPLY)
00166     // Multiply mode - noop    
00167     return 0;
00168   else {
00169     // Factorization - Needed for both 'factor' and 'invert' modes
00170     int NumBlocks=NumMyBlocks();
00171     for(int i=0;i<NumBlocks;i++){
00172       int Nb=BlockSize(i);
00173       if(Nb==1) {
00174         // Optimize for Block Size 1        
00175         Values_[DataMap_->FirstPointInElement(i)]=1.0/Values_[DataMap_->FirstPointInElement(i)];
00176       }
00177       else if(Nb==2) {
00178         // Optimize for Block Size 2
00179         double * v=&Values_[DataMap_->FirstPointInElement(i)];          
00180         double d=1/(v[0]*v[3]-v[1]*v[2]);
00181         double v0old=v[0];
00182         v[0]=v[3]*d;
00183         v[1]=-v[1]*d;
00184         v[2]=-v[2]*d;
00185         v[3]=v0old*d;
00186       }
00187       else{
00188         // "Large" Block - Use LAPACK
00189         LAPACK.GETRF(Nb,Nb,&Values_[DataMap_->FirstPointInElement(i)],Nb,&Pivots_[Map().FirstPointInElement(i)],&info);
00190         if(info) EPETRA_CHK_ERR(-2);
00191       }
00192     }
00193     
00194     if(ApplyMode_==AM_INVERT){
00195       // Invert - Use the factorization and invert the blocks in-place
00196       int lwork=3*DataMap_->MaxMyElementSize();
00197       std::vector<double> work(lwork);
00198       for(int i=0;i<NumBlocks;i++){
00199         int Nb=BlockSize(i);
00200         if(Nb==1 || Nb==2){
00201           // Optimize for Block Size 1 and 2
00202           // No need to do anything - factorization has already inverted the value
00203         }
00204         else{
00205           // "Large" Block - Use LAPACK
00206           LAPACK.GETRI(Nb,&Values_[DataMap_->FirstPointInElement(i)],Nb,&Pivots_[Map().FirstPointInElement(i)],&work[0],&lwork,&info);
00207           if(info) EPETRA_CHK_ERR(-3);
00208         }
00209       }
00210     }      
00211   }
00212   HasComputed_=true;
00213   return 0;
00214 }
00215 
00216 
00217 //=========================================================================
00218 int EpetraExt_BlockDiagMatrix::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
00219   int info;
00220   // Sanity Checks
00221   int NumVectors=X.NumVectors();  
00222   if(NumVectors!=Y.NumVectors())
00223     EPETRA_CHK_ERR(-1);
00224   if(!HasComputed_ && (ApplyMode_==AM_INVERT || ApplyMode_==AM_FACTOR))
00225     EPETRA_CHK_ERR(-2);
00226   
00227   //NTS: MultiVector's MyLength and [] Operators are  "points" level operators
00228   //not a "block/element" level operators.
00229 
00230   const int *vlist=DataMap_->FirstPointInElementList();
00231   const int *xlist=Map().FirstPointInElementList();
00232   const int *blocksize=Map().ElementSizeList();
00233   
00234   if(ApplyMode_==AM_MULTIPLY || ApplyMode_==AM_INVERT){
00235     // Multiply & Invert mode have the same apply
00236     int NumBlocks=NumMyBlocks();
00237     for(int i=0;i<NumBlocks;i++){
00238       int Nb=blocksize[i];
00239       int vidx0=vlist[i];
00240       int xidx0=xlist[i];
00241       for(int j=0;j<NumVectors;j++){  
00242   if(Nb==1) {
00243     // Optimize for size = 1
00244     Y[j][xidx0]=Values_[vidx0]*X[j][xidx0];
00245   }
00246   else if(Nb==2){
00247     // Optimize for size = 2
00248     Y[j][xidx0  ]=Values_[vidx0  ]*X[j][xidx0] + Values_[vidx0+2]*X[j][xidx0+1];
00249     Y[j][xidx0+1]=Values_[vidx0+1]*X[j][xidx0] + Values_[vidx0+3]*X[j][xidx0+1];
00250   }
00251   else{
00252     // "Large" Block - Use BLAS
00253     //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 
00254     GEMV('N',Nb,Nb,1.0,&Values_[vidx0],Nb,&X[j][xidx0],0.0,&Y[j][xidx0]);
00255   }
00256       }   
00257     }
00258   }
00259   else{
00260     // Factorization mode has a different apply
00261     int NumBlocks=NumMyBlocks();
00262     for(int i=0;i<NumBlocks;i++){
00263       int Nb=blocksize[i];
00264       int vidx0=vlist[i];
00265       int xidx0=xlist[i];      
00266       for(int j=0;j<NumVectors;j++){
00267   if(Nb==1) {
00268     // Optimize for size = 1 - use the inverse
00269     Y[j][xidx0]=Values_[vidx0]*X[j][xidx0];
00270   }
00271   else if(Nb==2){
00272     // Optimize for size = 2 - use the inverse
00273     Y[j][xidx0  ]=Values_[vidx0  ]*X[j][xidx0] + Values_[vidx0+2]*X[j][xidx0+1];
00274     Y[j][xidx0+1]=Values_[vidx0+1]*X[j][xidx0] + Values_[vidx0+3]*X[j][xidx0+1];
00275   }
00276   else{
00277     // "Large" Block - use LAPACK
00278     //    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 
00279     for(int k=0;k<Nb;k++) Y[j][xidx0+k]=X[j][xidx0+k];
00280     LAPACK.GETRS('N',Nb,1,&Values_[vidx0],Nb,&Pivots_[xidx0],&Y[j][xidx0],Nb,&info);
00281     if(info) EPETRA_CHK_ERR(info);
00282   }
00283       }
00284     }    
00285   }  
00286   return 0;
00287 }
00288 
00289 
00290 
00291 
00292 //=========================================================================
00293 // Print method
00294 void EpetraExt_BlockDiagMatrix::Print(ostream & os) const{
00295   int MyPID = DataMap_->Comm().MyPID();
00296   int NumProc = DataMap_->Comm().NumProc();
00297   
00298   for (int iproc=0; iproc < NumProc; iproc++) {
00299     if (MyPID==iproc) {
00300       int NumMyElements1 =DataMap_->NumMyElements();
00301       int MaxElementSize1 = DataMap_->MaxElementSize();
00302       int * MyGlobalElements1 = DataMap_->MyGlobalElements();
00303       int * FirstPointInElementList1;
00304       if (MaxElementSize1!=1) FirstPointInElementList1 = DataMap_->FirstPointInElementList();
00305 
00306       if (MyPID==0) {
00307   os.width(8);
00308   os <<  "     MyPID"; os << "    ";
00309   os.width(12);
00310   if (MaxElementSize1==1)
00311     os <<  "GID  ";
00312   else
00313     os <<  "     GID/Point";
00314   os.width(20);
00315   os <<  "Values ";
00316   os << endl;
00317       }
00318       for (int i=0; i < NumMyElements1; i++) {
00319   for (int ii=0; ii< DataMap_->ElementSize(i); ii++) {
00320     int iii;
00321     os.width(10);
00322     os <<  MyPID; os << "    ";
00323     os.width(10);
00324     if (MaxElementSize1==1) {
00325       os << MyGlobalElements1[i] << "    ";
00326       iii = i;
00327     }
00328     else {
00329       os <<  MyGlobalElements1[i]<< "/" << ii << "    ";
00330       iii = FirstPointInElementList1[i]+ii;
00331     }
00332           os.width(20);
00333           os <<  Values_[iii];
00334     os << endl;
00335   }
00336       }
00337       os << flush; 
00338     }
00339 
00340     // Do a few global ops to give I/O a chance to complete
00341     Map().Comm().Barrier();
00342     Map().Comm().Barrier();
00343     Map().Comm().Barrier();
00344   }
00345   return;
00346 }
00347 
00348 
00349 //=========================================================================
00350 // Allows the source and target (\e this) objects to be compared for compatibility, return nonzero if not.
00351 int EpetraExt_BlockDiagMatrix::CheckSizes(const Epetra_SrcDistObject& Source){
00352   return &Map() == &Source.Map();
00353 }
00354 
00355 
00356  //=========================================================================
00357 // Perform ID copies and permutations that are on processor.
00358 int EpetraExt_BlockDiagMatrix::CopyAndPermute(const Epetra_SrcDistObject& Source,
00359                                            int NumSameIDs, 
00360                                            int NumPermuteIDs,
00361                                            int * PermuteToLIDs,
00362                                            int * PermuteFromLIDs,
00363                                            const Epetra_OffsetIndex * Indexor){
00364   (void)Indexor;
00365 
00366   const EpetraExt_BlockDiagMatrix & A = dynamic_cast<const EpetraExt_BlockDiagMatrix &>(Source);
00367 
00368   double *From=A.Values();
00369   double *To = Values_;
00370 
00371   int * ToFirstPointInElementList = 0;
00372   int * FromFirstPointInElementList = 0;
00373   int * FromElementSizeList = 0;
00374   int MaxElementSize = DataMap().MaxElementSize();
00375   bool ConstantElementSize = DataMap().ConstantElementSize();
00376 
00377   if (!ConstantElementSize) {
00378     ToFirstPointInElementList =   DataMap().FirstPointInElementList();
00379     FromFirstPointInElementList = A.DataMap().FirstPointInElementList();
00380     FromElementSizeList = A.DataMap().ElementSizeList();
00381   }
00382   int j, jj, jjj, k;
00383   
00384   int NumSameEntries;
00385 
00386   bool Case1 = false;
00387   bool Case2 = false;
00388   // bool Case3 = false;
00389 
00390   if (MaxElementSize==1) {
00391     Case1 = true;
00392     NumSameEntries = NumSameIDs;
00393   }
00394   else if (ConstantElementSize) {
00395     Case2 = true;
00396     NumSameEntries = NumSameIDs * MaxElementSize;
00397   }
00398   else {
00399     // Case3 = true;
00400     NumSameEntries = FromFirstPointInElementList[NumSameIDs];
00401   }
00402 
00403   // Short circuit for the case where the source and target vector is the same.
00404   if (To==From) NumSameEntries = 0;
00405   
00406   // Do copy first
00407   if (NumSameIDs>0)
00408     if (To!=From) {
00409   for (j=0; j<NumSameEntries; j++)
00410     To[j] = From[j];
00411     }
00412   // Do local permutation next
00413   if (NumPermuteIDs>0) {
00414   
00415     // Point entry case
00416     if (Case1) {
00417       
00418       for (j=0; j<NumPermuteIDs; j++) 
00419   To[PermuteToLIDs[j]] = From[PermuteFromLIDs[j]];
00420     }
00421     // constant element size case
00422     else if (Case2) {
00423       
00424       for (j=0; j<NumPermuteIDs; j++) {
00425   jj = MaxElementSize*PermuteToLIDs[j];
00426   jjj = MaxElementSize*PermuteFromLIDs[j];
00427     for (k=0; k<MaxElementSize; k++)
00428       To[jj+k] = From[jjj+k];
00429       }
00430     }
00431     
00432     // variable element size case
00433     else {
00434       
00435       for (j=0; j<NumPermuteIDs; j++) {
00436   jj = ToFirstPointInElementList[PermuteToLIDs[j]];
00437   jjj = FromFirstPointInElementList[PermuteFromLIDs[j]];
00438   int ElementSize = FromElementSizeList[PermuteFromLIDs[j]];
00439     for (k=0; k<ElementSize; k++)
00440       To[jj+k] = From[jjj+k];
00441       }
00442     }
00443   }
00444   return(0);
00445 }
00446 
00447 //=========================================================================
00448 // Perform any packing or preparation required for call to DoTransfer().
00449 int EpetraExt_BlockDiagMatrix::PackAndPrepare(const Epetra_SrcDistObject& Source,
00450                                            int NumExportIDs,
00451                                            int* ExportLIDs,
00452                                            int& LenExports,
00453                                            char*& Exports,
00454                                            int& SizeOfPacket,
00455                                            int* Sizes,
00456                                            bool & VarSizes,
00457                                            Epetra_Distributor& Distor){
00458   (void)Sizes;
00459   (void)VarSizes;
00460   (void)Distor;
00461   const EpetraExt_BlockDiagMatrix & A = dynamic_cast<const EpetraExt_BlockDiagMatrix &>(Source);
00462 
00463   int j, jj, k;
00464 
00465   double *From=A.Values();
00466   int MaxElementSize = DataMap().MaxElementSize();
00467   bool ConstantElementSize = DataMap().ConstantElementSize();
00468 
00469   int * FromFirstPointInElementList = 0;
00470   int * FromElementSizeList = 0;
00471 
00472   if (!ConstantElementSize) {
00473     FromFirstPointInElementList = A.DataMap().FirstPointInElementList();
00474     FromElementSizeList = A.DataMap().ElementSizeList();
00475   }
00476 
00477   SizeOfPacket = MaxElementSize * (int)sizeof(double); 
00478 
00479   if(NumExportIDs*SizeOfPacket>LenExports) {
00480     if (LenExports>0) delete [] Exports;
00481     LenExports = NumExportIDs*SizeOfPacket;
00482     Exports = new char[LenExports];
00483   }
00484 
00485   double * ptr;
00486 
00487   if (NumExportIDs>0) {
00488     ptr = (double *) Exports;
00489     
00490     // Point entry case
00491     if (MaxElementSize==1) for (j=0; j<NumExportIDs; j++) *ptr++ = From[ExportLIDs[j]];
00492 
00493     // constant element size case
00494     else if (ConstantElementSize) {      
00495       for (j=0; j<NumExportIDs; j++) {
00496   jj = MaxElementSize*ExportLIDs[j];
00497     for (k=0; k<MaxElementSize; k++)
00498       *ptr++ = From[jj+k];
00499       }
00500     }
00501     
00502     // variable element size case
00503     else {     
00504       SizeOfPacket = MaxElementSize;
00505       for (j=0; j<NumExportIDs; j++) {
00506   ptr = (double *) Exports + j*SizeOfPacket;
00507   jj = FromFirstPointInElementList[ExportLIDs[j]];
00508   int ElementSize = FromElementSizeList[ExportLIDs[j]];
00509     for (k=0; k<ElementSize; k++)
00510       *ptr++ = From[jj+k];
00511       }
00512     }
00513   }
00514 
00515   return(0);
00516 }
00517 
00518 
00519 //=========================================================================
00520 // Perform any unpacking and combining after call to DoTransfer().
00521 int EpetraExt_BlockDiagMatrix::UnpackAndCombine(const Epetra_SrcDistObject& Source, 
00522                                              int NumImportIDs,
00523                                              int* ImportLIDs, 
00524                                              int LenImports,
00525                                              char* Imports,
00526                                              int& SizeOfPacket, 
00527                                              Epetra_Distributor& Distor,
00528                                              Epetra_CombineMode CombineMode,
00529                                              const Epetra_OffsetIndex * Indexor){
00530   (void)Source;
00531   (void)LenImports;
00532   (void)Distor;
00533   (void)Indexor;
00534   int j, jj, k;
00535   
00536   if(    CombineMode != Add
00537       && CombineMode != Zero
00538       && CombineMode != Insert
00539       && CombineMode != Average
00540       && CombineMode != AbsMax )
00541     EPETRA_CHK_ERR(-1); //Unsupported CombinedMode, will default to Zero
00542 
00543   if (NumImportIDs<=0) return(0);
00544 
00545   double * To = Values_;
00546   int MaxElementSize = DataMap().MaxElementSize();
00547   bool ConstantElementSize = DataMap().ConstantElementSize();
00548 
00549   int * ToFirstPointInElementList = 0;
00550   int * ToElementSizeList = 0;
00551 
00552   if (!ConstantElementSize) {
00553     ToFirstPointInElementList = DataMap().FirstPointInElementList();
00554     ToElementSizeList = DataMap().ElementSizeList();
00555   }
00556   
00557   double * ptr;
00558   // Unpack it...
00559 
00560   ptr = (double *) Imports;
00561     
00562   // Point entry case
00563   if (MaxElementSize==1) {
00564       
00565       if (CombineMode==Add)
00566   for (j=0; j<NumImportIDs; j++) To[ImportLIDs[j]] += *ptr++; // Add to existing value
00567       else if(CombineMode==Insert)
00568   for (j=0; j<NumImportIDs; j++) To[ImportLIDs[j]] = *ptr++;
00569       else if(CombineMode==AbsMax)
00570         for (j=0; j<NumImportIDs; j++) {
00571     To[ImportLIDs[j]] = EPETRA_MAX( To[ImportLIDs[j]],std::abs(*ptr));
00572     ptr++;
00573   }
00574       // Note:  The following form of averaging is not a true average if more that one value is combined.
00575       //        This might be an issue in the future, but we leave this way for now.
00576       else if(CombineMode==Average)
00577   for (j=0; j<NumImportIDs; j++) {To[ImportLIDs[j]] += *ptr++; To[ImportLIDs[j]] /= 2;}
00578   }
00579 
00580   // constant element size case
00581 
00582   else if (ConstantElementSize) {
00583    
00584     if (CombineMode==Add) {
00585       for (j=0; j<NumImportIDs; j++) {
00586   jj = MaxElementSize*ImportLIDs[j];
00587     for (k=0; k<MaxElementSize; k++)
00588       To[jj+k] += *ptr++; // Add to existing value
00589       }
00590     }
00591     else if(CombineMode==Insert) {
00592       for (j=0; j<NumImportIDs; j++) {
00593   jj = MaxElementSize*ImportLIDs[j];
00594     for (k=0; k<MaxElementSize; k++)
00595       To[jj+k] = *ptr++;
00596       }
00597     }
00598     else if(CombineMode==AbsMax) {
00599       for (j=0; j<NumImportIDs; j++) {
00600   jj = MaxElementSize*ImportLIDs[j];
00601   for (k=0; k<MaxElementSize; k++) {
00602       To[jj+k] = EPETRA_MAX( To[jj+k], std::abs(*ptr));
00603       ptr++;
00604   }
00605       }
00606     }
00607     // Note:  The following form of averaging is not a true average if more that one value is combined.
00608     //        This might be an issue in the future, but we leave this way for now.
00609     else if(CombineMode==Average) {
00610       for (j=0; j<NumImportIDs; j++) {
00611   jj = MaxElementSize*ImportLIDs[j];
00612     for (k=0; k<MaxElementSize; k++)
00613       { To[jj+k] += *ptr++; To[jj+k] /= 2;}
00614       }
00615     }
00616   }
00617     
00618   // variable element size case
00619 
00620   else {
00621       
00622     SizeOfPacket = MaxElementSize;
00623 
00624     if (CombineMode==Add) {
00625       for (j=0; j<NumImportIDs; j++) {
00626   ptr = (double *) Imports + j*SizeOfPacket;
00627   jj = ToFirstPointInElementList[ImportLIDs[j]];
00628   int ElementSize = ToElementSizeList[ImportLIDs[j]];
00629     for (k=0; k<ElementSize; k++)
00630       To[jj+k] += *ptr++; // Add to existing value
00631       }
00632     }
00633     else  if(CombineMode==Insert){
00634       for (j=0; j<NumImportIDs; j++) {
00635   ptr = (double *) Imports + j*SizeOfPacket;
00636   jj = ToFirstPointInElementList[ImportLIDs[j]];
00637   int ElementSize = ToElementSizeList[ImportLIDs[j]];
00638     for (k=0; k<ElementSize; k++)
00639       To[jj+k] = *ptr++;
00640       }
00641     }
00642     else  if(CombineMode==AbsMax){
00643       for (j=0; j<NumImportIDs; j++) {
00644   ptr = (double *) Imports + j*SizeOfPacket;
00645   jj = ToFirstPointInElementList[ImportLIDs[j]];
00646   int ElementSize = ToElementSizeList[ImportLIDs[j]];
00647   for (k=0; k<ElementSize; k++) {
00648       To[jj+k] = EPETRA_MAX( To[jj+k], std::abs(*ptr));
00649       ptr++;
00650   }
00651       }
00652     }
00653     // Note:  The following form of averaging is not a true average if more that one value is combined.
00654     //        This might be an issue in the future, but we leave this way for now.
00655     else if(CombineMode==Average) {
00656       for (j=0; j<NumImportIDs; j++) {
00657   ptr = (double *) Imports + j*SizeOfPacket;
00658   jj = ToFirstPointInElementList[ImportLIDs[j]];
00659   int ElementSize = ToElementSizeList[ImportLIDs[j]];
00660     for (k=0; k<ElementSize; k++)
00661       { To[jj+k] += *ptr++; To[jj+k] /= 2;}
00662       }
00663     }
00664   }
00665   
00666   return(0);
00667 }
00668   
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines