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         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> &other)
  58 {
  59         for(unsigned i=0; i<M; ++i)
  60                 for(unsigned j=0; j<N; ++j)
  61                         element(i, j) = other(i, j);
  62 }
  63
  64 template<typename T, unsigned M, unsigned N>
  65 inline Matrix<T, M, N> Matrix<T, M, N>::from_columns(const Vector<T, M> *v)
  66 {
  67         Matrix<T, M, N> m;
  68         for(unsigned i=0; i<M; ++i)
  69                 for(unsigned j=0; j<N; ++j)
  70                         m(i, j) = v[j][i];
  71         return m;
  72 }
  73
  74 template<typename T, unsigned M, unsigned N>
  75 inline Matrix<T, M, N> Matrix<T, M, N>::from_rows(const Vector<T, N> *v)
  76 {
  77         Matrix<T, M, N> m;
  78         for(unsigned i=0; i<M; ++i)
  79                 for(unsigned j=0; j<N; ++j)
  80                         m(i, j) = v[i][j];
  81         return m;
  82 }
  83
  84 template<typename T, unsigned M, unsigned N>
  85 inline Matrix<T, M, N> &Matrix<T, M, N>::operator*=(T s)
  86 {
  87         for(unsigned i=0; i<M*N; ++i)
  88                 data[i] *= s;
  89         return *this;
  90 }
  91
  92 template<typename T, unsigned M, unsigned N>
  93 inline Matrix<T, M, N> operator*(const Matrix<T, M, N> &m, T s)
  94 {
  95         Matrix<T, M, N> r(m);
  96         return r *= s;
  97 }
  98
  99 template<typename T, unsigned M, unsigned N>
 100 inline Matrix<T, M, N> operator*(T s, const Matrix<T, M, N> &m)
 101 {
 102         return m*s;
 103 }
 104
 105 template<typename T, unsigned M, unsigned P, unsigned N>
 106 inline Matrix<T, M, N> operator*(const Matrix<T, M, P> &m1, const Matrix<T, P, N> &m2)
 107 {
 108         Matrix<T, M, N> r;
 109         for(unsigned i=0; i<M; ++i)
 110                 for(unsigned j=0; j<N; ++j)
 111                         for(unsigned k=0; k<P; ++k)
 112                                 r(i, j) += m1(i, k)*m2(k, j);
 113         return r;
 114 }
 115
 116 template<typename T, unsigned M, unsigned N>
 117 inline Vector<T, M> operator*(const Matrix<T, M, N> &m, const Vector<T, N> &v)
 118 {
 119         Vector<T, M> r;
 120         for(unsigned i=0; i<M; ++i)
 121                 for(unsigned j=0; j<N; ++j)
 122                         r[i] += m(i, j)*v[j];
 123         return r;
 124 }
 125
 126 template<typename T, unsigned M, unsigned N>
 127 inline Vector<T, N> operator*(const Vector<T, M> &v, const Matrix<T, M, N> &m)
 128 {
 129         Vector<T, N> r;
 130         for(unsigned j=0; j<N; ++j)
 131                 for(unsigned i=0; i<M; ++i)
 132                         r[j] += v[i]*m(i, j);
 133         return r;
 134 }
 135
 136 template<typename T, unsigned M, unsigned N>
 137 inline Matrix<T, M, N> &Matrix<T, M, N>::operator/=(T s)
 138 {
 139         for(unsigned i=0; i<M*N; ++i)
 140                 data[i] /= s;
 141         return *this;
 142 }
 143
 144 template<typename T, unsigned M, unsigned N>
 145 inline Matrix<T, M, N> operator/(const Matrix<T, M, N> &m, T s)
 146 {
 147         Matrix<T, M, N> r(m);
 148         return r /= s;
 149 }
 150
 151 template<typename T, unsigned M, unsigned N>
 152 inline Matrix<T, M, N> &Matrix<T, M, N>::operator+=(const Matrix<T, M, N> &m)
 153 {
 154         for(unsigned i=0; i<M*N; ++i)
 155                 data[i] += m.data[i];
 156         return *this;
 157 }
 158
 159 template<typename T, unsigned M, unsigned N>
 160 inline Matrix<T, M, N> operator+(const Matrix<T, M, N> &m1, const Matrix<T, M, N> &m2)
 161 {
 162         Matrix<T, M, N> r(m1);
 163         return r += m2;
 164 }
 165
 166 template<typename T, unsigned M, unsigned N>
 167 inline Matrix<T, M, N> &Matrix<T, M, N>::operator-=(const Matrix<T, M, N> &m)
 168 {
 169         for(unsigned i=0; i<M*N; ++i)
 170                 data[i] -= m.data[i];
 171         return *this;
 172 }
 173
 174 template<typename T, unsigned M, unsigned N>
 175 inline Matrix<T, M, N> operator-(const Matrix<T, M, N> &m1, const Matrix<T, M, N> &m2)
 176 {
 177         Matrix<T, M, N> r(m1);
 178         return r -= m2;
 179 }
 180
 181 template<typename T, unsigned M, unsigned N>
 182 inline bool operator==(const Matrix<T, M, N> &a, const Matrix<T, M, N> &b)
 183 {
 184         for(unsigned j=0; j<N; ++j)
 185                 for(unsigned i=0; i<M; ++i)
 186                         if(a(i, j)!=b(i, j))
 187                                 return false;
 188         return true;
 189 }
 190
 191 template<typename T, unsigned M, unsigned N>
 192 inline Matrix<T, M, N> &Matrix<T, M, N>::exchange_rows(unsigned i, unsigned j)
 193 {
 194         using std::swap;
 195         for(unsigned k=0; k<N; ++k)
 196                 swap(element(i, k), element(j, k));
 197         return *this;
 198 }
 199
 200 template<typename T, unsigned M, unsigned N>
 201 inline Matrix<T, M, N> &Matrix<T, M, N>::multiply_row(unsigned i, T s)
 202 {
 203         for(unsigned k=0; k<N; ++k)
 204                 element(i, k) *= s;
 205         return *this;
 206 }
 207
 208 template<typename T, unsigned M, unsigned N>
 209 inline Matrix<T, M, N> &Matrix<T, M, N>::add_row(unsigned i, unsigned j, T s)
 210 {
 211         for(unsigned k=0; k<N; ++k)
 212                 element(j, k) += element(i, k)*s;
 213         return *this;
 214 }
 215
 216 template<typename T, unsigned M, unsigned N>
 217 inline Matrix<T, N, M> transpose(const Matrix<T, M, N> &m)
 218 {
 219         Matrix<T, N, M> r;
 220         for(unsigned j=0; j<N; ++j)
 221                 for(unsigned i=0; i<M; ++i)
 222                         r(j, i) = m(i, j);
 223         return r;
 224 }
 225
 226 } // namespace LinAl
 227 } // namespace Msp
 228
 229 #endif