EpetraExt Package Browser (Single Doxygen Collection)  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 
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
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
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
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
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 //==============================================================================
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())
136  else
137 #endif
138 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
141  else
142 #endif
143  throw "EpetraExt::BlockCrsMatrix::BlockCrsMatrix: Error, Global indices unknown.";
144 }
145 
146 //==============================================================================
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 )
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
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 //==============================================================================
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 //==============================================================================
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 
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 ...
void TBlockReplaceGlobalValues(const int_type BaseRow, int NumIndices, double *Values, const int_type *Indices, const int_type Row, const int_type Col)
virtual const Epetra_Map & RowMatrixRowMap() const =0
bool GlobalIndicesLongLong() const
std::vector< std::vector< int > > RowStencil_int_
int ExtractGlobalRowView(int_type Row, int &NumEntries, double *&values, int_type *&Indices) const
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.
const Epetra_Map & RowMatrixRowMap() const
virtual int SumIntoGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
void BlockReplaceGlobalValues(const int BaseRow, int NumIndices, double *Values, const int *Indices, const int Row, const int Col)
std::vector< long long > RowIndices_LL_
const Epetra_BlockMap & ColMap() const
void ExtractBlock(Epetra_CrsMatrix &BaseMatrix, const int Row, const int Col)
bool GlobalIndicesInt() const
static long long CalculateOffset64(const Epetra_BlockMap &BaseMap)
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
void BlockExtractGlobalRowView(const int BaseRow, int &NumEntries, double *&Values, const int Row, const int Col)
int NumMyElements() const
virtual int MaxNumEntries() const =0
long long GID64(int LID) const
std::vector< std::vector< long long > > RowStencil_LL_
int MaxNumEntries() const
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...
void TBlockExtractGlobalRowView(const int_type BaseRow, int &NumEntries, double *&Values, const int_type Row, const int_type Col)
void TSumIntoGlobalBlock(double alpha, const Epetra_RowMatrix &BaseMatrix, const int_type Row, const int_type Col)
const Epetra_BlockMap & RowMap() const
int LID(int GID) const
BlockCrsMatrix(const Epetra_CrsGraph &BaseGraph, const std::vector< int > &RowStencil, int RowIndex, const Epetra_Comm &GlobalComm)
BlockCrsMatrix constuctor with one block row per processor.
const Epetra_Map & RowMatrixColMap() const
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...
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const =0
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...
virtual const Epetra_Map & RowMatrixColMap() const =0
void TExtractBlock(Epetra_CrsMatrix &BaseMatrix, const int_type Row, const int_type Col)
virtual int ReplaceGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
void TLoadBlock(const Epetra_RowMatrix &BaseMatrix, const int_type Row, const int_type Col)
void TSumIntoBlock(double alpha, const Epetra_RowMatrix &BaseMatrix, const int_type Row, const int_type Col)
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.
void TBlockSumIntoGlobalValues(const int_type BaseRow, int NumIndices, double *Values, const int_type *Indices, const int_type Row, const int_type Col)