EpetraExt Development
Loading...
Searching...
No Matches
EpetraExt_BlockCrsMatrix.cpp
Go to the documentation of this file.
1//@HEADER
2// ***********************************************************************
3//
4// EpetraExt: Epetra Extended - Linear Algebra Services Package
5// Copyright (2011) Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41
42#include "EpetraExt_BlockCrsMatrix.h"
43#include "EpetraExt_BlockUtility.h"
44#include "Epetra_Map.h"
45
46namespace EpetraExt {
47
48using std::vector;
49
50//==============================================================================
51#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
52BlockCrsMatrix::BlockCrsMatrix(
53 const Epetra_CrsGraph & BaseGraph,
54 const vector<int> & RowStencil,
55 int rowIndex,
56 const Epetra_Comm & GlobalComm )
57 : Epetra_CrsMatrix( Copy, *(BlockUtility::GenerateBlockGraph( BaseGraph, vector< vector<int> >(1,RowStencil), vector<int>(1,rowIndex), GlobalComm )) ),
58 BaseGraph_( BaseGraph ),
59 RowStencil_int_( vector< vector<int> >(1,RowStencil) ),
60 RowIndices_int_( vector<int>(1,rowIndex) ),
61 ROffset_(BlockUtility::CalculateOffset64(BaseGraph.RowMap())),
62 COffset_(BlockUtility::CalculateOffset64(BaseGraph.ColMap()))
63{
64}
65#endif
66
67#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
68BlockCrsMatrix::BlockCrsMatrix(
69 const Epetra_CrsGraph & BaseGraph,
70 const vector<long long> & RowStencil,
71 long long rowIndex,
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) ),
77 ROffset_(BlockUtility::CalculateOffset64(BaseGraph.RowMap())),
78 COffset_(BlockUtility::CalculateOffset64(BaseGraph.ColMap()))
79{
80}
81#endif
82
83//==============================================================================
84#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
85BlockCrsMatrix::BlockCrsMatrix(
86 const Epetra_CrsGraph & BaseGraph,
87 const vector< vector<int> > & RowStencil,
88 const vector<int> & RowIndices,
89 const Epetra_Comm & GlobalComm )
90 : Epetra_CrsMatrix( Copy, *(BlockUtility::GenerateBlockGraph( BaseGraph, RowStencil, RowIndices, GlobalComm )) ),
91 BaseGraph_( BaseGraph ),
92 RowStencil_int_( RowStencil ),
93 RowIndices_int_( RowIndices ),
94 ROffset_(BlockUtility::CalculateOffset64(BaseGraph.RowMap())),
95 COffset_(BlockUtility::CalculateOffset64(BaseGraph.ColMap()))
96{
97}
98#endif
99#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
100BlockCrsMatrix::BlockCrsMatrix(
101 const Epetra_CrsGraph & BaseGraph,
102 const vector< vector<long long> > & RowStencil,
103 const vector<long long> & RowIndices,
104 const Epetra_Comm & GlobalComm )
105 : Epetra_CrsMatrix( Copy, *(BlockUtility::GenerateBlockGraph( BaseGraph, RowStencil, RowIndices, GlobalComm )) ),
106 BaseGraph_( BaseGraph ),
107 RowStencil_LL_( RowStencil ),
108 RowIndices_LL_( RowIndices ),
109 ROffset_(BlockUtility::CalculateOffset64(BaseGraph.RowMap())),
110 COffset_(BlockUtility::CalculateOffset64(BaseGraph.ColMap()))
111{
112}
113#endif
114
115//==============================================================================
116BlockCrsMatrix::BlockCrsMatrix(
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
123 RowStencil_int_( ),
124 RowIndices_int_( ),
125#endif
126#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
127 RowStencil_LL_( ),
128 RowIndices_LL_( ),
129#endif
130 ROffset_(BlockUtility::CalculateOffset64(BaseGraph.RowMap())),
131 COffset_(BlockUtility::CalculateOffset64(BaseGraph.ColMap()))
132{
133#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
134 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt() && LocalBlockGraph.RowMap().GlobalIndicesInt())
135 BlockUtility::GenerateRowStencil(LocalBlockGraph, RowIndices_int_, RowStencil_int_);
136 else
137#endif
138#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
139 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong() && LocalBlockGraph.RowMap().GlobalIndicesLongLong())
140 BlockUtility::GenerateRowStencil(LocalBlockGraph, RowIndices_LL_, RowStencil_LL_);
141 else
142#endif
143 throw "EpetraExt::BlockCrsMatrix::BlockCrsMatrix: Error, Global indices unknown.";
144}
145
146//==============================================================================
147#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
148BlockCrsMatrix::BlockCrsMatrix(
149 const Epetra_RowMatrix & BaseMatrix,
150 const vector< vector<int> > & RowStencil,
151 const vector<int> & RowIndices,
152 const Epetra_Comm & GlobalComm )
153 : Epetra_CrsMatrix( Copy, *(BlockUtility::GenerateBlockGraph( BaseMatrix, RowStencil, RowIndices, GlobalComm )) ),
154 BaseGraph_( Copy, BaseMatrix.RowMatrixRowMap(), 1 ), //Junk to satisfy constructor
155 RowStencil_int_( RowStencil ),
156 RowIndices_int_( RowIndices ),
157 ROffset_(BlockUtility::CalculateOffset64(BaseMatrix.RowMatrixRowMap())),
158 COffset_(BlockUtility::CalculateOffset64(BaseMatrix.RowMatrixColMap()))
159{
160}
161#endif
162
163#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
164BlockCrsMatrix::BlockCrsMatrix(
165 const Epetra_RowMatrix & BaseMatrix,
166 const vector< vector<long long> > & RowStencil,
167 const vector<long long> & RowIndices,
168 const Epetra_Comm & GlobalComm )
169 : Epetra_CrsMatrix( Copy, *(BlockUtility::GenerateBlockGraph( BaseMatrix, RowStencil, RowIndices, GlobalComm )) ),
170 BaseGraph_( Copy, BaseMatrix.RowMatrixRowMap(), 1 ), //Junk to satisfy constructor
171 RowStencil_LL_( RowStencil ),
172 RowIndices_LL_( RowIndices ),
173 ROffset_(BlockUtility::CalculateOffset64(BaseMatrix.RowMatrixRowMap())),
174 COffset_(BlockUtility::CalculateOffset64(BaseMatrix.RowMatrixColMap()))
175{
176}
177#endif
178
179//==============================================================================
180BlockCrsMatrix::BlockCrsMatrix( const BlockCrsMatrix & Matrix )
181 : Epetra_CrsMatrix( dynamic_cast<const Epetra_CrsMatrix &>( Matrix ) ),
182 BaseGraph_( Matrix.BaseGraph_ ),
183#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
184 RowStencil_int_( Matrix.RowStencil_int_ ),
185 RowIndices_int_( Matrix.RowIndices_int_ ),
186#endif
187#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
188 RowStencil_LL_( Matrix.RowStencil_LL_ ),
189 RowIndices_LL_( Matrix.RowIndices_LL_ ),
190#endif
191 ROffset_( Matrix.ROffset_ ),
192 COffset_( Matrix.COffset_ )
193{
194}
195
196//==============================================================================
200
201//==============================================================================
202template<typename int_type>
203void BlockCrsMatrix::TLoadBlock(const Epetra_RowMatrix & BaseMatrix, const int_type Row, const int_type Col)
204{
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_;
209
210// const Epetra_CrsGraph & BaseGraph = BaseMatrix.Graph();
211 const Epetra_BlockMap & BaseMap = BaseMatrix.RowMatrixRowMap();
212 const Epetra_BlockMap & BaseColMap = BaseMatrix.RowMatrixColMap();
213
214 // This routine copies entries of a BaseMatrix into big BlockCrsMatrix
215 // It performs the following operation on the global IDs row-by-row
216 // this->val[i+rowOffset][j+ColOffset] = BaseMatrix.val[i][j]
217
218 int MaxIndices = BaseMatrix.MaxNumEntries();
219 vector<int> Indices_local(MaxIndices);
220 vector<int_type> Indices_global(MaxIndices);
221 vector<double> vals(MaxIndices);
222 int NumIndices;
223 int ierr=0;
224
225 for (int i=0; i<BaseMap.NumMyElements(); i++) {
226 BaseMatrix.ExtractMyRowCopy( i, MaxIndices, NumIndices, &vals[0], &Indices_local[0] );
227
228 // Convert to BlockMatrix Global numbering scheme
229 for( int l = 0; l < NumIndices; ++l )
230 Indices_global[l] = ColOffset + (int_type) BaseColMap.GID64(Indices_local[l]);
231
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;
236
237 }
238}
239
240#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
241void BlockCrsMatrix::LoadBlock(const Epetra_RowMatrix & BaseMatrix, const int Row, const int Col)
242{
243 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt() && BaseMatrix.RowMatrixRowMap().GlobalIndicesInt())
244 return TLoadBlock<int>(BaseMatrix, Row, Col);
245 else
246 throw "EpetraExt::BlockCrsMatrix::LoadBlock: Global indices not int";
247}
248#endif
249
250#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
251void BlockCrsMatrix::LoadBlock(const Epetra_RowMatrix & BaseMatrix, const long long Row, const long long Col)
252{
253 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong() && BaseMatrix.RowMatrixRowMap().GlobalIndicesLongLong())
254 return TLoadBlock<long long>(BaseMatrix, Row, Col);
255 else
256 throw "EpetraExt::BlockCrsMatrix::LoadBlock: Global indices not long long";
257}
258#endif
259
260//==============================================================================
261template<typename int_type>
262void BlockCrsMatrix::TSumIntoBlock(double alpha, const Epetra_RowMatrix & BaseMatrix, const int_type Row, const int_type Col)
263{
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_;
268
269// const Epetra_CrsGraph & BaseGraph = BaseMatrix.Graph();
270 const Epetra_BlockMap & BaseMap = BaseMatrix.RowMatrixRowMap();
271 const Epetra_BlockMap & BaseColMap = BaseMatrix.RowMatrixColMap();
272
273 // This routine copies entries of a BaseMatrix into big BlockCrsMatrix
274 // It performs the following operation on the global IDs row-by-row
275 // this->val[i+rowOffset][j+ColOffset] = BaseMatrix.val[i][j]
276
277 int MaxIndices = BaseMatrix.MaxNumEntries();
278 vector<int> Indices_local(MaxIndices);
279 vector<int_type> Indices_global(MaxIndices);
280 vector<double> vals(MaxIndices);
281 int NumIndices;
282 int ierr=0;
283
284 for (int i=0; i<BaseMap.NumMyElements(); i++) {
285 BaseMatrix.ExtractMyRowCopy( i, MaxIndices, NumIndices, &vals[0], &Indices_local[0] );
286
287 // Convert to BlockMatrix Global numbering scheme
288 for( int l = 0; l < NumIndices; ++l ) {
289 Indices_global[l] = ColOffset + (int_type) BaseColMap.GID64(Indices_local[l]);
290 vals[l] *= alpha;
291 }
292
293 int_type BaseRow = (int_type) BaseMap.GID64(i);
294 ierr = this->SumIntoGlobalValues(BaseRow + RowOffset, NumIndices, &vals[0], &Indices_global[0]);
295 if (ierr != 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;
299 }
300 }
301}
302
303#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
304void BlockCrsMatrix::SumIntoBlock(double alpha, const Epetra_RowMatrix & BaseMatrix, const int Row, const int Col)
305{
306 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt() && BaseMatrix.RowMatrixRowMap().GlobalIndicesInt())
307 return TSumIntoBlock<int>(alpha, BaseMatrix, Row, Col);
308 else
309 throw "EpetraExt::BlockCrsMatrix::SumIntoBlock: Global indices not int";
310}
311#endif
312
313#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
314void BlockCrsMatrix::SumIntoBlock(double alpha, const Epetra_RowMatrix & BaseMatrix, const long long Row, const long long Col)
315{
316 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong() && BaseMatrix.RowMatrixRowMap().GlobalIndicesLongLong())
317 return TSumIntoBlock<long long>(alpha, BaseMatrix, Row, Col);
318 else
319 throw "EpetraExt::BlockCrsMatrix::SumIntoBlock: Global indices not long long";
320}
321#endif
322
323//==============================================================================
324template<typename int_type>
325void BlockCrsMatrix::TSumIntoGlobalBlock(double alpha, const Epetra_RowMatrix & BaseMatrix, const int_type Row, const int_type Col)
326{
327 int_type RowOffset = Row * ROffset_;
328 int_type ColOffset = Col * COffset_;
329
330// const Epetra_CrsGraph & BaseGraph = BaseMatrix.Graph();
331 const Epetra_BlockMap & BaseMap = BaseMatrix.RowMatrixRowMap();
332 const Epetra_BlockMap & BaseColMap = BaseMatrix.RowMatrixColMap();
333
334 // This routine copies entries of a BaseMatrix into big BlockCrsMatrix
335 // It performs the following operation on the global IDs row-by-row
336 // this->val[i+rowOffset][j+ColOffset] = BaseMatrix.val[i][j]
337
338 int MaxIndices = BaseMatrix.MaxNumEntries();
339 vector<int> Indices_local(MaxIndices);
340 vector<int_type> Indices_global(MaxIndices);
341 vector<double> vals(MaxIndices);
342 int NumIndices;
343 int ierr=0;
344
345 for (int i=0; i<BaseMap.NumMyElements(); i++) {
346 BaseMatrix.ExtractMyRowCopy( i, MaxIndices, NumIndices, &vals[0], &Indices_local[0] );
347
348 // Convert to BlockMatrix Global numbering scheme
349 for( int l = 0; l < NumIndices; ++l ) {
350 Indices_global[l] = ColOffset + (int_type) BaseColMap.GID64(Indices_local[l]);
351 vals[l] *= alpha;
352 }
353
354 int_type BaseRow = (int_type) BaseMap.GID64(i);
355 ierr = this->SumIntoGlobalValues(BaseRow + RowOffset, NumIndices, &vals[0], &Indices_global[0]);
356 if (ierr != 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;
361 }
362 }
363}
364
365#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
366void BlockCrsMatrix::SumIntoGlobalBlock(double alpha, const Epetra_RowMatrix & BaseMatrix, const int Row, const int Col)
367{
368 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt() && BaseMatrix.RowMatrixRowMap().GlobalIndicesInt())
369 return TSumIntoGlobalBlock<int>(alpha, BaseMatrix, Row, Col);
370 else
371 throw "EpetraExt::BlockCrsMatrix::SumIntoGlobalBlock: Global indices not int";
372}
373#endif
374
375#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
376void BlockCrsMatrix::SumIntoGlobalBlock(double alpha, const Epetra_RowMatrix & BaseMatrix, const long long Row, const long long Col)
377{
378 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong() && BaseMatrix.RowMatrixRowMap().GlobalIndicesLongLong())
379 return TSumIntoGlobalBlock<long long>(alpha, BaseMatrix, Row, Col);
380 else
381 throw "EpetraExt::BlockCrsMatrix::SumIntoGlobalBlock: Global indices not long long";
382}
383#endif
384
385
386//==============================================================================
387template<typename int_type>
388void BlockCrsMatrix::TBlockSumIntoGlobalValues(const int_type BaseRow, int NumIndices,
389 double* vals, const int_type* Indices, const int_type Row, const int_type Col)
390//All arguments could be const, except some were not set as const in CrsMatrix
391{
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_;
396
397 // Convert to BlockMatrix Global numbering scheme
398 vector<int_type> OffsetIndices(NumIndices);
399 for( int l = 0; l < NumIndices; ++l ) OffsetIndices[l] = Indices[l] + ColOffset;
400
401 int ierr = this->SumIntoGlobalValues(BaseRow + RowOffset, NumIndices,
402 vals, &OffsetIndices[0]);
403
404 if (ierr != 0) {
405 std::cout << "WARNING BlockCrsMatrix::BlockSumIntoGlobalValues err = "
406 << ierr << std::endl << "\t Row " << BaseRow + RowOffset
407 << " Col start" << Indices[0] << std::endl;
408 }
409}
410
411#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
412void BlockCrsMatrix::BlockSumIntoGlobalValues(const int BaseRow, int NumIndices,
413 double* vals, const int* Indices, const int Row, const int Col)
414{
415 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt())
416 return TBlockSumIntoGlobalValues<int>(BaseRow, NumIndices, vals, Indices, Row, Col);
417 else
418 throw "EpetraExt::BlockCrsMatrix::BlockSumIntoGlobalValues: Global indices not int";
419}
420#endif
421
422#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
423void BlockCrsMatrix::BlockSumIntoGlobalValues(const long long BaseRow, int NumIndices,
424 double* vals, const long long* Indices, const long long Row, const long long Col)
425{
426 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong() && Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong())
427 return TBlockSumIntoGlobalValues<long long>(BaseRow, NumIndices, vals, Indices, Row, Col);
428 else
429 throw "EpetraExt::BlockCrsMatrix::BlockSumIntoGlobalValues: Global indices not long long";
430}
431#endif
432
433//==============================================================================
434template<typename int_type>
435void BlockCrsMatrix::TBlockReplaceGlobalValues(const int_type BaseRow, int NumIndices,
436 double* vals, const int_type* Indices, const int_type Row, const int_type Col)
437//All arguments could be const, except some were not set as const in CrsMatrix
438{
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_;
443
444 // Convert to BlockMatrix Global numbering scheme
445 vector<int_type> OffsetIndices(NumIndices);
446 for( int l = 0; l < NumIndices; ++l ) OffsetIndices[l] = Indices[l] + ColOffset;
447
448 int ierr = this->ReplaceGlobalValues(BaseRow + RowOffset, NumIndices,
449 vals, &OffsetIndices[0]);
450
451 if (ierr != 0) {
452 std::cout << "WARNING BlockCrsMatrix::BlockReplaceGlobalValues err = "
453 << ierr << "\n\t Row " << BaseRow + RowOffset << "Col start"
454 << Indices[0] << std::endl;
455 }
456}
457
458#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
459void BlockCrsMatrix::BlockReplaceGlobalValues(const int BaseRow, int NumIndices,
460 double* vals, const int* Indices, const int Row, const int Col)
461{
462 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt())
463 return TBlockReplaceGlobalValues<int>(BaseRow, NumIndices, vals, Indices, Row, Col);
464 else
465 throw "EpetraExt::BlockCrsMatrix::BlockReplaceGlobalValues: Global indices not int";
466}
467#endif
468
469#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
470void BlockCrsMatrix::BlockReplaceGlobalValues(const long long BaseRow, int NumIndices,
471 double* vals, const long long* Indices, const long long Row, const long long Col)
472{
473 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong())
474 return TBlockReplaceGlobalValues<long long>(BaseRow, NumIndices, vals, Indices, Row, Col);
475 else
476 throw "EpetraExt::BlockCrsMatrix::BlockReplaceGlobalValues: Global indices not long long";
477}
478#endif
479
480//==============================================================================
481template<typename int_type>
482void BlockCrsMatrix::TBlockExtractGlobalRowView(const int_type BaseRow,
483 int& NumEntries,
484 double*& vals,
485 const int_type Row,
486 const int_type Col)
487//All arguments could be const, except some were not set as const in CrsMatrix
488{
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_;
493
494 // Get the whole row
495 int ierr = this->ExtractGlobalRowView(BaseRow + RowOffset, NumEntries,
496 vals);
497
498 // Adjust for just this block column
499 vals += ColOffset;
500 NumEntries -= ColOffset;
501
502 if (ierr != 0) {
503 std::cout << "WARNING BlockCrsMatrix::BlockExtractGlobalRowView err = "
504 << ierr << "\n\t Row " << BaseRow + RowOffset
505 << " Col " << Col+ColOffset << std::endl;
506 }
507}
508
509#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
510void BlockCrsMatrix::BlockExtractGlobalRowView(const int BaseRow, int& NumEntries,
511 double*& vals, const int Row, const int Col)
512{
513 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt())
514 return TBlockExtractGlobalRowView<int>(BaseRow, NumEntries, vals, Row, Col);
515 else
516 throw "EpetraExt::BlockCrsMatrix::BlockExtractGlobalRowView: Global indices not int";
517}
518#endif
519
520#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
521void BlockCrsMatrix::BlockExtractGlobalRowView(const long long BaseRow, int& NumEntries,
522 double*& vals, const long long Row, const long long Col)
523{
524 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong())
525 return TBlockExtractGlobalRowView<long long>(BaseRow, NumEntries, vals, Row, Col);
526 else
527 throw "EpetraExt::BlockCrsMatrix::BlockExtractGlobalRowView: Global indices not long long";
528}
529#endif
530
531//==============================================================================
532
533template<typename int_type>
534void BlockCrsMatrix::TExtractBlock(Epetra_CrsMatrix & BaseMatrix, const int_type Row, const int_type Col)
535{
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_;
540
541// const Epetra_CrsGraph & BaseGraph = BaseMatrix.Graph();
542 const Epetra_BlockMap & BaseMap = BaseMatrix.RowMatrixRowMap();
543 //const Epetra_BlockMap & BaseColMap = BaseMatrix.RowMatrixColMap();
544
545 // This routine extracts entries of a BaseMatrix from a big BlockCrsMatrix
546 // It performs the following operation on the global IDs row-by-row
547 // BaseMatrix.val[i][j] = this->val[i+rowOffset][j+ColOffset]
548
549 int MaxIndices = BaseMatrix.MaxNumEntries();
550 vector<int_type> Indices(MaxIndices);
551 vector<double> vals(MaxIndices);
552 int NumIndices;
553 int_type indx,icol;
554 double* BlkValues;
555 int *BlkIndices;
556 int BlkNumIndices;
557 int ierr=0;
558 (void) ierr; // Forestall compiler warning for unused variable.
559
560 for (int i=0; i<BaseMap.NumMyElements(); i++) {
561
562 // Get pointers to values and indices of whole block matrix row
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);
566
567 NumIndices = 0;
568 // Grab columns with global indices in correct range for this block
569 for( int l = 0; l < BlkNumIndices; ++l ) {
570 icol = (int_type) this->RowMatrixColMap().GID64(BlkIndices[l]);
571 indx = icol - ColOffset;
572 if (indx >= 0 && indx < COffset_) {
573 Indices[NumIndices] = indx;
574 vals[NumIndices] = BlkValues[l];
575 NumIndices++;
576 }
577 }
578
579 //Load this row into base matrix
580 BaseMatrix.ReplaceGlobalValues(BaseRow, NumIndices, &vals[0], &Indices[0]);
581 }
582}
583
584#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
585void BlockCrsMatrix::ExtractBlock(Epetra_CrsMatrix & BaseMatrix, const int Row, const int Col)
586{
587 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt() && BaseMatrix.RowMatrixRowMap().GlobalIndicesInt())
588 return TExtractBlock<int>(BaseMatrix, Row, Col);
589 else
590 throw "EpetraExt::BlockCrsMatrix::ExtractBlock: Global indices not int";
591}
592#endif
593
594#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
595void BlockCrsMatrix::ExtractBlock(Epetra_CrsMatrix & BaseMatrix, const long long Row, const long long Col)
596{
597 if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong() && BaseMatrix.RowMatrixRowMap().GlobalIndicesLongLong())
598 return TExtractBlock<long long>(BaseMatrix, Row, Col);
599 else
600 throw "EpetraExt::BlockCrsMatrix::ExtractBlock: Global indices not long long";
601}
602#endif
603
604} //namespace EpetraExt
605
EpetraExt_BlockCrsMatrix.h
EpetraExt_BlockUtility.h
Copy
Copy
EpetraExt::BlockCrsMatrix::SumIntoGlobalBlock
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...
Definition EpetraExt_BlockCrsMatrix.cpp:366
EpetraExt::BlockCrsMatrix::RowStencil_int_
std::vector< std::vector< int > > RowStencil_int_
Definition EpetraExt_BlockCrsMatrix.h:225
EpetraExt::BlockCrsMatrix::SumIntoBlock
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...
Definition EpetraExt_BlockCrsMatrix.cpp:304
EpetraExt::BlockCrsMatrix::BlockExtractGlobalRowView
void BlockExtractGlobalRowView(const int BaseRow, int &NumEntries, double *&Values, const int Row, const int Col)
Definition EpetraExt_BlockCrsMatrix.cpp:510
EpetraExt::BlockCrsMatrix::ExtractBlock
void ExtractBlock(Epetra_CrsMatrix &BaseMatrix, const int Row, const int Col)
Definition EpetraExt_BlockCrsMatrix.cpp:585
EpetraExt::BlockCrsMatrix::RowStencil_LL_
std::vector< std::vector< long long > > RowStencil_LL_
Definition EpetraExt_BlockCrsMatrix.h:230
EpetraExt::BlockCrsMatrix::BlockSumIntoGlobalValues
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...
Definition EpetraExt_BlockCrsMatrix.cpp:412
EpetraExt::BlockCrsMatrix::BlockReplaceGlobalValues
void BlockReplaceGlobalValues(const int BaseRow, int NumIndices, double *Values, const int *Indices, const int Row, const int Col)
Definition EpetraExt_BlockCrsMatrix.cpp:459
EpetraExt::BlockCrsMatrix::RowIndices_LL_
std::vector< long long > RowIndices_LL_
Definition EpetraExt_BlockCrsMatrix.h:232
EpetraExt::BlockCrsMatrix::LoadBlock
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 ...
Definition EpetraExt_BlockCrsMatrix.cpp:241
EpetraExt::BlockCrsMatrix::BlockCrsMatrix
BlockCrsMatrix(const Epetra_CrsGraph &BaseGraph, const std::vector< int > &RowStencil, int RowIndex, const Epetra_Comm &GlobalComm)
BlockCrsMatrix constuctor with one block row per processor.
EpetraExt::BlockUtility::GenerateRowStencil
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.
Definition EpetraExt_BlockUtility.cpp:498
EpetraExt::BlockUtility::GenerateBlockGraph
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.
EpetraExt::BlockUtility::CalculateOffset64
static long long CalculateOffset64(const Epetra_BlockMap &BaseMap)
Definition EpetraExt_BlockUtility.cpp:546
Epetra_BlockMap
Epetra_BlockMap::LID
int LID(int GID) const
Epetra_BlockMap::GlobalIndicesInt
bool GlobalIndicesInt() const
Epetra_BlockMap::NumMyElements
int NumMyElements() const
Epetra_BlockMap::GlobalIndicesLongLong
bool GlobalIndicesLongLong() const
Epetra_Comm
Epetra_CrsGraph
Epetra_CrsGraph::RowMap
const Epetra_BlockMap & RowMap() const
Epetra_CrsGraph::ColMap
const Epetra_BlockMap & ColMap() const
Epetra_CrsMatrix
Epetra_CrsMatrix::ExtractGlobalRowView
int ExtractGlobalRowView(int GlobalRow, int &NumEntries, double *&Values, int *&Indices) const
Epetra_CrsMatrix::SumIntoGlobalValues
virtual int SumIntoGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Epetra_CrsMatrix::ExtractMyRowView
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
Epetra_CrsMatrix::ReplaceGlobalValues
virtual int ReplaceGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Epetra_CrsMatrix::MaxNumEntries
int MaxNumEntries() const
Epetra_CrsMatrix::RowMatrixColMap
const Epetra_Map & RowMatrixColMap() const
Epetra_CrsMatrix::RowMatrixRowMap
const Epetra_Map & RowMatrixRowMap() const
Epetra_RowMatrix
Epetra_RowMatrix::RowMatrixColMap
virtual const Epetra_Map & RowMatrixColMap() const=0
Epetra_RowMatrix::RowMatrixRowMap
virtual const Epetra_Map & RowMatrixRowMap() const=0
Epetra_RowMatrix::MaxNumEntries
virtual int MaxNumEntries() const=0
Epetra_RowMatrix::ExtractMyRowCopy
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
EpetraExt
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.
Definition EpetraExt_BlockCrsMatrix.cpp:46
Generated by doxygen 1.9.8