triangular_matrix.hpp
Go to the documentation of this file.
1 /*!
2  * \brief Implements a matrix with triangular storage layout
3  *
4  * \author O. Krause
5  * \date 2015
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_TRIANGULAR_MATRIX_HPP
29 #define SHARK_LINALG_BLAS_TRIANGULAR_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 template<class T, class Orientation, class TriangularType>
44 class triangular_matrix:public matrix_container<triangular_matrix<T,Orientation,TriangularType> > {
46  typedef std::vector<T> array_type;
47 public:
48  typedef typename array_type::size_type size_type;
49  typedef typename array_type::difference_type difference_type;
50  typedef typename array_type::value_type value_type;
51  typedef value_type scalar_type;
52  typedef value_type const_reference;
53  typedef value_type reference;
54  typedef const T* const_pointer;
55  typedef T* pointer;
56 
57  typedef std::size_t index_type;
58  typedef index_type const* const_index_pointer;
59  typedef index_type index_pointer;
60 
63  typedef packed_tag storage_category;
64  typedef elementwise_tag evaluation_category;
65  typedef packed<Orientation,TriangularType> orientation;
66 
67  // Construction and destruction
68 
69  /// Default triangular_matrix constructor. Make a dense matrix of size (0,0)
70  triangular_matrix():m_size(0){}
71 
72  /** Packed matrix constructor with defined size
73  * \param size number of rows and columns
74  */
75  triangular_matrix(size_type size):m_size(size),m_data(size * (size+1)/2) {}
76 
77  /** Packed matrix constructor with defined size and an initial value for all triangular matrix elements
78  * \param size number of rows and columns
79  * \param init initial value of the non-zero elements
80  */
81  triangular_matrix(size_type size, scalar_type init):m_size(size),m_data(size * (size+1)/2,init) {}
82 
83  /** Copy-constructor of a dense matrix
84  * \param m is a dense matrix
85  */
86  triangular_matrix(const triangular_matrix& m):m_size(m.m_size), m_data(m.m_data) {}
87 
88  /** Copy-constructor of a dense matrix from a matrix expression
89  * \param e is a matrix expression which has to be triangular
90  */
91  template<class E>
93  :m_size(e().size1()), m_data(m_size * (m_size+1)/2)
94  {
95  assign(*this, e);
96  }
97 
98  // ---------
99  // Low level interface
100  // ---------
101 
102  ///\brief Returns the number of rows of the matrix.
103  size_type size1() const {
104  return m_size;
105  }
106  ///\brief Returns the number of columns of the matrix.
107  size_type size2() const {
108  return m_size;
109  }
110 
111  ///\brief Returns the pointer to the beginning of the matrix storage
112  ///
113  /// Grants low-level access to the matrix internals. Element order depends on whether the matrix is
114  /// row_major or column_major and upper or lower triangular
115  const_pointer storage()const{
116  if(m_data.empty())
117  return 0;
118  return &m_data[0];
119  }
120 
121  ///\brief Returns the pointer to the beginning of the matrix storage
122  ///
123  /// Grants low-level access to the matrix internals. Element order depends on whether the matrix is row_major or column_major.
124  /// to access element (i,j) use storage()[i*stride1()+j*stride2()].
125  pointer storage(){
126  if(m_data.empty())
127  return 0;
128  return &m_data[0];
129  }
130 
131  ///\brief Number of nonzero-elements stores in the matrix.
132  size_type nnz() const {
133  return m_data.size();
134  }
135 
136  // ---------
137  // High level interface
138  // ---------
139 
140  // Resizing
141  /** Resize a matrix to new dimensions. If resizing is performed, the data is not preserved.
142  * \param size the new number of rows and columns
143  */
144  void resize(size_type size) {
145  m_data.resize(size*(size+1)/2);
146  m_size = size;
147  }
148 
149  void resize(size_type size1, size_type size2) {
150  SIZE_CHECK(size1 == size2);
151  resize(size1);
152  (void)size2;
153  }
154 
155  void clear(){
156  std::fill(m_data.begin(), m_data.end(), value_type/*zero*/());
157  }
158 
159  // Element access read only
160  const_reference operator()(index_type i, index_type j) const {
161  if(!orientation::non_zero(i,j))
162  return value_type();
163  return m_data [orientation::element(i,j,size1())];
164  }
165 
166  // separate write access
167  void set_element(std::size_t i,std::size_t j, value_type t){
168  SIZE_CHECK(orientation::non_zero(i,j));
169  m_data [orientation::element(i,j,size1())] = t;
170  }
171 
172  bool non_zero(std::size_t i,std::size_t j)const{
173  return orientation::non_zero(i,j);
174  }
175 
176  /*! @note "pass by value" the key idea to enable move semantics */
178  swap(m);
179  return *this;
180  }
181  template<class C> // Container assignment without temporary
183  SIZE_CHECK(m().size1()==m().size2());
184  resize(m().size1());
185  assign(*this, m);
186  return *this;
187  }
188  template<class E>
190  self_type temporary(e);
191  swap(temporary);
192  return *this;
193  }
194 
195  // Swapping
197  std::swap(m_size, m.m_size);
198  m_data.swap(m.m_data);
199  }
200  friend void swap(triangular_matrix& m1, triangular_matrix& m2) {
201  m1.swap(m2);
202  }
203 
204 
205  ///////////iterators
206  template<class TIter>
208  public random_access_iterator_base<
209  major1_iterator<TIter>,
210  typename boost::remove_const<T>::type,
211  packed_random_access_iterator_tag
212  >{
213  private:
214  difference_type offset(difference_type n)const{
215  difference_type k = m_index;
216  if(n >= 0){
217  return (n*(2*k+n+1))/2;
218  }else{
219  k+= n;
220  n*=-1;
221  return -(n*(2*k+n+1))/2;
222  }
223  }
224  public:
225  typedef typename boost::remove_const<TIter>::type value_type;
226  typedef TIter& reference;
227  typedef TIter* pointer;
228 
229  // Construction
231  major1_iterator(pointer arrayBegin, size_type index, difference_type /*size*/)
232  :m_pos(arrayBegin), m_index(index){}
233 
234  template<class U>
236  :m_pos(iter.m_pos), m_index(iter.m_index){}
237 
238  template<class U>
240  m_pos = iter.m_pos;
241  m_index = iter.m_index;
242  return *this;
243  }
244 
245  // Arithmetic
247  ++m_index;
248  m_pos += m_index;
249  return *this;
250  }
252  m_pos -= m_index;
253  --m_index;
254  return *this;
255  }
256  major1_iterator& operator += (difference_type n) {
257  m_pos += offset(n);
258  m_index += n;
259  return *this;
260  }
261  major1_iterator& operator -= (difference_type n) {
262  m_pos += offset(-n);
263  m_index -= n;
264  return *this;
265  }
266  template<class U>
267  difference_type operator - (major1_iterator<U> const& it) const {
268  return m_index - it.m_index;
269  }
270 
271  // Dereference
272  reference operator*() const {
273  return *m_pos;
274  }
275  reference operator [](difference_type n) const {
276  return m_pos[offset(n)];
277  }
278 
279  // Index
280  size_type index() const {
281  return m_index;
282  }
283 
284  // Comparison
285  template<class U>
286  bool operator == (major1_iterator<U> const& it) const {
287  return m_index == it.m_index;
288  }
289  template<class U>
290  bool operator < (major1_iterator<U> const& it) const {
291  return m_index < it.m_index;
292  }
293 
294  private:
295  pointer m_pos;
296  difference_type m_index;
297  template<class> friend class major1_iterator;
298  };
299 
300  template<class TIter>
302  public random_access_iterator_base<
303  major2_iterator<TIter>,
304  typename boost::remove_const<T>::type,
305  packed_random_access_iterator_tag
306  >{
307  private:
308  difference_type offset(difference_type n)const{
309  difference_type k = m_size-m_index-1;
310  if(n >= 0){
311  return (2*k*n-n*n+n)/2;
312  }else{
313  n*=-1;
314  k+= n;
315  return -(2*k*n-n*n+n)/2;
316  }
317  }
318  public:
319  typedef typename boost::remove_const<TIter>::type value_type;
320  typedef TIter& reference;
321  typedef TIter* pointer;
322 
323  // Construction
325  major2_iterator(pointer arrayBegin, size_type index, difference_type size)
326  :m_pos(arrayBegin), m_index(index), m_size(size){}
327 
328  template<class U>
330  :m_pos(iter.m_pos), m_index(iter.m_index), m_size(iter.m_size){}
331 
332  template<class U>
334  m_pos = iter.m_pos;
335  m_index = iter.m_index;
336  m_size = iter.m_size;
337  return *this;
338  }
339 
340  // Arithmetic
342  ++m_index;
343  m_pos += m_size-m_index;
344  return *this;
345  }
347  m_pos -= m_size-m_index;
348  --m_index;
349  return *this;
350  }
351  major2_iterator& operator += (difference_type n) {
352  m_pos += offset(n);
353  m_index += n;
354  return *this;
355  }
356  major2_iterator& operator -= (difference_type n) {
357  m_pos += offset(-n);
358  m_index -= n;
359  return *this;
360  }
361  template<class U>
362  difference_type operator - (major2_iterator<U> const& it) const {
363  return m_index - it.m_index;
364  }
365 
366  // Dereference
367  reference operator*() const {
368  return *m_pos;
369  }
370  reference operator [](difference_type n) const {
371  return m_pos[offset(n)];
372  }
373 
374  // Index
375  size_type index() const {
376  return m_index;
377  }
378 
379  // Comparison
380  template<class U>
381  bool operator == (major2_iterator<U> const& it) const {
382  return m_index == it.m_index;
383  }
384  template<class U>
385  bool operator < (major2_iterator<U> const& it) const {
386  return m_index < it.m_index;
387  }
388 
389  private:
390  pointer m_pos;
391  difference_type m_index;
392  difference_type m_size;
393  template<class> friend class major2_iterator;
394  };
395 
396  typedef typename boost::mpl::if_<
397  boost::is_same<Orientation,row_major>,
398  dense_storage_iterator<value_type>,
399  typename boost::mpl::if_c<
400  TriangularType::is_upper,
403  >::type
404  >::type row_iterator;
405  typedef typename boost::mpl::if_<
406  boost::is_same<Orientation,row_major>,
407  typename boost::mpl::if_c<
408  TriangularType::is_upper,
410  major1_iterator<value_type>
411  >::type,
412  dense_storage_iterator<value_type>
414 
415  typedef typename boost::mpl::if_<
416  boost::is_same<Orientation,row_major>,
417  dense_storage_iterator<value_type const>,
418  typename boost::mpl::if_c<
419  TriangularType::is_upper,
422  >::type
424  typedef typename boost::mpl::if_<
425  boost::is_same<Orientation,row_major>,
426  typename boost::mpl::if_c<
427  TriangularType::is_upper,
429  major1_iterator<value_type const>
430  >::type,
431  dense_storage_iterator<value_type const>
433 
434 public:
435 
436  const_row_iterator row_begin(index_type i) const {
437  std::size_t index = TriangularType::is_upper?i:0;
438  return const_row_iterator(
439  storage()+orientation::element(i,index,size1())
440  ,index
441  ,orientation::stride2(size1(),size2())//1 if row_major, size2() otherwise
442  );
443  }
444  const_row_iterator row_end(index_type i) const {
445  std::size_t index = TriangularType::is_upper?size2():i+1;
446  return const_row_iterator(
447  storage() + orientation::element(i, index, size1())
448  ,index
449  ,orientation::stride2(size1(),size2())//1 if row_major, size2() otherwise
450  );
451  }
452  row_iterator row_begin(index_type i){
453  std::size_t index = TriangularType::is_upper?i:0;
454  return row_iterator(
455  storage() + orientation::element(i, index, size1())
456  ,index
457  ,orientation::stride2(size1(),size2())//1 if row_major, size2() otherwise
458  );
459  }
460  row_iterator row_end(index_type i){
461  std::size_t index = TriangularType::is_upper?size2():i+1;
462  return row_iterator(
463  storage() + orientation::element(i, index, size1())
464  ,index
465  ,orientation::stride2(size1(),size2())//1 if row_major, size2() otherwise
466  );
467  }
468 
469  const_column_iterator column_begin(index_type i) const {
470  std::size_t index = TriangularType::is_upper?0:i;
471  return const_column_iterator(
472  storage() + orientation::element(index, i, size1())
473  ,index
474  ,orientation::stride1(size1(),size2())//size1() if row_major, 1 otherwise
475  );
476  }
477  const_column_iterator column_end(index_type i) const {
478  std::size_t index = TriangularType::is_upper?i+1:size2();
479  return const_column_iterator(
480  storage() + orientation::element(index, i, size1())
481  ,index
482  ,orientation::stride1(size1(),size2())//size1() if row_major, 1 otherwise
483  );
484  }
485  column_iterator column_begin(index_type i){
486  std::size_t index = TriangularType::is_upper?0:i;
487  return column_iterator(
488  storage() + orientation::element(index, i, size1())
489  ,index
490  ,orientation::stride1(size1(),size2())//size1() if row_major, 1 otherwise
491  );
492  }
493  column_iterator column_end(index_type i){
494  std::size_t index = TriangularType::is_upper?i+1:size2();
495  return column_iterator(
496  storage() + orientation::element(index, i, size1())
497  ,index
498  ,orientation::stride1(size1(),size2())//size1() if row_major, 1 otherwise
499  );
500  }
501 
502  // Serialization
503  template<class Archive>
504  void serialize(Archive& ar, const unsigned int /* file_version */) {
505 
506  // we need to copy to a collection_size_type to get a portable
507  // and efficient boost::serialization
508  boost::serialization::collection_size_type s(m_size);
509 
510  // serialize the sizes
511  ar & boost::serialization::make_nvp("size",s);
512 
513  // copy the values back if loading
514  if (Archive::is_loading::value) {
515  m_size = s;
516  }
517  ar & boost::serialization::make_nvp("data",m_data);
518  }
519 
520 private:
521  size_type m_size;
522  array_type m_data;
523 };
524 
525 template<class T, class Orientation, class TriangularType>
526 struct const_expression<triangular_matrix<T,Orientation, TriangularType> >{
528 };
529 template<class T, class Orientation, class TriangularType>
530 struct const_expression<triangular_matrix<T,Orientation, TriangularType> const>{
532 };
533 
534 }
535 }
536 
537 #endif