43#include "Epetra_Map.h" 
   44#include "Epetra_Comm.h" 
   50template<
typename int_type>
 
   53        const int_type * RowIndices,
 
   58  int_type BaseIndex = (int_type) BaseMap.IndexBase64();
 
   60    Offset = BlockUtility::TCalculateOffset<int_type>(BaseMap);
 
   64  int TotalSize = NumBlockRows * Size;
 
   65  vector<int_type> GIDs(Size);
 
   68  vector<int_type> GlobalGIDs( TotalSize );
 
   69  for( 
int i = 0; i < NumBlockRows; ++i )
 
   71    for( 
int j = 0; j < Size; ++j )
 
   72      GlobalGIDs[i*Size+j] = GIDs[j] + RowIndices[i] * Offset;
 
   76  int_type TotalSize_int_type = TotalSize;
 
   77  GlobalComm.
SumAll( &TotalSize_int_type, &GlobalSize, 1 );
 
   80    new Epetra_Map( GlobalSize, TotalSize, &GlobalGIDs[0], BaseIndex,
 
   87#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 
   90        const int * RowIndices,
 
   96    return TGenerateBlockMap<int>(BaseMap, RowIndices, NumBlockRows, GlobalComm, Offset);
 
   98    throw "EpetraExt::BlockUtility::GenerateBlockMap: Global Indices not int.";
 
 
  102#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 
  105        const long long * RowIndices,
 
  111    return TGenerateBlockMap<long long>(BaseMap, RowIndices, NumBlockRows, GlobalComm, Offset);
 
  113    throw "EpetraExt::BlockUtility::GenerateBlockMap: Global Indices not long long.";
 
 
  117#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 
  120        const std::vector<int>& RowIndices,
 
 
  129#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 
  132        const std::vector<long long>& RowIndices,
 
 
  148#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 
  157#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 
  160                          BlockMap.MyGlobalElements64(),
 
  166    throw "EpetraExt::BlockUtility::GenerateBlockMap: Error Global Indices unknown.";
 
 
  170template<
typename int_type>
 
  173        const vector< vector<int_type> > & RowStencil,
 
  174        const vector<int_type> & RowIndices,
 
  179  int_type BaseIndex = (int_type) BaseMap.IndexBase64();
 
  180  int_type Offset = BlockUtility::TCalculateOffset<int_type>(BaseMap);
 
  183  int NumBlockRows = RowIndices.size();
 
  185  int TotalSize = NumBlockRows * Size;
 
  186  vector<int_type> GIDs(Size);
 
  189  vector<int_type> GlobalGIDs( TotalSize );
 
  190  for( 
int i = 0; i < NumBlockRows; ++i )
 
  192    for( 
int j = 0; j < Size; ++j )
 
  193      GlobalGIDs[i*Size+j] = GIDs[j] + RowIndices[i] * Offset;
 
  197  int_type TotalSize_int_type = TotalSize;
 
  198  GlobalComm.
SumAll( &TotalSize_int_type, &GlobalSize, 1 );
 
  200  Epetra_Map GlobalMap( GlobalSize, TotalSize, &GlobalGIDs[0], BaseIndex, GlobalComm );
 
  203  vector<int_type> Indices(MaxIndices);
 
  210  for( 
int i = 0; i < NumBlockRows; ++i )
 
  212    int StencilSize = RowStencil[i].size();
 
  213    for( 
int j = 0; j < Size; ++j )
 
  215      int_type BaseRow = (int_type) BaseMap.GID64(j);
 
  216      int_type GlobalRow = (int_type) GlobalMap.GID64(j+i*Size);
 
  219      for( 
int k = 0; k < StencilSize; ++k )
 
  221        int_type ColOffset = (RowIndices[i]+RowStencil[i][k]) * Offset;
 
  222        if( k > 0 ) ColOffset -= (RowIndices[i]+RowStencil[i][k-1]) * Offset;
 
  224        for( 
int l = 0; l < NumIndices; ++l )
 
  225          Indices[l] += ColOffset;
 
 
  237#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 
  240        const vector< vector<int> > & RowStencil,
 
  241        const vector<int> & RowIndices,
 
  245    return TGenerateBlockGraph<int>(BaseGraph, RowStencil, RowIndices, GlobalComm);
 
  247    throw "EpetraExt::BlockUtility::GenerateBlockGraph: Error Global Indices not int.";
 
  251#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 
  254        const vector< vector<long long> > & RowStencil,
 
  255        const vector<long long> & RowIndices,
 
  259    return TGenerateBlockGraph<long long>(BaseGraph, RowStencil, RowIndices, GlobalComm);
 
  261    throw "EpetraExt::BlockUtility::GenerateBlockGraph: Error Global Indices not long long.";
 
  266template<
typename int_type>
 
  269        const vector< vector<int_type> > & RowStencil,
 
  270        const vector<int_type> & RowIndices,
 
  276  int_type BaseIndex = (int_type) BaseMap.IndexBase64();
 
  277  int_type Offset = BlockUtility::TCalculateOffset<int_type>(BaseMap);
 
  280  int NumBlockRows = RowIndices.size();
 
  282  int TotalSize = NumBlockRows * Size;
 
  283  vector<int_type> GIDs(Size);
 
  286  vector<int_type> GlobalGIDs( TotalSize );
 
  287  for( 
int i = 0; i < NumBlockRows; ++i )
 
  289    for( 
int j = 0; j < Size; ++j )
 
  290      GlobalGIDs[i*Size+j] = GIDs[j] + RowIndices[i] * Offset;
 
  294  int_type TotalSize_int_type = TotalSize;
 
  295  GlobalComm.
SumAll( &TotalSize_int_type, &GlobalSize, 1 );
 
  297  Epetra_Map GlobalMap( GlobalSize, TotalSize, &GlobalGIDs[0], BaseIndex, GlobalComm );
 
  300  vector<int> Indices_local(MaxIndices);
 
  301  vector<int_type> Indices_global(MaxIndices);
 
  302  vector<double> Values(MaxIndices);
 
  309  for( 
int i = 0; i < NumBlockRows; ++i )
 
  311    int StencilSize = RowStencil[i].size();
 
  312    for( 
int j = 0; j < Size; ++j )
 
  314      int_type GlobalRow = (int_type) GlobalMap.GID64(j+i*Size);
 
  316      BaseMatrix.
ExtractMyRowCopy( j, MaxIndices, NumIndices, &Values[0], &Indices_local[0] );
 
  317      for( 
int l = 0; l < NumIndices; ++l ) Indices_global[l] = (int_type) BaseColMap.GID64(Indices_local[l]);
 
  319      for( 
int k = 0; k < StencilSize; ++k )
 
  321        int_type ColOffset = (RowIndices[i]+RowStencil[i][k]) * Offset;
 
  322        if( k > 0 ) ColOffset -= (RowIndices[i]+RowStencil[i][k-1]) * Offset;
 
  324        for( 
int l = 0; l < NumIndices; ++l )
 
  325          Indices_global[l] += ColOffset;
 
 
  337#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 
  340        const vector< vector<int> > & RowStencil,
 
  341        const vector<int> & RowIndices,
 
  345    return TGenerateBlockGraph<int>(BaseMatrix, RowStencil, RowIndices, GlobalComm);
 
  347    throw "EpetraExt::BlockUtility::GenerateBlockGraph: Error Global Indices not int.";
 
  351#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 
  354        const vector< vector<long long> > & RowStencil,
 
  355        const vector<long long> & RowIndices,
 
  359    return TGenerateBlockGraph<long long>(BaseMatrix, RowStencil, RowIndices, GlobalComm);
 
  361    throw "EpetraExt::BlockUtility::GenerateBlockGraph: Error Global Indices not long long.";
 
  366template<
typename int_type>
 
  374  int_type ROffset = BlockUtility::TCalculateOffset<int_type>(BaseRowMap);
 
  376  int_type COffset = BlockUtility::TCalculateOffset<int_type>(BaseColMap);
 
  383  vector<int_type> RowIndices(NumBlockRows);
 
  393  vector<int_type> Indices(MaxIndices);
 
  399  int NumBlockIndices, NumBaseIndices;
 
  400  int *BlockIndices, *BaseIndices;
 
  401  for( 
int i = 0; i < NumBlockRows; ++i )
 
  405    for( 
int j = 0; j < Size; ++j )
 
  407      int_type GlobalRow = (int_type) GlobalRowMap->GID64(j+i*Size);
 
  410      for( 
int k = 0; k < NumBlockIndices; ++k )
 
  412        int_type ColOffset = (int_type) BlockColMap.GID64(BlockIndices[k]) * COffset;
 
  414        for( 
int l = 0; l < NumBaseIndices; ++l )
 
  415          Indices[l] = (int_type) BaseGraph.GCID64(BaseIndices[l]) + ColOffset;
 
  432  GlobalGraph->
FillComplete(*GlobalDomainMap, *GlobalRangeMap);
 
  434  delete GlobalDomainMap;
 
  435  delete GlobalRangeMap;
 
  446#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 
  448    return TGenerateBlockGraph<int>(BaseGraph, LocalBlockGraph, GlobalComm);
 
  451#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 
  453      return TGenerateBlockGraph<long long>(BaseGraph, LocalBlockGraph, GlobalComm);
 
  456    throw "EpetraExt::BlockUtility::GenerateBlockGraph: Error Global Indices unknown.";
 
 
  460template<
typename int_type>
 
  461void BlockUtility::TGenerateRowStencil(
const Epetra_CrsGraph& LocalBlockGraph,
 
  462                                      std::vector<int_type> RowIndices,
 
  463                                      std::vector< std::vector<int_type> >& RowStencil)
 
  466  int NumMyRows = LocalBlockGraph.
NumMyRows();
 
  467  RowIndices.resize(NumMyRows);
 
  472  RowStencil.resize(NumMyRows);
 
  474    for (
int i=0; i<NumMyRows; i++) {
 
  475      int_type Row = RowIndices[i];
 
  477      RowStencil[i].resize(NumCols);
 
  480      for (
int k=0; k<NumCols; k++)
 
  481        RowStencil[i][k] -= Row;
 
  485    for (
int i=0; i<NumMyRows; i++) {
 
  487    std::vector<int> RowStencil_local(NumCols);
 
  488      RowStencil[i].resize(NumCols);
 
  490                                       &RowStencil_local[0]);
 
  491      for (
int k=0; k<NumCols; k++)
 
  492        RowStencil[i][k] = (int_type) ((int) (LocalBlockGraph.GCID64(RowStencil_local[k]) - RowIndices[i]));
 
  497#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 
  499                                      std::vector<int> RowIndices,
 
  500                                      std::vector< std::vector<int> >& RowStencil)
 
  503    BlockUtility::TGenerateRowStencil<int>(LocalBlockGraph, RowIndices, RowStencil);
 
  505    throw "EpetraExt::BlockUtility::GenerateRowStencil: Global Indices not int.";
 
 
  509#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 
  511                                      std::vector<long long> RowIndices,
 
  512                                      std::vector< std::vector<long long> >& RowStencil)
 
  515    BlockUtility::TGenerateRowStencil<long long>(LocalBlockGraph, RowIndices, RowStencil);
 
  517    throw "EpetraExt::BlockUtility::GenerateRowStencil: Global Indices not long long.";
 
 
  523template<
