49 #include <boost/lexical_cast.hpp>
55 template<
typename SparseStorageType>
58 sparseStoragePtrVector):
61 m_rowoffset(sparseStoragePtrVector.num_elements()+1, 0.0),
63 m_submatrix(sparseStoragePtrVector.num_elements(),sparseStoragePtrVector)
65 for (
int i = 0; i < sparseStoragePtrVector.num_elements(); i++)
67 const IndexType rows = sparseStoragePtrVector[i]->GetRows();
69 m_cols += sparseStoragePtrVector[i]->GetColumns();
74 template<
typename SparseStorageType>
78 m_rowoffset(src.m_rowoffset),
79 m_mulCallsCounter(src.m_mulCallsCounter),
80 m_submatrix(src.m_submatrix)
84 template<
typename SparseStorageType>
89 template<
typename SparseStorageType>
95 template<
typename SparseStorageType>
102 template<
typename SparseStorageType>
105 return m_submatrix[i]->GetRows();
109 template<
typename SparseStorageType>
112 return m_submatrix[i]->GetColumns();
115 template<
typename SparseStorageType>
118 return m_submatrix.num_elements();
121 template<
typename SparseStorageType>
125 for (
int i = 0; i < m_submatrix.num_elements(); i++)
127 nnz += m_submatrix[i]->GetNumNonZeroEntries();
133 template<
typename SparseStorageType>
136 return m_submatrix[i]->GetNumNonZeroEntries();
139 template<
typename SparseStorageType>
142 return m_submatrix[i]->GetFillInRatio();
145 template<
typename SparseStorageType>
150 for (
int i = 0; i < m_submatrix.num_elements(); i++)
152 stored += m_submatrix[i]->GetNumStoredDoubles();
153 nnz += m_submatrix[i]->GetNumNonZeroEntries();
159 template<
typename SparseStorageType>
162 return m_submatrix[block]->GetValue(loc_row, loc_column);
165 template<
typename SparseStorageType>
166 typename boost::call_traits<typename SparseStorageType::DataType>::const_reference
179 static DataType defaultReturnValue = 0;
181 signed int local_row = glob_row;
182 signed int local_col = glob_column;
184 while ((local_row >= m_submatrix[i]->GetRows()) &&
187 local_row -= m_submatrix[i]->
GetRows();
188 local_col -= m_submatrix[i]->GetColumns();
192 if (local_col < 0)
return defaultReturnValue;
194 return m_submatrix[i]->GetValue(local_row, local_col);
198 template<
typename SparseStorageType>
203 for (
int i = 0; i < m_submatrix.num_elements(); ++i)
205 m_submatrix[i]->Multiply(&in[m_rowoffset[i]], &out[m_rowoffset[i]]);
210 template<
typename SparseStorageType>
214 for (
int i = 0; i < m_submatrix.num_elements(); ++i)
216 m_submatrix[i]->Multiply(&in[m_rowoffset[i]], &out[m_rowoffset[i]]);
221 template<
typename SparseStorageType>
227 m_submatrix[blockNum]->Multiply(in + m_rowoffset[blockNum], out + m_rowoffset[blockNum]);
231 template<
typename SparseStorageType>
236 sizeof(
unsigned long) +
238 sizeof(
IndexType)*m_rowoffset.capacity() +
241 for (
int i = 0; i < m_submatrix.num_elements(); i++)
243 bytes += m_submatrix[i]->GetMemoryUsage(
244 m_submatrix[i]->GetNumNonZeroEntries(),
245 m_submatrix[i]->GetRows()
251 template<
typename SparseStorageType>
254 return m_submatrix[i]->GetMemoryUsage(
255 m_submatrix[i]->GetNumNonZeroEntries(),
256 m_submatrix[i]->GetRows()
260 template<
typename SparseStorageType>
263 return m_mulCallsCounter;
266 template<
typename SparseStorageType>
270 for (
int i = 0; i < m_submatrix.num_elements(); i++)
272 avgRowDensity += (
DataType) m_submatrix[i]->GetNumNonZeroEntries() /
273 (
DataType) m_submatrix[i]->GetRows();
275 return avgRowDensity / m_submatrix.num_elements();
278 template<
typename SparseStorageType>
281 return (
DataType) m_submatrix[i]->GetNumNonZeroEntries() /
282 (
DataType) m_submatrix[i]->GetRows();
285 template<
typename SparseStorageType>
289 for (
int i = 0; i < m_submatrix.num_elements(); i++)
291 typename SparseStorageType::const_iterator entry = m_submatrix[i]->begin();
292 for (; entry != m_submatrix[i]->end(); ++entry)
294 bandwidth = (std::max)(static_cast<int>(bandwidth),
295 2*abs( static_cast<int>(entry->first.first - entry->first.second)+1));
301 template<
typename SparseStorageType>
305 typename SparseStorageType::const_iterator entry = m_submatrix[i]->begin();
306 for (; entry != m_submatrix[i]->end(); ++entry)
308 bandwidth = (std::max)(static_cast<int>(bandwidth),
309 2*abs( static_cast<int>(entry->first.first - entry->first.second)+1));
314 template<
typename SparseStorageType>
318 typename SparseStorageType::const_iterator entry = m_submatrix[i]->begin();
319 for (; entry != m_submatrix[i]->end(); entry++)
323 (*coo)[std::make_pair(loc_row, loc_col) ] = entry->second;
328 template<
typename SparseStorageType>
334 for (
IndexType i = 0; i < m_submatrix.num_elements(); i++)
336 typename SparseStorageType::const_iterator entry = m_submatrix[i]->begin();
337 for (; entry != m_submatrix[i]->end(); entry++)
341 (*coo)[std::make_pair(loc_row + row_offset, loc_col + col_offset) ] = entry->second;
343 row_offset += m_submatrix[i]->GetRows();
344 col_offset += m_submatrix[i]->GetColumns();
356 template<
typename SparseStorageType>
362 const int matRows = m_submatrix[subMatrixIdx]->GetRows();
363 const int gridRows = matRows / blockSize + (matRows % blockSize > 0);
364 const int gridCols = gridRows;
366 std::vector< std::vector<int> > grid (gridRows);
367 for (
int row = 0; row < gridRows; row++)
369 grid[row].resize(gridCols,0.0);
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)
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]++;
385 for (
int row = 0; row < gridRows; row++)
387 for (
int col = 0; col < gridCols; col++)
389 out << grid[row][col] <<
" ";
395 template<
typename SparseStorageType>
400 const int matRows = GetRows();
401 const int gridRows = matRows / blockSize + (matRows % blockSize > 0);
402 const int gridCols = gridRows;
404 std::vector< std::vector<int> > grid (gridRows);
405 for (
int row = 0; row < gridRows; row++)
407 grid[row].resize(gridCols,0.0);
412 for (
int i = 0; i < m_submatrix.num_elements(); i++)
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)
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]++;
424 row_offset += m_submatrix[i]->GetRows();
425 col_offset += m_submatrix[i]->GetColumns();
428 for (
int row = 0; row < gridRows; row++)
430 for (
int col = 0; col < gridCols; col++)
432 out << grid[row][col] <<
" ";
444 template<
typename SparseStorageType>
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++)
456 grid[row].resize(blockSize,0.0);
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)
465 const IndexType row = entry->first.first;
466 const IndexType col = entry->first.second;
468 if (blk_row != row / blockSize )
continue;
469 if (blk_col != col / blockSize )
continue;
470 grid[row % blockSize][col % blockSize]++;
473 for (
int row = 0; row < blockSize; row++)
475 for (
int col = 0; col < blockSize; col++)
477 out << grid[row][col] <<
" ";
483 template<
typename SparseStorageType>
490 blockSize = (std::min)(blockSize, GetRows());
491 std::vector< std::vector<int> > grid (blockSize);
492 for (
int row = 0; row < blockSize; row++)
494 grid[row].resize(blockSize,0.0);
499 for (
int i = 0; i < m_submatrix.num_elements(); i++)
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)
507 const IndexType row = entry->first.first + row_offset;
508 const IndexType col = entry->first.second + col_offset;
510 if (blk_row != row / blockSize )
continue;
511 if (blk_col != col / blockSize )
continue;
512 grid[row % blockSize][col % blockSize]++;
514 row_offset += m_submatrix[i]->GetRows();
515 col_offset += m_submatrix[i]->GetColumns();
518 for (
int row = 0; row < blockSize; row++)
520 for (
int col = 0; col < blockSize; col++)
522 out << grid[row][col] <<
" ";
IndexType GetRows() const
std::map< CoordType, NekDouble > COOMatType
boost::shared_ptr< SparseStorageType > SparseStorageSharedPtr
boost::shared_ptr< COOMatType > COOMatTypeSharedPtr
COOMatTypeSharedPtr GetCooStorage()
size_t GetMemoryFootprint()
Array< OneD, IndexType > IndexVector
void MultiplySubMatrix(const IndexType blockNum, DataType *in, DataType *out)
DataType GetAvgRowDensity()
DataType GetFillInRatio() const
SparseStorageType::DataType DataType
IndexType GetNumberOfMatrixBlocks() const
unsigned long GetMulCallsCounter() const
boost::call_traits< DataType >::const_reference operator()(const IndexType row, const IndexType column) const
IndexType GetColumns() const
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...
void writeBlockSparsityPatternTo(std::ostream &out, const IndexType blk_row=0, const IndexType blk_col=0, IndexType blockSize=64)
void Multiply(const DataVectorType &in, DataVectorType &out)
void writeSparsityPatternTo(std::ostream &out, IndexType blockSize=64)
IndexType GetNumNonZeroEntries()
void writeSubmatrixSparsityPatternTo(std::ostream &out, const IndexType subMatrixIdx, IndexType blockSize=64)
~NekSparseDiagBlkMatrix()
NekSparseDiagBlkMatrix(const SparseStorageSharedPtrVector &sparseStoragePtrVector)