source/linal/matrix.h

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