1 #ifndef RIVET_MATH_MATRIXDIAG
2 #define RIVET_MATH_MATRIXDIAG
4 #include "Rivet/Math/MathHeader.hh"
5 #include "Rivet/Math/MatrixN.hh"
7 #include "gsl/gsl_vector.h"
8 #include "gsl/gsl_matrix.h"
9 #include "gsl/gsl_eigen.h"
54 typedef pair<double, Vector<N> > EigenPair;
55 typedef vector<EigenPair> EigenPairs;
58 assert(_eigenPairs.size() == N);
60 for (
size_t i = 0; i < N; ++i) {
61 ret.
set(i, _eigenPairs[i].first);
66 Matrix<N> getDiagMatrix()
const {
67 return Matrix<N>::mkDiag(getDiagVector());
70 EigenPairs getEigenPairs()
const {
74 vector<double> getEigenValues()
const {
75 assert(_eigenPairs.size() == N);
77 for (
size_t i = 0; i < N; ++i) {
78 ret.push_back(_eigenPairs[i].first);
83 vector<Vector<N> > getEigenVectors()
const {
84 assert(_eigenPairs.size() == N);
85 vector<Vector<N> > ret;
86 for (
size_t i = 0; i < N; ++i) {
87 ret.push_back(_eigenPairs[i].second);
93 EigenPairs _eigenPairs;
100 public std::binary_function<const typename EigenSystem<N>::EigenPair&,
101 const typename EigenSystem<N>::EigenPair&, bool> {
102 bool operator()(
const typename EigenSystem<N>::EigenPair& a,
103 const typename EigenSystem<N>::EigenPair& b) {
104 return a.first < b.first;
116 gsl_matrix* A = gsl_matrix_alloc(N, N);
117 for (
size_t i = 0; i < N; ++i) {
118 for (
size_t j = 0; j < N; ++j) {
119 gsl_matrix_set(A, i, j, m.get(i, j));
124 gsl_matrix* vecs = gsl_matrix_alloc(N, N);
125 gsl_vector* vals = gsl_vector_alloc(N);
126 gsl_eigen_symmv_workspace* workspace = gsl_eigen_symmv_alloc(N);
127 gsl_eigen_symmv(A, vals, vecs, workspace);
128 gsl_eigen_symmv_sort(vals, vecs, GSL_EIGEN_SORT_VAL_DESC);
131 typename EigenSystem<N>::EigenPairs eigensolns;
132 for (
size_t i = 0; i < N; ++i) {
133 typename EigenSystem<N>::EigenPair ep;
134 ep.first = gsl_vector_get(vals, i);
136 for (
size_t j = 0; j < N; ++j) {
137 ev.
set(j, gsl_matrix_get(vecs, j, i));
140 eigensolns.push_back(ep);
144 gsl_eigen_symmv_free(workspace);
146 gsl_matrix_free(vecs);
147 gsl_vector_free(vals);
150 esys._eigenPairs = eigensolns;
156 inline const string toString(
const typename EigenSystem<N>::EigenPair& e) {
159 ss << e->first <<
" -> " << e->second;
166 inline ostream& operator<<(std::ostream& out, const typename EigenSystem<N>::EigenPair& e) {
Definition: MC_JetAnalysis.hh:9
EigenSystem< N > diagonalize(const Matrix< N > &m)
Definition: MatrixDiag.hh:112
Handy object containing results of a diagonalization.
Definition: MatrixDiag.hh:42
Vector< N > & set(const size_t index, const double value)
Set indexed value.
Definition: VectorN.hh:52
General -dimensional mathematical matrix object.
Definition: MatrixN.hh:14
A minimal base class for -dimensional vectors.
Definition: VectorN.hh:13
Comparison functor for "eigen-pairs".
Definition: MatrixDiag.hh:99
std::string toString(const AnalysisInfo &ai)
String representation.
Definition: AnalysisInfo.cc:234