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