1 #ifndef MSP_LINAL_MATRIX_H_
2 #define MSP_LINAL_MATRIX_H_
13 A general mathematical matrix with M rows and N columns.
15 template<typename T, unsigned M, unsigned N>
19 typedef T ElementType;
28 Matrix(const Matrix<U, M, N> &);
30 static Matrix identity();
31 static Matrix from_columns(const Vector<T, M> *);
32 static Matrix from_rows(const Vector<T, N> *);
34 unsigned rows() const { return M; }
35 unsigned columns() const { return N; }
37 T &element(unsigned i, unsigned j) { return data[i+M*j]; }
38 const T &element(unsigned i, unsigned j) const { return data[i+M*j]; }
39 T &operator()(unsigned i, unsigned j) { return element(i, j); }
40 const T &operator()(unsigned i, unsigned j) const { return element(i, j); }
42 Vector<T, M> column(unsigned i) const { return Vector<T, M>(data+M*i); }
43 Vector<T, N> row(unsigned i) const { return Vector<T, N>(data+i, M); }
45 template<unsigned P, unsigned Q>
46 Matrix<T, P, Q> select(const Vector<unsigned, P> &, const Vector<unsigned, Q> &) const;
48 template<unsigned P, unsigned Q>
49 Matrix<T, P, Q> block(unsigned, unsigned) const;
51 Matrix &operator*=(T);
52 Matrix &operator*=(const Matrix &);
53 Matrix &operator/=(T);
54 Matrix &operator+=(const Matrix &);
55 Matrix &operator-=(const Matrix &);
59 Matrix &exchange_columns(unsigned, unsigned);
60 Matrix &multiply_column(unsigned, T);
61 Matrix &add_column(unsigned, unsigned, T);
64 template<typename T, unsigned M, unsigned N>
65 inline Matrix<T, M, N>::Matrix()
67 std::fill(data, data+M*N, T());
70 template<typename T, unsigned M, unsigned N>
71 inline Matrix<T, M, N>::Matrix(const T *d)
73 std::copy(d, d+M*N, data);
76 template<typename T, unsigned M, unsigned N>
78 inline Matrix<T, M, N>::Matrix(const Matrix<U, M, N> &other)
80 for(unsigned i=0; i<M; ++i)
81 for(unsigned j=0; j<N; ++j)
82 element(i, j) = other(i, j);
85 template<typename T, unsigned M, unsigned N>
86 inline Matrix<T, M, N> Matrix<T, M, N>::identity()
88 static_assert(M==N, "An identity matrix must be square");
90 for(unsigned i=0; i<M; ++i)
95 template<typename T, unsigned M, unsigned N>
96 inline Matrix<T, M, N> Matrix<T, M, N>::from_columns(const Vector<T, M> *v)
99 for(unsigned i=0; i<M; ++i)
100 for(unsigned j=0; j<N; ++j)
105 template<typename T, unsigned M, unsigned N>
106 inline Matrix<T, M, N> Matrix<T, M, N>::from_rows(const Vector<T, N> *v)
109 for(unsigned i=0; i<M; ++i)
110 for(unsigned j=0; j<N; ++j)
115 template<typename T, unsigned M, unsigned N>
116 template<unsigned P, unsigned Q>
117 inline Matrix<T, P, Q> Matrix<T, M, N>::select(const Vector<unsigned, P> &row_indices, const Vector<unsigned, Q> &col_indices) const
120 for(unsigned j=0; j<P; ++j)
121 for(unsigned i=0; i<Q; ++i)
122 r(j, i) = element(row_indices[j], col_indices[i]);
126 template<typename T, unsigned M, unsigned N>
127 template<unsigned P, unsigned Q>
128 inline Matrix<T, P, Q> Matrix<T, M, N>::block(unsigned y, unsigned x) const
131 for(unsigned j=0; j<P; ++j)
132 for(unsigned i=0; i<Q; ++i)
133 r(j, i) = element(y+j, x+i);
137 template<typename T, unsigned M, unsigned N>
138 inline Matrix<T, M, N> &Matrix<T, M, N>::operator*=(T s)
140 for(unsigned i=0; i<M*N; ++i)
145 template<typename T, unsigned M, unsigned N>
146 inline Matrix<T, M, N> &Matrix<T, M, N>::operator*=(const Matrix<T, M, N> &m)
148 static_assert(M==N, "Multiplication-assignment is only possible on square matrices");
149 return *this = *this*m;
152 template<typename T, unsigned M, unsigned N>
153 inline Matrix<T, M, N> operator*(const Matrix<T, M, N> &m, T s)
155 Matrix<T, M, N> r(m);
159 template<typename T, unsigned M, unsigned N>
160 inline Matrix<T, M, N> operator*(T s, const Matrix<T, M, N> &m)
165 template<typename T, unsigned M, unsigned P, unsigned N>
166 inline Matrix<T, M, N> operator*(const Matrix<T, M, P> &m1, const Matrix<T, P, N> &m2)
169 for(unsigned i=0; i<M; ++i)
170 for(unsigned j=0; j<N; ++j)
171 for(unsigned k=0; k<P; ++k)
172 r(i, j) += m1(i, k)*m2(k, j);
176 template<typename T, unsigned M, unsigned N>
177 inline Vector<T, M> operator*(const Matrix<T, M, N> &m, const Vector<T, N> &v)
180 for(unsigned i=0; i<M; ++i)
181 for(unsigned j=0; j<N; ++j)
182 r[i] += m(i, j)*v[j];
186 template<typename T, unsigned M, unsigned N>
187 inline Vector<T, N> operator*(const Vector<T, M> &v, const Matrix<T, M, N> &m)
190 for(unsigned j=0; j<N; ++j)
191 for(unsigned i=0; i<M; ++i)
192 r[j] += v[i]*m(i, j);
196 template<typename T, unsigned M, unsigned N>
197 inline Matrix<T, M, N> &Matrix<T, M, N>::operator/=(T s)
199 for(unsigned i=0; i<M*N; ++i)
204 template<typename T, unsigned M, unsigned N>
205 inline Matrix<T, M, N> operator/(const Matrix<T, M, N> &m, T s)
207 Matrix<T, M, N> r(m);
211 template<typename T, unsigned M, unsigned N>
212 inline Matrix<T, M, N> &Matrix<T, M, N>::operator+=(const Matrix<T, M, N> &m)
214 for(unsigned i=0; i<M*N; ++i)
215 data[i] += m.data[i];
219 template<typename T, unsigned M, unsigned N>
220 inline Matrix<T, M, N> operator+(const Matrix<T, M, N> &m1, const Matrix<T, M, N> &m2)
222 Matrix<T, M, N> r(m1);
226 template<typename T, unsigned M, unsigned N>
227 inline Matrix<T, M, N> &Matrix<T, M, N>::operator-=(const Matrix<T, M, N> &m)
229 for(unsigned i=0; i<M*N; ++i)
230 data[i] -= m.data[i];
234 template<typename T, unsigned M, unsigned N>
235 inline Matrix<T, M, N> operator-(const Matrix<T, M, N> &m1, const Matrix<T, M, N> &m2)
237 Matrix<T, M, N> r(m1);
241 template<typename T, unsigned M, unsigned N>
242 inline Matrix<T, M, N>& Matrix<T, M, N>::invert()
244 static_assert(M==N, "Inversion is only possible on square matrices");
245 Matrix<T, M, N> r = identity();
246 gauss_jordan(*this, r);
250 template<typename T, unsigned S>
251 inline Matrix<T, S, S> invert(const Matrix<T, S, S> &m)
253 Matrix<T, S, S> temp = m;
254 Matrix<T, S, S> r = Matrix<T, S, S>::identity();
255 return gauss_jordan(temp, r);
258 template<typename T, unsigned M, unsigned N>
259 inline bool operator==(const Matrix<T, M, N> &a, const Matrix<T, M, N> &b)
261 for(unsigned j=0; j<N; ++j)
262 for(unsigned i=0; i<M; ++i)
268 template<typename T, unsigned M, unsigned N>
269 inline Matrix<T, M, N> &Matrix<T, M, N>::exchange_columns(unsigned i, unsigned j)
272 for(unsigned k=0; k<M; ++k)
273 swap(element(k, i), element(k, j));
277 template<typename T, unsigned M, unsigned N>
278 inline Matrix<T, M, N> &Matrix<T, M, N>::multiply_column(unsigned i, T s)
280 for(unsigned k=0; k<M; ++k)
285 template<typename T, unsigned M, unsigned N>
286 inline Matrix<T, M, N> &Matrix<T, M, N>::add_column(unsigned i, unsigned j, T s)
288 for(unsigned k=0; k<M; ++k)
289 element(k, j) += element(k, i)*s;
293 template<typename T, unsigned M, unsigned N>
294 inline Matrix<T, N, M> transpose(const Matrix<T, M, N> &m)
297 for(unsigned j=0; j<N; ++j)
298 for(unsigned i=0; i<M; ++i)
303 template<typename T, unsigned M, unsigned N>
304 inline std::ostream &operator<<(std::ostream &s, const Matrix<T, M, N> &m)
306 s << "Matrix" << M << 'x' << N << '(';
307 for(unsigned i=0; i<N; ++i)
312 for(unsigned j=0; j<M; ++j)