EpetraExt_PointToBlockDiagPermute.cpp

Go to the documentation of this file.
00001 #include "EpetraExt_PointToBlockDiagPermute.h"
00002 #include "EpetraExt_BlockDiagMatrix.h"
00003 #include "Epetra_MultiVector.h"
00004 #include "Epetra_Import.h"
00005 #include "Epetra_Export.h"
00006 #include "Epetra_Comm.h"
00007 
00008 #include <stdio.h>
00009 #include <fstream>
00010 
00011 #define MAX(x,y) ((x)>(y)?(x):(y))
00012 
00013 //=========================================================================
00014 EpetraExt_PointToBlockDiagPermute::EpetraExt_PointToBlockDiagPermute(const Epetra_CrsMatrix& MAT)
00015   :Epetra_DistObject(MAT.RowMap()),
00016    Matrix_(&MAT),
00017    PurelyLocalMode_(true),
00018    NumBlocks_(0),
00019    Blockstart_(0),
00020    Blockids_(0),
00021    BDMap_(0),
00022    CompatibleMap_(0),
00023    BDMat_(0),
00024    Importer_(0),
00025    Exporter_(0),
00026    ImportVector_(0),
00027    ExportVector_(0)
00028 {
00029 
00030 }
00031   
00032 //=========================================================================  
00033 // Destructor
00034 EpetraExt_PointToBlockDiagPermute::~EpetraExt_PointToBlockDiagPermute()
00035 {
00036   if(BDMap_) delete BDMap_;
00037   if(CompatibleMap_) delete CompatibleMap_;
00038   if(BDMat_) delete BDMat_;
00039   if(Importer_) delete Importer_;
00040   if(Exporter_) delete Exporter_;
00041   if(ImportVector_) delete ImportVector_;
00042   if(ExportVector_) delete ExportVector_;
00043 }
00044 
00045 
00046 
00047 //=========================================================================  
00048 // Set list
00049 int EpetraExt_PointToBlockDiagPermute::SetParameters(Teuchos::ParameterList & List){
00050   List_=List;
00051 
00052   // Local vs. global ids & mode
00053   NumBlocks_=List_.get("number of local blocks",0);
00054   if(NumBlocks_ <= 0) EPETRA_CHK_ERR(-1);
00055   
00056   Blockstart_=List_.get("block start index",(int*)0);
00057   if(!NumBlocks_) EPETRA_CHK_ERR(-2);
00058   
00059   Blockids_=List_.get("block entry lids",(int*)0);
00060   if(Blockids_)
00061     PurelyLocalMode_=true;
00062   else{
00063     Blockids_=List_.get("block entry gids",(int*)0);
00064     PurelyLocalMode_=false;
00065   }
00066   if(!Blockids_) EPETRA_CHK_ERR(-3);
00067   
00068   return 0;
00069 }
00070 
00071 
00072 //=========================================================================  
00073 // Extracts the block-diagonal, builds maps, etc.
00074 int EpetraExt_PointToBlockDiagPermute::Compute(){
00075   int rv=ExtractBlockDiagonal();
00076   return rv;
00077 }
00078 
00079 //=========================================================================  
00080 // Returns the result of a Epetra_Operator applied to a Epetra_MultiVector X in Y.
00081 int EpetraExt_PointToBlockDiagPermute::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
00082   return -1;
00083 
00084 }
00085 
00086 //=========================================================================  
00087 // Returns the result of a Epetra_Operator inverse applied to an Epetra_MultiVector X in Y.
00088 int EpetraExt_PointToBlockDiagPermute::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
00089   // Stuff borrowed from Epetra_CrsMatrix
00090   int NumVectors = X.NumVectors();
00091   if (NumVectors!=Y.NumVectors()) {
00092     EPETRA_CHK_ERR(-2); // Need same number of vectors in each MV
00093   }
00094 
00095   const Epetra_MultiVector *Xp=&X;
00096   Epetra_MultiVector *Yp=&Y;
00097 
00098   // Allocate temp workspace if X==Y and there are no imports or exports
00099   Epetra_MultiVector * Xcopy = 0;
00100   if (&X==&Y && Importer_==0 && Exporter_==0) {
00101     Xcopy = new Epetra_MultiVector(X);
00102     Xp=Xcopy;
00103   }
00104   
00105   UpdateImportVector(NumVectors); // Make sure Import and Export Vectors are compatible
00106   UpdateExportVector(NumVectors);
00107 
00108   // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
00109   if (Importer_){
00110     EPETRA_CHK_ERR(ImportVector_->Import(X, *Importer_, Insert));
00111     Xp=ImportVector_;
00112   }
00113   
00114   // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
00115   if (Exporter_) {
00116     Yp=ExportVector_;
00117   }
00118   
00119   // Do the matvec 
00120   BDMat_->ApplyInverse(*Xp,*Yp);
00121 
00122   // Export if needed
00123   if (Exporter_) {
00124     Y.PutScalar(0.0);  // Make sure target is zero
00125     Y.Export(*ExportVector_, *Exporter_, Add); // Fill Y with Values from export vector
00126   }
00127   
00128   // Cleanup
00129   if(Xcopy) {
00130     delete Xcopy;
00131     EPETRA_CHK_ERR(1); // Return positive code to alert the user about needing extra copy of X
00132     return 1;
00133   }
00134 
00135   return 0;
00136 }
00137 
00138 //=========================================================================  
00139 // Print method
00140 void EpetraExt_PointToBlockDiagPermute::Print(ostream& os) const{
00141   if(Importer_) cout<<*Importer_<<endl;
00142   if(Exporter_) cout<<*Exporter_<<endl;
00143   if(BDMat_) cout<<*BDMat_<<endl;
00144 }
00145 
00146 //=========================================================================  
00147 // Pulls the block diagonal of the matrix and then builds the BDMat_
00148 int EpetraExt_PointToBlockDiagPermute::ExtractBlockDiagonal(){  
00149   int i,j;
00150   std::vector<int> ExtRows;
00151   int* LocalColIDS=0;
00152   int ExtSize=0,ExtCols=0,MainCols=0;
00153   int Nrows=Matrix_->NumMyRows();
00154   int *l2b,*block_offset,*l2blockid;
00155   const Epetra_Map &RowMap=Matrix_->RowMap();
00156   int *bsize=new int[NumBlocks_]; 
00157   int index,col,row_in_block,col_in_block,length,*colind;
00158   double *values;
00159 
00160   bool verbose=(bool)(List_.get("output",0) > 0);
00161   
00162   // Compute block size lists
00163   for(i=0;i<NumBlocks_;i++) 
00164     bsize[i]=Blockstart_[i+1]-Blockstart_[i];
00165   
00166   // Use the ScanSum function to compute a prefix sum of the number of points
00167   int MyMinGID;
00168   Matrix_->Comm().ScanSum(&NumBlocks_,&MyMinGID, 1);
00169   MyMinGID-=NumBlocks_;
00170   int *MyBlockGIDs=new int[NumBlocks_];
00171   for(i=0;i<NumBlocks_;i++)
00172     MyBlockGIDs[i]=MyMinGID+i;
00173 
00174   BDMap_=new Epetra_BlockMap(-1,NumBlocks_,MyBlockGIDs,bsize,0,Matrix_->Comm());
00175   
00176   // Allocate the Epetra_DistBlockMatrix
00177   BDMat_=new EpetraExt_BlockDiagMatrix(*BDMap_,true);  
00178   
00179   // First check to see if we can switch back to PurelyLocalMode_ if we're not
00180   // in it
00181   if(!PurelyLocalMode_){
00182     // Find the non-local row IDs
00183     ExtSize=MAX(0,Blockstart_[NumBlocks_]-RowMap.NumMyElements());// estimate
00184     ExtRows.reserve(ExtSize);
00185     
00186     for(i=0;i<Blockstart_[NumBlocks_];i++){
00187       if(RowMap.LID(Blockids_[i])==-1)
00188         ExtRows.push_back(Blockids_[i]);
00189     }
00190 
00191     // Switch to PurelyLocalMode_ if we never need GIDs
00192     int size_sum;
00193     ExtSize=ExtRows.size();   
00194     RowMap.Comm().SumAll(&ExtSize,&size_sum,1);
00195     if(size_sum==0){
00196       if(verbose && !Matrix_->Comm().MyPID()) printf("EpetraExt_PointToBlockDiagPermute: Switching to purely local mode\n");
00197       PurelyLocalMode_=true;
00198       for(i=0;i<Blockstart_[NumBlocks_];i++){
00199         Blockids_[i]=RowMap.LID(Blockids_[i]);
00200       }     
00201     }
00202   }
00203   
00204   if(PurelyLocalMode_){
00205     /*******************************************************/
00206     // Allocations
00207     l2b=new int[Nrows];
00208     block_offset=new int[Nrows];
00209 
00210     // Build the local-id-to-block-id list, block offset list    
00211     for(i=0;i<NumBlocks_;i++) {
00212       for(j=Blockstart_[i];j<Blockstart_[i+1];j++){
00213         block_offset[Blockids_[j]]=j-Blockstart_[i];
00214         l2b[Blockids_[j]]=i;
00215       }
00216     }
00217     
00218     // Copy the data to the EpetraExt_BlockDiagMatrix        
00219     for(i=0;i<Nrows;i++) {
00220       int block_num = l2b[i];
00221       if(block_num>=0 && block_num<NumBlocks_) {
00222         row_in_block=block_offset[i];
00223         Matrix_->ExtractMyRowView(i,length,values,colind);
00224         int Nnz=0;
00225         for (j = 0; j < length; j++) {
00226           col = colind[j];
00227           if(col < Nrows && l2b[col]==block_num){
00228             Nnz++;
00229             col_in_block = block_offset[col];
00230             index = col_in_block * bsize[block_num] + row_in_block;
00231             (*BDMat_)[block_num][index]=values[j];
00232           }
00233         }
00234         /* Handle the case of a zero row. */
00235         /* By just putting a 1 on the diagonal. */
00236         if (Nnz == 0) {
00237           index = row_in_block * bsize[block_num] + row_in_block;
00238           (*BDMat_)[block_num][index]=1.0;
00239         }
00240       }
00241     }
00242 
00243     // Build the compatible map for import/export
00244     l2blockid=new int[Nrows];
00245     for(i=0;i<Nrows;i++) 
00246       l2blockid[Blockstart_[l2b[i]]+block_offset[i]]=Matrix_->RowMap().GID(i);
00247 
00248     // Build the Compatible Map, Import/Export Objects
00249     CompatibleMap_=new Epetra_Map(-1,Nrows,l2blockid,0,Matrix_->Comm());
00250     
00251   }
00252   else{
00253     /*******************************************************/      
00254     // Do the import to grab matrix entries
00255     // Allocate temporaries for import
00256     Epetra_Map TmpMap(-1,ExtSize, &ExtRows[0],0,Matrix_->Comm());; 
00257     Epetra_CrsMatrix TmpMatrix(Copy,TmpMap,0);
00258     Epetra_Import TmpImporter(TmpMap,RowMap);
00259 
00260     TmpMatrix.Import(*Matrix_,TmpImporter,Insert);
00261     TmpMatrix.FillComplete();
00262       
00263     ExtCols=TmpMatrix.NumMyCols();
00264     MainCols=Matrix_->NumMyCols();
00265         
00266 
00267     // Build the column reidex - main matrix
00268     LocalColIDS=new int[MainCols+ExtCols];      
00269     for(i=0;i<MainCols;i++){
00270       int GID=Matrix_->GCID(i);
00271       int MainLID=RowMap.LID(GID);
00272       if(MainLID!=-1) LocalColIDS[i]=MainLID;
00273       else{
00274         int ExtLID=TmpMatrix.LRID(GID);
00275         if(ExtLID!=-1) LocalColIDS[i]=Nrows+ExtLID;
00276         else LocalColIDS[i]=-1;
00277       }
00278     }
00279 
00280     // Build the column reidex - ext matrix
00281     for(i=0;i<ExtCols;i++){
00282       int GID=TmpMatrix.GCID(i);
00283       int MainLID=RowMap.LID(GID);
00284       if(MainLID!=-1) LocalColIDS[MainCols+i]=MainLID;
00285       else{
00286         int ExtLID=TmpMatrix.LRID(GID);          
00287         if(ExtLID!=-1) LocalColIDS[MainCols+i]=Nrows+ExtLID;
00288         else LocalColIDS[MainCols+i]=-1;
00289       }
00290     }
00291 
00292     // Allocations
00293     l2b=new int[Nrows+ExtSize];
00294     block_offset=new int[Nrows+ExtSize];  
00295   
00296     // Build l2b/block_offset with the expanded local index
00297     //NTS: Uses the same ordering of operation as the above routine, which is why
00298     //it works.
00299     for(i=0;i<Nrows+ExtSize;i++) block_offset[i]=l2b[i]=-1;    
00300     int ext_idx=0;
00301     for(i=0;i<NumBlocks_;i++) {
00302       for(j=Blockstart_[i];j<Blockstart_[i+1];j++){
00303         int LID=RowMap.LID(Blockids_[j]);
00304         if(LID==-1) {LID=Nrows+ext_idx;ext_idx++;}
00305         block_offset[LID]=j-Blockstart_[i];
00306         l2b[LID]=i;          
00307       }
00308     }
00309     
00310     // Copy the data to the EpetraExt_BlockDiagMatrix from Matrix_
00311     for(i=0;i<Nrows;i++) {
00312       int block_num = l2b[i];
00313       if(block_num>=0 && block_num<NumBlocks_) {
00314         row_in_block=block_offset[i];
00315         Matrix_->ExtractMyRowView(i,length,values,colind);
00316         int Nnz=0;
00317         for (j = 0; j < length; j++) {
00318           col = LocalColIDS[colind[j]];
00319           if(col!=-1 && l2b[col]==block_num){
00320             Nnz++;
00321             col_in_block = block_offset[col];
00322             index = col_in_block * bsize[block_num] + row_in_block;
00323             (*BDMat_)[block_num][index]=values[j];
00324           }
00325         }
00326         /* Handle the case of a zero row. */
00327         /* By just putting a 1 on the diagonal. */
00328         if (Nnz == 0) {
00329           index = row_in_block * bsize[block_num] + row_in_block;
00330           (*BDMat_)[block_num][index]=values[j];
00331         }
00332       }
00333     }
00334     
00335     // Copy the data to the EpetraExt_BlockDiagMarix from TmpMatrix
00336     for(i=0;i<ExtSize;i++) {
00337       int block_num = l2b[Nrows+i];
00338       if(block_num>=0 && block_num<NumBlocks_) {
00339         row_in_block=block_offset[Nrows+i];
00340         TmpMatrix.ExtractMyRowView(i,length,values,colind);
00341         int Nnz=0;
00342         for (j = 0; j < length; j++) {
00343           col = LocalColIDS[MainCols+colind[j]];
00344           if(col!=-1 && l2b[col]==block_num){
00345             Nnz++;
00346             col_in_block = block_offset[col];
00347             index = col_in_block * bsize[block_num] + row_in_block;
00348             (*BDMat_)[block_num][index]=values[j];
00349           }
00350         }
00351         /* Handle the case of a zero row. */
00352         /* By just putting a 1 on the diagonal. */
00353         if (Nnz == 0) {
00354           index = row_in_block * bsize[block_num] + row_in_block;
00355           (*BDMat_)[block_num][index]=1.0;
00356         }
00357       }      
00358     }
00359     
00360     // Build the compatible map for import/export
00361     l2blockid=new int[Blockstart_[NumBlocks_]];
00362     for(i=0;i<Nrows;i++){
00363       int bkid=l2b[i];
00364       if(bkid>-1) l2blockid[Blockstart_[bkid]+block_offset[i]]=RowMap.GID(i);
00365     }
00366     // NTS: This is easier - if we imported it, we know we need it
00367     for(i=0;i<ExtSize;i++)
00368       l2blockid[Blockstart_[l2b[Nrows+i]]+block_offset[Nrows+i]]=TmpMatrix.GRID(i);
00369     
00370     // Build the Compatible Map, Import/Export Objects
00371     CompatibleMap_=new Epetra_Map(-1,Blockstart_[NumBlocks_],l2blockid,0,Matrix_->Comm());
00372         
00373   }//end else
00374 
00375   // Set BDMat_'s parameter list and compute!
00376   Teuchos::ParameterList dummy,inList;
00377   inList=List_.get("blockdiagmatrix: list",dummy);
00378   BDMat_->SetParameters(inList);  
00379   BDMat_->Compute();
00380   
00381   // Build importer/exporter
00382   if(!CompatibleMap_->SameAs(Matrix_->DomainMap()))
00383     Importer_ = new Epetra_Import(*CompatibleMap_,Matrix_->DomainMap());
00384   if(!CompatibleMap_->SameAs(Matrix_->RangeMap()))
00385     Exporter_ = new Epetra_Export(*CompatibleMap_,Matrix_->RangeMap());
00386   
00387   // Cleanup
00388   delete [] LocalColIDS;
00389   delete [] block_offset;
00390   delete [] l2b;
00391   delete [] l2blockid;
00392   delete [] bsize;
00393   delete [] MyBlockGIDs;
00394 
00395   return 0;
00396 }
00397 
00398 
00399 
00400 //=======================================================================================================
00401 void EpetraExt_PointToBlockDiagPermute::UpdateImportVector(int NumVectors) const {    
00402   if(Importer_ != 0) {
00403     if(ImportVector_ != 0) {
00404       if(ImportVector_->NumVectors() != NumVectors) { 
00405   delete ImportVector_; 
00406   ImportVector_= 0;
00407       }
00408     }
00409     if(ImportVector_ == 0) 
00410       ImportVector_ = new Epetra_MultiVector(*CompatibleMap_,NumVectors); // Create import vector if needed
00411   }
00412   return;
00413 }
00414 //=======================================================================================================
00415 void EpetraExt_PointToBlockDiagPermute::UpdateExportVector(int NumVectors) const {    
00416   if(Exporter_ != 0) {
00417     if(ExportVector_ != 0) {
00418       if(ExportVector_->NumVectors() != NumVectors) { 
00419   delete ExportVector_; 
00420   ExportVector_= 0;
00421       }
00422     }
00423     if(ExportVector_ == 0) 
00424       ExportVector_ = new Epetra_MultiVector(*CompatibleMap_,NumVectors); // Create Export vector if needed
00425   }
00426   return;
00427 }
00428 
00429 
00430 
00431 
00432 //=========================================================================  
00433 int EpetraExt_PointToBlockDiagPermute::Import(const Epetra_SrcDistObject& A, const Epetra_Import& Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor){
00434  return -1;
00435 }
00436 
00437 //=========================================================================  
00438 int EpetraExt_PointToBlockDiagPermute::Import(const Epetra_SrcDistObject& A, const Epetra_Export& Exporter, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor){
00439  return -1;
00440 }
00441 
00442 //=========================================================================  
00443 int EpetraExt_PointToBlockDiagPermute::Export(const Epetra_SrcDistObject& A, const Epetra_Import & Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor){
00444  return -1;
00445 }
00446 
00447 //=========================================================================  
00448 int EpetraExt_PointToBlockDiagPermute::Export(const Epetra_SrcDistObject& A, const Epetra_Export& Exporter, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor){
00449  return -1;
00450 }
00451 
00452 //=========================================================================  
00453 // Allows the source and target (\e this) objects to be compared for compatibility, return nonzero if not.
00454 int EpetraExt_PointToBlockDiagPermute::CheckSizes(const Epetra_SrcDistObject& Source){
00455   return -1;
00456 }
00457 
00458 //=========================================================================  
00459 // Perform ID copies and permutations that are on processor.
00460 int EpetraExt_PointToBlockDiagPermute::CopyAndPermute(const Epetra_SrcDistObject& Source,
00461                    int NumSameIDs, 
00462                    int NumPermuteIDs,
00463                    int * PermuteToLIDs,
00464                    int * PermuteFromLIDs,
00465                    const Epetra_OffsetIndex * Indexor){
00466   return -1;
00467 }
00468 
00469 //=========================================================================  
00470 // Perform any packing or preparation required for call to DoTransfer().
00471 int EpetraExt_PointToBlockDiagPermute::PackAndPrepare(const Epetra_SrcDistObject& Source,
00472                    int NumExportIDs,
00473                    int* ExportLIDs,
00474                    int& LenExports,
00475                    char*& Exports,
00476                    int& SizeOfPacket,
00477                    int* Sizes,
00478                    bool & VarSizes,
00479                    Epetra_Distributor& Distor){
00480   return -1;
00481 }
00482 
00483 //=========================================================================  
00484 // Perform any unpacking and combining after call to DoTransfer().
00485 int EpetraExt_PointToBlockDiagPermute::UnpackAndCombine(const Epetra_SrcDistObject& Source, 
00486                      int NumImportIDs,
00487                      int* ImportLIDs, 
00488                      int LenImports,
00489                      char* Imports,
00490                      int& SizeOfPacket, 
00491                      Epetra_Distributor& Distor,
00492                      Epetra_CombineMode CombineMode,
00493                      const Epetra_OffsetIndex * Indexor){
00494   return -1;
00495 }
00496 
00497 

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