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]);
56 static Matrix<N> mkIdentity() {
58 for (
size_t i = 0; i < N; ++i) {
71 Matrix(
const Matrix<N>& other) {
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));
100 Matrix<N>& setRow(
const size_t row,
const Vector<N>& r) {
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));
115 Matrix<N>& setColumn(
const size_t col,
const Vector<N>& c) {
116 for (
size_t i = 0; i < N; ++i) {
117 _matrix(i,col) = c.get(i);
122 Matrix<N> transpose()
const {
123 Matrix<N> tmp = *
this;
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;
227 Matrix<N>& operator*=(
const Matrix<N>& m) {
228 _matrix = _matrix * m._matrix;
232 Matrix<N>& operator*=(
const double a) {
237 Matrix<N>& operator/=(
const double a) {
242 Matrix<N>& operator+=(
const Matrix<N>& m) {
243 _matrix += m._matrix;
247 Matrix<N>& operator-=(
const Matrix<N>& m) {
248 _matrix -= m._matrix;
253 typedef Eigen::Matrix<double,N> EMatrix;
262 inline Matrix<N> add(
const Matrix<N>& a,
const Matrix<N>& b) {
264 result._matrix = a._matrix + b._matrix;
269 inline Matrix<N> subtract(
const Matrix<N>& a,
const Matrix<N>& b) {
274 inline Matrix<N> operator+(
const Matrix<N> a,
const Matrix<N>& b) {
279 inline Matrix<N> operator-(
const Matrix<N> a,
const Matrix<N>& b) {
280 return subtract(a, b);
284 inline Matrix<N> multiply(
const double a,
const Matrix<N>& m) {
286 rtn._matrix = a * m._matrix;
291 inline Matrix<N> multiply(
const Matrix<N>& m,
const double a) {
292 return multiply(a, m);
296 inline Matrix<N> divide(
const Matrix<N>& m,
const double a) {
297 return multiply(1/a, m);
301 inline Matrix<N> operator*(
const double a,
const Matrix<N>& m) {
302 return multiply(a, m);
306 inline Matrix<N> operator*(
const Matrix<N>& m,
const double a) {
307 return multiply(a, m);
311 inline Matrix<N> multiply(
const Matrix<N>& a,
const Matrix<N>& b) {
313 tmp._matrix = a._matrix * b._matrix;
318 inline Matrix<N> operator*(
const Matrix<N>& a,
const Matrix<N>& b) {
319 return multiply(a, b);
324 inline Vector<N> multiply(
const Matrix<N>& a,
const Vector<N>& b) {
326 tmp._vec = a._matrix * b._vec;
331 inline Vector<N> operator*(
const Matrix<N>& a,
const Vector<N>& b) {
332 return multiply(a, b);
336 inline Matrix<N> transpose(
const Matrix<N>& m) {
344 return m.transpose();
348 inline Matrix<N> inverse(
const Matrix<N>& m) {
353 inline double det(
const Matrix<N>& m) {
354 return m.determinant();
358 inline double trace(
const Matrix<N>& m) {
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);