EpetraExt Development
EpetraExt_PointToBlockDiagPermute.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 "Epetra_ConfigDefs.h"
00045 #include "EpetraExt_PointToBlockDiagPermute.h"
00046 #include "EpetraExt_BlockDiagMatrix.h"
00047 #include "Epetra_MultiVector.h"
00048 #include "Epetra_Import.h"
00049 #include "Epetra_Export.h"
00050 #include "Epetra_Comm.h"
00051 
00052 #include <stdio.h>
00053 #include <fstream>
00054 
00055 #define MAX(x,y) ((x)>(y)?(x):(y))
00056 
00057 //=========================================================================
00058 EpetraExt_PointToBlockDiagPermute::EpetraExt_PointToBlockDiagPermute(const Epetra_CrsMatrix& MAT)
00059   :Epetra_DistObject(MAT.RowMap()),
00060    Matrix_(&MAT),
00061    PurelyLocalMode_(true),
00062    ContiguousBlockMode_(false),
00063    ContiguousBlockSize_(0),
00064    NumBlocks_(0),
00065    Blockstart_(0),
00066    Blockids_int_(0),
00067 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00068    Blockids_LL_(0),
00069 #endif
00070    BDMap_(0),
00071    CompatibleMap_(0),
00072    BDMat_(0),
00073    Importer_(0),
00074    Exporter_(0),
00075    ImportVector_(0),
00076    ExportVector_(0)
00077 {
00078 
00079 }
00080   
00081 //=========================================================================  
00082 // Destructor
00083 EpetraExt_PointToBlockDiagPermute::~EpetraExt_PointToBlockDiagPermute()
00084 {
00085   if(BDMap_) delete BDMap_;
00086   if(CompatibleMap_) delete CompatibleMap_;
00087   if(BDMat_) delete BDMat_;
00088   if(Importer_) delete Importer_;
00089   if(Exporter_) delete Exporter_;
00090   if(ImportVector_) delete ImportVector_;
00091   if(ExportVector_) delete ExportVector_;
00092 }
00093 
00094 //=========================================================================  
00095 // Set list
00096 template<typename int_type>
00097 int EpetraExt_PointToBlockDiagPermute::TSetParameters(Teuchos::ParameterList & List){
00098   List_=List;
00099 
00100   // Check for contiguous blocking first
00101   ContiguousBlockSize_=List_.get("contiguous block size",0);
00102   if(ContiguousBlockSize_!=0){
00103     ContiguousBlockMode_=true;
00104     PurelyLocalMode_=false;
00105   }  
00106   
00107   // Local vs. global ids & mode
00108   NumBlocks_=List_.get("number of local blocks",0);  
00109   Blockstart_=List_.get("block start index",(int*)0);    
00110   Blockids_ref<int>()=List_.get("block entry lids",(int*)0);
00111   if(Blockids_const_ptr<int>())
00112     PurelyLocalMode_=true;
00113   else{
00114     Blockids_ref<int_type>()=List_.get("block entry gids",(int_type*)0);
00115     if(Blockids_const_ptr<int_type>()) {
00116       PurelyLocalMode_=false;
00117       if(sizeof(int) != sizeof(int_type)) {
00118         Blockids_ref<int>() = new int[Matrix_->RowMap().NumMyElements()];
00119       }
00120     }
00121   }
00122   
00123   // Sanity checks
00124   if(ContiguousBlockMode_){
00125     // Can't use contiguous at the same time as the other modes
00126     if(NumBlocks_ || Blockstart_ || Blockids_const_ptr<int>() || Blockids_const_ptr<int_type>()) EPETRA_CHK_ERR(-4);
00127   }
00128   else {
00129     if(NumBlocks_ <= 0) EPETRA_CHK_ERR(-1);
00130     if(!Blockstart_) EPETRA_CHK_ERR(-2);
00131     if(!Blockids_const_ptr<int>() && !Blockids_const_ptr<int_type>()) EPETRA_CHK_ERR(-3);
00132   }
00133   
00134   return 0;
00135 }
00136 
00137 int EpetraExt_PointToBlockDiagPermute::SetParameters(Teuchos::ParameterList & List){
00138 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00139   if(Matrix_->RowMap().GlobalIndicesInt()) {
00140   return TSetParameters<int>(List);
00141   }
00142   else
00143 #endif
00144 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00145   if(Matrix_->RowMap().GlobalIndicesLongLong()) {
00146   return TSetParameters<long long>(List);
00147   }
00148   else
00149 #endif
00150     throw "EpetraExt_PointToBlockDiagPermute::SetParameters: ERROR, GlobalIndices type unknown.";
00151 }
00152 
00153 //=========================================================================  
00154 // Extracts the block-diagonal, builds maps, etc.
00155 int EpetraExt_PointToBlockDiagPermute::Compute(){
00156   int rv=ExtractBlockDiagonal();
00157   return rv;
00158 }
00159 
00160 //=========================================================================  
00161 // Returns the result of a Epetra_Operator applied to a Epetra_MultiVector X in Y.
00162 int EpetraExt_PointToBlockDiagPermute::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
00163   return -1;
00164 
00165 }
00166 
00167 //=========================================================================  
00168 // Returns the result of a Epetra_Operator inverse applied to an Epetra_MultiVector X in Y.
00169 int EpetraExt_PointToBlockDiagPermute::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
00170   // Stuff borrowed from Epetra_CrsMatrix
00171   int NumVectors = X.NumVectors();
00172   if (NumVectors!=Y.NumVectors()) {
00173     EPETRA_CHK_ERR(-2); // Need same number of vectors in each MV
00174   }
00175 
00176   const Epetra_MultiVector *Xp=&X;
00177   Epetra_MultiVector *Yp=&Y;
00178 
00179   // Allocate temp workspace if X==Y and there are no imports or exports
00180   Epetra_MultiVector * Xcopy = 0;
00181   if (&X==&Y && Importer_==0 && Exporter_==0) {
00182     Xcopy = new Epetra_MultiVector(X);
00183     Xp=Xcopy;
00184   }
00185   
00186   UpdateImportVector(NumVectors); // Make sure Import and Export Vectors are compatible
00187   UpdateExportVector(NumVectors);
00188 
00189   // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
00190   if (Importer_){
00191     EPETRA_CHK_ERR(ImportVector_->Import(X, *Importer_, Insert));
00192     Xp=ImportVector_;
00193   }
00194   
00195   // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
00196   if (Exporter_) {
00197     Yp=ExportVector_;
00198   }
00199   
00200   // Do the matvec 
00201   BDMat_->ApplyInverse(*Xp,*Yp);
00202 
00203   // Export if needed
00204   if (Exporter_) {
00205     Y.PutScalar(0.0);  // Make sure target is zero
00206     Y.Export(*ExportVector_, *Exporter_, Add); // Fill Y with Values from export vector
00207   }
00208   
00209   // Cleanup
00210   if(Xcopy) {
00211     delete Xcopy;
00212     EPETRA_CHK_ERR(1); // Return positive code to alert the user about needing extra copy of X
00213     return 1;
00214   }
00215 
00216   return 0;
00217 }
00218 
00219 //=========================================================================  
00220 // Print method
00221 void EpetraExt_PointToBlockDiagPermute::Print(std::ostream& os) const{
00222   if(Importer_) std::cout<<*Importer_<<std::endl;
00223   if(Exporter_) std::cout<<*Exporter_<<std::endl;
00224   if(BDMat_) std::cout<<*BDMat_<<std::endl;
00225 }
00226 
00227 //=========================================================================  
00228 // Pulls the block diagonal of the matrix and then builds the BDMat_
00229 template<typename int_type>
00230 int EpetraExt_PointToBlockDiagPermute::TExtractBlockDiagonal(){  
00231   int i,j;
00232   std::vector<int_type> ExtRows;
00233   int* LocalColIDS=0;
00234   int ExtSize=0,ExtCols=0,MainCols=0;
00235   int Nrows=Matrix_->NumMyRows();
00236   int *l2b,*block_offset;
00237   int_type *l2blockid;
00238   const Epetra_Map &RowMap=Matrix_->RowMap();
00239   int index,col,row_in_block,col_in_block,length,*colind;
00240   double *values;
00241 
00242   bool verbose=(bool)(List_.get("output",0) > 0);
00243   
00244   // Contiguous Setup
00245   SetupContiguousMode();
00246 
00247   // Compute block size lists
00248   int *bsize=new int[NumBlocks_]; 
00249   for(i=0;i<NumBlocks_;i++) 
00250     bsize[i]=Blockstart_[i+1]-Blockstart_[i];
00251   
00252   // Use the ScanSum function to compute a prefix sum of the number of points
00253   int_type MyMinGID;
00254   int_type NumBlocks_int_type = NumBlocks_;
00255   Matrix_->Comm().ScanSum(&NumBlocks_int_type,&MyMinGID, 1);
00256   MyMinGID-=NumBlocks_;
00257   int_type *MyBlockGIDs=new int_type[NumBlocks_];
00258   for(i=0;i<NumBlocks_;i++)
00259     MyBlockGIDs[i]=MyMinGID+i;
00260 
00261   BDMap_=new Epetra_BlockMap((int_type) -1,NumBlocks_,MyBlockGIDs,bsize,0,Matrix_->Comm());
00262   
00263   // Allocate the Epetra_DistBlockMatrix
00264   BDMat_=new EpetraExt_BlockDiagMatrix(*BDMap_,true);  
00265   
00266   // First check to see if we can switch back to PurelyLocalMode_ if we're not
00267   // in it
00268   if(!PurelyLocalMode_){
00269     // Find the non-local row IDs
00270     ExtSize=MAX(0,Blockstart_[NumBlocks_]-RowMap.NumMyElements());// estimate
00271     ExtRows.reserve(ExtSize);
00272     
00273     for(i=0;i<Blockstart_[NumBlocks_];i++){
00274       if(RowMap.LID(Blockids_const_ptr<int_type>()[i])==-1)
00275         ExtRows.push_back(Blockids_const_ptr<int_type>()[i]);
00276     }
00277 
00278     // Switch to PurelyLocalMode_ if we never need GIDs
00279     int size_sum;
00280     ExtSize=ExtRows.size();   
00281     RowMap.Comm().SumAll(&ExtSize,&size_sum,1);
00282     if(size_sum==0){
00283       if(verbose && !Matrix_->Comm().MyPID()) printf("EpetraExt_PointToBlockDiagPermute: Switching to purely local mode\n");
00284       PurelyLocalMode_=true;
00285       for(i=0;i<Blockstart_[NumBlocks_];i++){
00286         Blockids_ref<int>()[i]=RowMap.LID(Blockids_const_ptr<int_type>()[i]);
00287       }     
00288     }
00289   }
00290   
00291   if(PurelyLocalMode_){
00292     /*******************************************************/
00293     // Allocations
00294     l2b=new int[Nrows];
00295     block_offset=new int[Nrows];
00296 
00297     // Build the local-id-to-block-id list, block offset list    
00298     for(i=0;i<NumBlocks_;i++) {
00299       for(j=Blockstart_[i];j<Blockstart_[i+1];j++){
00300         block_offset[Blockids_const_ptr<int>()[j]]=j-Blockstart_[i];
00301         l2b[Blockids_const_ptr<int>()[j]]=i;
00302       }
00303     }
00304     
00305     // Copy the data to the EpetraExt_BlockDiagMatrix        
00306     for(i=0;i<Nrows;i++) {
00307       int block_num = l2b[i];
00308       if(block_num>=0 && block_num<NumBlocks_) {
00309         row_in_block=block_offset[i];
00310         Matrix_->ExtractMyRowView(i,length,values,colind);
00311         int Nnz=0;
00312         for (j = 0; j < length; j++) {
00313           col = colind[j];
00314           if(col < Nrows && l2b[col]==block_num){
00315             Nnz++;
00316             col_in_block = block_offset[col];
00317             index = col_in_block * bsize[block_num] + row_in_block;
00318             (*BDMat_)[block_num][index]=values[j];
00319           }
00320         }
00321         /* Handle the case of a zero row. */
00322         /* By just putting a 1 on the diagonal. */
00323         if (Nnz == 0) {
00324           index = row_in_block * bsize[block_num] + row_in_block;
00325           (*BDMat_)[block_num][index]=1.0;
00326         }
00327       }
00328     }
00329 
00330     // Build the compatible map for import/export
00331     l2blockid=new int_type[Nrows];
00332     for(i=0;i<Nrows;i++) 
00333       l2blockid[Blockstart_[l2b[i]]+block_offset[i]]= (int_type) Matrix_->RowMap().GID64(i);
00334 
00335     // Build the Compatible Map, Import/Export Objects
00336     CompatibleMap_=new Epetra_Map(-1,Nrows,l2blockid,0,Matrix_->Comm());
00337     
00338   }
00339   else{
00340     /*******************************************************/      
00341     // Do the import to grab matrix entries
00342     // Allocate temporaries for import
00343     int_type* ExtRowsPtr = ExtRows.size()>0 ? &ExtRows[0] : NULL;
00344     Epetra_Map TmpMap((int_type) -1,ExtSize, ExtRowsPtr,0,Matrix_->Comm());; 
00345     Epetra_CrsMatrix TmpMatrix(Copy,TmpMap,0);
00346     Epetra_Import TmpImporter(TmpMap,RowMap);
00347 
00348     TmpMatrix.Import(*Matrix_,TmpImporter,Insert);
00349     TmpMatrix.FillComplete();
00350       
00351     ExtCols=TmpMatrix.NumMyCols();
00352     MainCols=Matrix_->NumMyCols();
00353         
00354 
00355     // Build the column reidex - main matrix
00356     LocalColIDS=new int[MainCols+ExtCols];      
00357     for(i=0;i<MainCols;i++){
00358       int_type GID= (int_type) Matrix_->GCID64(i);
00359       int MainLID=RowMap.LID(GID);
00360       if(MainLID!=-1) LocalColIDS[i]=MainLID;
00361       else{
00362         int ExtLID=TmpMatrix.LRID(GID);
00363         if(ExtLID!=-1) LocalColIDS[i]=Nrows+ExtLID;
00364         else LocalColIDS[i]=-1;
00365       }
00366     }
00367 
00368     // Build the column reidex - ext matrix
00369     for(i=0;i<ExtCols;i++){
00370       int_type GID= (int_type) TmpMatrix.GCID64(i);
00371       int MainLID=RowMap.LID(GID);
00372       if(MainLID!=-1) LocalColIDS[MainCols+i]=MainLID;
00373       else{
00374         int ExtLID=TmpMatrix.LRID(GID);          
00375         if(ExtLID!=-1) LocalColIDS[MainCols+i]=Nrows+ExtLID;
00376         else LocalColIDS[MainCols+i]=-1;
00377       }
00378     }
00379 
00380     // Allocations
00381     l2b=new int[Nrows+ExtSize];
00382     block_offset=new int[Nrows+ExtSize];  
00383   
00384     // Build l2b/block_offset with the expanded local index
00385     //NTS: Uses the same ordering of operation as the above routine, which is why
00386     //it works.
00387     for(i=0;i<Nrows+ExtSize;i++) block_offset[i]=l2b[i]=-1;    
00388     int ext_idx=0;
00389     for(i=0;i<NumBlocks_;i++) {
00390       for(j=Blockstart_[i];j<Blockstart_[i+1];j++){
00391         int LID=RowMap.LID(Blockids_const_ptr<int_type>()[j]);
00392         if(LID==-1) {LID=Nrows+ext_idx;ext_idx++;}
00393         block_offset[LID]=j-Blockstart_[i];
00394         l2b[LID]=i;          
00395       }
00396     }
00397     
00398     // Copy the data to the EpetraExt_BlockDiagMatrix from Matrix_
00399     for(i=0;i<Nrows;i++) {
00400       int block_num = l2b[i];
00401       if(block_num>=0 && block_num<NumBlocks_) {
00402         row_in_block=block_offset[i];
00403         Matrix_->ExtractMyRowView(i,length,values,colind);
00404         int Nnz=0;
00405         for (j = 0; j < length; j++) {
00406           col = LocalColIDS[colind[j]];
00407           if(col!=-1 && l2b[col]==block_num){
00408             Nnz++;
00409             col_in_block = block_offset[col];
00410             index = col_in_block * bsize[block_num] + row_in_block;
00411             (*BDMat_)[block_num][index]=values[j];
00412           }
00413         }
00414         /* Handle the case of a zero row. */
00415         /* By just putting a 1 on the diagonal. */
00416         if (Nnz == 0) {
00417           index = row_in_block * bsize[block_num] + row_in_block;
00418           (*BDMat_)[block_num][index]=values[j];
00419         }
00420       }
00421     }
00422     
00423     // Copy the data to the EpetraExt_BlockDiagMatrix from TmpMatrix
00424     for(i=0;i<ExtSize;i++) {
00425       int block_num = l2b[Nrows+i];
00426       if(block_num>=0 && block_num<NumBlocks_) {
00427         row_in_block=block_offset[Nrows+i];
00428         TmpMatrix.ExtractMyRowView(i,length,values,colind);
00429         int Nnz=0;
00430         for (j = 0; j < length; j++) {
00431           col = LocalColIDS[MainCols+colind[j]];
00432           if(col!=-1 && l2b[col]==block_num){
00433             Nnz++;
00434             col_in_block = block_offset[col];
00435             index = col_in_block * bsize[block_num] + row_in_block;
00436             (*BDMat_)[block_num][index]=values[j];
00437           }
00438         }
00439         /* Handle the case of a zero row. */
00440         /* By just putting a 1 on the diagonal. */
00441         if (Nnz == 0) {
00442           index = row_in_block * bsize[block_num] + row_in_block;
00443           (*BDMat_)[block_num][index]=1.0;
00444         }
00445       }      
00446     }
00447     
00448     // Build the compatible map for import/export
00449     l2blockid=new int_type[Blockstart_[NumBlocks_]];
00450     for(i=0;i<Nrows;i++){
00451       int bkid=l2b[i];
00452       if(bkid>-1) l2blockid[Blockstart_[bkid]+block_offset[i]]= (int_type) RowMap.GID64(i);
00453     }
00454     // NTS: This is easier - if we imported it, we know we need it
00455     for(i=0;i<ExtSize;i++)
00456       l2blockid[Blockstart_[l2b[Nrows+i]]+block_offset[Nrows+i]]= (int_type) TmpMatrix.GRID64(i);
00457     
00458     // Build the Compatible Map, Import/Export Objects
00459     CompatibleMap_=new Epetra_Map(-1,Blockstart_[NumBlocks_],l2blockid,0,Matrix_->Comm());
00460         
00461   }//end else
00462 
00463   // Set BDMat_'s parameter list and compute!
00464   Teuchos::ParameterList dummy,inList;
00465   inList=List_.get("blockdiagmatrix: list",dummy);
00466   BDMat_->SetParameters(inList);  
00467   BDMat_->Compute();
00468   
00469   // Build importer/exporter
00470   if(!CompatibleMap_->SameAs(Matrix_->DomainMap()))
00471     Importer_ = new Epetra_Import(*CompatibleMap_,Matrix_->DomainMap());
00472   if(!CompatibleMap_->SameAs(Matrix_->RangeMap()))
00473     Exporter_ = new Epetra_Export(*CompatibleMap_,Matrix_->RangeMap());
00474 
00475   // Cleanup
00476   delete [] LocalColIDS;
00477   delete [] block_offset;
00478   delete [] l2b;
00479   delete [] l2blockid;
00480   delete [] bsize;
00481   delete [] MyBlockGIDs;
00482 
00483   // Contiguous Cleanup
00484   CleanupContiguousMode();
00485 
00486   return 0;
00487 }
00488 
00489 int EpetraExt_PointToBlockDiagPermute::ExtractBlockDiagonal(){
00490 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00491   if(Matrix_->RowMap().GlobalIndicesInt()) {
00492   return TExtractBlockDiagonal<int>();
00493   }
00494   else
00495 #endif
00496 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00497   if(Matrix_->RowMap().GlobalIndicesLongLong()) {
00498   return TExtractBlockDiagonal<long long>();
00499   }
00500   else
00501 #endif
00502     throw "EpetraExt_PointToBlockDiagPermute::ExtractBlockDiagonal: ERROR, GlobalIndices type unknown.";
00503 }
00504 
00505 //=======================================================================================================
00506 template<typename int_type>
00507 int EpetraExt_PointToBlockDiagPermute::TSetupContiguousMode(){
00508   if(!ContiguousBlockMode_) return 0;
00509   // NTS: In case of processor-crossing blocks, the lowest PID always gets the block;
00510   const Epetra_Map &RowMap=Matrix_->RowMap();
00511   
00512   int_type MinMyGID= (int_type) RowMap.MinMyGID64(); 
00513   int_type MaxMyGID= (int_type) RowMap.MaxMyGID64(); 
00514   int_type Base= (int_type) Matrix_->IndexBase64();
00515 
00516   // Find the GID that begins my first block
00517   int_type MyFirstBlockGID=ContiguousBlockSize_*(int_type)ceil(((double)(MinMyGID - Base))/ContiguousBlockSize_)+Base;
00518   NumBlocks_=(int)ceil((double)((MaxMyGID-MyFirstBlockGID+1.0)) / ContiguousBlockSize_);
00519 
00520   // Allocate memory
00521   Blockstart_=new int[NumBlocks_+1];
00522   Blockids_ref<int_type>()=new int_type[NumBlocks_*ContiguousBlockSize_];
00523   if(sizeof(int) != sizeof(int_type))
00524     Blockids_ref<int>()=new int[NumBlocks_*ContiguousBlockSize_];
00525   Blockstart_[NumBlocks_]=NumBlocks_*ContiguousBlockSize_;
00526 
00527   // Fill the arrays
00528   for(int i=0,ct=0;i<NumBlocks_;i++){
00529     Blockstart_[i]=ct;
00530     for(int j=0;j<ContiguousBlockSize_;j++,ct++){
00531       Blockids_ref<int_type>()[ct]=MyFirstBlockGID+ct;
00532     }
00533   }
00534   
00535   return 0;
00536 }
00537 
00538 int EpetraExt_PointToBlockDiagPermute::SetupContiguousMode(){
00539 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00540   if(Matrix_->RowMap().GlobalIndicesInt()) {
00541   return TSetupContiguousMode<int>();
00542   }
00543   else
00544 #endif
00545 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00546   if(Matrix_->RowMap().GlobalIndicesLongLong()) {
00547   return TSetupContiguousMode<long long>();
00548   }
00549   else
00550 #endif
00551     throw "EpetraExt_PointToBlockDiagPermute::SetupContiguousMode: ERROR, GlobalIndices type unknown.";
00552 }
00553 
00554 //=======================================================================================================
00555 int EpetraExt_PointToBlockDiagPermute::CleanupContiguousMode(){
00556   if(!ContiguousBlockMode_) return 0;
00557   NumBlocks_=0;
00558   if(Blockstart_) {delete [] Blockstart_; Blockstart_=0;}
00559   if(Blockids_int_)   {delete [] Blockids_int_; Blockids_int_=0;}
00560 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00561   if(Blockids_LL_)   {delete [] Blockids_LL_; Blockids_LL_=0;}
00562 #endif
00563   return 0;
00564 }
00565 
00566 
00567 //=======================================================================================================
00568 // Creates an Epetra_CrsMatrix from the BlockDiagMatrix.  This is generally only useful if you want to do a matrix-matrix multiply.
00569 template<typename int_type>
00570 Epetra_FECrsMatrix * EpetraExt_PointToBlockDiagPermute::TCreateFECrsMatrix(){
00571   Epetra_FECrsMatrix * NewMat=new Epetra_FECrsMatrix(Copy,Matrix_->RowMap(),0);
00572   
00573   const Epetra_BlockMap &BlockMap=BDMat_->BlockMap();
00574   const Epetra_BlockMap &DataMap=BDMat_->DataMap();
00575   const int *vlist=DataMap.FirstPointInElementList();
00576   const int *xlist=BlockMap.FirstPointInElementList();
00577   const int *blocksize=BlockMap.ElementSizeList();
00578   const double *values=BDMat_->Values();
00579   int NumBlocks=BDMat_->NumMyBlocks();
00580 
00581   // Maximum size vector for stashing GIDs
00582   std::vector<int_type> GIDs;
00583   GIDs.resize(BlockMap.MaxMyElementSize());
00584 
00585 
00586   for(int i=0;i<NumBlocks;i++){
00587     int Nb=blocksize[i];
00588     int vidx0=vlist[i];
00589     int xidx0=xlist[i];
00590     // Get global indices
00591     for(int j=0;j<Nb;j++)
00592       GIDs[j]= (int_type) CompatibleMap_->GID64(xidx0+j);
00593     
00594     // Remember: We're using column-major storage for LAPACK's benefit    
00595     int ierr=NewMat->InsertGlobalValues(Nb,&GIDs[0],&values[vidx0],Epetra_FECrsMatrix::COLUMN_MAJOR);
00596     if(ierr < 0) throw "CreateFECrsMatrix: ERROR in InsertGlobalValues";
00597   }   
00598   NewMat->GlobalAssemble();
00599   return NewMat;
00600 }
00601 
00602 Epetra_FECrsMatrix * EpetraExt_PointToBlockDiagPermute::CreateFECrsMatrix(){
00603 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00604   if(Matrix_->RowMap().GlobalIndicesInt()) {
00605   return TCreateFECrsMatrix<int>();
00606   }
00607   else
00608 #endif
00609 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00610   if(Matrix_->RowMap().GlobalIndicesLongLong()) {
00611   return TCreateFECrsMatrix<long long>();
00612   }
00613   else
00614 #endif
00615     throw "EpetraExt_PointToBlockDiagPermute::CreateFECrsMatrix: ERROR, GlobalIndices type unknown.";
00616 }
00617 
00618 //=======================================================================================================
00619 void EpetraExt_PointToBlockDiagPermute::UpdateImportVector(int NumVectors) const {    
00620   if(Importer_ != 0) {
00621     if(ImportVector_ != 0) {
00622       if(ImportVector_->NumVectors() != NumVectors) { 
00623   delete ImportVector_; 
00624   ImportVector_= 0;
00625       }
00626     }
00627     if(ImportVector_ == 0) 
00628       ImportVector_ = new Epetra_MultiVector(*CompatibleMap_,NumVectors); // Create import vector if needed
00629   }
00630   return;
00631 }
00632 //=======================================================================================================
00633 void EpetraExt_PointToBlockDiagPermute::UpdateExportVector(int NumVectors) const {    
00634   if(Exporter_ != 0) {
00635     if(ExportVector_ != 0) {
00636       if(ExportVector_->NumVectors() != NumVectors) { 
00637   delete ExportVector_; 
00638   ExportVector_= 0;
00639       }
00640     }
00641     if(ExportVector_ == 0) 
00642       ExportVector_ = new Epetra_MultiVector(*CompatibleMap_,NumVectors); // Create Export vector if needed
00643   }
00644   return;
00645 }
00646 
00647 
00648 
00649 
00650 //=========================================================================  
00651 int EpetraExt_PointToBlockDiagPermute::Import(const Epetra_SrcDistObject& A, const Epetra_Import& Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor){
00652  return -1;
00653 }
00654 
00655 //=========================================================================  
00656 int EpetraExt_PointToBlockDiagPermute::Import(const Epetra_SrcDistObject& A, const Epetra_Export& Exporter, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor){
00657  return -1;
00658 }
00659 
00660 //=========================================================================  
00661 int EpetraExt_PointToBlockDiagPermute::Export(const Epetra_SrcDistObject& A, const Epetra_Import & Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor){
00662  return -1;
00663 }
00664 
00665 //=========================================================================  
00666 int EpetraExt_PointToBlockDiagPermute::Export(const Epetra_SrcDistObject& A, const Epetra_Export& Exporter, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor){
00667  return -1;
00668 }
00669 
00670 //=========================================================================  
00671 // Allows the source and target (\e this) objects to be compared for compatibility, return nonzero if not.
00672 int EpetraExt_PointToBlockDiagPermute::CheckSizes(const Epetra_SrcDistObject& Source){
00673   return -1;
00674 }
00675 
00676 //=========================================================================  
00677 // Perform ID copies and permutations that are on processor.
00678 int EpetraExt_PointToBlockDiagPermute::CopyAndPermute(const Epetra_SrcDistObject& Source,
00679                    int NumSameIDs, 
00680                    int NumPermuteIDs,
00681                    int * PermuteToLIDs,
00682                    int * PermuteFromLIDs,
00683                    const Epetra_OffsetIndex * Indexor){
00684   return -1;
00685 }
00686 
00687 //=========================================================================  
00688 // Perform any packing or preparation required for call to DoTransfer().
00689 int EpetraExt_PointToBlockDiagPermute::PackAndPrepare(const Epetra_SrcDistObject& Source,
00690                    int NumExportIDs,
00691                    int* ExportLIDs,
00692                    int& LenExports,
00693                    char*& Exports,
00694                    int& SizeOfPacket,
00695                    int* Sizes,
00696                    bool & VarSizes,
00697                    Epetra_Distributor& Distor){
00698   return -1;
00699 }
00700 
00701 //=========================================================================  
00702 // Perform any unpacking and combining after call to DoTransfer().
00703 int EpetraExt_PointToBlockDiagPermute::UnpackAndCombine(const Epetra_SrcDistObject& Source, 
00704                      int NumImportIDs,
00705                      int* ImportLIDs, 
00706                      int LenImports,
00707                      char* Imports,
00708                      int& SizeOfPacket, 
00709                      Epetra_Distributor& Distor,
00710                      Epetra_CombineMode CombineMode,
00711                      const Epetra_OffsetIndex * Indexor){
00712   return -1;
00713 }
00714 
00715 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines