gambit is hosted by Hepforge, IPPP Durham
 GAMBIT  v1.5.0-2191-ga4742ac a Global And Modular Bsm Inference Tool
cholesky.hpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
10 //
16
17 #ifndef CHOLESKY_HPP
18 #define CHOLESKY_HPP
19
20 #include <vector>
21 #include <cmath>
22 #include <iostream>
23 namespace Gambit
24 {
25  // Cholesky decomposition class
26  class Cholesky
27  {
28  private:
29  std::vector<std::vector<double>> el;
30
31  public:
33
34  Cholesky(const int num) : el(num, std::vector<double>(num)) {}
35
36  bool EnterMat(std::vector<std::vector<double>> &a)
37  {
38  el = a;
39  int num = el.size();
40  double sum = 0;
41  int i, j, k;
42
43  for (i = 0; i < num; i++)
44  {
45  for (j = i; j < num; j++)
46  {
47  for(sum = el[i][j], k = i - 1; k >= 0; k--)
48  sum -= el[i][k]*el[j][k];
49  if(i == j)
50  {
51  if(sum <= 0.0)
52  {
53  return false;
54  }
55  el[i][i] = std::sqrt(sum);
56  }
57  else
58  el[j][i] = sum/el[i][i];
59  }
60  }
61
62  for (i = 0; i < num; i++)
63  for (j = 0; j < i; j++)
64  el[j][i] = 0.0;
65
66  return true;
67  }
68
69  void ElMult (std::vector<double> &y) const
70  {
71  int i, j;
72  int num = el.size();
73  std::vector<double> b(num);
74  for(i = 0; i < num; i++)
75  {
76  b[i] = 0.0;
77  for (j = 0; j <= i; j++)
78  {
79  b[i] += el[i][j]*y[j];
80  }
81  }
82
83  y = b;
84  }
85
91  std::vector<double> invElMult(const std::vector<double> &y) const
92  {
93  std::vector<double> x(y.size(), 0.);
94  for (int i = 0, n = x.size(); i < n; i++)
95  {
96  double sum_ = 0;
97  for (int j = 0; j < i; j++)
98  {
99  sum_ += el[i][j] * x[j];
100  }
101  x[i] = (y[i] - sum_) / el[i][i];
102  }
103  return x;
104  }
105
106  double Square(const std::vector<double> &y, const std::vector<double> &y0)
107  {
108  int i, j, num = y.size();
109  double sum;
110  std::vector<double> x(num);
111
112  for (i = 0; i < num; i++)
113  {
114  for (sum = (y[i]-y0[i]), j=0; j < i; j++)
115  sum -= el[i][j]*x[j];
116
117  x[i]=sum/el[i][i];
118  }
119
120  sum = 0.0;
121  for (i = 0; i < num; i++)
122  sum += x[i]*x[i];
123
124  return sum;
125  }
126
127  double DetSqrt()
128  {
129  double temp = 1.0;
130  int num = el.size();
131  for (int i = 0; i < num; i++)
132  temp *= el[i][i];
133  return temp;
134  }
135  };
136 }
137
138 #endif
DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry double
bool EnterMat(std::vector< std::vector< double >> &a)
Definition: cholesky.hpp:36
double Square(const std::vector< double > &y, const std::vector< double > &y0)
Definition: cholesky.hpp:106
STL namespace.
START_MODEL b
Definition: demo.hpp:270
Cholesky(const int num)
Definition: cholesky.hpp:34
std::vector< std::vector< double > > el
Definition: cholesky.hpp:29
double DetSqrt()
Definition: cholesky.hpp:127
std::vector< double > invElMult(const std::vector< double > &y) const
x = L^-1 y where L is the lower-diagonal Cholesky matrix
Definition: cholesky.hpp:91
void ElMult(std::vector< double > &y) const
Definition: cholesky.hpp:69
TODO: see if we can use this one:
Definition: Analysis.hpp:33