EpetraExt Package Browser (Single Doxygen Collection) Development
EpetraExt_BlockUtility.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2011) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041 
00042 #include "EpetraExt_BlockUtility.h"
00043 #include "Epetra_Map.h"
00044 #include "Epetra_Comm.h"
00045 
00046 namespace EpetraExt {
00047 
00048 using std::vector;
00049 
00050 //==============================================================================
00051 Epetra_Map * BlockUtility::GenerateBlockMap(
00052   const Epetra_BlockMap & BaseMap,
00053         const int * RowIndices,
00054   int NumBlockRows,
00055         const Epetra_Comm & GlobalComm,
00056   int Offset) 
00057 {
00058   int BaseIndex = BaseMap.IndexBase();
00059   if (Offset == 0)
00060     Offset = BlockUtility::CalculateOffset(BaseMap);
00061 
00062   //Get Base Global IDs
00063   int Size = BaseMap.NumMyElements();
00064   int TotalSize = NumBlockRows * Size;
00065   vector<int> GIDs(Size);
00066   BaseMap.MyGlobalElements( &GIDs[0] );
00067 
00068   vector<int> GlobalGIDs( TotalSize );
00069   for( int i = 0; i < NumBlockRows; ++i )
00070   {
00071     for( int j = 0; j < Size; ++j )
00072       GlobalGIDs[i*Size+j] = GIDs[j] + RowIndices[i] * Offset;
00073   }
00074 
00075   int GlobalSize;
00076   GlobalComm.SumAll( &TotalSize, &GlobalSize, 1 );
00077 
00078   Epetra_Map *GlobalMap = 
00079     new Epetra_Map( GlobalSize, TotalSize, &GlobalGIDs[0], BaseIndex, 
00080         GlobalComm );
00081 
00082   return GlobalMap;
00083 }
00084 
00085 //==============================================================================
00086 Epetra_Map * BlockUtility::GenerateBlockMap(
00087   const Epetra_BlockMap & BaseMap,
00088         const std::vector<int>& RowIndices,
00089   const Epetra_Comm & GlobalComm,
00090   int Offset ) 
00091 {
00092   return GenerateBlockMap(BaseMap, &RowIndices[0], RowIndices.size(), 
00093         GlobalComm, Offset);
00094 }
00095 
00096 //==============================================================================
00097 Epetra_Map * BlockUtility::GenerateBlockMap(
00098   const Epetra_BlockMap & BaseMap,
00099         const Epetra_BlockMap& BlockMap,
00100   const Epetra_Comm & GlobalComm,
00101   int Offset) 
00102 {
00103   return GenerateBlockMap(BaseMap, 
00104         BlockMap.MyGlobalElements(), 
00105         BlockMap.NumMyElements(), 
00106         GlobalComm,
00107         Offset);
00108 }
00109 
00110 //==============================================================================
00111 Epetra_CrsGraph * BlockUtility::GenerateBlockGraph(
00112         const Epetra_CrsGraph & BaseGraph,
00113         const vector< vector<int> > & RowStencil,
00114         const vector<int> & RowIndices,
00115         const Epetra_Comm & GlobalComm ) 
00116 {
00117 
00118   const Epetra_BlockMap & BaseMap = BaseGraph.RowMap();
00119   int BaseIndex = BaseMap.IndexBase();
00120   int Offset = BlockUtility::CalculateOffset(BaseMap);
00121 
00122   //Get Base Global IDs
00123   int NumBlockRows = RowIndices.size();
00124   int Size = BaseMap.NumMyElements();
00125   int TotalSize = NumBlockRows * Size;
00126   vector<int> GIDs(Size);
00127   BaseMap.MyGlobalElements( &GIDs[0] );
00128 
00129   vector<int> GlobalGIDs( TotalSize );
00130   for( int i = 0; i < NumBlockRows; ++i )
00131   {
00132     for( int j = 0; j < Size; ++j )
00133       GlobalGIDs[i*Size+j] = GIDs[j] + RowIndices[i] * Offset;
00134   }
00135 
00136   int GlobalSize;
00137   GlobalComm.SumAll( &TotalSize, &GlobalSize, 1 );
00138 
00139   Epetra_Map GlobalMap( GlobalSize, TotalSize, &GlobalGIDs[0], BaseIndex, GlobalComm );
00140 
00141   int MaxIndices = BaseGraph.MaxNumIndices();
00142   vector<int> Indices(MaxIndices);
00143   int NumIndices;
00144 
00145   Epetra_CrsGraph * GlobalGraph = new Epetra_CrsGraph( Copy, 
00146                                dynamic_cast<Epetra_BlockMap&>(GlobalMap),
00147                                0 );
00148 
00149   for( int i = 0; i < NumBlockRows; ++i )
00150   {
00151     int StencilSize = RowStencil[i].size();
00152     for( int j = 0; j < Size; ++j )
00153     {
00154       int BaseRow = BaseMap.GID(j);
00155       int GlobalRow = GlobalMap.GID(j+i*Size);
00156 
00157       BaseGraph.ExtractGlobalRowCopy( BaseRow, MaxIndices, NumIndices, &Indices[0] );
00158       for( int k = 0; k < StencilSize; ++k )
00159       {
00160         int ColOffset = (RowIndices[i]+RowStencil[i][k]) * Offset;
00161         if( k > 0 ) ColOffset -= (RowIndices[i]+RowStencil[i][k-1]) * Offset;
00162 
00163         for( int l = 0; l < NumIndices; ++l )
00164           Indices[l] += ColOffset;
00165 
00166         GlobalGraph->InsertGlobalIndices( GlobalRow, NumIndices, &Indices[0] );
00167       }
00168     }
00169   }
00170 
00171   GlobalGraph->FillComplete();
00172 
00173   return GlobalGraph;
00174 }
00175 
00176 //==============================================================================
00177 Epetra_CrsGraph * BlockUtility::GenerateBlockGraph(
00178         const Epetra_RowMatrix & BaseMatrix,
00179         const vector< vector<int> > & RowStencil,
00180         const vector<int> & RowIndices,
00181         const Epetra_Comm & GlobalComm ) 
00182 {
00183 
00184   const Epetra_BlockMap & BaseMap = BaseMatrix.RowMatrixRowMap();
00185   const Epetra_BlockMap & BaseColMap = BaseMatrix.RowMatrixColMap();
00186   int BaseIndex = BaseMap.IndexBase();
00187   int Offset = BlockUtility::CalculateOffset(BaseMap);
00188 
00189   //Get Base Global IDs
00190   int NumBlockRows = RowIndices.size();
00191   int Size = BaseMap.NumMyElements();
00192   int TotalSize = NumBlockRows * Size;
00193   vector<int> GIDs(Size);
00194   BaseMap.MyGlobalElements( &GIDs[0] );
00195 
00196   vector<int> GlobalGIDs( TotalSize );
00197   for( int i = 0; i < NumBlockRows; ++i )
00198   {
00199     for( int j = 0; j < Size; ++j )
00200       GlobalGIDs[i*Size+j] = GIDs[j] + RowIndices[i] * Offset;
00201   }
00202 
00203   int GlobalSize;
00204   GlobalComm.SumAll( &TotalSize, &GlobalSize, 1 );
00205 
00206   Epetra_Map GlobalMap( GlobalSize, TotalSize, &GlobalGIDs[0], BaseIndex, GlobalComm );
00207 
00208   int MaxIndices = BaseMatrix.MaxNumEntries();
00209   vector<int> Indices(MaxIndices);
00210   vector<double> Values(MaxIndices); 
00211   int NumIndices;
00212 
00213   Epetra_CrsGraph * GlobalGraph = new Epetra_CrsGraph( Copy, 
00214                                dynamic_cast<Epetra_BlockMap&>(GlobalMap),
00215                                0 );
00216 
00217   for( int i = 0; i < NumBlockRows; ++i )
00218   {
00219     int StencilSize = RowStencil[i].size();
00220     for( int j = 0; j < Size; ++j )
00221     {
00222       int GlobalRow = GlobalMap.GID(j+i*Size);
00223 
00224       BaseMatrix.ExtractMyRowCopy( j, MaxIndices, NumIndices, &Values[0], &Indices[0] );
00225       for( int l = 0; l < NumIndices; ++l ) Indices[l] = BaseColMap.GID(Indices[l]);
00226 
00227       for( int k = 0; k < StencilSize; ++k )
00228       {
00229         int ColOffset = (RowIndices[i]+RowStencil[i][k]) * Offset;
00230         if( k > 0 ) ColOffset -= (RowIndices[i]+RowStencil[i][k-1]) * Offset;
00231 
00232         for( int l = 0; l < NumIndices; ++l )
00233           Indices[l] += ColOffset;
00234 
00235         GlobalGraph->InsertGlobalIndices( GlobalRow, NumIndices, &Indices[0] );
00236       }
00237     }
00238   }
00239 
00240   GlobalGraph->FillComplete();
00241 
00242   return GlobalGraph;
00243 }
00244 
00245 //==============================================================================
00246 Epetra_CrsGraph * BlockUtility::GenerateBlockGraph(
00247         const Epetra_CrsGraph & BaseGraph,
00248         const Epetra_CrsGraph & LocalBlockGraph,
00249         const Epetra_Comm & GlobalComm ) 
00250 {
00251   const Epetra_BlockMap & BaseRowMap = BaseGraph.RowMap();
00252   const Epetra_BlockMap & BaseColMap = BaseGraph.ColMap();
00253   int ROffset = BlockUtility::CalculateOffset(BaseRowMap);
00254   int COffset = BlockUtility::CalculateOffset(BaseColMap);
00255 
00256   //Get Base Global IDs
00257   const Epetra_BlockMap & BlockRowMap = LocalBlockGraph.RowMap();
00258   const Epetra_BlockMap & BlockColMap = LocalBlockGraph.ColMap();
00259   
00260   int NumBlockRows = BlockRowMap.NumMyElements();
00261   vector<int> RowIndices(NumBlockRows);
00262   BlockRowMap.MyGlobalElements(&RowIndices[0]);
00263 
00264   int Size = BaseRowMap.NumMyElements();
00265   
00266   Epetra_Map *GlobalRowMap = 
00267     GenerateBlockMap(BaseRowMap, BlockRowMap, GlobalComm);
00268   
00269 
00270   int MaxIndices = BaseGraph.MaxNumIndices();
00271   vector<int> Indices(MaxIndices);
00272 
00273   Epetra_CrsGraph * GlobalGraph = new Epetra_CrsGraph( Copy, 
00274                                dynamic_cast<Epetra_BlockMap&>(*GlobalRowMap),
00275                                0 );
00276 
00277   int NumBlockIndices, NumBaseIndices;
00278   int *BlockIndices, *BaseIndices;
00279   for( int i = 0; i < NumBlockRows; ++i )
00280   {
00281     LocalBlockGraph.ExtractMyRowView(i, NumBlockIndices, BlockIndices);
00282     
00283     for( int j = 0; j < Size; ++j )
00284     {
00285       int GlobalRow = GlobalRowMap->GID(j+i*Size);
00286 
00287       BaseGraph.ExtractMyRowView( j, NumBaseIndices, BaseIndices );
00288       for( int k = 0; k < NumBlockIndices; ++k )
00289       {
00290         int ColOffset = BlockColMap.GID(BlockIndices[k]) * COffset;
00291 
00292         for( int l = 0; l < NumBaseIndices; ++l )
00293           Indices[l] = BaseGraph.GCID(BaseIndices[l]) + ColOffset;
00294 
00295         GlobalGraph->InsertGlobalIndices( GlobalRow, NumBaseIndices, &Indices[0] );
00296       }
00297     }
00298   }
00299 
00300   const Epetra_BlockMap & BaseDomainMap = BaseGraph.DomainMap();
00301   const Epetra_BlockMap & BaseRangeMap = BaseGraph.RangeMap();
00302   const Epetra_BlockMap & BlockDomainMap = LocalBlockGraph.DomainMap();
00303   const Epetra_BlockMap & BlockRangeMap = LocalBlockGraph.RangeMap();
00304 
00305   Epetra_Map *GlobalDomainMap = 
00306     GenerateBlockMap(BaseDomainMap, BlockDomainMap, GlobalComm);
00307   Epetra_Map *GlobalRangeMap = 
00308     GenerateBlockMap(BaseRangeMap, BlockRangeMap, GlobalComm);
00309 
00310   GlobalGraph->FillComplete(*GlobalDomainMap, *GlobalRangeMap);
00311 
00312   delete GlobalDomainMap;
00313   delete GlobalRangeMap;
00314   delete GlobalRowMap;
00315 
00316   return GlobalGraph;
00317 }
00318 
00319 //==============================================================================
00320 void BlockUtility::GenerateRowStencil(const Epetra_CrsGraph& LocalBlockGraph, 
00321               std::vector<int> RowIndices, 
00322               std::vector< std::vector<int> >& RowStencil)
00323 {
00324   // Get row indices
00325   int NumMyRows = LocalBlockGraph.NumMyRows();
00326   RowIndices.resize(NumMyRows);
00327   const Epetra_BlockMap& RowMap = LocalBlockGraph.RowMap();
00328   RowMap.MyGlobalElements(&RowIndices[0]);
00329 
00330   // Get stencil
00331   RowStencil.resize(NumMyRows);
00332   if (LocalBlockGraph.IndicesAreGlobal()) {
00333     for (int i=0; i<NumMyRows; i++) {
00334       int Row = RowIndices[i];
00335       int NumCols = LocalBlockGraph.NumGlobalIndices(Row);
00336       RowStencil[i].resize(NumCols);
00337       LocalBlockGraph.ExtractGlobalRowCopy(Row, NumCols, NumCols, 
00338              &RowStencil[i][0]);
00339       for (int k=0; k<NumCols; k++)
00340   RowStencil[i][k] -= Row;
00341     }
00342   }
00343   else {
00344     for (int i=0; i<NumMyRows; i++) {
00345       int NumCols = LocalBlockGraph.NumMyIndices(i);
00346       RowStencil[i].resize(NumCols);
00347       LocalBlockGraph.ExtractMyRowCopy(i, NumCols, NumCols, 
00348                &RowStencil[i][0]);
00349       for (int k=0; k<NumCols; k++)
00350   RowStencil[i][k] = LocalBlockGraph.GCID(RowStencil[i][k]) - 
00351     RowIndices[i];
00352     }
00353   }
00354 }
00355 
00356 //==============================================================================
00357 int BlockUtility::CalculateOffset(const Epetra_BlockMap & BaseMap)
00358 {
00359   int MaxGID = BaseMap.MaxAllGID();
00360 
00361 //   int Offset = 1;
00362 //   while( Offset <= MaxGID ) Offset *= 10;
00363 
00364 //   return Offset;
00365 
00366   return MaxGID+1;
00367 }
00368 
00369 
00370 } //namespace EpetraExt
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines