/usr/include/Rivet/Math/MatrixDiag.hh is in librivet-dev 1.8.3-1.1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 | #ifndef RIVET_MATH_MATRIXDIAG
#define RIVET_MATH_MATRIXDIAG
#include "Rivet/Math/MathHeader.hh"
#include "Rivet/Math/MatrixN.hh"
#include "gsl/gsl_vector.h"
#include "gsl/gsl_matrix.h"
#include "gsl/gsl_eigen.h"
namespace Rivet {
// // GSL forward declarations (avoids need for GSL header files)
// extern "C" {
// struct gsl_vector;
// gsl_vector* gsl_vector_alloc(size_t);
// double gsl_vector_get(gsl_vector*, size_t);
// void gsl_vector_set(gsl_vector*, size_t, double);
// void gsl_vector_free(gsl_vector*);
// struct gsl_matrix;
// gsl_matrix* gsl_matrix_alloc(size_t, size_t);
// double gsl_matrix_get(gsl_matrix*, size_t, size_t);
// void gsl_matrix_set(gsl_matrix*, size_t, size_t, double);
// void gsl_matrix_free(gsl_matrix*);
// struct gsl_eigen_symmv_workspace;
// gsl_eigen_symmv_workspace* gsl_eigen_symmv_alloc(size_t);
// void gsl_eigen_symmv(gsl_matrix*, gsl_vector*, gsl_matrix*, gsl_eigen_symmv_workspace*);
// void gsl_eigen_symmv_free(gsl_eigen_symmv_workspace*);
// typedef enum {
// GSL_EIGEN_SORT_VAL_ASC,
// GSL_EIGEN_SORT_VAL_DESC,
// GSL_EIGEN_SORT_ABS_ASC,
// GSL_EIGEN_SORT_ABS_DESC
// }
// gsl_eigen_sort_t;
// int gsl_eigen_symmv_sort(gsl_vector * eval, gsl_matrix * evec, gsl_eigen_sort_t sort_type);
// }
template <size_t N>
class EigenSystem;
template <size_t N>
EigenSystem<N> diagonalize(const Matrix<N>& m);
/// Handy object containing results of a diagonalization.
template <size_t N>
class EigenSystem {
template <size_t M>
friend EigenSystem<M> diagonalize(const Matrix<M>&);
public:
typedef pair<double, Vector<N> > EigenPair;
typedef vector<EigenPair> EigenPairs;
Vector<N> getDiagVector() const {
assert(_eigenPairs.size() == N);
Vector<N> ret;
for (size_t i = 0; i < N; ++i) {
ret.set(i, _eigenPairs[i].first);
}
return ret;
}
Matrix<N> getDiagMatrix() const {
return Matrix<N>::mkDiag(getDiagVector());
}
EigenPairs getEigenPairs() const {
return _eigenPairs;
}
vector<double> getEigenValues() const {
assert(_eigenPairs.size() == N);
vector<double> ret;
for (size_t i = 0; i < N; ++i) {
ret.push_back(_eigenPairs[i].first);
}
return ret;
}
vector<Vector<N> > getEigenVectors() const {
assert(_eigenPairs.size() == N);
vector<Vector<N> > ret;
for (size_t i = 0; i < N; ++i) {
ret.push_back(_eigenPairs[i].second);
}
return ret;
}
private:
EigenPairs _eigenPairs;
};
/// Comparison functor for "eigen-pairs".
template <size_t N>
struct EigenPairCmp :
public std::binary_function<const typename EigenSystem<N>::EigenPair&,
const typename EigenSystem<N>::EigenPair&, bool> {
bool operator()(const typename EigenSystem<N>::EigenPair& a,
const typename EigenSystem<N>::EigenPair& b) {
return a.first < b.first;
}
};
/// Diagonalize an NxN matrix, returning a collection of pairs of
/// eigenvalues and eigenvectors, ordered decreasing in eigenvalue.
template <size_t N>
EigenSystem<N> diagonalize(const Matrix<N>& m) {
EigenSystem<N> esys;
// Make a GSL matrix.
gsl_matrix* A = gsl_matrix_alloc(N, N);
for (size_t i = 0; i < N; ++i) {
for (size_t j = 0; j < N; ++j) {
gsl_matrix_set(A, i, j, m.get(i, j));
}
}
// Use GSL diagonalization.
gsl_matrix* vecs = gsl_matrix_alloc(N, N);
gsl_vector* vals = gsl_vector_alloc(N);
gsl_eigen_symmv_workspace* workspace = gsl_eigen_symmv_alloc(N);
gsl_eigen_symmv(A, vals, vecs, workspace);
gsl_eigen_symmv_sort(vals, vecs, GSL_EIGEN_SORT_VAL_DESC);
// Build the vector of "eigen-pairs".
typename EigenSystem<N>::EigenPairs eigensolns;
for (size_t i = 0; i < N; ++i) {
typename EigenSystem<N>::EigenPair ep;
ep.first = gsl_vector_get(vals, i);
Vector<N> ev;
for (size_t j = 0; j < N; ++j) {
ev.set(j, gsl_matrix_get(vecs, j, i));
}
ep.second = ev;
eigensolns.push_back(ep);
}
// Free GSL memory.
gsl_eigen_symmv_free(workspace);
gsl_matrix_free(A);
gsl_matrix_free(vecs);
gsl_vector_free(vals);
// Populate the returned object.
esys._eigenPairs = eigensolns;
return esys;
}
template <size_t N>
inline const string toString(const typename EigenSystem<N>::EigenPair& e) {
ostringstream ss;
//for (typename EigenSystem<N>::EigenPairs::const_iterator i = e.begin(); i != e.end(); ++i) {
ss << e->first << " -> " << e->second;
// if (i+1 != e.end()) ss << endl;
//}
return ss.str();
}
template <size_t N>
inline ostream& operator<<(std::ostream& out, const typename EigenSystem<N>::EigenPair& e) {
out << toString(e);
return out;
}
}
#endif
|