typename int_type>
 
  524int_type BlockUtility::TCalculateOffset(
const Epetra_BlockMap & BaseMap)
 
  526  int_type MaxGID = (int_type) BaseMap.MaxAllGID64();
 
  536#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 
  540    return TCalculateOffset<int>(BaseMap);
 
  542    throw "EpetraExt::BlockUtility::GenerateBlockMap: Global Indices not int.";
 
 
  548  return TCalculateOffset<long long>(BaseMap);
 
 
static void GenerateRowStencil(const Epetra_CrsGraph &LocalBlockGraph, std::vector< int > RowIndices, std::vector< std::vector< int > > &RowStencil)
Generate stencil arrays from a local block graph.
 
static Epetra_CrsGraph * GenerateBlockGraph(const Epetra_CrsGraph &BaseGraph, const std::vector< std::vector< int > > &RowStencil, const std::vector< int > &RowIndices, const Epetra_Comm &GlobalComm)
BlockCrsMatrix constuctor.
 
static int CalculateOffset(const Epetra_BlockMap &BaseMap)
Routine for calculating Offset for creating unique global IDs for Block representation.
 
static Epetra_Map * GenerateBlockMap(const Epetra_BlockMap &BaseMap, const int *RowIndices, int num_indices, const Epetra_Comm &GlobalComm, int Offset=0)
 
static long long CalculateOffset64(const Epetra_BlockMap &BaseMap)
 
int MyGlobalElements(int *MyGlobalElementList) const
 
bool GlobalIndicesInt() const
 
int NumMyElements() const
 
bool GlobalIndicesLongLong() const
 
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const=0
 
int InsertGlobalIndices(int GlobalRow, int NumIndices, int *Indices)
 
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
 
int MaxNumIndices() const
 
const Epetra_BlockMap & DomainMap() const
 
const Epetra_BlockMap & RowMap() const
 
bool IndicesAreGlobal() const
 
int ExtractMyRowCopy(int LocalRow, int LenOfIndices, int &NumIndices, int *Indices) const
 
int NumMyIndices(int Row) const
 
int ExtractGlobalRowCopy(int GlobalRow, int LenOfIndices, int &NumIndices, int *Indices) const
 
int NumGlobalIndices(long long Row) const
 
const Epetra_BlockMap & RangeMap() const
 
const Epetra_BlockMap & ColMap() const
 
virtual const Epetra_Map & RowMatrixColMap() const=0
 
virtual const Epetra_Map & RowMatrixRowMap() const=0
 
virtual int MaxNumEntries() const=0
 
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
 
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.