50template<
typename int_type>
53 const int_type * RowIndices,
58 int_type BaseIndex = (int_type) BaseMap.IndexBase64();
60 Offset = BlockUtility::TCalculateOffset<int_type>(BaseMap);
63 int Size = BaseMap.NumMyElements();
64 int TotalSize = NumBlockRows * Size;
65 vector<int_type> GIDs(Size);
66 BaseMap.MyGlobalElements( &GIDs[0] );
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,
95 if(BaseMap.GlobalIndicesInt())
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,
110 if(BaseMap.GlobalIndicesLongLong())
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
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();
184 int Size = BaseMap.NumMyElements();
185 int TotalSize = NumBlockRows * Size;
186 vector<int_type> GIDs(Size);
187 BaseMap.MyGlobalElements( &GIDs[0] );
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 );
202 int MaxIndices = BaseGraph.MaxNumIndices();
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);
218 BaseGraph.ExtractGlobalRowCopy( BaseRow, MaxIndices, NumIndices, &Indices[0] );
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;
227 GlobalGraph->InsertGlobalIndices( GlobalRow, NumIndices, &Indices[0] );
232 GlobalGraph->FillComplete();
237#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
240 const vector< vector<int> > & RowStencil,
241 const vector<int> & RowIndices,
244 if(BaseGraph.RowMap().GlobalIndicesInt())
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,
258 if(BaseGraph.RowMap().GlobalIndicesLongLong())
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();
281 int Size = BaseMap.NumMyElements();
282 int TotalSize = NumBlockRows * Size;
283 vector<int_type> GIDs(Size);
284 BaseMap.MyGlobalElements( &GIDs[0] );
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 );
299 int MaxIndices = BaseMatrix.MaxNumEntries();
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;
327 GlobalGraph->InsertGlobalIndices( GlobalRow, NumIndices, &Indices_global[0] );
332 GlobalGraph->FillComplete();
337#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
340 const vector< vector<int> > & RowStencil,
341 const vector<int> & RowIndices,
344 if(BaseMatrix.RowMatrixRowMap().GlobalIndicesInt())
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,
358 if(BaseMatrix.RowMatrixRowMap().GlobalIndicesLongLong())
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);
384 BlockRowMap.MyGlobalElements(&RowIndices[0]);
386 int Size = BaseRowMap.NumMyElements();
392 int MaxIndices = BaseGraph.MaxNumIndices();
393 vector<int_type> Indices(MaxIndices);
399 int NumBlockIndices, NumBaseIndices;
400 int *BlockIndices, *BaseIndices;
401 for(
int i = 0; i < NumBlockRows; ++i )
403 LocalBlockGraph.ExtractMyRowView(i, NumBlockIndices, BlockIndices);
405 for(
int j = 0; j < Size; ++j )
407 int_type GlobalRow = (int_type) GlobalRowMap->GID64(j+i*Size);
409 BaseGraph.ExtractMyRowView( j, NumBaseIndices, BaseIndices );
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;
417 GlobalGraph->InsertGlobalIndices( GlobalRow, NumBaseIndices, &Indices[0] );
432 GlobalGraph->FillComplete(*GlobalDomainMap, *GlobalRangeMap);
434 delete GlobalDomainMap;
435 delete GlobalRangeMap;
446#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
447 if(BaseGraph.RowMap().GlobalIndicesInt() && LocalBlockGraph.RowMap().GlobalIndicesInt())
448 return TGenerateBlockGraph<int>(BaseGraph, LocalBlockGraph, GlobalComm);
451#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
452 if(BaseGraph.RowMap().GlobalIndicesLongLong() && LocalBlockGraph.RowMap().GlobalIndicesLongLong())
453 return TGenerateBlockGraph<long long>(BaseGraph, LocalBlockGraph, GlobalComm);
456 throw "EpetraExt::BlockUtility::GenerateBlockGraph: Error Global Indices unknown.";
460template<
typename int_type>
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);
473 if (LocalBlockGraph.IndicesAreGlobal()) {
474 for (
int i=0; i<NumMyRows; i++) {
475 int_type Row = RowIndices[i];
476 int NumCols = LocalBlockGraph.NumGlobalIndices(Row);
477 RowStencil[i].resize(NumCols);
478 LocalBlockGraph.ExtractGlobalRowCopy(Row, NumCols, NumCols,
480 for (
int k=0; k<NumCols; k++)
481 RowStencil[i][k] -= Row;
485 for (
int i=0; i<NumMyRows; i++) {
486 int NumCols = LocalBlockGraph.NumMyIndices(i);
487 std::vector<int> RowStencil_local(NumCols);
488 RowStencil[i].resize(NumCols);
489 LocalBlockGraph.ExtractMyRowCopy(i, NumCols, 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)
502 if(LocalBlockGraph.RowMap().GlobalIndicesInt())
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)
514 if(LocalBlockGraph.RowMap().GlobalIndicesLongLong())
515 BlockUtility::TGenerateRowStencil<long long>(LocalBlockGraph, RowIndices, RowStencil);
517 throw "EpetraExt::BlockUtility::GenerateRowStencil: Global Indices not long long.";
523template<
typename int_type>
526 int_type MaxGID = (int_type) BaseMap.MaxAllGID64();
536#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
539 if(BaseMap.GlobalIndicesInt())
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)
static int_type TCalculateOffset(const Epetra_BlockMap &BaseMap)
static void TGenerateRowStencil(const Epetra_CrsGraph &LocalBlockGraph, std::vector< int_type > RowIndices, std::vector< std::vector< int_type > > &RowStencil)
static Epetra_Map * TGenerateBlockMap(const Epetra_BlockMap &BaseMap, const int_type *RowIndices, int num_indices, const Epetra_Comm &GlobalComm, int_type Offset=0)
static Epetra_CrsGraph * TGenerateBlockGraph(const Epetra_CrsGraph &BaseGraph, const std::vector< std::vector< int_type > > &RowStencil, const std::vector< int_type > &RowIndices, const Epetra_Comm &GlobalComm)
long long IndexBase64() const
int MyGlobalElements(int *MyGlobalElementList) const
long long * MyGlobalElements64() const
bool GlobalIndicesInt() const
int NumMyElements() const
bool GlobalIndicesLongLong() const
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.