gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
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 
86  double Square(const std::vector<double> &y, const std::vector<double> &y0)
87  {
88  int i, j, num = y.size();
89  double sum;
90  std::vector<double> x(num);
91 
92  for (i = 0; i < num; i++)
93  {
94  for (sum = (y[i]-y0[i]), j=0; j < i; j++)
95  sum -= el[i][j]*x[j];
96 
97  x[i]=sum/el[i][i];
98  }
99 
100  sum = 0.0;
101  for (i = 0; i < num; i++)
102  sum += x[i]*x[i];
103 
104  return sum;
105  }
106 
107  double DetSqrt()
108  {
109  double temp = 1.0;
110  int num = el.size();
111  for (int i = 0; i < num; i++)
112  temp *= el[i][i];
113  return temp;
114  }
115  };
116 }
117 
118 #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:86
STL namespace.
START_MODEL b
Definition: demo.hpp:235
Cholesky(const int num)
Definition: cholesky.hpp:34
std::vector< std::vector< double > > el
Definition: cholesky.hpp:29
double DetSqrt()
Definition: cholesky.hpp:107
void ElMult(std::vector< double > &y) const
Definition: cholesky.hpp:69
TODO: see if we can use this one:
Definition: Analysis.hpp:33