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