source/linal/dynamicmatrix.h

   1 #ifndef MSP_LINAL_DYNAMICMATRIX_H_
   2 #define MSP_LINAL_DYNAMICMATRIX_H_
   3
   4 #include <algorithm>
   5 #include "dynamicvector.h"
   6 #include "matrixops.h"
   7
   8 namespace Msp {
   9 namespace LinAl {
  10
  11 /**
  12 A general mathematical matrix.  Size is specified at runtime.  May be slower
  13 than a fixed-size matrix.  There are no compile-time diagnostics for mismatched
  14 matrix sizes.
  15 */
  16 template<typename T>
  17 class DynamicMatrix
  18 {
  19 public:
  20         typedef T ElementType;
  21
  22 private:
  23         unsigned rows_;
  24         unsigned columns_;
  25         T *data;
  26
  27 public:
  28         DynamicMatrix(unsigned, unsigned);
  29         DynamicMatrix(unsigned, unsigned, const T *);
  30         DynamicMatrix(const DynamicMatrix &);
  31         DynamicMatrix &operator=(const DynamicMatrix &);
  32         ~DynamicMatrix();
  33
  34         unsigned rows() const { return rows_; }
  35         unsigned columns() const { return columns_; }
  36
  37         T &element(unsigned, unsigned);
  38         const T &element(unsigned, unsigned) const;
  39         T &operator()(unsigned i, unsigned j) { return element(i, j); }
  40         const T &operator()(unsigned i, unsigned j) const { return element(i, j); }
  41
  42         DynamicMatrix &operator*=(T);
  43         DynamicMatrix &operator/=(T);
  44         DynamicMatrix &operator+=(const DynamicMatrix &);
  45         DynamicMatrix &operator-=(const DynamicMatrix &);
  46
  47         DynamicMatrix &exchange_rows(unsigned, unsigned);
  48         DynamicMatrix &multiply_row(unsigned, T);
  49         DynamicMatrix &add_row(unsigned, unsigned, T);
  50
  51         DynamicMatrix &invert();
  52 };
  53
  54
  55 template<typename T>
  56 inline DynamicMatrix<T>::DynamicMatrix(unsigned r, unsigned c):
  57         rows_(r),
  58         columns_(c),
  59         data(new T[rows_*columns_])
  60 {
  61         std::fill(data, data+rows_*columns_, T());
  62 }
  63
  64 template<typename T>
  65 inline DynamicMatrix<T>::DynamicMatrix(unsigned r, unsigned c, const T *d):
  66         rows_(r),
  67         columns_(c),
  68         data(new T[rows_*columns_])
  69 {
  70         std::copy(d, d+rows_*columns_, data);
  71 }
  72
  73 template<typename T>
  74 inline DynamicMatrix<T>::DynamicMatrix(const DynamicMatrix &other):
  75         rows_(other.rows()),
  76         columns_(other.columns()),
  77         data(new T[rows_*columns_])
  78 {
  79         std::copy(other.data, other.data+rows_*columns_, data);
  80 }
  81
  82 template<typename T>
  83 inline DynamicMatrix<T> &DynamicMatrix<T>::operator=(const DynamicMatrix &other)
  84 {
  85         if(rows_!=other.rows() || columns_!=other.columns())
  86         {
  87                 delete[] data;
  88                 rows_ = other.rows();
  89                 columns_ = other.columns();
  90                 data = new T[rows_*columns_];
  91         }
  92
  93         std::copy(other.data, other.data+rows_*columns_, data);
  94
  95         return *this;
  96 }
  97
  98 template<typename T>
  99 inline DynamicMatrix<T>::~DynamicMatrix()
 100 {
 101         delete[] data;
 102 }
 103
 104 template<typename T>
 105 inline T &DynamicMatrix<T>::element(unsigned i, unsigned j)
 106 {
 107         if(i>=rows_ || j>=columns_)
 108                 throw std::out_of_range("DynamicMatrix::element");
 109
 110         return data[i+rows_*j];
 111 }
 112
 113 template<typename T>
 114 inline const T &DynamicMatrix<T>::element(unsigned i, unsigned j) const
 115 {
 116         if(i>=rows_ || j>=columns_)
 117                 throw std::out_of_range("DynamicMatrix::element");
 118
 119         return data[i+rows_*j];
 120 }
 121
 122 template<typename T>
 123 inline DynamicMatrix<T> &DynamicMatrix<T>::operator*=(T s)
 124 {
 125         for(unsigned i=0; i<rows_*columns_; ++i)
 126                 data[i] *= s;
 127         return *this;
 128 }
 129
 130 template<typename T>
 131 inline DynamicMatrix<T> operator*(const DynamicMatrix<T> &m, T s)
 132 {
 133         DynamicMatrix<T> r(m);
 134         return r *= s;
 135 }
 136
 137 template<typename T>
 138 inline DynamicMatrix<T> operator*(T s, const DynamicMatrix<T> &m)
 139 {
 140         return m*s;
 141 }
 142
 143 template<typename T>
 144 inline DynamicMatrix<T> operator*(const DynamicMatrix<T> &m1, const DynamicMatrix<T> &m2)
 145 {
 146         if(m1.columns()!=m2.rows())
 147                 throw size_mismatch("matrix*matrix");
 148
 149         DynamicMatrix<T> r(m1.rows(), m2.columns());
 150         for(unsigned i=0; i<m1.rows(); ++i)
 151                 for(unsigned j=0; j<m2.columns(); ++j)
 152                         for(unsigned k=0; k<m1.columns(); ++k)
 153                                 r(i, j) += m1(i, k)*m2(k, j);
 154
 155         return r;
 156 }
 157
 158 template<typename T>
 159 inline DynamicVector<T> operator*(const DynamicMatrix<T> &m, const DynamicVector<T> &v)
 160 {
 161         if(m.columns()!=v.size())
 162                 throw size_mismatch("matrix*vector");
 163
 164         DynamicVector<T> r(m.rows());
 165         for(unsigned i=0; i<m.rows(); ++i)
 166                 for(unsigned j=0; j<m.columns(); ++j)
 167                         r[i] += m(i, j)*v[j];
 168
 169         return r;
 170 }
 171
 172 template<typename T>
 173 inline DynamicVector<T> operator*(const DynamicVector<T> &v, const DynamicMatrix<T> &m)
 174 {
 175         if(v.size()!=m.rows())
 176                 throw size_mismatch("vector*matrix");
 177
 178         DynamicVector<T> r(m.columns());
 179         for(unsigned j=0; j<m.columns(); ++j)
 180                 for(unsigned i=0; i<m.rows(); ++i)
 181                         r[i] += v[i]*m(i, j);
 182
 183         return r;
 184 }
 185
 186 template<typename T>
 187 inline DynamicMatrix<T> &DynamicMatrix<T>::operator/=(T s)
 188 {
 189         for(unsigned i=0; i<rows_*columns_; ++i)
 190                 data[i] /= s;
 191         return *this;
 192 }
 193
 194 template<typename T>
 195 inline DynamicMatrix<T> operator/(const DynamicMatrix<T> &m, T s)
 196 {
 197         DynamicMatrix<T> r(m);
 198         return r /= s;
 199 }
 200
 201 template<typename T>
 202 inline DynamicMatrix<T> &DynamicMatrix<T>::operator+=(const DynamicMatrix<T> &m)
 203 {
 204         if(rows_!=m.rows_ || columns_!=m.columns_)
 205                 throw size_mismatch("matrix+matrix");
 206
 207         for(unsigned i=0; i<rows_*columns_; ++i)
 208                 data[i] += m.data[i];
 209
 210         return *this;
 211 }
 212
 213 template<typename T>
 214 inline DynamicMatrix<T> operator+(const DynamicMatrix<T> &m1, const DynamicMatrix<T> &m2)
 215 {
 216         DynamicMatrix<T> r(m1);
 217         return r += m2;
 218 }
 219
 220 template<typename T>
 221 inline DynamicMatrix<T> &DynamicMatrix<T>::operator-=(const DynamicMatrix<T> &m)
 222 {
 223         if(rows_!=m.rows_ || columns_!=m.columns_)
 224                 throw size_mismatch("matrix-matrix");
 225
 226         for(unsigned i=0; i<rows_*columns_; ++i)
 227                 data[i] += m.data[i];
 228
 229         return *this;
 230 }
 231
 232 template<typename T>
 233 inline DynamicMatrix<T> operator-(const DynamicMatrix<T> &m1, const DynamicMatrix<T> &m2)
 234 {
 235         DynamicMatrix<T> r(m1);
 236         return r += m2;
 237 }
 238
 239 template<typename T>
 240 inline bool operator==(const DynamicMatrix<T> &m1, const DynamicMatrix<T> &m2)
 241 {
 242         if(m1.rows()!=m2.rows() || m1.columns()!=m2.columns())
 243                 throw size_mismatch("matrix==matrix");
 244
 245         for(unsigned i=0; i<m1.rows(); ++i)
 246                 for(unsigned j=0; j<m1.columns(); ++j)
 247                         if(m1(i, j)!=m2(i, j))
 248                                 return false;
 249
 250         return true;
 251 }
 252
 253 template<typename T>
 254 inline DynamicMatrix<T> &DynamicMatrix<T>::exchange_rows(unsigned i, unsigned j)
 255 {
 256         if(i>=rows_ || j>=rows_)
 257                 throw std::out_of_range("DynamicMatrix::exchange_rows");
 258
 259         using std::swap;
 260         for(unsigned k=0; k<columns_; ++k)
 261                 swap(element(i, k), element(j, k));
 262
 263         return *this;
 264 }
 265
 266 template<typename T>
 267 inline DynamicMatrix<T> &DynamicMatrix<T>::multiply_row(unsigned i, T s)
 268 {
 269         if(i>=rows_)
 270                 throw std::out_of_range("DynamicMatrix::multiply_row");
 271
 272         for(unsigned k=0; k<columns_; ++k)
 273                 element(i, k) *= s;
 274
 275         return *this;
 276 }
 277
 278 template<typename T>
 279 inline DynamicMatrix<T> &DynamicMatrix<T>::add_row(unsigned i, unsigned j, T s)
 280 {
 281         if(i>=rows_ || j>=rows_)
 282                 throw std::out_of_range("DynamicMatrix::exchange_rows");
 283
 284         for(unsigned k=0; k<columns_; ++k)
 285                 element(j, k) += element(i, k)*s;
 286
 287         return *this;
 288 }
 289
 290 template<typename T>
 291 inline DynamicMatrix<T> transpose(const DynamicMatrix<T> &m)
 292 {
 293         DynamicMatrix<T> r(m.columns(), m.rows());
 294         for(unsigned i=0; i<m.rows(); ++i)
 295                 for(unsigned j=0; j<m.columns(); ++j)
 296                         r(j, i) = m(i, j);
 297         return r;
 298 }
 299
 300 template<typename T>
 301 inline DynamicMatrix<T> &DynamicMatrix<T>::invert()
 302 {
 303         if(rows_!=columns_)
 304                 throw size_mismatch("matrix invert");
 305
 306         DynamicMatrix<T> r(rows_, columns_);
 307         for(unsigned i=0; i<rows_; ++i)
 308                 r(i, i) = T(1);
 309
 310         return gauss_jordan(*this, r);
 311 }
 312
 313 template<typename T>
 314 inline DynamicMatrix<T> invert(const DynamicMatrix<T> &m)
 315 {
 316         DynamicMatrix<T> r = m;
 317         return r.invert();
 318 }
 319
 320 } // namespace LinAl
 321 } // namespace Msp
 322
 323 #endif