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