]> git.tdb.fi Git - libs/math.git/blob - source/linal/matrix.h
Another big batch of stuff
[libs/math.git] / source / linal / matrix.h
1 #ifndef MSP_LINAL_MATRIX_H_
2 #define MSP_LINAL_MATRIX_H_
3
4 #include <algorithm>
5 #include "vector.h"
6
7 namespace Msp {
8 namespace LinAl {
9
10 /**
11 A general mathematical matrix.
12 */
13 template<typename T, unsigned M, unsigned N>
14 class Matrix
15 {
16 private:
17         T data[M*N];
18
19 public:
20         Matrix();
21         Matrix(const T *);
22         template<typename U>
23         Matrix(const Matrix<U, M, N> &);
24
25         static Matrix from_columns(const Vector<T, M> *);
26         static Matrix from_rows(const Vector<T, N> *);
27
28         T &element(unsigned, unsigned);
29         const T &element(unsigned, unsigned) const;
30         T &operator()(unsigned, unsigned);
31         const T &operator()(unsigned, unsigned) const;
32
33         Matrix &operator*=(T);
34         Matrix &operator/=(T);
35         Matrix &operator+=(const Matrix &);
36         Matrix &operator-=(const Matrix &);
37
38         Matrix &exchange_rows(unsigned, unsigned);
39         Matrix &multiply_row(unsigned, T);
40         Matrix &add_row(unsigned, unsigned, T);
41 };
42
43 template<typename T, unsigned M, unsigned N>
44 inline Matrix<T, M, N>::Matrix()
45 {
46         std::fill(data, data+M*N, T());
47 }
48
49 template<typename T, unsigned M, unsigned N>
50 inline Matrix<T, M, N>::Matrix(const T *d)
51 {
52         std::copy(d, d+M*N, data);
53 }
54
55 template<typename T, unsigned M, unsigned N>
56 template<typename U>
57 inline Matrix<T, M, N>::Matrix(const Matrix<U, M, N> &m)
58 {
59         std::copy(m.data, m.data+M*N, data);
60 }
61
62 template<typename T, unsigned M, unsigned N>
63 inline Matrix<T, M, N> Matrix<T, M, N>::from_columns(const Vector<T, M> *v)
64 {
65         Matrix<T, M, N> m;
66         for(unsigned i=0; i<M; ++i)
67                 for(unsigned j=0; j<N; ++j)
68                         m(i, j) = v[j][i];
69 }
70
71 template<typename T, unsigned M, unsigned N>
72 inline Matrix<T, M, N> Matrix<T, M, N>::from_rows(const Vector<T, N> *v)
73 {
74         Matrix<T, M, N> m;
75         for(unsigned i=0; i<M; ++i)
76                 for(unsigned j=0; j<N; ++j)
77                         m(i, j) = v[i][j];
78 }
79
80 template<typename T, unsigned M, unsigned N>
81 inline T &Matrix<T, M, N>::element(unsigned i, unsigned j)
82 {
83         return data[i+M*j];
84 }
85
86 template<typename T, unsigned M, unsigned N>
87 inline const T &Matrix<T, M, N>::element(unsigned i, unsigned j) const
88 {
89         return data[i+M*j];
90 }
91
92 template<typename T, unsigned M, unsigned N>
93 inline T &Matrix<T, M, N>::operator()(unsigned i, unsigned j)
94 {
95         return element(i, j);
96 }
97
98 template<typename T, unsigned M, unsigned N>
99 inline const T &Matrix<T, M, N>::operator()(unsigned i, unsigned j) const
100 {
101         return element(i, j);
102 }
103
104 template<typename T, unsigned M, unsigned N>
105 inline Matrix<T, M, N> &Matrix<T, M, N>::operator*=(T s)
106 {
107         for(unsigned i=0; i<M*N; ++i)
108                 data[i] *= s;
109         return *this;
110 }
111
112 template<typename T, unsigned M, unsigned N>
113 inline Matrix<T, M, N> operator*(const Matrix<T, M, N> &m, T s)
114 {
115         Matrix<T, M, N> r(m);
116         return r *= s;
117 }
118
119 template<typename T, unsigned M, unsigned N>
120 inline Matrix<T, M, N> operator*(T s, const Matrix<T, M, N> &m)
121 {
122         return m*s;
123 }
124
125 template<typename T, unsigned M, unsigned P, unsigned N>
126 inline Matrix<T, M, N> operator*(const Matrix<T, M, P> &m1, const Matrix<T, P, N> &m2)
127 {
128         Matrix<T, M, N> r;
129         for(unsigned i=0; i<M; ++i)
130                 for(unsigned j=0; j<N; ++j)
131                         for(unsigned k=0; k<P; ++k)
132                                 r(i, j) += m1(i, k)*m2(k, j);
133         return r;
134 }
135
136 template<typename T, unsigned M, unsigned N>
137 inline Vector<T, M> operator*(const Matrix<T, M, N> &m, const Vector<T, N> &v)
138 {
139         Vector<T, M> r;
140         for(unsigned i=0; i<M; ++i)
141                 for(unsigned j=0; j<N; ++j)
142                         r[i] += m(i, j)*v[j];
143         return r;
144 }
145
146 template<typename T, unsigned M, unsigned N>
147 inline Vector<T, N> operator*(const Vector<T, M> &v, const Matrix<T, M, N> &m)
148 {
149         Vector<T, N> r;
150         for(unsigned j=0; j<N; ++j)
151                 for(unsigned i=0; i<M; ++i)
152                         r[j] += v[i]*m(i, j);
153         return r;
154 }
155
156 template<typename T, unsigned M, unsigned N>
157 inline Matrix<T, M, N> &Matrix<T, M, N>::operator/=(T s)
158 {
159         for(unsigned i=0; i<M*N; ++i)
160                 data[i] /= s;
161         return *this;
162 }
163
164 template<typename T, unsigned M, unsigned N>
165 inline Matrix<T, M, N> operator/(const Matrix<T, M, N> &m, T s)
166 {
167         Matrix<T, M, N> r(m);
168         return r /= s;
169 }
170
171 template<typename T, unsigned M, unsigned N>
172 inline Matrix<T, M, N> &Matrix<T, M, N>::operator+=(const Matrix<T, M, N> &m)
173 {
174         for(unsigned i=0; i<M*N; ++i)
175                 data[i] += m.data[i];
176         return *this;
177 }
178
179 template<typename T, unsigned M, unsigned N>
180 inline Matrix<T, M, N> operator+(const Matrix<T, M, N> &m1, const Matrix<T, M, N> &m2)
181 {
182         Matrix<T, M, N> r(m1);
183         return r += m2;
184 }
185
186 template<typename T, unsigned M, unsigned N>
187 inline Matrix<T, M, N> &Matrix<T, M, N>::operator-=(const Matrix<T, M, N> &m)
188 {
189         for(unsigned i=0; i<M*N; ++i)
190                 data[i] -= m.data[i];
191         return *this;
192 }
193
194 template<typename T, unsigned M, unsigned N>
195 inline Matrix<T, M, N> operator-(const Matrix<T, M, N> &m1, const Matrix<T, M, N> &m2)
196 {
197         Matrix<T, M, N> r(m1);
198         return r -= m2;
199 }
200
201 template<typename T, unsigned M, unsigned N>
202 inline bool operator==(const Matrix<T, M, N> &a, const Matrix<T, M, N> &b)
203 {
204         for(unsigned j=0; j<N; ++j)
205                 for(unsigned i=0; i<M; ++i)
206                         if(a(i, j)!=b(i, j))
207                                 return false;
208         return true;
209 }
210
211 template<typename T, unsigned M, unsigned N>
212 inline Matrix<T, M, N> &Matrix<T, M, N>::exchange_rows(unsigned i, unsigned j)
213 {
214         for(unsigned k=0; k<N; ++k)
215                 std::swap(element(i, k), element(j, k));
216         return *this;
217 }
218
219 template<typename T, unsigned M, unsigned N>
220 inline Matrix<T, M, N> &Matrix<T, M, N>::multiply_row(unsigned i, T s)
221 {
222         for(unsigned k=0; k<N; ++k)
223                 element(i, k) *= s;
224         return *this;
225 }
226
227 template<typename T, unsigned M, unsigned N>
228 inline Matrix<T, M, N> &Matrix<T, M, N>::add_row(unsigned i, unsigned j, T s)
229 {
230         for(unsigned k=0; k<N; ++k)
231                 element(j, k) += element(i, k)*s;
232         return *this;
233 }
234
235 template<typename T, unsigned M, unsigned N>
236 inline Matrix<T, N, M> transpose(const Matrix<T, M, N> &m)
237 {
238         Matrix<T, N, M> r;
239         for(unsigned j=0; j<N; ++j)
240                 for(unsigned i=0; i<M; ++i)
241                         r(j, i) = m(i, j);
242         return r;
243 }
244
245 } // namespace LinAl
246 } // namespace Msp
247
248 #endif