rivet is hosted by Hepforge, IPPP Durham
Rivet 4.0.0
MatrixDiag.hh
1#ifndef RIVET_MATH_MATRIXDIAG
2#define RIVET_MATH_MATRIXDIAG
3
4#include "Rivet/Math/MathConstants.hh"
5#include "Rivet/Math/MatrixN.hh"
6
7// #include "gsl/gsl_vector.h"
8// #include "gsl/gsl_matrix.h"
9// #include "gsl/gsl_eigen.h"
10
11namespace Rivet {
12
13
14template <size_t N>
15class EigenSystem;
16template <size_t N>
17EigenSystem<N> diagonalize(const Matrix<N>& m);
18
19
21template <size_t N>
22class EigenSystem {
23 template <size_t M>
24 friend EigenSystem<M> diagonalize(const Matrix<M>&);
25
26public:
27
28 typedef pair<double, Vector<N> > EigenPair;
29 typedef vector<EigenPair> EigenPairs;
30
31 Vector<N> getDiagVector() const {
32 assert(_eigenPairs.size() == N);
33 Vector<N> ret;
34 for (size_t i = 0; i < N; ++i) {
35 ret.set(i, _eigenPairs[i].first);
36 }
37 return ret;
38 }
39
40 Matrix<N> getDiagMatrix() const {
41 return Matrix<N>::mkDiag(getDiagVector());
42 }
43
44 EigenPairs getEigenPairs() const {
45 return _eigenPairs;
46 }
47
48 vector<double> getEigenValues() const {
49 assert(_eigenPairs.size() == N);
50 vector<double> ret;
51 for (size_t i = 0; i < N; ++i) {
52 ret.push_back(_eigenPairs[i].first);
53 }
54 return ret;
55 }
56
57 vector<Vector<N> > getEigenVectors() const {
58 assert(_eigenPairs.size() == N);
59 vector<Vector<N> > ret;
60 for (size_t i = 0; i < N; ++i) {
61 ret.push_back(_eigenPairs[i].second);
62 }
63 return ret;
64 }
65
66private:
67
68 EigenPairs _eigenPairs;
69
70};
71
72
74template <size_t N>
75struct EigenPairCmp :
76 public std::binary_function<const typename EigenSystem<N>::EigenPair&,
77 const typename EigenSystem<N>::EigenPair&, bool> {
78 bool operator()(const typename EigenSystem<N>::EigenPair& a,
79 const typename EigenSystem<N>::EigenPair& b) {
80 return a.first < b.first;
81 }
82};
83
84
85// /// Diagonalize an NxN matrix, returning a collection of pairs of
86// /// eigenvalues and eigenvectors, ordered decreasing in eigenvalue.
87// template <size_t N>
88// EigenSystem<N> diagonalize(const Matrix<N>& m) {
89// EigenSystem<N> esys;
90
91// // Make a GSL matrix.
92// gsl_matrix* A = gsl_matrix_alloc(N, N);
93// for (size_t i = 0; i < N; ++i) {
94// for (size_t j = 0; j < N; ++j) {
95// gsl_matrix_set(A, i, j, m.get(i, j));
96// }
97// }
98
99// // Use GSL diagonalization.
100// gsl_matrix* vecs = gsl_matrix_alloc(N, N);
101// gsl_vector* vals = gsl_vector_alloc(N);
102// gsl_eigen_symmv_workspace* workspace = gsl_eigen_symmv_alloc(N);
103// gsl_eigen_symmv(A, vals, vecs, workspace);
104// gsl_eigen_symmv_sort(vals, vecs, GSL_EIGEN_SORT_VAL_DESC);
105
106// // Build the vector of "eigen-pairs".
107// typename EigenSystem<N>::EigenPairs eigensolns;
108// for (size_t i = 0; i < N; ++i) {
109// typename EigenSystem<N>::EigenPair ep;
110// ep.first = gsl_vector_get(vals, i);
111// Vector<N> ev;
112// for (size_t j = 0; j < N; ++j) {
113// ev.set(j, gsl_matrix_get(vecs, j, i));
114// }
115// ep.second = ev;
116// eigensolns.push_back(ep);
117// }
118
119// // Free GSL memory.
120// gsl_eigen_symmv_free(workspace);
121// gsl_matrix_free(A);
122// gsl_matrix_free(vecs);
123// gsl_vector_free(vals);
124
125// // Populate the returned object.
126// esys._eigenPairs = eigensolns;
127// return esys;
128// }
129
130
131template <size_t N>
132inline const string toString(const typename EigenSystem<N>::EigenPair& e) {
133 ostringstream ss;
134 //for (typename EigenSystem<N>::EigenPairs::const_iterator i = e.begin(); i != e.end(); ++i) {
135 ss << e->first << " -> " << e->second;
136 // if (i+1 != e.end()) ss << endl;
137 //}
138 return ss.str();
139}
140
141template <size_t N>
142inline ostream& operator<<(std::ostream& out, const typename EigenSystem<N>::EigenPair& e) {
143 out << toString(e);
144 return out;
145}
146
147
148}
149
150#endif
Definition MC_CENT_PPB_Projections.hh:10
std::ostream & operator<<(std::ostream &os, const AnalysisInfo &ai)
Stream an AnalysisInfo as a text description.
Definition AnalysisInfo.hh:362
std::string toString(const AnalysisInfo &ai)
String representation.