Nektar++
SparseDiagBlkMatrix.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: SparseDiagBlkMatrix.hpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Diagonal block sparse matrix class templated by underlying
32 // sparse
33 // storage format
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
37 #include <algorithm>
38 #include <fstream>
39 #include <map>
40 #include <utility>
41 #include <vector>
42 
48 
49 namespace Nektar
50 {
51 
52 template <typename SparseStorageType>
54  const SparseStorageSharedPtrVector &sparseStoragePtrVector)
55  : m_rows(0), m_cols(0), m_rowoffset(sparseStoragePtrVector.size() + 1, 0.0),
56  m_mulCallsCounter(0),
57  m_submatrix(sparseStoragePtrVector.size(), sparseStoragePtrVector)
58 {
59  for (int i = 0; i < sparseStoragePtrVector.size(); i++)
60  {
61  const IndexType rows = sparseStoragePtrVector[i]->GetRows();
62  m_rows += rows;
63  m_cols += sparseStoragePtrVector[i]->GetColumns();
64  m_rowoffset[i + 1] = m_rowoffset[i] + rows;
65  }
66 }
67 
68 template <typename SparseStorageType>
70  const NekSparseDiagBlkMatrix &src)
71  : m_rows(src.m_rows), m_cols(src.m_cols), m_rowoffset(src.m_rowoffset),
72  m_mulCallsCounter(src.m_mulCallsCounter), m_submatrix(src.m_submatrix)
73 {
74 }
75 
76 template <typename SparseStorageType>
78 {
79 }
80 
81 template <typename SparseStorageType>
83 {
84  return m_rows;
85 }
86 
87 template <typename SparseStorageType>
89 {
90  return m_cols;
91 }
92 
93 /// number of rows at i-th submatrix
94 template <typename SparseStorageType>
96 {
97  return m_submatrix[i]->GetRows();
98 }
99 
100 /// number of columns at i-th submatrix
101 template <typename SparseStorageType>
103 {
104  return m_submatrix[i]->GetColumns();
105 }
106 
107 template <typename SparseStorageType>
109  const
110 {
111  return m_submatrix.size();
112 }
113 
114 template <typename SparseStorageType>
116 {
117  IndexType nnz = 0;
118  for (int i = 0; i < m_submatrix.size(); i++)
119  {
120  nnz += m_submatrix[i]->GetNumNonZeroEntries();
121  }
122  return nnz;
123 }
124 
125 // nnz of i-th CSR matrix
126 template <typename SparseStorageType>
128  int i) const
129 {
130  return m_submatrix[i]->GetNumNonZeroEntries();
131 }
132 
133 template <typename SparseStorageType>
135  SparseStorageType>::GetFillInRatio(int i) const
136 {
137  return m_submatrix[i]->GetFillInRatio();
138 }
139 
140 template <typename SparseStorageType>
142  SparseStorageType>::GetFillInRatio() const
143 {
144  IndexType stored = 0;
145  IndexType nnz = 0;
146  for (int i = 0; i < m_submatrix.size(); i++)
147  {
148  stored += m_submatrix[i]->GetNumStoredDoubles();
149  nnz += m_submatrix[i]->GetNumNonZeroEntries();
150  }
151  return (DataType)stored / (DataType)nnz;
152 }
153 
154 template <typename SparseStorageType>
155 typename boost::call_traits<
156  typename SparseStorageType::DataType>::const_reference
158  IndexType block, IndexType loc_row, IndexType loc_column) const
159 {
160  return m_submatrix[block]->GetValue(loc_row, loc_column);
161 }
162 
163 template <typename SparseStorageType>
164 typename boost::call_traits<
165  typename SparseStorageType::DataType>::const_reference
167  IndexType glob_row, IndexType glob_column) const
168 {
169  IndexType i = 0;
170  static DataType defaultReturnValue = 0;
171 
172  signed int local_row = glob_row;
173  signed int local_col = glob_column;
174 
175  while ((local_row >= m_submatrix[i]->GetRows()) && (local_col >= 0))
176  {
177  local_row -= m_submatrix[i]->GetRows();
178  local_col -= m_submatrix[i]->GetColumns();
179  i++;
180  }
181  /// \todo double check, might be a bug when local_col > local column
182  if (local_col < 0)
183  return defaultReturnValue;
184 
185  return m_submatrix[i]->GetValue(local_row, local_col);
186 }
187 
188 template <typename SparseStorageType>
190  const DataVectorType &in, DataVectorType &out)
191 {
192  for (int i = 0; i < m_submatrix.size(); ++i)
193  {
194  m_submatrix[i]->Multiply(&in[m_rowoffset[i]], &out[m_rowoffset[i]]);
195  }
196  m_mulCallsCounter++;
197 }
198 
199 template <typename SparseStorageType>
201  DataType *out)
202 {
203  for (int i = 0; i < m_submatrix.size(); ++i)
204  {
205  m_submatrix[i]->Multiply(&in[m_rowoffset[i]], &out[m_rowoffset[i]]);
206  }
207  m_mulCallsCounter++;
208 }
209 
210 template <typename SparseStorageType>
212  const IndexType blockNum, DataType *in, DataType *out)
213 {
214  m_submatrix[blockNum]->Multiply(in + m_rowoffset[blockNum],
215  out + m_rowoffset[blockNum]);
216  // m_mulCallsCounter++;
217 }
218 
219 template <typename SparseStorageType>
221 {
222  size_t bytes = sizeof(IndexType) * 2 + // sizes
223  sizeof(unsigned long) + // mulCallsCounter
224  sizeof(IndexVector) +
225  sizeof(IndexType) * m_rowoffset.capacity() +
227  sizeof(SparseStorageSharedPtr) * m_submatrix.capacity();
228  for (int i = 0; i < m_submatrix.size(); i++)
229  {
230  bytes += m_submatrix[i]->GetMemoryUsage();
231  }
232  return bytes;
233 }
234 
235 template <typename SparseStorageType>
237  IndexType i) const
238 {
239  return m_submatrix[i]->GetMemoryUsage();
240 }
241 
242 template <typename SparseStorageType>
244  const
245 {
246  return m_mulCallsCounter;
247 }
248 
249 template <typename SparseStorageType>
250 typename SparseStorageType::DataType NekSparseDiagBlkMatrix<
251  SparseStorageType>::GetAvgRowDensity()
252 {
253  DataType avgRowDensity = 0.0;
254  for (int i = 0; i < m_submatrix.size(); i++)
255  {
256  avgRowDensity += (DataType)m_submatrix[i]->GetNumNonZeroEntries() /
257  (DataType)m_submatrix[i]->GetRows();
258  }
259  return avgRowDensity / m_submatrix.size();
260 }
261 
262 template <typename SparseStorageType>
263 typename SparseStorageType::DataType NekSparseDiagBlkMatrix<
264  SparseStorageType>::GetAvgRowDensity(IndexType i) const
265 {
266  return (DataType)m_submatrix[i]->GetNumNonZeroEntries() /
267  (DataType)m_submatrix[i]->GetRows();
268 }
269 
270 template <typename SparseStorageType>
272 {
273  IndexType bandwidth = 0;
274  for (int i = 0; i < m_submatrix.size(); i++)
275  {
276  typename SparseStorageType::const_iterator entry =
277  m_submatrix[i]->begin();
278  for (; entry != m_submatrix[i]->end(); ++entry)
279  {
280  bandwidth =
281  (std::max)(static_cast<int>(bandwidth),
282  2 * abs(static_cast<int>(entry->first.first -
283  entry->first.second) +
284  1));
285  }
286  }
287  return bandwidth;
288 }
289 
290 template <typename SparseStorageType>
292 {
293  IndexType bandwidth = 0;
294  typename SparseStorageType::const_iterator entry = m_submatrix[i]->begin();
295  for (; entry != m_submatrix[i]->end(); ++entry)
296  {
297  bandwidth = (std::max)(
298  static_cast<int>(bandwidth),
299  2 * abs(static_cast<int>(entry->first.first - entry->first.second) +
300  1));
301  }
302  return bandwidth;
303 }
304 
305 template <typename SparseStorageType>
307  IndexType i)
308 {
309  COOMatTypeSharedPtr coo(new COOMatType());
310  typename SparseStorageType::const_iterator entry = m_submatrix[i]->begin();
311  for (; entry != m_submatrix[i]->end(); entry++)
312  {
313  IndexType loc_row = entry->first.first;
314  IndexType loc_col = entry->first.second;
315  (*coo)[std::make_pair(loc_row, loc_col)] = entry->second;
316  }
317  return coo;
318 }
319 
320 template <typename SparseStorageType>
322 {
323  COOMatTypeSharedPtr coo(new COOMatType());
324  IndexType row_offset = 0;
325  IndexType col_offset = 0;
326  for (IndexType i = 0; i < m_submatrix.size(); i++)
327  {
328  typename SparseStorageType::const_iterator entry =
329  m_submatrix[i]->begin();
330  for (; entry != m_submatrix[i]->end(); entry++)
331  {
332  IndexType loc_row = entry->first.first;
333  IndexType loc_col = entry->first.second;
334  (*coo)[std::make_pair(loc_row + row_offset, loc_col + col_offset)] =
335  entry->second;
336  }
337  row_offset += m_submatrix[i]->GetRows();
338  col_offset += m_submatrix[i]->GetColumns();
339  }
340  return coo;
341 }
342 
343 // Generates a matrix each element of which is a number of
344 // nonzero entries to block-matrix associated with *this object.
345 // E.g. for unit 6x6 matrix and block size = 2 it will generate
346 // 3x3 matrix with elements 2 along the diagonal.
347 //
348 // The output can be visualised with Python's matplotlib's spy().
349 //
350 template <typename SparseStorageType>
352  std::ostream &out, const IndexType subMatrixIdx, const IndexType blockSize)
353 {
354  const int matRows = m_submatrix[subMatrixIdx]->GetRows();
355  const int gridRows = matRows / blockSize + (matRows % blockSize > 0);
356  const int gridCols = gridRows;
357 
358  std::vector<std::vector<int>> grid(gridRows);
359  for (int row = 0; row < gridRows; row++)
360  {
361  grid[row].resize(gridCols, 0.0);
362  }
363 
364  typename SparseStorageType::const_iterator entry =
365  m_submatrix[subMatrixIdx]->begin();
366  typename SparseStorageType::const_iterator stop =
367  m_submatrix[subMatrixIdx]->end();
368  for (; entry != stop; ++entry)
369  {
370  const IndexType row = entry->first.first;
371  const IndexType col = entry->first.second;
372  const int gridRow = row / blockSize;
373  const int gridCol = col / blockSize;
374  grid[gridRow][gridCol]++;
375  }
376 
377  for (int row = 0; row < gridRows; row++)
378  {
379  for (int col = 0; col < gridCols; col++)
380  {
381  out << grid[row][col] << " ";
382  }
383  out << std::endl;
384  }
385 }
386 
387 template <typename SparseStorageType>
389  std::ostream &out, const IndexType blockSize)
390 {
391  const int matRows = GetRows();
392  const int gridRows = matRows / blockSize + (matRows % blockSize > 0);
393  const int gridCols = gridRows;
394 
395  std::vector<std::vector<int>> grid(gridRows);
396  for (int row = 0; row < gridRows; row++)
397  {
398  grid[row].resize(gridCols, 0.0);
399  }
400 
401  IndexType row_offset = 0;
402  IndexType col_offset = 0;
403  for (int i = 0; i < m_submatrix.size(); i++)
404  {
405  typename SparseStorageType::const_iterator entry =
406  m_submatrix[i]->begin();
407  typename SparseStorageType::const_iterator stop = m_submatrix[i]->end();
408  for (; entry != stop; ++entry)
409  {
410  const IndexType row = entry->first.first + row_offset;
411  const IndexType col = entry->first.second + col_offset;
412  const int gridRow = row / blockSize;
413  const int gridCol = col / blockSize;
414  grid[gridRow][gridCol]++;
415  }
416  row_offset += m_submatrix[i]->GetRows();
417  col_offset += m_submatrix[i]->GetColumns();
418  }
419 
420  for (int row = 0; row < gridRows; row++)
421  {
422  for (int col = 0; col < gridCols; col++)
423  {
424  out << grid[row][col] << " ";
425  }
426  out << std::endl;
427  }
428 }
429 
430 /// Complementary routine to the previous. It generates
431 /// exact non-zero pattern of a given block matrix entry
432 /// to *this object. E.g., for a 6x6 matrix defined as
433 /// A(i,i):=i, A(i,j):=0, block size = 2 and
434 /// blk_row=blk_col=2 it produces 2x2 matrix with 5 and 6
435 /// along the diagonal.
436 template <typename SparseStorageType>
438  writeSubmatrixBlockSparsityPatternTo(std::ostream &out,
439  const IndexType subMatrixIdx,
440  const IndexType blk_row,
441  const IndexType blk_col,
442  IndexType blockSize)
443 {
444  blockSize = (std::min)(blockSize, m_submatrix[subMatrixIdx]->GetRows());
445  std::vector<std::vector<int>> grid(blockSize);
446  for (int row = 0; row < blockSize; row++)
447  {
448  grid[row].resize(blockSize, 0.0);
449  }
450 
451  typename SparseStorageType::const_iterator entry =
452  m_submatrix[subMatrixIdx]->begin();
453  typename SparseStorageType::const_iterator stop =
454  m_submatrix[subMatrixIdx]->end();
455  for (; entry != stop; ++entry)
456  {
457  const IndexType row = entry->first.first;
458  const IndexType col = entry->first.second;
459 
460  if (blk_row != row / blockSize)
461  continue;
462  if (blk_col != col / blockSize)
463  continue;
464  grid[row % blockSize][col % blockSize]++;
465  }
466 
467  for (int row = 0; row < blockSize; row++)
468  {
469  for (int col = 0; col < blockSize; col++)
470  {
471  out << grid[row][col] << " ";
472  }
473  out << std::endl;
474  }
475 }
476 
477 template <typename SparseStorageType>
479  std::ostream &out, const IndexType blk_row, const IndexType blk_col,
480  IndexType blockSize)
481 {
482  blockSize = (std::min)(blockSize, GetRows());
483  std::vector<std::vector<int>> grid(blockSize);
484  for (int row = 0; row < blockSize; row++)
485  {
486  grid[row].resize(blockSize, 0.0);
487  }
488 
489  IndexType row_offset = 0;
490  IndexType col_offset = 0;
491  for (int i = 0; i < m_submatrix.size(); i++)
492  {
493  typename SparseStorageType::const_iterator entry =
494  m_submatrix[i]->begin();
495  typename SparseStorageType::const_iterator stop = m_submatrix[i]->end();
496  for (; entry != stop; ++entry)
497  {
498  const IndexType row = entry->first.first + row_offset;
499  const IndexType col = entry->first.second + col_offset;
500 
501  if (blk_row != row / blockSize)
502  continue;
503  if (blk_col != col / blockSize)
504  continue;
505  grid[row % blockSize][col % blockSize]++;
506  }
507  row_offset += m_submatrix[i]->GetRows();
508  col_offset += m_submatrix[i]->GetColumns();
509  }
510 
511  for (int row = 0; row < blockSize; row++)
512  {
513  for (int col = 0; col < blockSize; col++)
514  {
515  out << grid[row][col] << " ";
516  }
517  out << std::endl;
518  }
519 }
520 
521 // explicit instantiation
523 
524 } // namespace Nektar
SparseStorageType::DataType DataType
boost::call_traits< DataType >::const_reference operator()(const IndexType row, const IndexType column) const
std::shared_ptr< SparseStorageType > SparseStorageSharedPtr
void MultiplySubMatrix(const IndexType blockNum, DataType *in, DataType *out)
void writeSubmatrixBlockSparsityPatternTo(std::ostream &out, const IndexType subMatrixIdx, const IndexType blk_row=0, const IndexType blk_col=0, IndexType blockSize=64)
Complementary routine to the previous. It generates exact non-zero pattern of a given block matrix en...
unsigned long GetMulCallsCounter() const
void writeSubmatrixSparsityPatternTo(std::ostream &out, const IndexType subMatrixIdx, IndexType blockSize=64)
void writeSparsityPatternTo(std::ostream &out, IndexType blockSize=64)
void Multiply(const DataVectorType &in, DataVectorType &out)
void writeBlockSparsityPatternTo(std::ostream &out, const IndexType blk_row=0, const IndexType blk_col=0, IndexType blockSize=64)
NekSparseDiagBlkMatrix(const SparseStorageSharedPtrVector &sparseStoragePtrVector)
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
Array< OneD, IndexType > IndexVector
std::shared_ptr< COOMatType > COOMatTypeSharedPtr
unsigned int IndexType
std::map< CoordType, NekDouble > COOMatType
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:295