matrix.hpp
Go to the documentation of this file.
1 /*!
2  * \brief Dense Matrix class
3  *
4  * \author O. Krause
5  * \date 2013
6  *
7  *
8  * \par Copyright 1995-2015 Shark Development Team
9  *
10  * <BR><HR>
11  * This file is part of Shark.
12  * <http://image.diku.dk/shark/>
13  *
14  * Shark is free software: you can redistribute it and/or modify
15  * it under the terms of the GNU Lesser General Public License as published
16  * by the Free Software Foundation, either version 3 of the License, or
17  * (at your option) any later version.
18  *
19  * Shark is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  * GNU Lesser General Public License for more details.
23  *
24  * You should have received a copy of the GNU Lesser General Public License
25  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
26  *
27  */
28 #ifndef SHARK_LINALG_BLAS_MATRIX_HPP
29 #define SHARK_LINALG_BLAS_MATRIX_HPP
30 
31 #include "matrix_proxy.hpp"
32 #include "vector_proxy.hpp"
34 
35 #include <boost/serialization/collection_size_type.hpp>
36 #include <boost/serialization/array.hpp>
37 #include <boost/serialization/nvp.hpp>
38 #include <boost/serialization/vector.hpp>
39 
40 namespace shark {
41 namespace blas {
42 
43 /** \brief A dense matrix of values of type \c T.
44  *
45  * For a \f$(m \times n)\f$-dimensional matrix and \f$ 0 \leq i < m, 0 \leq j < n\f$, every element \f$ m_{i,j} \f$ is mapped to
46  * the \f$(i.n + j)\f$-th element of the container for row major orientation or the \f$ (i + j.m) \f$-th element of
47  * the container for column major orientation. In a dense matrix all elements are represented in memory in a
48  * contiguous chunk of memory by definition.
49  *
50  * Orientation can also be specified, otherwise a \c row_major is used.
51  *
52  * \tparam T the type of object stored in the matrix (like double, float, complex, etc...)
53  * \tparam L the storage organization. It can be either \c row_major or \c column_major. Default is \c row_major
54  */
55 template<class T, class L=row_major>
56 class matrix:public matrix_container<matrix<T, L> > {
57  typedef matrix<T, L> self_type;
58  typedef std::vector<T> array_type;
59 public:
60  typedef typename array_type::size_type size_type;
61  typedef typename array_type::difference_type difference_type;
62  typedef typename array_type::value_type value_type;
63  typedef value_type scalar_type;
64  typedef typename array_type::const_reference const_reference;
65  typedef typename array_type::reference reference;
66  typedef const T* const_pointer;
67  typedef T* pointer;
68 
69  typedef std::size_t index_type;
70  typedef index_type const* const_index_pointer;
71  typedef index_type index_pointer;
72 
73  typedef const matrix_reference<const self_type> const_closure_type;
74  typedef matrix_reference<self_type> closure_type;
75  typedef dense_tag storage_category;
76  typedef elementwise_tag evaluation_category;
77  typedef L orientation;
78 
79  // Construction and destruction
80 
81  /// Default dense matrix constructor. Make a dense matrix of size (0,0)
82  matrix():m_size1(0), m_size2(0){}
83 
84  /** Dense matrix constructor with defined size
85  * \param size1 number of rows
86  * \param size2 number of columns
87  */
88  matrix(size_type size1, size_type size2)
89  :m_size1(size1), m_size2(size2), m_data(size1 * size2) {}
90 
91  /** Dense matrix constructor with defined size a initial value for all the matrix elements
92  * \param size1 number of rows
93  * \param size2 number of columns
94  * \param init initial value assigned to all elements
95  */
96  matrix(size_type size1, size_type size2, const value_type& init):
97  m_size1(size1), m_size2(size2), m_data(size1 * size2, init) {}
98 
99  /** Dense matrix constructor with defined size and an initial data array
100  * \param size1 number of rows
101  * \param size2 number of columns
102  * \param data array to copy into the matrix. Must have the same dimension as the matrix
103  */
104  matrix(size_type size1, size_type size2, const array_type& data):
105  m_size1(size1), m_size2(size2), m_data(data) {}
106 
107  /** Copy-constructor of a dense matrix
108  * \param m is a dense matrix
109  */
110  matrix(const matrix& m):
111  m_size1(m.m_size1), m_size2(m.m_size2), m_data(m.m_data) {}
112 
113  /** Copy-constructor of a dense matrix from a matrix expression
114  * \param e is a matrix expression
115  */
116  template<class E>
118  m_size1(e().size1()), m_size2(e().size2()), m_data(m_size1 * m_size2) {
119  assign(*this,e);
120  }
121 
122  // ---------
123  // Dense low level interface
124  // ---------
125 
126  ///\brief Returns the number of rows of the matrix.
127  size_type size1() const {
128  return m_size1;
129  }
130  ///\brief Returns the number of columns of the matrix.
131  size_type size2() const {
132  return m_size2;
133  }
134 
135  ///\brief Returns the stride in memory between two rows.
136  difference_type stride1()const{
137  return orientation::stride1(m_size1,m_size2);
138  }
139  ///\brief Returns the stride in memory between two columns.
140  difference_type stride2()const{
141  return orientation::stride2(m_size1,m_size2);
142  }
143 
144  ///\brief Returns the pointer to the beginning of the matrix storage
145  ///
146  /// Grants low-level access to the matrix internals. Element order depends on whether the matrix is row_major or column_major.
147  /// to access element (i,j) use storage()[i*stride1()+j*stride2()].
148  const_pointer storage()const{
149  return &m_data[0];
150  }
151 
152  ///\brief Returns the pointer to the beginning of the matrix storage
153  ///
154  /// Grants low-level access to the matrix internals. Element order depends on whether the matrix is row_major or column_major.
155  /// to access element (i,j) use storage()[i*stride1()+j*stride2()].
156  pointer storage(){
157  return &m_data[0];
158  }
159 
160  // ---------
161  // High level interface
162  // ---------
163 
164  // Resizing
165  /** Resize a matrix to new dimensions. If resizing is performed, the data is not preserved.
166  * \param size1 the new number of rows
167  * \param size2 the new number of colums
168  */
169  void resize(size_type size1, size_type size2) {
170  m_data.resize(size1* size2);
171  m_size1 = size1;
172  m_size2 = size2;
173  }
174 
175  void clear(){
176  std::fill(m_data.begin(), m_data.end(), value_type/*zero*/());
177  }
178 
179  // Element access
180  const_reference operator()(index_type i, index_type j) const {
181  return m_data [orientation::element(i, m_size1, j, m_size2)];
182  }
183  reference operator()(index_type i, index_type j) {
184  return m_data [orientation::element(i, m_size1, j, m_size2)];
185  }
186 
187  void set_element(size_type i, size_type j,value_type t){
188  m_data [orientation::element(i, m_size1, j, m_size2)] = t;
189  }
190 
191  // Assignment
192  /*! @note "pass by value" the key idea to enable move semantics */
194  swap(m);
195  return *this;
196  }
197  template<class C> // Container assignment without temporary
199  resize(m().size1(), m().size2());
200  assign(*this, m);
201  return *this;
202  }
203  template<class E>
205  self_type temporary(e);
206  swap(temporary);
207  return *this;
208  }
209 
210  // Swapping
211  void swap(matrix& m) {
212  std::swap(m_size1, m.m_size1);
213  std::swap(m_size2, m.m_size2);
214  m_data.swap(m.m_data);
215  }
216  friend void swap(matrix& m1, matrix& m2) {
217  m1.swap(m2);
218  }
219 
220  friend void swap_rows(matrix& a, index_type i, matrix& b, index_type j){
221  SIZE_CHECK(i < a.size1());
222  SIZE_CHECK(j < b.size1());
223  SIZE_CHECK(a.size2() == b.size2());
224  for(std::size_t k = 0; k != a.size2(); ++k){
225  std::swap(a(i,k),b(j,k));
226  }
227  }
228 
229  friend void swap_rows(matrix& a, index_type i, index_type j) {
230  if(i == j) return;
231  swap_rows(a,i,a,j);
232  }
233 
234 
235  friend void swap_columns(matrix& a, index_type i, matrix& b, index_type j){
236  SIZE_CHECK(i < a.size2());
237  SIZE_CHECK(j < b.size2());
238  SIZE_CHECK(a.size1() == b.size1());
239  for(std::size_t k = 0; k != a.size1(); ++k){
240  std::swap(a(k,i),b(k,j));
241  }
242  }
243 
244  friend void swap_columns(matrix& a, index_type i, index_type j) {
245  if(i == j) return;
246  swap_columns(a,i,a,j);
247  }
248 
249  //Iterators
250  typedef dense_storage_iterator<value_type> row_iterator;
251  typedef dense_storage_iterator<value_type> column_iterator;
252  typedef dense_storage_iterator<value_type const> const_row_iterator;
253  typedef dense_storage_iterator<value_type const> const_column_iterator;
254 
255  const_row_iterator row_begin(index_type i) const {
256  return const_row_iterator(&m_data[0] + i*stride1(),0,stride2());
257  }
258  const_row_iterator row_end(index_type i) const {
259  return const_row_iterator(&m_data[0] + i*stride1()+stride2()*size2(),size2(),stride2());
260  }
261  row_iterator row_begin(index_type i){
262  return row_iterator(&m_data[0] + i*stride1(),0,stride2());
263  }
264  row_iterator row_end(index_type i){
265  return row_iterator(&m_data[0] + i*stride1()+stride2()*size2(),size2(),stride2());
266  }
267 
268  const_row_iterator column_begin(std::size_t j) const {
269  return const_column_iterator(&m_data[0]+j*stride2(),0,stride1());
270  }
271  const_column_iterator column_end(std::size_t j) const {
272  return const_column_iterator(&m_data[0]+j*stride2()+ stride1()*size1(),size1(),stride1());
273  }
274  column_iterator column_begin(std::size_t j){
275  return column_iterator(&m_data[0]+j*stride2(),0,stride1());
276  }
277  column_iterator column_end(std::size_t j){
278  return column_iterator(&m_data[0]+j*stride2()+ stride1()*size1(),size1(),stride1());
279  }
280 
281  typedef typename blas::major_iterator<self_type>::type major_iterator;
282 
283  //sparse interface
284  major_iterator set_element(major_iterator pos, index_type index, value_type value) {
285  RANGE_CHECK(pos.index() == index);
286  *pos=value;
287  return pos;
288  }
289 
290  major_iterator clear_element(major_iterator elem) {
291  *elem = value_type();
292  return elem+1;
293  }
294 
295  major_iterator clear_range(major_iterator start, major_iterator end) {
296  std::fill(start,end,value_type());
297  return end;
298  }
299 
300  void reserve(size_type non_zeros) {}
301  void reserve_row(std::size_t, std::size_t){}
302  void reserve_column(std::size_t, std::size_t){}
303 
304  // Serialization
305  template<class Archive>
306  void serialize(Archive& ar, const unsigned int /* file_version */) {
307 
308  // we need to copy to a collection_size_type to get a portable
309  // and efficient boost::serialization
310  boost::serialization::collection_size_type s1(m_size1);
311  boost::serialization::collection_size_type s2(m_size2);
312 
313  // serialize the sizes
314  ar& boost::serialization::make_nvp("size1",s1)
315  & boost::serialization::make_nvp("size2",s2);
316 
317  // copy the values back if loading
318  if (Archive::is_loading::value) {
319  m_size1 = s1;
320  m_size2 = s2;
321  }
322  ar& boost::serialization::make_nvp("data",m_data);
323  }
324 
325 private:
326  size_type m_size1;
327  size_type m_size2;
328  array_type m_data;
329 };
330 template<class T, class L>
331 struct matrix_temporary_type<T,L,dense_random_access_iterator_tag>{
332  typedef matrix<T,L> type;
333 };
334 
335 template<class T>
336 struct matrix_temporary_type<T,unknown_orientation,dense_random_access_iterator_tag>{
337  typedef matrix<T,row_major> type;
338 };
339 
340 /** \brief An diagonal matrix with values stored inside a diagonal vector
341  *
342  * the matrix stores a Vector representing the diagonal.
343  */
344 template<class VectorType>
345 class diagonal_matrix: public matrix_container<diagonal_matrix< VectorType > > {
347 public:
348  typedef std::size_t size_type;
349  typedef std::ptrdiff_t difference_type;
350  typedef typename VectorType::value_type value_type;
351  typedef typename VectorType::scalar_type scalar_type;
352  typedef typename VectorType::const_reference const_reference;
353  typedef typename VectorType::reference reference;
354  typedef typename VectorType::pointer pointer;
355  typedef typename VectorType::const_pointer const_pointer;
356 
357  typedef std::size_t index_type;
358  typedef index_type const* const_index_pointer;
359  typedef index_type index_pointer;
360 
361  typedef const matrix_reference<const self_type> const_closure_type;
362  typedef matrix_reference<self_type> closure_type;
363  typedef sparse_tag storage_category;
364  typedef elementwise_tag evaluation_category;
365  typedef row_major orientation;
366 
367  // Construction and destruction
368  diagonal_matrix():m_zero(){}
369  diagonal_matrix(VectorType const& diagonal):m_zero(),m_diagonal(diagonal){}
370 
371  // Accessors
372  size_type size1() const {
373  return m_diagonal.size();
374  }
375  size_type size2() const {
376  return m_diagonal.size();
377  }
378 
379  // Element access
380  const_reference operator()(index_type i, index_type j) const {
381  if (i == j)
382  return m_diagonal(i);
383  else
384  return m_zero;
385  }
386 
387  void set_element(size_type i, size_type j,value_type t){
388  RANGE_CHECK(i == j);
389  m_diagonal(i) = t;
390  }
391 
392  // Assignment
394  m_diagonal = m.m_diagonal;
395  return *this;
396  }
397 
398  // Swapping
400  swap(m_diagonal,m.m_diagonal);
401  }
402  friend void swap(diagonal_matrix& m1, diagonal_matrix& m2) {
403  m1.swap(m2);
404  }
405 
406  //Iterators
407 
408  class const_row_iterator:public bidirectional_iterator_base<const_row_iterator, value_type> {
409  public:
410  typedef typename diagonal_matrix::value_type value_type;
411  typedef typename diagonal_matrix::difference_type difference_type;
413  typedef value_type const* pointer;
414 
415  // Construction and destruction
417  const_row_iterator(index_type index, value_type value, bool isEnd)
418  :m_index(index),m_value(value),m_isEnd(isEnd){}
419 
420  // Arithmetic
421  const_row_iterator& operator ++ () {
422  m_isEnd = true;
423  return *this;
424  }
425  const_row_iterator& operator -- () {
426  m_isEnd = false;
427  return *this;
428  }
429 
430  // Dereference
431  const_reference operator*() const {
432  return m_value;
433  }
434 
435  // Indices
436  index_type index() const{
437  return m_index;
438  }
439 
440  // Assignment
442  m_index = it.m_index;
443  return *this;
444  }
445 
446  // Comparison
447  bool operator == (const_row_iterator const& it) const {
448  RANGE_CHECK(m_index == it.m_index);
449  return m_isEnd == it.m_isEnd;
450  }
451 
452  private:
453  index_type m_index;
454  value_type m_value;
455  bool m_isEnd;
456  };
459  typedef const_column_iterator column_iterator;
460 
461 
462  const_row_iterator row_begin(index_type i) const {
463  return row_iterator(i, m_diagonal(i),false);
464  }
465  const_row_iterator row_end(index_type i) const {
466  return const_row_iterator(i, m_zero,true);
467  }
468  const_column_iterator column_begin(index_type i) const {
469  return column_iterator(i, m_diagonal(i),false);
470  }
471  const_column_iterator column_end(index_type i) const {
472  return const_column_iterator(i, m_zero,true);
473  }
474 
475 private:
476  value_type const m_zero;
477  VectorType m_diagonal;
478 };
479 
480 
481 template<class T, class Orientation>
482 struct const_expression<matrix<T,Orientation> >{
483  typedef matrix<T,Orientation> const type;
484 };
485 template<class T, class Orientation>
486 struct const_expression<matrix<T,Orientation> const>{
487  typedef matrix<T,Orientation> const type;
488 };
489 
490 }
491 }
492 
493 #endif