43 #include "Epetra_Map.h" 44 #include "Epetra_Comm.h" 50 template<
typename int_type>
52 const Epetra_BlockMap & BaseMap,
53 const int_type * RowIndices,
55 const Epetra_Comm & GlobalComm,
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 );
79 Epetra_Map *GlobalMap =
80 new Epetra_Map( GlobalSize, TotalSize, &GlobalGIDs[0], BaseIndex,
87 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 89 const Epetra_BlockMap & BaseMap,
90 const int * RowIndices,
92 const Epetra_Comm & GlobalComm,
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 104 const Epetra_BlockMap & BaseMap,
105 const long long * RowIndices,
107 const Epetra_Comm & GlobalComm,
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 119 const Epetra_BlockMap & BaseMap,
120 const std::vector<int>& RowIndices,
121 const Epetra_Comm & GlobalComm,
129 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 131 const Epetra_BlockMap & BaseMap,
132 const std::vector<long long>& RowIndices,
133 const Epetra_Comm & GlobalComm,
143 const Epetra_BlockMap & BaseMap,
144 const Epetra_BlockMap& BlockMap,
145 const Epetra_Comm & GlobalComm,
148 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 149 if(BaseMap.GlobalIndicesInt() && BlockMap.GlobalIndicesInt())
151 BlockMap.MyGlobalElements(),
152 BlockMap.NumMyElements(),
157 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 158 if(BaseMap.GlobalIndicesLongLong() && BlockMap.GlobalIndicesLongLong())
160 BlockMap.MyGlobalElements64(),
161 BlockMap.NumMyElements(),
166 throw "EpetraExt::BlockUtility::GenerateBlockMap: Error Global Indices unknown.";
170 template<
typename int_type>
172 const Epetra_CrsGraph & BaseGraph,
173 const vector< vector<int_type> > & RowStencil,
174 const vector<int_type> & RowIndices,
175 const Epetra_Comm & GlobalComm )
178 const Epetra_BlockMap & BaseMap = BaseGraph.RowMap();
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);
206 Epetra_CrsGraph * GlobalGraph =
new Epetra_CrsGraph(
Copy,
207 dynamic_cast<Epetra_BlockMap&>(GlobalMap),
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 239 const Epetra_CrsGraph & BaseGraph,
240 const vector< vector<int> > & RowStencil,
241 const vector<int> & RowIndices,
242 const Epetra_Comm & GlobalComm )
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 253 const Epetra_CrsGraph & BaseGraph,
254 const vector< vector<long long> > & RowStencil,
255 const vector<long long> & RowIndices,
256 const Epetra_Comm & GlobalComm )
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.";
266 template<
typename int_type>
268 const Epetra_RowMatrix & BaseMatrix,
269 const vector< vector<int_type> > & RowStencil,
270 const vector<int_type> & RowIndices,
271 const Epetra_Comm & GlobalComm )
274 const Epetra_BlockMap & BaseMap = BaseMatrix.RowMatrixRowMap();
275 const Epetra_BlockMap & BaseColMap = BaseMatrix.RowMatrixColMap();
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);
305 Epetra_CrsGraph * GlobalGraph =
new Epetra_CrsGraph(
Copy,
306 dynamic_cast<Epetra_BlockMap&>(GlobalMap),
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 339 const Epetra_RowMatrix & BaseMatrix,
340 const vector< vector<int> > & RowStencil,
341 const vector<int> & RowIndices,
342 const Epetra_Comm & GlobalComm )
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 353 const Epetra_RowMatrix & BaseMatrix,
354 const vector< vector<long long> > & RowStencil,
355 const vector<long long> & RowIndices,
356 const Epetra_Comm & GlobalComm )
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.";
366 template<
typename int_type>
368 const Epetra_CrsGraph & BaseGraph,
369 const Epetra_CrsGraph & LocalBlockGraph,
370 const Epetra_Comm & GlobalComm )
372 const Epetra_BlockMap & BaseRowMap = BaseGraph.RowMap();
373 const Epetra_BlockMap & BaseColMap = BaseGraph.ColMap();
374 int_type ROffset = BlockUtility::TCalculateOffset<int_type>(BaseRowMap);
376 int_type COffset = BlockUtility::TCalculateOffset<int_type>(BaseColMap);
379 const Epetra_BlockMap & BlockRowMap = LocalBlockGraph.RowMap();
380 const Epetra_BlockMap & BlockColMap = LocalBlockGraph.ColMap();
382 int NumBlockRows = BlockRowMap.NumMyElements();
383 vector<int_type> RowIndices(NumBlockRows);
384 BlockRowMap.MyGlobalElements(&RowIndices[0]);
386 int Size = BaseRowMap.NumMyElements();
388 Epetra_Map *GlobalRowMap =
392 int MaxIndices = BaseGraph.MaxNumIndices();
393 vector<int_type> Indices(MaxIndices);
395 Epetra_CrsGraph * GlobalGraph =
new Epetra_CrsGraph(
Copy,
396 dynamic_cast<Epetra_BlockMap&>(*GlobalRowMap),
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] );
422 const Epetra_BlockMap & BaseDomainMap = BaseGraph.DomainMap();
423 const Epetra_BlockMap & BaseRangeMap = BaseGraph.RangeMap();
424 const Epetra_BlockMap & BlockDomainMap = LocalBlockGraph.DomainMap();
425 const Epetra_BlockMap & BlockRangeMap = LocalBlockGraph.RangeMap();
427 Epetra_Map *GlobalDomainMap =
429 Epetra_Map *GlobalRangeMap =
432 GlobalGraph->FillComplete(*GlobalDomainMap, *GlobalRangeMap);
434 delete GlobalDomainMap;
435 delete GlobalRangeMap;
442 const Epetra_CrsGraph & BaseGraph,
443 const Epetra_CrsGraph & LocalBlockGraph,
444 const Epetra_Comm & GlobalComm )
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.";
460 template<
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);
468 const Epetra_BlockMap& RowMap = LocalBlockGraph.RowMap();
469 RowMap.MyGlobalElements(&RowIndices[0]);
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.";
523 template<
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 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)
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.
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.
static long long CalculateOffset64(const Epetra_BlockMap &BaseMap)
static Epetra_Map * TGenerateBlockMap(const Epetra_BlockMap &BaseMap, const int_type *RowIndices, int num_indices, const Epetra_Comm &GlobalComm, int_type Offset=0)
static int CalculateOffset(const Epetra_BlockMap &BaseMap)
Routine for calculating Offset for creating unique global IDs for Block representation.
static int_type TCalculateOffset(const Epetra_BlockMap &BaseMap)
static Epetra_Map * GenerateBlockMap(const Epetra_BlockMap &BaseMap, const int *RowIndices, int num_indices, const Epetra_Comm &GlobalComm, int Offset=0)
static void TGenerateRowStencil(const Epetra_CrsGraph &LocalBlockGraph, std::vector< int_type > RowIndices, std::vector< std::vector< int_type > > &RowStencil)
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.