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