Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
SparseMatrix.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: SparseMatrix.cpp
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: Generic 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 
48 
49 #include <boost/lexical_cast.hpp>
50 
51 using std::min;
52 using std::max;
53 
54 namespace Nektar
55 {
56 
57  template<typename SparseStorageType>
59  m_mulCallsCounter(0),
60  m_sparseStorage(sparseStoragePtr)
61  {
62  }
63 
64  template<typename SparseStorageType>
66  m_mulCallsCounter(src.m_mulCallsCounter),
67  m_sparseStorage(src.m_sparseStorage)
68  {
69  }
70 
71  template<typename SparseStorageType>
73  {
74  }
75 
76 
77  template<typename SparseStorageType>
79  {
80  return m_sparseStorage->GetRows();
81  }
82 
83  template<typename SparseStorageType>
85  {
86  return m_sparseStorage->GetColumns();
87  }
88 
89  template<typename SparseStorageType>
91  {
92  return m_sparseStorage->GetNumNonZeroEntries();
93  }
94 
95 
96  template<typename SparseStorageType>
98  {
99  return m_sparseStorage->GetFillInRatio();
100  }
101 
102 
103 
104  template<typename SparseStorageType>
105  typename boost::call_traits<typename SparseStorageType::DataType>::const_reference NekSparseMatrix<SparseStorageType>::operator()(unsigned int row, unsigned int column) const
106  {
107  return m_sparseStorage->GetValue(row, column);
108  }
109 
110  template<typename SparseStorageType>
111  typename SparseStorageType::const_iterator NekSparseMatrix<SparseStorageType>::begin() const
112  {
113  return m_sparseStorage->begin();
114  }
115 
116  template<typename SparseStorageType>
117  typename SparseStorageType::const_iterator NekSparseMatrix<SparseStorageType>::end() const
118  {
119  return m_sparseStorage->end();
120  }
121 
122 
123  template<typename SparseStorageType>
125  DataVectorType &out)
126  {
127  m_sparseStorage->Multiply(in,out);
128  m_mulCallsCounter++;
129  }
130 
131  template<typename SparseStorageType>
133  DataType* out)
134  {
135  m_sparseStorage->Multiply(in,out);
136  m_mulCallsCounter++;
137  }
138 
139  template<typename SparseStorageType>
141  {
142  return m_sparseStorage->GetMemoryUsage(
143  m_sparseStorage->GetNumNonZeroEntries(),
144  m_sparseStorage->GetRows()
145  ) +
146  sizeof(SparseStorageSharedPtr) +
147  sizeof(unsigned long); // mulCallsCounter
148  }
149 
150  template<typename SparseStorageType>
152  {
153  return m_mulCallsCounter;
154  }
155 
156  template<typename SparseStorageType>
157  const typename SparseStorageType::DataType NekSparseMatrix<SparseStorageType>::GetAvgRowDensity() const
158  {
159  return (DataType) m_sparseStorage->GetNumNonZeroEntries() /
160  (DataType) m_sparseStorage->GetRows();
161  }
162 
163  template<typename SparseStorageType>
165  {
166  int bandwidth = 0;
167 
168  typename SparseStorageType::const_iterator entry = m_sparseStorage->begin();
169 
170  for (; entry != m_sparseStorage->end(); ++entry)
171  {
172  bandwidth = (std::max)(bandwidth, 2*std::abs((int)(entry->first.first - entry->first.second+1)));
173  }
174  return bandwidth;
175  }
176 
177  template<typename SparseStorageType>
179  {
180  COOMatTypeSharedPtr coo (new COOMatType());
181 
182  typename SparseStorageType::const_iterator entry = m_sparseStorage->begin();
183  for (; entry != m_sparseStorage->end(); entry++)
184  {
185  (*coo)[std::make_pair(entry->first.first, entry->first.second) ] = entry->second;
186  }
187  return coo;
188  }
189 
190  // Generates a matrix each element of which is a number of
191  // nonzero entries to block-matrix associated with *this object.
192  // E.g. for unit 6x6 matrix and block size = 2 it will generate
193  // 3x3 matrix with elements 2 along the diagonal.
194  //
195  // The output can be visualised with Python's matplotlib's spy().
196  //
197  template<typename SparseStorageType>
199  {
200  const int matRows = m_sparseStorage->GetRows();
201  const int gridRows = matRows / blockSize + (matRows % blockSize > 0);
202  const int gridCols = gridRows;
203 
204  std::vector< std::vector<int> > grid (gridRows);
205  for (int row = 0; row < gridRows; row++)
206  {
207  grid[row].resize(gridCols,0);
208  }
209 
210  typename SparseStorageType::const_iterator entry = m_sparseStorage->begin();
211  for (; entry != m_sparseStorage->end(); entry++)
212  {
213  const IndexType row = entry->first.first;
214  const IndexType col = entry->first.second;
215  const int gridRow = row / blockSize;
216  const int gridCol = col / blockSize;
217  grid[gridRow][gridCol]++;
218  }
219 
220  for (int row = 0; row < gridRows; row++)
221  {
222  for (int col = 0; col < gridCols; col++)
223  {
224  out << grid[row][col] << " ";
225  }
226  out << std::endl;
227  }
228  }
229 
230  /// Complementary routine to the previous. It generates
231  /// exact non-zero pattern of a given block matrix entry
232  /// to *this object. E.g., for a 6x6 matrix defined as
233  /// A(i,i):=i, A(i,j):=0, block size = 2 and
234  /// blk_row=blk_col=2 it produces 2x2 matrix with 5 and 6
235  /// along the diagonal.
236  template<typename SparseStorageType>
238  const IndexType blk_row, const IndexType blk_col, IndexType blockSize)
239  {
240  blockSize = (std::min)(blockSize, m_sparseStorage->GetRows());
241  std::vector< std::vector<int> > grid (blockSize);
242  for (int row = 0; row < blockSize; row++)
243  {
244  grid[row].resize(blockSize,0);
245  }
246 
247  typename SparseStorageType::const_iterator entry = m_sparseStorage->begin();
248  for (; entry != m_sparseStorage->end(); entry++)
249  {
250  const IndexType row = entry->first.first;
251  const IndexType col = entry->first.second;
252 
253  if (blk_row != row / blockSize ) continue;
254  if (blk_col != col / blockSize ) continue;
255  grid[row % blockSize][col % blockSize]++;
256  }
257 
258  for (int row = 0; row < blockSize; row++)
259  {
260  for (int col = 0; col < blockSize; col++)
261  {
262  out << grid[row][col] << " ";
263  }
264  out << std::endl;
265  }
266  }
267 
268 
269 
270  // explicit instantiation
272 
273 
274 } // namespace
275 
boost::shared_ptr< SparseStorageType > SparseStorageSharedPtr
COOMatTypeSharedPtr GetCooStorage()
std::map< CoordType, NekDouble > COOMatType
boost::shared_ptr< COOMatType > COOMatTypeSharedPtr
boost::call_traits< DataType >::const_reference operator()(const IndexType row, const IndexType column) const
SparseStorageType::const_iterator begin() const
const size_t GetMemoryFootprint() const
SparseStorageType::DataType DataType
void Multiply(const DataVectorType &in, DataVectorType &out)
const IndexType GetColumns() const
const IndexType GetBandwidth()
const IndexType GetNumNonZeroEntries() const
const unsigned long GetMulCallsCounter() const
void writeBlockSparsityPatternTo(std::ostream &out, 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...
const IndexType GetRows() const
unsigned int IndexType
SparseStorageType::const_iterator end() const
const DataType GetAvgRowDensity() const
NekSparseMatrix(const SparseStorageSharedPtr &sparseStoragePtr)
void writeSparsityPatternTo(std::ostream &out, IndexType blockSize=64)
const DataType GetFillInRatio() const