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   if(X.NumVectors()!=Y.NumVectors())
00179     EPETRA_CHK_ERR(-1);
00180   if(!HasComputed_ && (ApplyMode_==AM_INVERT || ApplyMode_==AM_FACTOR))
00181     EPETRA_CHK_ERR(-2);
00182   
00183   //NTS: Add multivector support
00184   assert(X.NumVectors()==1);
00185 
00186   //NTS: MultiVector's MyLength and [] Operators are  "points" level operators
00187   //not a "block/element" level operators.
00188 
00189   const int *vlist=DataMap_->FirstPointInElementList();
00190   const int *xlist=Map().FirstPointInElementList();
00191   const int *blocksize=Map().ElementSizeList();
00192   
00193   if(ApplyMode_==AM_MULTIPLY || ApplyMode_==AM_INVERT){
00194     // Multiply & Invert mode have the same apply
00195     int NumBlocks=NumMyBlocks();
00196     for(int i=0;i<NumBlocks;i++){
00197       int Nb=blocksize[i];
00198       int vidx0=vlist[i];
00199       int xidx0=xlist[i];
00200 
00201       if(Nb==1) {
00202         // Optimize for size = 1
00203         Y[0][xidx0]=Values_[vidx0]*X[0][xidx0];
00204       }
00205       else if(Nb==2){
00206         // Optimize for size = 2
00207         Y[0][xidx0  ]=Values_[vidx0  ]*X[0][xidx0] + Values_[vidx0+2]*X[0][xidx0+1];
00208         Y[0][xidx0+1]=Values_[vidx0+1]*X[0][xidx0] + Values_[vidx0+3]*X[0][xidx0+1];
00209       }
00210       else{
00211         // "Large" Block - Use BLAS
00212         //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 
00213         GEMV('N',Nb,Nb,1.0,&Values_[vidx0],Nb,&X[0][xidx0],0.0,&Y[0][xidx0]);
00214       }
00215             
00216     }
00217   }
00218   else{
00219     // Factorization mode has a different apply
00220     int NumBlocks=NumMyBlocks();
00221     for(int i=0;i<NumBlocks;i++){
00222       int Nb=blocksize[i];
00223       int vidx0=vlist[i];
00224       int xidx0=xlist[i];      
00225       if(Nb==1) {
00226         // Optimize for size = 1 - use the inverse
00227         Y[0][xidx0]=Values_[vidx0]*X[0][xidx0];
00228       }
00229       else if(Nb==2){
00230         // Optimize for size = 2 - use the inverse
00231         Y[0][xidx0  ]=Values_[vidx0  ]*X[0][xidx0] + Values_[vidx0+2]*X[0][xidx0+1];
00232         Y[0][xidx0+1]=Values_[vidx0+1]*X[0][xidx0] + Values_[vidx0+3]*X[0][xidx0+1];
00233       }
00234       else{
00235         // "Large" Block - use LAPACK
00236         //    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 
00237         for(int j=0;j<Nb;j++) Y[0][xidx0+j]=X[0][xidx0+j];
00238         LAPACK.GETRS('N',Nb,1,&Values_[vidx0],Nb,&Pivots_[xidx0],&Y[0][xidx0],Nb,&info);
00239         if(info) EPETRA_CHK_ERR(info);
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   exit(-1);
00323 
00324   // NEEDS A FIX
00325 
00326   const EpetraExt_BlockDiagMatrix & A = dynamic_cast<const EpetraExt_BlockDiagMatrix &>(Source);
00327 
00328   double *From=A.Values();
00329   double *To = Values_;
00330 
00331   int * ToFirstPointInElementList = 0;
00332   int * FromFirstPointInElementList = 0;
00333   int * FromElementSizeList = 0;
00334   int MaxElementSize = Map().MaxElementSize();
00335   bool ConstantElementSize = Map().ConstantElementSize();
00336 
00337   if (!ConstantElementSize) {
00338     ToFirstPointInElementList =   Map().FirstPointInElementList();
00339     FromFirstPointInElementList = A.Map().FirstPointInElementList();
00340     FromElementSizeList = A.Map().ElementSizeList();
00341   }
00342   int j, jj, jjj, k;
00343   
00344   int NumSameEntries;
00345 
00346   bool Case1 = false;
00347   bool Case2 = false;
00348   // bool Case3 = false;
00349 
00350   if (MaxElementSize==1) {
00351     Case1 = true;
00352     NumSameEntries = NumSameIDs;
00353   }
00354   else if (ConstantElementSize) {
00355     Case2 = true;
00356     NumSameEntries = NumSameIDs * MaxElementSize;
00357   }
00358   else {
00359     // Case3 = true;
00360     NumSameEntries = FromFirstPointInElementList[NumSameIDs];
00361   }
00362 
00363   // Short circuit for the case where the source and target vector is the same.
00364   if (To==From) NumSameEntries = 0;
00365   
00366   // Do copy first
00367   if (NumSameIDs>0)
00368     if (To!=From) {
00369   for (j=0; j<NumSameEntries; j++)
00370     To[j] = From[j];
00371     }
00372   // Do local permutation next
00373   if (NumPermuteIDs>0) {
00374   
00375     // Point entry case
00376     if (Case1) {
00377       
00378       for (j=0; j<NumPermuteIDs; j++) 
00379   To[PermuteToLIDs[j]] = From[PermuteFromLIDs[j]];
00380     }
00381     // constant element size case
00382     else if (Case2) {
00383       
00384       for (j=0; j<NumPermuteIDs; j++) {
00385   jj = MaxElementSize*PermuteToLIDs[j];
00386   jjj = MaxElementSize*PermuteFromLIDs[j];
00387     for (k=0; k<MaxElementSize; k++)
00388       To[jj+k] = From[jjj+k];
00389       }
00390     }
00391     
00392     // variable element size case
00393     else {
00394       
00395       for (j=0; j<NumPermuteIDs; j++) {
00396   jj = ToFirstPointInElementList[PermuteToLIDs[j]];
00397   jjj = FromFirstPointInElementList[PermuteFromLIDs[j]];
00398   int ElementSize = FromElementSizeList[PermuteFromLIDs[j]];
00399     for (k=0; k<ElementSize; k++)
00400       To[jj+k] = From[jjj+k];
00401       }
00402     }
00403   }
00404   return(0);
00405 }
00406 
00407 //=========================================================================
00408 // Perform any packing or preparation required for call to DoTransfer().
00409 int EpetraExt_BlockDiagMatrix::PackAndPrepare(const Epetra_SrcDistObject& Source,
00410                                            int NumExportIDs,
00411                                            int* ExportLIDs,
00412                                            int& LenExports,
00413                                            char*& Exports,
00414                                            int& SizeOfPacket,
00415                                            int* Sizes,
00416                                            bool & VarSizes,
00417                                            Epetra_Distributor& Distor){
00418   (void)Sizes;
00419   (void)VarSizes;
00420   (void)Distor;
00421   const EpetraExt_BlockDiagMatrix & A = dynamic_cast<const EpetraExt_BlockDiagMatrix &>(Source);
00422 
00423   exit(-1);//FIX
00424 
00425   int j, jj, k;
00426 
00427   double *From=A.Values();
00428   int MaxElementSize = Map().MaxElementSize();
00429   bool ConstantElementSize = Map().ConstantElementSize();
00430 
00431   int * FromFirstPointInElementList = 0;
00432   int * FromElementSizeList = 0;
00433 
00434   if (!ConstantElementSize) {
00435     FromFirstPointInElementList = A.Map().FirstPointInElementList();
00436     FromElementSizeList = A.Map().ElementSizeList();
00437   }
00438 
00439   SizeOfPacket = MaxElementSize * (int)sizeof(double); 
00440 
00441   if(NumExportIDs*SizeOfPacket>LenExports) {
00442     if (LenExports>0) delete [] Exports;
00443     LenExports = NumExportIDs*SizeOfPacket;
00444     Exports = new char[LenExports];
00445   }
00446 
00447   double * ptr;
00448 
00449   if (NumExportIDs>0) {
00450     ptr = (double *) Exports;
00451     
00452     // Point entry case
00453     if (MaxElementSize==1) for (j=0; j<NumExportIDs; j++) *ptr++ = From[ExportLIDs[j]];
00454 
00455     // constant element size case
00456     else if (ConstantElementSize) {      
00457       for (j=0; j<NumExportIDs; j++) {
00458   jj = MaxElementSize*ExportLIDs[j];
00459     for (k=0; k<MaxElementSize; k++)
00460       *ptr++ = From[jj+k];
00461       }
00462     }
00463     
00464     // variable element size case
00465     else {     
00466       SizeOfPacket = MaxElementSize;
00467       for (j=0; j<NumExportIDs; j++) {
00468   ptr = (double *) Exports + j*SizeOfPacket;
00469   jj = FromFirstPointInElementList[ExportLIDs[j]];
00470   int ElementSize = FromElementSizeList[ExportLIDs[j]];
00471     for (k=0; k<ElementSize; k++)
00472       *ptr++ = From[jj+k];
00473       }
00474     }
00475   }
00476 
00477   return(0);
00478 }
00479 
00480 
00481 //=========================================================================
00482 // Perform any unpacking and combining after call to DoTransfer().
00483 int EpetraExt_BlockDiagMatrix::UnpackAndCombine(const Epetra_SrcDistObject& Source, 
00484                                              int NumImportIDs,
00485                                              int* ImportLIDs, 
00486                                              int LenImports,
00487                                              char* Imports,
00488                                              int& SizeOfPacket, 
00489                                              Epetra_Distributor& Distor,
00490                                              Epetra_CombineMode CombineMode,
00491                                              const Epetra_OffsetIndex * Indexor){
00492   (void)Source;
00493   (void)LenImports;
00494   (void)Distor;
00495   (void)Indexor;
00496   int j, jj, k;
00497 
00498   exit(-1);//FIX
00499   
00500   if(    CombineMode != Add
00501       && CombineMode != Zero
00502       && CombineMode != Insert
00503       && CombineMode != Average
00504       && CombineMode != AbsMax )
00505     EPETRA_CHK_ERR(-1); //Unsupported CombinedMode, will default to Zero
00506 
00507   if (NumImportIDs<=0) return(0);
00508 
00509   double * To = Values_;
00510   int MaxElementSize = Map().MaxElementSize();
00511   bool ConstantElementSize = Map().ConstantElementSize();
00512 
00513   int * ToFirstPointInElementList = 0;
00514   int * ToElementSizeList = 0;
00515 
00516   if (!ConstantElementSize) {
00517     ToFirstPointInElementList = Map().FirstPointInElementList();
00518     ToElementSizeList = Map().ElementSizeList();
00519   }
00520   
00521   double * ptr;
00522   // Unpack it...
00523 
00524   ptr = (double *) Imports;
00525     
00526   // Point entry case
00527   if (MaxElementSize==1) {
00528       
00529       if (CombineMode==Add)
00530   for (j=0; j<NumImportIDs; j++) To[ImportLIDs[j]] += *ptr++; // Add to existing value
00531       else if(CombineMode==Insert)
00532   for (j=0; j<NumImportIDs; j++) To[ImportLIDs[j]] = *ptr++;
00533       else if(CombineMode==AbsMax)
00534         for (j=0; j<NumImportIDs; j++) {
00535     To[ImportLIDs[j]] = EPETRA_MAX( To[ImportLIDs[j]],std::abs(*ptr));
00536     ptr++;
00537   }
00538       // Note:  The following form of averaging is not a true average if more that one value is combined.
00539       //        This might be an issue in the future, but we leave this way for now.
00540       else if(CombineMode==Average)
00541   for (j=0; j<NumImportIDs; j++) {To[ImportLIDs[j]] += *ptr++; To[ImportLIDs[j]] /= 2;}
00542   }
00543 
00544   // constant element size case
00545 
00546   else if (ConstantElementSize) {
00547    
00548     if (CombineMode==Add) {
00549       for (j=0; j<NumImportIDs; j++) {
00550   jj = MaxElementSize*ImportLIDs[j];
00551     for (k=0; k<MaxElementSize; k++)
00552       To[jj+k] += *ptr++; // Add to existing value
00553       }
00554     }
00555     else if(CombineMode==Insert) {
00556       for (j=0; j<NumImportIDs; j++) {
00557   jj = MaxElementSize*ImportLIDs[j];
00558     for (k=0; k<MaxElementSize; k++)
00559       To[jj+k] = *ptr++;
00560       }
00561     }
00562     else if(CombineMode==AbsMax) {
00563       for (j=0; j<NumImportIDs; j++) {
00564   jj = MaxElementSize*ImportLIDs[j];
00565   for (k=0; k<MaxElementSize; k++) {
00566       To[jj+k] = EPETRA_MAX( To[jj+k], std::abs(*ptr));
00567       ptr++;
00568   }
00569       }
00570     }
00571     // Note:  The following form of averaging is not a true average if more that one value is combined.
00572     //        This might be an issue in the future, but we leave this way for now.
00573     else if(CombineMode==Average) {
00574       for (j=0; j<NumImportIDs; j++) {
00575   jj = MaxElementSize*ImportLIDs[j];
00576     for (k=0; k<MaxElementSize; k++)
00577       { To[jj+k] += *ptr++; To[jj+k] /= 2;}
00578       }
00579     }
00580   }
00581     
00582   // variable element size case
00583 
00584   else {
00585       
00586     SizeOfPacket = MaxElementSize;
00587 
00588     if (CombineMode==Add) {
00589       for (j=0; j<NumImportIDs; j++) {
00590   ptr = (double *) Imports + j*SizeOfPacket;
00591   jj = ToFirstPointInElementList[ImportLIDs[j]];
00592   int ElementSize = ToElementSizeList[ImportLIDs[j]];
00593     for (k=0; k<ElementSize; k++)
00594       To[jj+k] += *ptr++; // Add to existing value
00595       }
00596     }
00597     else  if(CombineMode==Insert){
00598       for (j=0; j<NumImportIDs; j++) {
00599   ptr = (double *) Imports + j*SizeOfPacket;
00600   jj = ToFirstPointInElementList[ImportLIDs[j]];
00601   int ElementSize = ToElementSizeList[ImportLIDs[j]];
00602     for (k=0; k<ElementSize; k++)
00603       To[jj+k] = *ptr++;
00604       }
00605     }
00606     else  if(CombineMode==AbsMax){
00607       for (j=0; j<NumImportIDs; j++) {
00608   ptr = (double *) Imports + j*SizeOfPacket;
00609   jj = ToFirstPointInElementList[ImportLIDs[j]];
00610   int ElementSize = ToElementSizeList[ImportLIDs[j]];
00611   for (k=0; k<ElementSize; k++) {
00612       To[jj+k] = EPETRA_MAX( To[jj+k], std::abs(*ptr));
00613       ptr++;
00614   }
00615       }
00616     }
00617     // Note:  The following form of averaging is not a true average if more that one value is combined.
00618     //        This might be an issue in the future, but we leave this way for now.
00619     else if(CombineMode==Average) {
00620       for (j=0; j<NumImportIDs; j++) {
00621   ptr = (double *) Imports + j*SizeOfPacket;
00622   jj = ToFirstPointInElementList[ImportLIDs[j]];
00623   int ElementSize = ToElementSizeList[ImportLIDs[j]];
00624     for (k=0; k<ElementSize; k++)
00625       { To[jj+k] += *ptr++; To[jj+k] /= 2;}
00626       }
00627     }
00628   }
00629   
00630   return(0);
00631 }
00632   

Generated on Wed May 12 21:24:45 2010 for EpetraExt by  doxygen 1.4.7