source/linal/dynamicmatrix.h

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