EpetraExt  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
46 namespace EpetraExt {
47 
48 using std::vector;
49 
50 //==============================================================================
51 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
52 BlockCrsMatrix::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
68 BlockCrsMatrix::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
85 BlockCrsMatrix::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
100 BlockCrsMatrix::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 //==============================================================================
116 BlockCrsMatrix::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
148 BlockCrsMatrix::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
164 BlockCrsMatrix::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 //==============================================================================
180 BlockCrsMatrix::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 //==============================================================================
197 BlockCrsMatrix::~BlockCrsMatrix()
198 {
199 }
200 
201 //==============================================================================
202 template<typename int_type>
203 void 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
241 void 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
251 void 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 //==============================================================================
261 template<typename int_type>
262 void 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
304 void 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
314 void 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 //==============================================================================
324 template<typename int_type>
325 void 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
366 void 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
376 void 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 //==============================================================================
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)
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
412 void 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
423 void 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 //==============================================================================
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)
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
459 void 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
470 void 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 //==============================================================================
481 template<typename int_type>
482 void 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
510 void 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
521 void 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 
533 template<typename int_type>
534 void 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
585 void 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
595 void 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::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
Epetra_CrsMatrix::ReplaceGlobalValues
virtual int ReplaceGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Epetra_RowMatrix::RowMatrixRowMap
virtual const Epetra_Map & RowMatrixRowMap() const =0
Epetra_CrsMatrix::SumIntoGlobalValues
virtual int SumIntoGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Copy
Epetra_BlockMap::GlobalIndicesLongLong
bool GlobalIndicesLongLong() const
EpetraExt_BlockUtility.h
EpetraExt::BlockCrsMatrix::RowStencil_int_
std::vector< std::vector< int > > RowStencil_int_
Definition: EpetraExt_BlockCrsMatrix.h:219
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
Epetra_CrsMatrix::RowMatrixRowMap
const Epetra_Map & RowMatrixRowMap() const
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:226
Epetra_CrsGraph::ColMap
const Epetra_BlockMap & ColMap() const
EpetraExt::BlockCrsMatrix::ExtractBlock
void ExtractBlock(Epetra_CrsMatrix &BaseMatrix, const int Row, const int Col)
Definition: EpetraExt_BlockCrsMatrix.cpp:585
Epetra_BlockMap::GlobalIndicesInt
bool GlobalIndicesInt() const
EpetraExt::BlockUtility::CalculateOffset64
static long long CalculateOffset64(const Epetra_BlockMap &BaseMap)
Definition: EpetraExt_BlockUtility.cpp:546
Epetra_CrsMatrix::ExtractMyRowView
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
EpetraExt::BlockCrsMatrix::BlockExtractGlobalRowView
void BlockExtractGlobalRowView(const int BaseRow, int &NumEntries, double *&Values, const int Row, const int Col)
Definition: EpetraExt_BlockCrsMatrix.cpp:510
Epetra_BlockMap::NumMyElements
int NumMyElements() const
Epetra_RowMatrix::MaxNumEntries
virtual int MaxNumEntries() const =0
EpetraExt::BlockCrsMatrix::RowStencil_LL_
std::vector< std::vector< long long > > RowStencil_LL_
Definition: EpetraExt_BlockCrsMatrix.h:224
Epetra_Comm
Epetra_CrsMatrix::ExtractGlobalRowView
int ExtractGlobalRowView(int GlobalRow, int &NumEntries, double *&Values, int *&Indices) const
Epetra_CrsMatrix::MaxNumEntries
int MaxNumEntries() const
EpetraExt_BlockCrsMatrix.h
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
Epetra_BlockMap
Epetra_CrsGraph::RowMap
const Epetra_BlockMap & RowMap() const
Epetra_BlockMap::LID
int LID(int GID) const
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.
Epetra_CrsMatrix::RowMatrixColMap
const Epetra_Map & RowMatrixColMap() const
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
Epetra_CrsMatrix
Epetra_RowMatrix::ExtractMyRowCopy
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const =0
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
Epetra_RowMatrix::RowMatrixColMap
virtual const Epetra_Map & RowMatrixColMap() const =0
Epetra_RowMatrix
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.
Epetra_CrsGraph
Generated by   doxygen 1.8.5