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