This file is indexed.

/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