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