atlas  0.6
matrix.h
Go to the documentation of this file.
1 /*
2  This is matrix.h
3
4  Copyright (C) 2004,2005 Fokko du Cloux
5  Copyright (C) 2006-2016 Marc van Leeuwen
6  part of the Atlas of Lie Groups and Representations
7
8  For license information see the LICENSE file
9 */
10
11 #ifndef MATRIX_H /* guard against multiple inclusions */
12 #define MATRIX_H
13
14 #include <vector>
15 #include <stdexcept>
16
17 #include "matrix_fwd.h"
18
19 #include "tags.h"
20
21 // extra defs for windows compilation -spc
22 #ifdef WIN32
23 #include "constants.h"
24 #endif
25
26 namespace atlas {
27
28 /******** function definitions ***********************************************/
29
30 namespace matrix {
31
32 template<typename C>
33  std::vector<Vector<C> > standard_basis(size_t n);
34
35 template<typename C>
36  void initBasis(std::vector<Vector<C> >& v, size_t n)
37  { v=standard_basis<C>(n); }
38
39 template<typename C>
40  PID_Matrix<C>& operator+= (PID_Matrix<C>& A, C c); // |A=A+c|, avoiding copy
41
42 // non-destructive form takes value parameter, which maybe selects rvalue copy
43 template<typename C>
44  PID_Matrix<C> operator+ (PID_Matrix<C> A, C c) // add scalar matrix
45  { return A += c; }
46
47 template<typename C>
48  PID_Matrix<C>& operator-= (PID_Matrix<C>& A, C c) { return A += -c; }
49
50 template<typename C>
51  PID_Matrix<C> operator- (PID_Matrix<C> A, C c) { return A += -c; }
52
53 template<typename C>
54  PID_Matrix<C> operator- (C c, PID_Matrix<C> A) { return A.negate() += c; }
55
56 }
57
58 /******** type definitions ***************************************************/
59
60 namespace matrix {
61
62 template <typename C> class Matrix;
63
64 template<typename C>
65  class Vector : public std::vector<C>
66 {
67  typedef std::vector<C> base;
68 public:
69  Vector () : base() {}
70  explicit Vector (size_t n) : base(n) {} // entries remain undefined!
71  explicit Vector (const base& b) : base(b) {} // std::vector -> Vector
72  Vector (size_t n,C c) : base(n,c) {}
73  template<typename I>
74  Vector (I b, I e) : base(b,e) {} // construction from iterators
75
76  Vector& operator+= (const Vector&);
77  Vector& operator-= (const Vector&);
78  Vector& operator*= (C);
79  Vector& negate (); // negates argument in place
80  Vector& negate_add (const Vector& y); // reversed -=, so *this = y - *this
81
82  // the following two methods do not take an end iterator; count is |size()|
83  template<typename I> Vector& add(I b,C c); // add |Vector(b,b+size())*c|
84  template<typename I> Vector& subtract( I b,C c) { return add(b,-c); }
85
86
87  template<typename C1> C1 dot (const Vector<C1>& v) const;
88  bool isZero() const;
89
90 #if __GNUC__ < 4 || \
91  __GNUC__ == 4 && ( __GNUC_MINOR__ < 8 || \
92  __GNUC_MINOR__ == 8 && __GNUC_PATCHLEVEL__ < 1)
93  Vector operator+ (const Vector& v) const
94  { Vector result(*this); return result +=v; }
95  Vector operator- (const Vector&v) const
96  { Vector result(*this); return result -=v; }
97  Vector operator* (C c) const
98  { Vector result(*this); return result *=c; }
100  { Vector result(*this); return result.negate(); }
101 #else
102  Vector operator+ (const Vector& v) const &
103  { Vector result(*this); return result +=v; }
104  Vector operator- (const Vector&v) const &
105  { Vector result(*this); return result -=v; }
106  Vector operator* (C c) const &
107  { Vector result(*this); return result *=c; }
108  Vector operator- () const &
109  { Vector result(*this); return result.negate(); }
110
111  Vector operator+ (const Vector& v) &&
112  { Vector result(std::move(*this)); return result +=v; }
113  Vector operator- (const Vector&v) &&
114  { Vector result(std::move(*this)); return result -=v; }
115  Vector operator* (C c) &&
116  { Vector result(std::move(*this)); return result *=c; }
117  Vector operator- () &&
118  { Vector result(std::move(*this)); return result.negate(); }
119 #endif
120
121  template<typename C1>
122  Vector<C1> scaled (C1 c) const; // like operator*, but forces type to C1
123
124  Matrix<C> row_matrix() const; // vector as 1-row matrix
125  Matrix<C> column_matrix() const; // vector as 1-column matrix
126 }; // Vector
127
128 // these are external functions, to allow instantiating |Vector<Pol>|
129 template<typename C>
130  Vector<C>& operator/= (Vector<C>& v, C c); // this division must be exact
131 template<typename C>
132  Vector<C>& divide (Vector<C>& v, C c); // this integer division rounds down
133 template<typename C>
134  Vector<C>& operator%= (Vector<C>& v, C c);
135 template<typename C>
136  Vector<C> operator/ (Vector<C> v,C c) { return v /= c; }
137
138
139 template<typename C> class Matrix_base
140 {
141  size_t d_rows;
142  size_t d_columns;
143  protected: // derived classes may sometimes need to acces this
144  Vector<C> d_data; // vector of elements, concatenated by rows
145
146  public:
147
148 // constructors
149  Matrix_base(): d_rows(0),d_columns(0),d_data() {}
150
151  Matrix_base(size_t m, size_t n) : d_rows(m),d_columns(n),d_data(m*n) {}
152  Matrix_base(size_t m, size_t n, const C& c)
153  : d_rows(m), d_columns(n), d_data(m*n,c) {}
154
156  (const std::vector<Vector<C> >&, size_t n_rows); // with explicit #rows
157
158  template<typename I> // from sequence of columns obtained via iterator
159  Matrix_base(I begin, I end, size_t n_rows, tags::IteratorTag);
160
161  void swap(Matrix_base&);
162
163  // accessors
164  size_t numRows() const { return d_rows; }
165  size_t numColumns() const { return d_columns; }
166  size_t rowSize() const { return d_columns; }
167  size_t columnSize() const { return d_rows; }
168
169  const C& operator() (size_t i,size_t j)
170  const { return d_data[i*d_columns+j]; }
171
172  void get_row(Vector<C>&, size_t) const;
173  void get_column(Vector<C>&, size_t) const;
174
175  Vector<C> row(size_t i) const { return Vector<C>(at(i,0),at(i+1,0)); }
176  std::vector<Vector<C> > rows() const;
177  Vector<C> column(size_t j) const { Vector<C> c; get_column(c,j); return c; }
178  std::vector<Vector<C> > columns() const;
179
180  bool operator== (const Matrix_base<C>&) const;
181  bool operator!= (const Matrix_base<C>& m) const {return not(operator==(m)); }
182
183  bool is_zero() const; // whether all entries are zero
184  bool isEmpty() const { return d_data.size() == 0; }
185 // manipulators
186  C& operator() (size_t i, size_t j) { return d_data[i*d_columns+j]; }
187
188  void set_row(size_t,const Vector<C>&);
189  void set_column(size_t,const Vector<C>&);
190  void add_row(const Vector<C>&);
191  void add_column(const Vector<C>&);
192
193 // resize undefined; use |Matrix<C>(m,n).swap(M)| instead of |M.resize(m,n)|
194
195  void eraseColumn(size_t);
196  void eraseRow(size_t);
197  void reset() { d_data.assign(d_data.size(),C(0)); }
198
199  protected:
200  const C* at (size_t i,size_t j) const { return &operator()(i,j); }
201  C* at (size_t i,size_t j) { return &operator()(i,j); }
202 }; // |template<typename C> class Matrix_base|
203
204 template<typename C> class Matrix : public Matrix_base<C>
205 {
206  protected:
208
209  public:
210 // constructors
211  Matrix() : base() {}
212  Matrix(const base& M) : base(M) {}
213  Matrix(base&& M) : base(std::move(M)) {}
214
215 #if __GNUC__ > 4 || __GNUC__ == 4 && __GNUC_MINOR__ >= 8
216  // that is, if compiler version is sufficiently new
217  using base::base; // inherit all constructors
218 #else
219  Matrix(size_t m, size_t n) : base(m,n) {}
220  Matrix(size_t m, size_t n, const C& c) : base(m,n,c) {}
221
222  Matrix(const std::vector<Vector<C> >&cols, size_t n_rows)
223  : base(cols,n_rows) {}
224
225  template<typename I> // from sequence of columns obtained via iterator
226  Matrix(I begin, I end, size_t n_rows, tags::IteratorTag)
227  : base(begin,end,n_rows,tags::IteratorTag()) {}
228 #endif
229  explicit Matrix(size_t n); // square identity matrix
230
231 // accessors
232  // no non-destructive additive matrix operations; no need for them (yet?)
233
234 #if __GNUC__ < 4 || \
235  __GNUC__ == 4 && ( __GNUC_MINOR__ < 8 || \
236  __GNUC_MINOR__ == 8 && __GNUC_PATCHLEVEL__ < 1)
237  // that is, if compiler version is too old
238  Matrix transposed() const;
239  Matrix negative_transposed() const { return transposed().negate(); }
240 #else
241  Matrix transposed() const &;
242  Matrix transposed() &&
243  { return Matrix(std::move(*this)).transpose(); }
244  Matrix negative_transposed() const & { return transposed().negate(); }
245  Matrix negative_transposed() &&
246  { return Matrix(std::move(*this)).negate().transpose(); }
247 #endif
248
249  // there is no point in trying to do these matrix multiplications in-place
250  template<typename C1> Vector<C1> operator* (const Vector<C1>&) const;
251  template<typename C1> Vector<C1> right_prod(const Vector<C1>&) const;
252
253  template<typename C1> void apply_to(Vector<C1>& v) const
254  { v= operator*(v); }
255  template<typename C1> void right_mult(Vector<C1>& v) const
256  { v= right_prod(v); }
257
258  Matrix operator* (const Matrix&) const;
259
260 // manipulators
261
262  Matrix& operator+= (const Matrix&);
263  Matrix& operator-= (const Matrix&);
264  Matrix& operator*= (const Matrix& Q) { return *this = *this * Q; }
265  Matrix& leftMult (const Matrix& P) { return *this = P * *this; }
266
267  Matrix& negate() { base::d_data.negate(); return *this; }
268  Matrix& transpose();
269
270  // secondary manipulators
271
272  void rowOperation(size_t, size_t, const C&);
273  void columnOperation(size_t j, size_t k, const C& c); // |col(j) += c*col(k)|
274
275  void rowMultiply(size_t i, C f);
276  void columnMultiply(size_t j,C f);
277
278  void swapColumns(size_t, size_t);
279  void swapRows(size_t, size_t);
280
281 }; // |template<typename C> class Matrix|
282
283 // The following derived class is for integer types only
284 template<typename C> class PID_Matrix : public Matrix<C>
285 {
286  typedef Matrix<C> base;
287
288  public:
289  PID_Matrix() : base() {}
290  PID_Matrix(const base& M) : base(M) {}
291  PID_Matrix(base&& M) : base(std::move(M)) {}
292
293 // forward constructors to Matrix
294 #if __GNUC__ > 4 || __GNUC__ == 4 && __GNUC_MINOR__ >= 8
295  // that is, if compiler version is sufficiently new
296  using base::base; // inherit all constructors
297 #else
298  PID_Matrix(size_t m, size_t n) : base(m,n) {}
299  PID_Matrix(size_t m, size_t n, const C& c) : base(m,n,c) {}
300
301  explicit PID_Matrix(size_t n) : base(n) {} // square identity matrix
302  PID_Matrix(const std::vector<Vector<C> >&cols, size_t n_rows)
303  : base(cols,n_rows) {}
304
305  template<typename I> // from sequence of columns obtained via iterator
306  PID_Matrix(I begin, I end, size_t n_rows, tags::IteratorTag)
307  : base(begin,end,n_rows,tags::IteratorTag()) {}
308 #endif
309
310 // manipulators
312  { base::operator+=(M); return *this; }
314  { base::operator-=(M); return *this; }
315  PID_Matrix& operator*= (const Matrix<C>& Q)
316  { base::operator*=(Q); return *this; }
318  { base::leftMult(P); return *this; }
319
320  PID_Matrix& operator/= (const C& c) throw (std::runtime_error)
321  { if (c != C(1)) base::base::d_data /= c; return *this; }
322
323  PID_Matrix& negate() { base::negate(); return *this; }
324  PID_Matrix& transpose() { base::transpose(); return *this; }
325
326  PID_Matrix& invert();
327  PID_Matrix& invert(C& d);
328
329 // accessors
330 #if __GNUC__ < 4 || \
331  __GNUC__ == 4 && ( __GNUC_MINOR__ < 8 || \
332  __GNUC_MINOR__ == 8 && _GNUC_PATCHLEVEL__ < 1)
333  // that is, if compiler version is too old
334  PID_Matrix transposed() const { return PID_Matrix(base::transposed()); }
336  { return PID_Matrix(base::negative_transposed()); }
337  PID_Matrix inverse() const { return PID_Matrix(*this).invert(); }
338  PID_Matrix inverse(C& d) const { return PID_Matrix(*this).invert(d); }
339 #else
340  PID_Matrix transposed() const & { return PID_Matrix(base::transposed()); }
341  PID_Matrix transposed() &&
342  { return PID_Matrix(std::move(*this)).transpose(); }
343  PID_Matrix negative_transposed() const &
344  { return PID_Matrix(base::negative_transposed()); }
345  PID_Matrix negative_transposed() &&
346  { return PID_Matrix(std::move(*this)).negate().transpose(); }
347  PID_Matrix inverse() const & { return PID_Matrix(*this).invert(); }
348  PID_Matrix inverse() && { return PID_Matrix(std::move(*this)).invert(); }
349  PID_Matrix inverse(C& d) const & { return PID_Matrix(*this).invert(d); }
350  PID_Matrix inverse(C& d) && { return PID_Matrix(std::move(*this)).invert(d); }
351 #endif
352
353  using base::operator*;
354
356  { return PID_Matrix(base::operator*(Q)); }
357
358
359
360  PID_Matrix on_basis(const std::vector<Vector<C> >& basis) const;
361
362  // secondary accessors (for inversion algorithms)
363  bool divisible(C) const;
364  PID_Matrix block(size_t i0, size_t j0, size_t i1, size_t j1) const;
365
366 }; // |class PID_Matrix|
367
368 } // |namespace matrix|
369
370 } // |namespace atlas|
371
372 #endif
Definition: ctangle.c:136
Vector(const base &b)
Definition: matrix.h:71
PID_Matrix inverse() const
Definition: matrix.h:337
Vector< C > operator/(Vector< C > v, C c)
Definition: matrix.h:136
PID_Matrix negative_transposed() const
Definition: matrix.h:335
size_t d_columns
Definition: matrix.h:142
size_t rowSize() const
Definition: matrix.h:166
bool operator!=(const type_expr &x, const type_expr &y)
Definition: axis-types.h:374
PID_Matrix()
Definition: matrix.h:289
Matrix(const base &M)
Definition: matrix.h:212
Definition: Atlas.h:79
Matrix negative_transposed() const
Definition: matrix.h:239
Vector & negate()
Definition: matrix.cpp:138
PID_Matrix & negate()
Definition: matrix.h:323
PID_Matrix< C > operator-(PID_Matrix< C > A, C c)
Definition: matrix.h:51
std::vector< C > base
Definition: matrix.h:67
Matrix_base(size_t m, size_t n, const C &c)
Definition: matrix.h:152
void reset()
Definition: matrix.h:197
Vector< C > & divide(Vector< C > &v, C c)
Definition: matrix.cpp:111
Matrix_base()
Definition: matrix.h:149
Matrix_base(size_t m, size_t n)
Definition: matrix.h:151
Matrix(I begin, I end, size_t n_rows, tags::IteratorTag)
Definition: matrix.h:226
Vector< C > & operator%=(Vector< C > &v, C c)
Definition: matrix.cpp:127
Matrix()
Definition: matrix.h:211
Matrix & negate()
Definition: matrix.h:267
const C * at(size_t i, size_t j) const
Definition: matrix.h:200
Matrix & leftMult(const Matrix &P)
Definition: matrix.h:265
void swap(simple_list< T, Alloc > &x, simple_list< T, Alloc > &y)
Definition: sl_list.h:674
Definition: lists.cpp:13
Vector & subtract(I b, C c)
Definition: matrix.h:84
PID_Matrix & invert()
Definition: matrix.cpp:519
PID_Matrix< C > operator+(PID_Matrix< C > A, C c)
Definition: matrix.h:44
Vector()
Definition: matrix.h:69
size_t columnSize() const
Definition: matrix.h:167
Vector(size_t n, C c)
Definition: matrix.h:72
Matrix< C > base
Definition: matrix.h:286
PID_Matrix transposed() const
Definition: matrix.h:334
Vector< C > column(size_t j) const
Definition: matrix.h:177
Matrix(size_t m, size_t n, const C &c)
Definition: matrix.h:220
size_t d_rows
Definition: matrix.h:141
Matrix(base &&M)
Definition: matrix.h:213
Matrix & transpose()
Definition: matrix.cpp:503
RationalVector< C2 > operator*(const matrix::Matrix< C1 > &M, const RationalVector< C2 > &v)
Definition: ratvec.cpp:158
size_t numRows() const
Definition: matrix.h:164
PID_Matrix(size_t m, size_t n)
Definition: matrix.h:298
PID_Matrix< C > & operator+=(PID_Matrix< C > &A, C c)
Definition: matrix.cpp:733
Matrix_base< C > base
Definition: matrix.h:207
C * at(size_t i, size_t j)
Definition: matrix.h:201
void right_mult(Vector< C1 > &v) const
Definition: matrix.h:255
Definition: tags.h:32
void apply_to(Vector< C1 > &v) const
Definition: matrix.h:253
Vector< C > d_data
Definition: matrix.h:144
void initBasis(std::vector< Vector< C > > &v, size_t n)
Definition: matrix.h:36
Matrix(size_t m, size_t n)
Definition: matrix.h:219
size_t numColumns() const
Definition: matrix.h:165
PID_Matrix(size_t m, size_t n, const C &c)
Definition: matrix.h:299
PID_Matrix< C > & operator-=(PID_Matrix< C > &A, C c)
Definition: matrix.h:48
PID_Matrix(base &&M)
Definition: matrix.h:291
Vector< C > & operator/=(Vector< C > &v, C c)
Definition: matrix.cpp:99
Definition: Atlas.h:78
Vector(size_t n)
Definition: matrix.h:70
Vector< C > row(size_t i) const
Definition: matrix.h:175
PID_Matrix & transpose()
Definition: matrix.h:324
Matrix(const std::vector< Vector< C > > &cols, size_t n_rows)
Definition: matrix.h:222
std::vector< Vector< C > > standard_basis(size_t r)
Definition: matrix.cpp:722
PID_Matrix inverse(C &d) const
Definition: matrix.h:338
unsigned long n
Definition: axis.cpp:77
bool operator==(const type_expr &x, const type_expr &y)
Definition: axis-types.cpp:257
simple_list< T, Alloc >::const_iterator end(const simple_list< T, Alloc > &l)
Definition: sl_list.h:650
Definition: Atlas.h:38
Definition of dummy argument tags used for constructor overloading.
bool isEmpty() const
Definition: matrix.h:184
PID_Matrix(I begin, I end, size_t n_rows, tags::IteratorTag)
Definition: matrix.h:306
PID_Matrix(const base &M)
Definition: matrix.h:290
void basis(std::vector< matrix::Vector< int > > &b, const bitmap::BitMap &B, const FiniteAbelianGroup &A)
Definition: abelian.cpp:427
const expr & e
Definition: axis.cpp:95
Definition: common.h:36
PID_Matrix(size_t n)
Definition: matrix.h:301
Definition: Atlas.h:81
PID_Matrix(const std::vector< Vector< C > > &cols, size_t n_rows)
Definition: matrix.h:302
void transpose(Endomorphism &e, const FiniteAbelianGroup &A)
Definition: abelian.cpp:675
Vector(I b, I e)
Definition: matrix.h:74
Vertex v
Definition: graph.cpp:116
Definition: Atlas.h:80
PID_Matrix & leftMult(const Matrix< C > &P)
Definition: matrix.h:317