MatrixDiag.hh
#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);
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;
};
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;
}
};
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
Rivet::diagonalize
EigenSystem< N > diagonalize(const Matrix< N > &m)
Definition: MatrixDiag.hh:112
Rivet::EigenSystem
Handy object containing results of a diagonalization.
Definition: MatrixDiag.hh:42
Rivet::Vector::set
Vector< N > & set(const size_t index, const double value)
Set indexed value.
Definition: VectorN.hh:52
Rivet::Matrix
General -dimensional mathematical matrix object.
Definition: MatrixN.hh:14
Rivet::Vector
A minimal base class for -dimensional vectors.
Definition: VectorN.hh:13
Rivet::EigenPairCmp
Comparison functor for &quot;eigen-pairs&quot;.
Definition: MatrixDiag.hh:99
Rivet::toString
std::string toString(const AnalysisInfo &ai)
String representation.
Definition: AnalysisInfo.cc:234