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