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