1 #ifndef RIVET_MATH_MATRIXN 2 #define RIVET_MATH_MATRIXN 4 #include "Rivet/Math/MathHeader.hh" 5 #include "Rivet/Math/MathUtils.hh" 6 #include "Rivet/Math/Vectors.hh" 8 #include "Rivet/Math/eigen/matrix.h" 50 for (
size_t i = 0; i < N; ++i) {
51 rtn.set(i, i, diag[i]);
58 for (
size_t i = 0; i < N; ++i) {
72 _matrix = other._matrix;
75 Matrix&
set(
const size_t i,
const size_t j,
const double value) {
77 _matrix(i, j) = value;
79 throw std::runtime_error(
"Attempted set access outside matrix bounds.");
84 double get(
const size_t i,
const size_t j)
const {
88 throw std::runtime_error(
"Attempted get access outside matrix bounds.");
92 Vector<N> getRow(
const size_t row)
const {
94 for (
size_t i = 0; i < N; ++i) {
95 rtn.
set(i, _matrix(row,i));
101 for (
size_t i = 0; i < N; ++i) {
102 _matrix(row,i) = r.get(i);
107 Vector<N> getColumn(
const size_t col)
const {
109 for (
size_t i = 0; i < N; ++i) {
110 rtn.
set(i, _matrix(i,col));
116 for (
size_t i = 0; i < N; ++i) {
117 _matrix(i,col) = c.get(i);
124 tmp._matrix.replaceWithAdjoint();
136 tmp._matrix = _matrix.
inverse();
142 return _matrix.determinant();
148 for (
size_t i = 0; i < N; ++i) {
158 rtn._matrix = -_matrix;
168 bool isZero(
double tolerance=1E-5)
const {
169 for (
size_t i=0; i < N; ++i) {
170 for (
size_t j=0; j < N; ++j) {
179 for (
size_t i=0; i < N; ++i) {
180 for (
size_t j=i; j < N; ++j) {
181 if (!
Rivet::isZero(_matrix(i,j) - other._matrix(i,j)) )
return false;
189 return isEqual(this->transpose());
194 for (
size_t i=0; i < N; ++i) {
195 for (
size_t j=0; j < N; ++j) {
196 if (i == j)
continue;
203 bool operator==(
const Matrix<N>& a)
const {
204 return _matrix == a._matrix;
207 bool operator!=(
const Matrix<N>& a)
const {
208 return _matrix != a._matrix;
211 bool operator<(const Matrix<N>& a)
const {
212 return _matrix < a._matrix;
215 bool operator<=(const Matrix<N>& a)
const {
216 return _matrix <= a._matrix;
219 bool operator>(
const Matrix<N>& a)
const {
220 return _matrix > a._matrix;
223 bool operator>=(
const Matrix<N>& a)
const {
224 return _matrix >= a._matrix;
228 _matrix = _matrix * m._matrix;
243 _matrix += m._matrix;
248 _matrix -= m._matrix;
253 typedef Eigen::Matrix<double,N> EMatrix;
264 result._matrix = a._matrix + b._matrix;
280 return subtract(a, b);
286 rtn._matrix = a * m._matrix;
292 return multiply(a, m);
297 return multiply(1/a, m);
302 return multiply(a, m);
307 return multiply(a, m);
313 tmp._matrix = a._matrix * b._matrix;
319 return multiply(a, b);
326 tmp._vec = a._matrix * b._vec;
332 return multiply(a, b);
344 return m.transpose();
354 return m.determinant();
371 for (
size_t i = 0; i < m.
size(); ++i) {
373 for (
size_t j = 0; j < m.
size(); ++j) {
374 const double e = m.get(i, j);
386 inline ostream& operator<<(std::ostream& out, const Matrix<N>& m) {
398 for (
size_t i = 0; i < N; ++i) {
399 for (
size_t j = 0; j < N; ++j) {
400 const double a = ma.get(i, j);
401 const double b = mb.get(i, j);
412 return m.
isZero(tolerance);
Definition: MC_JetAnalysis.hh:9
Matrix< N > inverse() const
Calculate inverse.
Definition: MatrixN.hh:134
bool isSymm() const
Check for symmetry under transposition.
Definition: MatrixN.hh:188
double det() const
Calculate determinant.
Definition: MatrixN.hh:141
Vector< N > & set(const size_t index, const double value)
Set indexed value.
Definition: VectorN.hh:52
bool isZero(double tolerance=1E-5) const
Index-wise check for nullness, allowing for numerical precision.
Definition: MatrixN.hh:168
General -dimensional mathematical matrix object.
Definition: MatrixN.hh:14
bool isZero(double val, double tolerance=1E-8)
Definition: MathUtils.hh:17
A minimal base class for -dimensional vectors.
Definition: VectorN.hh:13
double trace() const
Calculate trace.
Definition: MatrixN.hh:146
Matrix< N > operator-() const
Negate.
Definition: MatrixN.hh:156
bool isEqual(Matrix< N > other) const
Check for index-wise equality, allowing for numerical precision.
Definition: MatrixN.hh:178
std::string toString(const AnalysisInfo &ai)
String representation.
Definition: AnalysisInfo.cc:234
size_t size() const
Get dimensionality.
Definition: MatrixN.hh:163
bool fuzzyEquals(double a, double b, double tolerance=1E-5)
Compare two floating point numbers for equality with a degree of fuzziness.
Definition: MathUtils.hh:34
bool isDiag() const
Check that all off-diagonal elements are zero, allowing for numerical precision.
Definition: MatrixN.hh:193