44 #include "Epetra_Map.h" 51 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 53 const Epetra_CrsGraph & BaseGraph,
54 const vector<int> & RowStencil,
56 const Epetra_Comm & GlobalComm )
58 BaseGraph_( BaseGraph ),
59 RowStencil_int_( vector< vector<int> >(1,RowStencil) ),
60 RowIndices_int_( vector<int>(1,rowIndex) ),
67 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 69 const Epetra_CrsGraph & BaseGraph,
70 const vector<long long> & RowStencil,
72 const Epetra_Comm & GlobalComm )
73 : Epetra_CrsMatrix( Copy, *(
BlockUtility::GenerateBlockGraph( BaseGraph, vector< vector<long long> >(1,RowStencil), vector<long long>(1,rowIndex), GlobalComm )) ),
74 BaseGraph_( BaseGraph ),
75 RowStencil_LL_( vector< vector<long long> >(1,RowStencil) ),
76 RowIndices_LL_( vector<long long>(1,rowIndex) ),
84 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 86 const Epetra_CrsGraph & BaseGraph,
87 const vector< vector<int> > & RowStencil,
88 const vector<int> & RowIndices,
89 const Epetra_Comm & GlobalComm )
91 BaseGraph_( BaseGraph ),
92 RowStencil_int_( RowStencil ),
93 RowIndices_int_( RowIndices ),
99 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 101 const Epetra_CrsGraph & BaseGraph,
102 const vector< vector<long long> > & RowStencil,
103 const vector<long long> & RowIndices,
104 const Epetra_Comm & GlobalComm )
106 BaseGraph_( BaseGraph ),
107 RowStencil_LL_( RowStencil ),
108 RowIndices_LL_( RowIndices ),
117 const Epetra_CrsGraph & BaseGraph,
118 const Epetra_CrsGraph & LocalBlockGraph,
119 const Epetra_Comm & GlobalComm )
120 : Epetra_CrsMatrix( Copy, *(
BlockUtility::GenerateBlockGraph( BaseGraph, LocalBlockGraph, GlobalComm )) ),
121 BaseGraph_( BaseGraph ),
122 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
126 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
130 ROffset_(
BlockUtility::CalculateOffset64(BaseGraph.RowMap())),
131 COffset_(
BlockUtility::CalculateOffset64(BaseGraph.ColMap()))
133 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 134 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt() && LocalBlockGraph.RowMap().GlobalIndicesInt())
138 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 139 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong() && LocalBlockGraph.RowMap().GlobalIndicesLongLong())
143 throw "EpetraExt::BlockCrsMatrix::BlockCrsMatrix: Error, Global indices unknown.";
147 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 149 const Epetra_RowMatrix & BaseMatrix,
150 const vector< vector<int> > & RowStencil,
151 const vector<int> & RowIndices,
152 const Epetra_Comm & GlobalComm )
154 BaseGraph_( Copy, BaseMatrix.RowMatrixRowMap(), 1 ),
163 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 165 const Epetra_RowMatrix & BaseMatrix,
166 const vector< vector<long long> > & RowStencil,
167 const vector<long long> & RowIndices,
168 const Epetra_Comm & GlobalComm )
170 BaseGraph_( Copy, BaseMatrix.RowMatrixRowMap(), 1 ),
181 : Epetra_CrsMatrix( dynamic_cast<const Epetra_CrsMatrix &>( Matrix ) ),
183 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
187 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
202 template<
typename int_type>
203 void BlockCrsMatrix::TLoadBlock(
const Epetra_RowMatrix & BaseMatrix,
const int_type Row,
const int_type Col)
205 std::vector<int_type>& RowIndices_ = TRowIndices<int_type>();
206 std::vector< std::vector<int_type> >& RowStencil_ = TRowStencil<int_type>();
207 int_type RowOffset = RowIndices_[(std::size_t)Row] *
ROffset_;
208 int_type ColOffset = (RowIndices_[(std::size_t)Row] + RowStencil_[(std::size_t)Row][(std::size_t)Col]) *
COffset_;
211 const Epetra_BlockMap & BaseMap = BaseMatrix.RowMatrixRowMap();
212 const Epetra_BlockMap & BaseColMap = BaseMatrix.RowMatrixColMap();
218 int MaxIndices = BaseMatrix.MaxNumEntries();
219 vector<int> Indices_local(MaxIndices);
220 vector<int_type> Indices_global(MaxIndices);
221 vector<double> vals(MaxIndices);
225 for (
int i=0; i<BaseMap.NumMyElements(); i++) {
226 BaseMatrix.ExtractMyRowCopy( i, MaxIndices, NumIndices, &vals[0], &Indices_local[0] );
229 for(
int l = 0; l < NumIndices; ++l )
230 Indices_global[l] = ColOffset + (int_type) BaseColMap.GID64(Indices_local[l]);
232 int_type BaseRow = (int_type) BaseMap.GID64(i);
233 ierr = this->ReplaceGlobalValues(BaseRow + RowOffset, NumIndices, &vals[0], &Indices_global[0]);
234 if (ierr != 0) std::cout <<
"WARNING BlockCrsMatrix::LoadBlock ReplaceGlobalValues err = " << ierr <<
235 "\n\t Row " << BaseRow + RowOffset <<
"Col start" << Indices_global[0] << std::endl;
240 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 243 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt() && BaseMatrix.RowMatrixRowMap().GlobalIndicesInt())
244 return TLoadBlock<int>(BaseMatrix, Row, Col);
246 throw "EpetraExt::BlockCrsMatrix::LoadBlock: Global indices not int";
250 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 253 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong() && BaseMatrix.RowMatrixRowMap().GlobalIndicesLongLong())
254 return TLoadBlock<long long>(BaseMatrix, Row, Col);
256 throw "EpetraExt::BlockCrsMatrix::LoadBlock: Global indices not long long";
261 template<
typename int_type>
262 void BlockCrsMatrix::TSumIntoBlock(
double alpha,
const Epetra_RowMatrix & BaseMatrix,
const int_type Row,
const int_type Col)
264 std::vector<int_type>& RowIndices_ = TRowIndices<int_type>();
265 std::vector< std::vector<int_type> >& RowStencil_ = TRowStencil<int_type>();
266 int_type RowOffset = RowIndices_[(std::size_t)Row] *
ROffset_;
267 int_type ColOffset = (RowIndices_[(std::size_t)Row] + RowStencil_[(std::size_t)Row][(std::size_t)Col]) *
COffset_;
270 const Epetra_BlockMap & BaseMap = BaseMatrix.RowMatrixRowMap();
271 const Epetra_BlockMap & BaseColMap = BaseMatrix.RowMatrixColMap();
277 int MaxIndices = BaseMatrix.MaxNumEntries();
278 vector<int> Indices_local(MaxIndices);
279 vector<int_type> Indices_global(MaxIndices);
280 vector<double> vals(MaxIndices);
284 for (
int i=0; i<BaseMap.NumMyElements(); i++) {
285 BaseMatrix.ExtractMyRowCopy( i, MaxIndices, NumIndices, &vals[0], &Indices_local[0] );
288 for(
int l = 0; l < NumIndices; ++l ) {
289 Indices_global[l] = ColOffset + (int_type) BaseColMap.GID64(Indices_local[l]);
293 int_type BaseRow = (int_type) BaseMap.GID64(i);
294 ierr = this->SumIntoGlobalValues(BaseRow + RowOffset, NumIndices, &vals[0], &Indices_global[0]);
296 std::cout <<
"WARNING BlockCrsMatrix::SumIntoBlock SumIntoGlobalValues " 297 "err = " << ierr << std::endl <<
"\t Row " << BaseRow + RowOffset <<
298 "Col start" << Indices_global[0] << std::endl;
303 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 306 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt() && BaseMatrix.RowMatrixRowMap().GlobalIndicesInt())
307 return TSumIntoBlock<int>(alpha, BaseMatrix, Row, Col);
309 throw "EpetraExt::BlockCrsMatrix::SumIntoBlock: Global indices not int";
313 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 316 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong() && BaseMatrix.RowMatrixRowMap().GlobalIndicesLongLong())
317 return TSumIntoBlock<long long>(alpha, BaseMatrix, Row, Col);
319 throw "EpetraExt::BlockCrsMatrix::SumIntoBlock: Global indices not long long";
324 template<
typename int_type>
325 void BlockCrsMatrix::TSumIntoGlobalBlock(
double alpha,
const Epetra_RowMatrix & BaseMatrix,
const int_type Row,
const int_type Col)
327 int_type RowOffset = Row *
ROffset_;
328 int_type ColOffset = Col *
COffset_;
331 const Epetra_BlockMap & BaseMap = BaseMatrix.RowMatrixRowMap();
332 const Epetra_BlockMap & BaseColMap = BaseMatrix.RowMatrixColMap();
338 int MaxIndices = BaseMatrix.MaxNumEntries();
339 vector<int> Indices_local(MaxIndices);
340 vector<int_type> Indices_global(MaxIndices);
341 vector<double> vals(MaxIndices);
345 for (
int i=0; i<BaseMap.NumMyElements(); i++) {
346 BaseMatrix.ExtractMyRowCopy( i, MaxIndices, NumIndices, &vals[0], &Indices_local[0] );
349 for(
int l = 0; l < NumIndices; ++l ) {
350 Indices_global[l] = ColOffset + (int_type) BaseColMap.GID64(Indices_local[l]);
354 int_type BaseRow = (int_type) BaseMap.GID64(i);
355 ierr = this->SumIntoGlobalValues(BaseRow + RowOffset, NumIndices, &vals[0], &Indices_global[0]);
357 std::cout <<
"WARNING BlockCrsMatrix::SumIntoBlock SumIntoGlobalValues " 358 <<
"err = " << ierr << std::endl
359 <<
"\t Row " << BaseRow + RowOffset
360 <<
" Col start" << Indices_global[0] << std::endl;
365 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 368 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt() && BaseMatrix.RowMatrixRowMap().GlobalIndicesInt())
369 return TSumIntoGlobalBlock<int>(alpha, BaseMatrix, Row, Col);
371 throw "EpetraExt::BlockCrsMatrix::SumIntoGlobalBlock: Global indices not int";
375 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 378 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong() && BaseMatrix.RowMatrixRowMap().GlobalIndicesLongLong())
379 return TSumIntoGlobalBlock<long long>(alpha, BaseMatrix, Row, Col);
381 throw "EpetraExt::BlockCrsMatrix::SumIntoGlobalBlock: Global indices not long long";
387 template<
typename int_type>
388 void BlockCrsMatrix::TBlockSumIntoGlobalValues(
const int_type BaseRow,
int NumIndices,
389 double* vals,
const int_type* Indices,
const int_type Row,
const int_type Col)
392 std::vector<int_type>& RowIndices_ = TRowIndices<int_type>();
393 std::vector< std::vector<int_type> >& RowStencil_ = TRowStencil<int_type>();
394 int_type RowOffset = RowIndices_[(std::size_t)Row] *
ROffset_;
395 int_type ColOffset = (RowIndices_[(std::size_t)Row] + RowStencil_[(std::size_t)Row][(std::size_t)Col]) *
COffset_;
398 vector<int_type> OffsetIndices(NumIndices);
399 for(
int l = 0; l < NumIndices; ++l ) OffsetIndices[l] = Indices[l] + ColOffset;
401 int ierr = this->SumIntoGlobalValues(BaseRow + RowOffset, NumIndices,
402 vals, &OffsetIndices[0]);
405 std::cout <<
"WARNING BlockCrsMatrix::BlockSumIntoGlobalValues err = " 406 << ierr << std::endl <<
"\t Row " << BaseRow + RowOffset
407 <<
" Col start" << Indices[0] << std::endl;
411 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 413 double* vals,
const int* Indices,
const int Row,
const int Col)
415 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt())
416 return TBlockSumIntoGlobalValues<int>(BaseRow, NumIndices, vals, Indices, Row, Col);
418 throw "EpetraExt::BlockCrsMatrix::BlockSumIntoGlobalValues: Global indices not int";
422 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 424 double* vals,
const long long* Indices,
const long long Row,
const long long Col)
426 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong() && Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong())
427 return TBlockSumIntoGlobalValues<long long>(BaseRow, NumIndices, vals, Indices, Row, Col);
429 throw "EpetraExt::BlockCrsMatrix::BlockSumIntoGlobalValues: Global indices not long long";
434 template<
typename int_type>
435 void BlockCrsMatrix::TBlockReplaceGlobalValues(
const int_type BaseRow,
int NumIndices,
436 double* vals,
const int_type* Indices,
const int_type Row,
const int_type Col)
439 std::vector<int_type>& RowIndices_ = TRowIndices<int_type>();
440 std::vector< std::vector<int_type> >& RowStencil_ = TRowStencil<int_type>();
441 int_type RowOffset = RowIndices_[(std::size_t)Row] *
ROffset_;
442 int_type ColOffset = (RowIndices_[(std::size_t)Row] + RowStencil_[(std::size_t)Row][(std::size_t)Col]) *
COffset_;
445 vector<int_type> OffsetIndices(NumIndices);
446 for(
int l = 0; l < NumIndices; ++l ) OffsetIndices[l] = Indices[l] + ColOffset;
448 int ierr = this->ReplaceGlobalValues(BaseRow + RowOffset, NumIndices,
449 vals, &OffsetIndices[0]);
452 std::cout <<
"WARNING BlockCrsMatrix::BlockReplaceGlobalValues err = " 453 << ierr <<
"\n\t Row " << BaseRow + RowOffset <<
"Col start" 454 << Indices[0] << std::endl;
458 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 460 double* vals,
const int* Indices,
const int Row,
const int Col)
462 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt())
463 return TBlockReplaceGlobalValues<int>(BaseRow, NumIndices, vals, Indices, Row, Col);
465 throw "EpetraExt::BlockCrsMatrix::BlockReplaceGlobalValues: Global indices not int";
469 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 471 double* vals,
const long long* Indices,
const long long Row,
const long long Col)
473 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong())
474 return TBlockReplaceGlobalValues<long long>(BaseRow, NumIndices, vals, Indices, Row, Col);
476 throw "EpetraExt::BlockCrsMatrix::BlockReplaceGlobalValues: Global indices not long long";
481 template<
typename int_type>
482 void BlockCrsMatrix::TBlockExtractGlobalRowView(
const int_type BaseRow,
489 std::vector<int_type>& RowIndices_ = TRowIndices<int_type>();
490 std::vector< std::vector<int_type> >& RowStencil_ = TRowStencil<int_type>();
491 int_type RowOffset = RowIndices_[(std::size_t)Row] *
ROffset_;
492 int_type ColOffset = (RowIndices_[(std::size_t)Row] + RowStencil_[(std::size_t)Row][(std::size_t)Col]) *
COffset_;
495 int ierr = this->ExtractGlobalRowView(BaseRow + RowOffset, NumEntries,
500 NumEntries -= ColOffset;
503 std::cout <<
"WARNING BlockCrsMatrix::BlockExtractGlobalRowView err = " 504 << ierr <<
"\n\t Row " << BaseRow + RowOffset
505 <<
" Col " << Col+ColOffset << std::endl;
509 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 511 double*& vals,
const int Row,
const int Col)
513 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt())
514 return TBlockExtractGlobalRowView<int>(BaseRow, NumEntries, vals, Row, Col);
516 throw "EpetraExt::BlockCrsMatrix::BlockExtractGlobalRowView: Global indices not int";
520 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 522 double*& vals,
const long long Row,
const long long Col)
524 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong())
525 return TBlockExtractGlobalRowView<long long>(BaseRow, NumEntries, vals, Row, Col);
527 throw "EpetraExt::BlockCrsMatrix::BlockExtractGlobalRowView: Global indices not long long";
533 template<
typename int_type>
534 void BlockCrsMatrix::TExtractBlock(Epetra_CrsMatrix & BaseMatrix,
const int_type Row,
const int_type Col)
536 std::vector<int_type>& RowIndices_ = TRowIndices<int_type>();
537 std::vector< std::vector<int_type> >& RowStencil_ = TRowStencil<int_type>();
538 int_type RowOffset = RowIndices_[(std::size_t)Row] *
ROffset_;
539 int_type ColOffset = (RowIndices_[(std::size_t)Row] + RowStencil_[(std::size_t)Row][(std::size_t)Col]) *
COffset_;
542 const Epetra_BlockMap & BaseMap = BaseMatrix.RowMatrixRowMap();
549 int MaxIndices = BaseMatrix.MaxNumEntries();
550 vector<int_type> Indices(MaxIndices);
551 vector<double> vals(MaxIndices);
560 for (
int i=0; i<BaseMap.NumMyElements(); i++) {
563 int_type BaseRow = (int_type) BaseMap.GID64(i);
564 int myBlkBaseRow = this->RowMatrixRowMap().LID(BaseRow + RowOffset);
565 ierr = this->ExtractMyRowView(myBlkBaseRow, BlkNumIndices, BlkValues, BlkIndices);
569 for(
int l = 0; l < BlkNumIndices; ++l ) {
570 icol = (int_type) this->RowMatrixColMap().GID64(BlkIndices[l]);
571 indx = icol - ColOffset;
573 Indices[NumIndices] = indx;
574 vals[NumIndices] = BlkValues[l];
580 BaseMatrix.ReplaceGlobalValues(BaseRow, NumIndices, &vals[0], &Indices[0]);
584 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 587 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt() && BaseMatrix.RowMatrixRowMap().GlobalIndicesInt())
588 return TExtractBlock<int>(BaseMatrix, Row, Col);
590 throw "EpetraExt::BlockCrsMatrix::ExtractBlock: Global indices not int";
594 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 597 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong() && BaseMatrix.RowMatrixRowMap().GlobalIndicesLongLong())
598 return TExtractBlock<long long>(BaseMatrix, Row, Col);
600 throw "EpetraExt::BlockCrsMatrix::ExtractBlock: Global indices not long long";
void LoadBlock(const Epetra_RowMatrix &BaseMatrix, const int Row, const int Col)
Routine for loading a base matrices values into the large Block Matrix The Row and Col arguments are ...
std::vector< std::vector< int > > RowStencil_int_
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.
void BlockReplaceGlobalValues(const int BaseRow, int NumIndices, double *Values, const int *Indices, const int Row, const int Col)
std::vector< long long > RowIndices_LL_
void ExtractBlock(Epetra_CrsMatrix &BaseMatrix, const int Row, const int Col)
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.
static long long CalculateOffset64(const Epetra_BlockMap &BaseMap)
void BlockExtractGlobalRowView(const int BaseRow, int &NumEntries, double *&Values, const int Row, const int Col)
std::vector< std::vector< long long > > RowStencil_LL_
std::vector< int > RowIndices_int_
void SumIntoGlobalBlock(double alpha, const Epetra_RowMatrix &BaseMatrix, const int Row, const int Col)
Routine for summing base matrices values into the large Block Matrix The Row and Col arguments are gl...
virtual ~BlockCrsMatrix()
Destructor.
BlockCrsMatrix(const Epetra_CrsGraph &BaseGraph, const std::vector< int > &RowStencil, int RowIndex, const Epetra_Comm &GlobalComm)
BlockCrsMatrix constuctor with one block row per processor.
void SumIntoBlock(double alpha, const Epetra_RowMatrix &BaseMatrix, const int Row, const int Col)
Routine for summing base matrices values into the large Block Matrix The Row and Col arguments are in...
Epetra_CrsGraph BaseGraph_
void BlockSumIntoGlobalValues(const int BaseRow, int NumIndices, double *Values, const int *Indices, const int Row, const int Col)
Sum Entries into Block matrix using base-matrix numbering plus block Row and Col The Row and Col argu...
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.