rivet is hosted by Hepforge, IPPP Durham
 The Rivet MC analysis system  2.5.2
MatrixN.hh
Go to the documentation of this file.
```00001 #ifndef RIVET_MATH_MATRIXN
00002 #define RIVET_MATH_MATRIXN
00003
00005 #include "Rivet/Math/MathUtils.hh"
00006 #include "Rivet/Math/Vectors.hh"
00007
00008 #include "Rivet/Math/eigen/matrix.h"
00009
00010 namespace Rivet {
00011
00012
00013   template <size_t N>
00014   class Matrix;
00015   typedef Matrix<4> Matrix4;
00016
00017   template <size_t N>
00018   Matrix<N> multiply(const Matrix<N>& a, const Matrix<N>& b);
00019   template <size_t N>
00020   Matrix<N> divide(const Matrix<N>&, const double);
00021   template <size_t N>
00022   Matrix<N> operator*(const Matrix<N>& a, const Matrix<N>& b);
00023
00024
00025   ///////////////////////////////////
00026
00027
00028   /// @brief General \f\$ N \f\$-dimensional mathematical matrix object.
00029   template <size_t N>
00030   class Matrix {
00031     template <size_t M>
00032     friend Matrix<M> add(const Matrix<M>&, const Matrix<M>&);
00033     template <size_t M>
00034     friend Matrix<M> multiply(const double, const Matrix<M>&);
00035     template <size_t M>
00036     friend Matrix<M> multiply(const Matrix<M>&, const Matrix<M>&);
00037     template <size_t M>
00038     friend Vector<M> multiply(const Matrix<M>&, const Vector<M>&);
00039     template <size_t M>
00040     friend Matrix<M> divide(const Matrix<M>&, const double);
00041
00042   public:
00043     static Matrix<N> mkZero() {
00044       Matrix<N> rtn;
00045       return rtn;
00046     }
00047
00048     static Matrix<N> mkDiag(Vector<N> diag) {
00049       Matrix<N> rtn;
00050       for (size_t i = 0; i < N; ++i) {
00051         rtn.set(i, i, diag[i]);
00052       }
00053       return rtn;
00054     }
00055
00056     static Matrix<N> mkIdentity() {
00057       Matrix<N> rtn;
00058       for (size_t i = 0; i < N; ++i) {
00059         rtn.set(i, i, 1);
00060       }
00061       return rtn;
00062     }
00063
00064
00065   public:
00066
00067     Matrix() {
00069     }
00070
00071     Matrix(const Matrix<N>& other) {
00072       _matrix = other._matrix;
00073     }
00074
00075     Matrix& set(const size_t i, const size_t j, const double value) {
00076       if (i < N && j < N) {
00077         _matrix(i, j) = value;
00078       } else {
00079         throw std::runtime_error("Attempted set access outside matrix bounds.");
00080       }
00081       return *this;
00082     }
00083
00084     double get(const size_t i, const size_t j) const {
00085       if (i < N && j < N) {
00086         return _matrix(i, j);
00087       } else {
00088         throw std::runtime_error("Attempted get access outside matrix bounds.");
00089       }
00090     }
00091
00092     Vector<N> getRow(const size_t row) const {
00093       Vector<N> rtn;
00094       for (size_t i = 0; i < N; ++i) {
00095         rtn.set(i, _matrix(row,i));
00096       }
00097       return rtn;
00098     }
00099
00100     Matrix<N>& setRow(const size_t row, const Vector<N>& r) {
00101       for (size_t i = 0; i < N; ++i) {
00102         _matrix(row,i) = r.get(i);
00103       }
00104       return *this;
00105     }
00106
00107     Vector<N> getColumn(const size_t col) const {
00108       Vector<N> rtn;
00109       for (size_t i = 0; i < N; ++i) {
00110         rtn.set(i, _matrix(i,col));
00111       }
00112       return rtn;
00113     }
00114
00115     Matrix<N>& setColumn(const size_t col, const Vector<N>& c) {
00116       for (size_t i = 0; i < N; ++i) {
00117         _matrix(i,col) = c.get(i);
00118       }
00119       return *this;
00120     }
00121
00122     Matrix<N> transpose() const {
00123       Matrix<N> tmp = *this;
00125       return tmp;
00126     }
00127
00128     // Matrix<N>& transposeInPlace() {
00130     //   return *this;
00131     // }
00132
00133     /// Calculate inverse
00134     Matrix<N> inverse() const {
00135       Matrix<N> tmp;
00136       tmp._matrix = _matrix.inverse();
00137       return tmp;
00138     }
00139
00140     /// Calculate determinant
00141     double det() const  {
00142       return _matrix.determinant();
00143     }
00144
00145     /// Calculate trace
00146     double trace() const {
00147       double tr = 0.0;
00148       for (size_t i = 0; i < N; ++i) {
00149         tr += _matrix(i,i);
00150       }
00151       return tr;
00152       // return _matrix.trace();
00153     }
00154
00155     /// Negate
00156     Matrix<N> operator-() const {
00157       Matrix<N> rtn;
00158       rtn._matrix = -_matrix;
00159       return rtn;
00160     }
00161
00162     /// Get dimensionality
00163     size_t size() const {
00164       return N;
00165     }
00166
00167     /// Index-wise check for nullness, allowing for numerical precision
00168     bool isZero(double tolerance=1E-5) const {
00169       for (size_t i=0; i < N; ++i) {
00170         for (size_t j=0; j < N; ++j) {
00171           if (! Rivet::isZero(_matrix(i,j), tolerance) ) return false;
00172         }
00173       }
00174       return true;
00175     }
00176
00177     /// Check for index-wise equality, allowing for numerical precision
00178     bool isEqual(Matrix<N> other) const {
00179       for (size_t i=0; i < N; ++i) {
00180         for (size_t j=i; j < N; ++j) {
00181           if (! Rivet::isZero(_matrix(i,j) - other._matrix(i,j)) ) return false;
00182         }
00183       }
00184       return true;
00185     }
00186
00187     /// Check for symmetry under transposition
00188     bool isSymm() const {
00189       return isEqual(this->transpose());
00190     }
00191
00192     /// Check that all off-diagonal elements are zero, allowing for numerical precision
00193     bool isDiag() const {
00194       for (size_t i=0; i < N; ++i) {
00195         for (size_t j=0; j < N; ++j) {
00196           if (i == j) continue;
00197           if (! Rivet::isZero(_matrix(i,j)) ) return false;
00198         }
00199       }
00200       return true;
00201     }
00202
00203     bool operator == (const Matrix<N>& a) const {
00204       return _matrix == a._matrix;
00205     }
00206
00207     bool operator != (const Matrix<N>& a) const {
00208       return _matrix != a._matrix;
00209     }
00210
00211     bool operator < (const Matrix<N>& a) const {
00212       return _matrix < a._matrix;
00213     }
00214
00215     bool operator <= (const Matrix<N>& a) const {
00216       return _matrix <= a._matrix;
00217     }
00218
00219     bool operator > (const Matrix<N>& a) const {
00220       return _matrix > a._matrix;
00221     }
00222
00223     bool operator >= (const Matrix<N>& a) const {
00224       return _matrix >= a._matrix;
00225     }
00226
00227     Matrix<N>& operator *= (const Matrix<N>& m) {
00228       _matrix = _matrix * m._matrix;
00229       return *this;
00230     }
00231
00232     Matrix<N>& operator *= (const double a) {
00233       _matrix *= a;
00234       return *this;
00235     }
00236
00237     Matrix<N>& operator /= (const double a) {
00238       _matrix /= a;
00239       return *this;
00240     }
00241
00242     Matrix<N>& operator += (const Matrix<N>& m) {
00243       _matrix += m._matrix;
00244       return *this;
00245     }
00246
00247     Matrix<N>& operator -= (const Matrix<N>& m) {
00248       _matrix -= m._matrix;
00249       return *this;
00250     }
00251
00252   protected:
00253
00254     typedef Eigen::Matrix<double,N> EMatrix;
00255     EMatrix _matrix;
00256
00257   };
00258
00259
00260   /////////////////////////////////
00261
00262
00263   template <size_t N>
00264   inline Matrix<N> add(const Matrix<N>& a, const Matrix<N>& b) {
00265     Matrix<N> result;
00266     result._matrix = a._matrix + b._matrix;
00267     return result;
00268   }
00269
00270   template <size_t N>
00271   inline Matrix<N> subtract(const Matrix<N>& a, const Matrix<N>& b) {
00273   }
00274
00275   template <size_t N>
00276   inline Matrix<N> operator + (const Matrix<N> a, const Matrix<N>& b) {
00278   }
00279
00280   template <size_t N>
00281   inline Matrix<N> operator - (const Matrix<N> a, const Matrix<N>& b) {
00282     return subtract(a, b);
00283   }
00284
00285   template <size_t N>
00286   inline Matrix<N> multiply(const double a, const Matrix<N>& m) {
00287     Matrix<N> rtn;
00288     rtn._matrix = a * m._matrix;
00289     return rtn;
00290   }
00291
00292   template <size_t N>
00293   inline Matrix<N> multiply(const Matrix<N>& m, const double a) {
00294     return multiply(a, m);
00295   }
00296
00297   template <size_t N>
00298   inline Matrix<N> divide(const Matrix<N>& m, const double a) {
00299     return multiply(1/a, m);
00300   }
00301
00302   template <size_t N>
00303   inline Matrix<N> operator * (const double a, const Matrix<N>& m) {
00304     return multiply(a, m);
00305   }
00306
00307   template <size_t N>
00308   inline Matrix<N> operator * (const Matrix<N>& m, const double a) {
00309     return multiply(a, m);
00310   }
00311
00312   template <size_t N>
00313   inline Matrix<N> multiply(const Matrix<N>& a, const Matrix<N>& b) {
00314     Matrix<N> tmp;
00315     tmp._matrix = a._matrix * b._matrix;
00316     return tmp;
00317   }
00318
00319   template <size_t N>
00320   inline Matrix<N> operator * (const Matrix<N>& a, const Matrix<N>& b) {
00321     return multiply(a, b);
00322   }
00323
00324
00325   template <size_t N>
00326   inline Vector<N> multiply(const Matrix<N>& a, const Vector<N>& b) {
00327     Vector<N> tmp;
00328     tmp._vec = a._matrix * b._vec;
00329     return tmp;
00330   }
00331
00332   template <size_t N>
00333   inline Vector<N> operator * (const Matrix<N>& a, const Vector<N>& b) {
00334     return multiply(a, b);
00335   }
00336
00337   template <size_t N>
00338   inline Matrix<N> transpose(const Matrix<N>& m) {
00339     // Matrix<N> tmp;
00340     // for (size_t i = 0; i < N; ++i) {
00341     //   for (size_t j = 0; j < N; ++j) {
00342     //     tmp.set(i, j, m.get(j, i));
00343     //   }
00344     // }
00345     // return tmp;
00346     return m.transpose();
00347   }
00348
00349   template <size_t N>
00350   inline Matrix<N> inverse(const Matrix<N>& m) {
00351     return m.inverse();
00352   }
00353
00354   template <size_t N>
00355   inline double det(const Matrix<N>& m) {
00356     return m.determinant();
00357   }
00358
00359   template <size_t N>
00360   inline double trace(const Matrix<N>& m) {
00361     return m.trace();
00362   }
00363
00364
00365   /////////////////////////////////
00366
00367
00368   /// Make string representation
00369   template <size_t N>
00370   inline string toString(const Matrix<N>& m) {
00371     ostringstream ss;
00372     ss << "[ ";
00373     for (size_t i = 0; i < m.size(); ++i) {
00374       ss << "( ";
00375       for (size_t j = 0; j < m.size(); ++j) {
00376         const double e = m.get(i, j);
00377         ss << (Rivet::isZero(e) ? 0.0 : e) << " ";
00378       }
00379       ss << ") ";
00380     }
00381     ss << "]";
00382     return ss.str();
00383   }
00384
00385
00386   /// Stream out string representation
00387   template <size_t N>
00388   inline ostream& operator << (std::ostream& out, const Matrix<N>& m) {
00389     out << toString(m);
00390     return out;
00391   }
00392
00393
00394   /////////////////////////////////////////////////
00395
00396
00397   /// Compare two matrices by index, allowing for numerical precision
00398   template <size_t N>
00399   inline bool fuzzyEquals(const Matrix<N>& ma, const Matrix<N>& mb, double tolerance=1E-5) {
00400     for (size_t i = 0; i < N; ++i) {
00401       for (size_t j = 0; j < N; ++j) {
00402         const double a = ma.get(i, j);
00403         const double b = mb.get(i, j);
00404         if (!Rivet::fuzzyEquals(a, b, tolerance)) return false;
00405       }
00406     }
00407     return true;
00408   }
00409
00410
00411   /// External form of numerically safe nullness check
00412   template <size_t N>
00413   inline bool isZero(const Matrix<N>& m, double tolerance=1E-5) {
00414     return m.isZero(tolerance);
00415   }
00416
00417
00418 }
00419
00420 #endif
```