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