source/matrix.h

   1 #ifndef MSP_LIEAL_MATMIX_H_
   2 #define MSP_LIEAL_MATMIX_H_
   3
   4 #include "vector.h"
   5
   6 namespace Msp {
   7 namespace LinAl {
   8
   9 /**
  10 A general mathematical matrix.
  11 */
  12 template<typename T, unsigned M, unsigned N>
  13 class Matrix
  14 {
  15 private:
  16         T data[M*N];
  17
  18 public:
  19         Matrix();
  20         Matrix(const T *);
  21         template<typename U>
  22         Matrix(const Matrix<U, M, N> &);
  23         static Matrix from_columns(const Vector<T, M> *);
  24         static Matrix from_rows(const Vector<T, N> *);
  25
  26         T &operator()(unsigned, unsigned);
  27         const T &operator()(unsigned, unsigned) const;
  28
  29         Matrix &operator*=(T);
  30         Matrix &operator/=(T);
  31         Matrix &operator+=(const Matrix &);
  32         Matrix &operator-=(const Matrix &);
  33
  34         Matrix<T, N, M> transpose() const;
  35 };
  36
  37 template<typename T, unsigned M, unsigned N>
  38 inline T &Matrix<T, M, N>::operator()(unsigned i, unsigned j)
  39 {
  40         return data[i+M*j];
  41 }
  42
  43 template<typename T, unsigned M, unsigned N>
  44 inline const T &Matrix<T, M, N>::operator()(unsigned i, unsigned j) const
  45 {
  46         return data[i+M*j];
  47 }
  48
  49 template<typename T, unsigned M, unsigned N>
  50 inline Matrix<T, M, N> &Matrix<T, M, N>::operator*=(T s)
  51 {
  52         for(unsigned i=0; i<M*N; ++i)
  53                 data[i] *= s;
  54         return *this;
  55 }
  56
  57 template<typename T, unsigned M, unsigned N>
  58 inline Matrix<T, M, N> operator*(const Matrix<T, M, N> &m, T s)
  59 {
  60         Matrix<T, M, N> r(m);
  61         return r *= s;
  62 }
  63
  64 template<typename T, unsigned M, unsigned N>
  65 inline Matrix<T, M, N> operator*(T s, const Matrix<T, M, N> &m)
  66 {
  67         return m*s;
  68 }
  69
  70 template<typename T, unsigned M, unsigned N>
  71 inline Matrix<T, M, N> &Matrix<T, M, N>::operator/=(T s)
  72 {
  73         for(unsigned i=0; i<M*N; ++i)
  74                 data[i] /= s;
  75         return *this;
  76 }
  77
  78 template<typename T, unsigned M, unsigned N>
  79 inline Matrix<T, M, N> operator/(const Matrix<T, M, N> &m, T s)
  80 {
  81         Matrix<T, M, N> r(m);
  82         return r /= s;
  83 }
  84
  85 template<typename T, unsigned M, unsigned N>
  86 inline Matrix<T, M, N> &Matrix<T, M, N>::operator+=(const Matrix<T, M, N> &m)
  87 {
  88         for(unsigned i=0; i<M*N; ++i)
  89                 data[i] += m.data[i];
  90         return *this;
  91 }
  92
  93 template<typename T, unsigned M, unsigned N>
  94 inline Matrix<T, M, N> operator+(const Matrix<T, M, N> &m1, const Matrix<T, M, N> &m2)
  95 {
  96         Matrix<T, M, N> r(m1);
  97         return r += m2;
  98 }
  99
 100 template<typename T, unsigned M, unsigned N>
 101 inline Matrix<T, M, N> &Matrix<T, M, N>::operator-=(const Matrix<T, M, N> &m)
 102 {
 103         for(unsigned i=0; i<M*N; ++i)
 104                 data[i] -= m.data[i];
 105         return *this;
 106 }
 107
 108 template<typename T, unsigned M, unsigned N>
 109 inline Matrix<T, M, N> operator-(const Matrix<T, M, N> &m1, const Matrix<T, M, N> &m2)
 110 {
 111         Matrix<T, M, N> r(m1);
 112         return r -= m2;
 113 }
 114
 115 template<typename T, unsigned M, unsigned P, unsigned N>
 116 Matrix<T, M, N> operator*(const Matrix<T, M, P> &m1, const Matrix<T, P, N> &m2)
 117 {
 118         Matrix<T, M, N> r;
 119         for(unsigned i=0; i<M; ++i)
 120                 for(unsigned j=0; j<N; ++j)
 121                         for(unsigned k=0; k<P; ++k)
 122                                 r.at(i, j) += m1(i, k)*m2(k, j);
 123         return r;
 124 }
 125
 126 template<typename T, unsigned M, unsigned N>
 127 Vector<T, M> operator*(const Matrix<T, M, N> &m, const Vector<T, N> &v)
 128 {
 129         Vector<T, M> r;
 130         for(unsigned i=0; i<M; ++i)
 131                 for(unsigned j=0; j<N; ++j)
 132                         r[i] += m(i, j)*v[j];
 133         return r;
 134 }
 135
 136 template<typename T, unsigned M, unsigned N>
 137 Vector<T, N> operator*(const Vector<T, M> &v, const Matrix<T, M, N> &m)
 138 {
 139         Vector<T, N> r;
 140         for(unsigned j=0; j<N; ++j)
 141                 for(unsigned i=0; i<M; ++i)
 142                         r[j] += v[i]*m.at(i, j);
 143         return r;
 144 }
 145
 146 } // namespace LinAl
 147 } // namespace Msp
 148
 149 #endif