Logo ROOT   6.13/01
Reference Guide
CholeskyDecomp.h
Go to the documentation of this file.
1 // @(#)root/smatrix:$Id$
2 // Author: M. Schiller 2009
3 
4 #ifndef ROOT_Math_CholeskyDecomp
5 #define ROOT_Math_CholeskyDecomp
6 
7 /** @file
8  * header file containing the templated implementation of matrix inversion
9  * routines for use with ROOT's SMatrix classes (symmetric positive
10  * definite case)
11  *
12  * @author Manuel Schiller
13  * @date Aug 29 2008
14  * initial release inside LHCb
15  * @date May 7 2009
16  * factored code to provide a nice Cholesky decomposition class, along
17  * with separate methods for solving a single linear system and to
18  * obtain the inverse matrix from the decomposition
19  * @date July 15th 2013
20  * provide a version of that class which works if the dimension of the
21  * problem is only known at run time
22  * @date September 30th 2013
23  * provide routines to access the result of the decomposition L and its
24  * inverse
25  */
26 
27 #include <cmath>
28 #include <algorithm>
29 
30 namespace ROOT {
31 
32  namespace Math {
33 
34 /// helpers for CholeskyDecomp
35 namespace CholeskyDecompHelpers {
36  // forward decls
37  template<class F, class M> struct _decomposerGenDim;
38  template<class F, unsigned N, class M> struct _decomposer;
39  template<class F, class M> struct _inverterGenDim;
40  template<class F, unsigned N, class M> struct _inverter;
41  template<class F, class V> struct _solverGenDim;
42  template<class F, unsigned N, class V> struct _solver;
43  template<typename G> class PackedArrayAdapter;
44 }
45 
46 /// class to compute the Cholesky decomposition of a matrix
47 /** class to compute the Cholesky decomposition of a symmetric
48  * positive definite matrix
49  *
50  * provides routines to check if the decomposition succeeded (i.e. if
51  * matrix is positive definite and non-singular), to solve a linear
52  * system for the given matrix and to obtain its inverse
53  *
54  * the actual functionality is implemented in templated helper
55  * classes which have specializations for dimensions N = 1 to 6
56  * to achieve a gain in speed for common matrix sizes
57  *
58  * usage example:
59  * @code
60  * // let m be a symmetric positive definite SMatrix (use type float
61  * // for internal computations, matrix size is 4x4)
62  * CholeskyDecomp<float, 4> decomp(m);
63  * // check if the decomposition succeeded
64  * if (!decomp) {
65  * std::cerr << "decomposition failed!" << std::endl;
66  * } else {
67  * // let rhs be a vector; we seek a vector x such that m * x = rhs
68  * decomp.Solve(rhs);
69  * // rhs now contains the solution we are looking for
70  *
71  * // obtain the inverse of m, put it into m itself
72  * decomp.Invert(m);
73  * }
74  * @endcode
75  */
76 template<class F, unsigned N> class CholeskyDecomp
77 {
78 private:
79  /// lower triangular matrix L
80  /** lower triangular matrix L, packed storage, with diagonal
81  * elements pre-inverted */
82  F fL[N * (N + 1) / 2];
83  /// flag indicating a successful decomposition
84  bool fOk;
85 public:
86  /// perform a Cholesky decomposition
87  /** perfrom a Cholesky decomposition of a symmetric positive
88  * definite matrix m
89  *
90  * this is the constructor to uses with an SMatrix (and objects
91  * that behave like an SMatrix in terms of using
92  * operator()(int i, int j) for access to elements)
93  */
94  template<class M> CholeskyDecomp(const M& m) :
95  fL( ), fOk(false)
96  {
98  fOk = _decomposer<F, N, M>()(fL, m);
99  }
100 
101  /// perform a Cholesky decomposition
102  /** perfrom a Cholesky decomposition of a symmetric positive
103  * definite matrix m
104  *
105  * this is the constructor to use in special applications where
106  * plain arrays are used
107  *
108  * NOTE: the matrix is given in packed representation, matrix
109  * element m(i,j) (j <= i) is supposed to be in array element
110  * (i * (i + 1)) / 2 + j
111  */
112  template<typename G> CholeskyDecomp(G* m) :
113  fL(), fOk(false)
114  {
117  fOk = _decomposer<F, N, PackedArrayAdapter<G> >()(
118  fL, PackedArrayAdapter<G>(m));
119  }
120 
121  /// returns true if decomposition was successful
122  /** @returns true if decomposition was successful */
123  bool ok() const { return fOk; }
124  /// returns true if decomposition was successful
125  /** @returns true if decomposition was successful */
126  operator bool() const { return fOk; }
127 
128  /** @brief solves a linear system for the given right hand side
129  *
130  * Note that you can use both SVector classes and plain arrays for
131  * rhs. (Make sure that the sizes match!). It will work with any vector
132  * implementing the operator [i]
133  *
134  * @returns if the decomposition was successful
135  */
136  template<class V> bool Solve(V& rhs) const
137  {
139  if (fOk) _solver<F,N,V>()(rhs, fL);
140  return fOk;
141  }
142 
143  /** @brief place the inverse into m
144  *
145  * This is the method to use with an SMatrix.
146  *
147  * @returns if the decomposition was successful
148  */
149  template<class M> bool Invert(M& m) const
150  {
152  if (fOk) _inverter<F,N,M>()(m, fL);
153  return fOk;
154  }
155 
156  /** @brief place the inverse into m
157  *
158  * This is the method to use with a plain array.
159  *
160  * @returns if the decomposition was successful
161  *
162  * NOTE: the matrix is given in packed representation, matrix
163  * element m(i,j) (j <= i) is supposed to be in array element
164  * (i * (i + 1)) / 2 + j
165  */
166  template<typename G> bool Invert(G* m) const
167  {
170  if (fOk) {
171  PackedArrayAdapter<G> adapted(m);
172  _inverter<F,N,PackedArrayAdapter<G> >()(adapted, fL);
173  }
174  return fOk;
175  }
176 
177  /** @brief obtain the decomposed matrix L
178  *
179  * This is the method to use with a plain array.
180  *
181  * @returns if the decomposition was successful
182  */
183  template<class M> bool getL(M& m) const
184  {
185  if (!fOk) return false;
186  for (unsigned i = 0; i < N; ++i) {
187  // zero upper half of matrix
188  for (unsigned j = i + 1; j < N; ++j)
189  m(i, j) = F(0);
190  // copy the rest
191  for (unsigned j = 0; j <= i; ++j)
192  m(i, j) = fL[i * (i + 1) / 2 + j];
193  // adjust the diagonal - we save 1/L(i, i) in that position, so
194  // convert to what caller expects
195  m(i, i) = F(1) / m(i, i);
196  }
197  return true;
198  }
199 
200  /** @brief obtain the decomposed matrix L
201  *
202  * @returns if the decomposition was successful
203  *
204  * NOTE: the matrix is given in packed representation, matrix
205  * element m(i,j) (j <= i) is supposed to be in array element
206  * (i * (i + 1)) / 2 + j
207  */
208  template<typename G> bool getL(G* m) const
209  {
210  if (!fOk) return false;
211  // copy L
212  for (unsigned i = 0; i < (N * (N + 1)) / 2; ++i)
213  m[i] = fL[i];
214  // adjust diagonal - we save 1/L(i, i) in that position, so convert to
215  // what caller expects
216  for (unsigned i = 0; i < N; ++i)
217  m[(i * (i + 1)) / 2 + i] = F(1) / fL[(i * (i + 1)) / 2 + i];
218  return true;
219  }
220 
221  /** @brief obtain the inverse of the decomposed matrix L
222  *
223  * This is the method to use with a plain array.
224  *
225  * @returns if the decomposition was successful
226  */
227  template<class M> bool getLi(M& m) const
228  {
229  if (!fOk) return false;
230  for (unsigned i = 0; i < N; ++i) {
231  // zero lower half of matrix
232  for (unsigned j = i + 1; j < N; ++j)
233  m(j, i) = F(0);
234  // copy the rest
235  for (unsigned j = 0; j <= i; ++j)
236  m(j, i) = fL[i * (i + 1) / 2 + j];
237  }
238  // invert the off-diagonal part of what we just copied
239  for (unsigned i = 1; i < N; ++i) {
240  for (unsigned j = 0; j < i; ++j) {
241  typename M::value_type tmp = F(0);
242  for (unsigned k = i; k-- > j;)
243  tmp -= m(k, i) * m(j, k);
244  m(j, i) = tmp * m(i, i);
245  }
246  }
247  return true;
248  }
249 
250  /** @brief obtain the inverse of the decomposed matrix L
251  *
252  * @returns if the decomposition was successful
253  *
254  * NOTE: the matrix is given in packed representation, matrix
255  * element m(j,i) (j <= i) is supposed to be in array element
256  * (i * (i + 1)) / 2 + j
257  */
258  template<typename G> bool getLi(G* m) const
259  {
260  if (!fOk) return false;
261  // copy L
262  for (unsigned i = 0; i < (N * (N + 1)) / 2; ++i)
263  m[i] = fL[i];
264  // invert the off-diagonal part of what we just copied
265  G* base1 = &m[1];
266  for (unsigned i = 1; i < N; base1 += ++i) {
267  for (unsigned j = 0; j < i; ++j) {
268  G tmp = F(0);
269  const G *base2 = &m[(i * (i - 1)) / 2];
270  for (unsigned k = i; k-- > j; base2 -= k)
271  tmp -= base1[k] * base2[j];
272  base1[j] = tmp * base1[i];
273  }
274  }
275  return true;
276  }
277 };
278 
279 /// class to compute the Cholesky decomposition of a matrix
280 /** class to compute the Cholesky decomposition of a symmetric
281  * positive definite matrix when the dimensionality of the problem is not known
282  * at compile time
283  *
284  * provides routines to check if the decomposition succeeded (i.e. if
285  * matrix is positive definite and non-singular), to solve a linear
286  * system for the given matrix and to obtain its inverse
287  *
288  * the actual functionality is implemented in templated helper
289  * classes which have specializations for dimensions N = 1 to 6
290  * to achieve a gain in speed for common matrix sizes
291  *
292  * usage example:
293  * @code
294  * // let m be a symmetric positive definite SMatrix (use type float
295  * // for internal computations, matrix size is 4x4)
296  * CholeskyDecompGenDim<float> decomp(4, m);
297  * // check if the decomposition succeeded
298  * if (!decomp) {
299  * std::cerr << "decomposition failed!" << std::endl;
300  * } else {
301  * // let rhs be a vector; we seek a vector x such that m * x = rhs
302  * decomp.Solve(rhs);
303  * // rhs now contains the solution we are looking for
304  *
305  * // obtain the inverse of m, put it into m itself
306  * decomp.Invert(m);
307  * }
308  * @endcode
309  */
310 template<class F> class CholeskyDecompGenDim
311 {
312 private:
313  /** @brief dimensionality
314  * dimensionality of the problem */
315  unsigned fN;
316  /// lower triangular matrix L
317  /** lower triangular matrix L, packed storage, with diagonal
318  * elements pre-inverted */
319  F *fL;
320  /// flag indicating a successful decomposition
321  bool fOk;
322 public:
323  /// perform a Cholesky decomposition
324  /** perfrom a Cholesky decomposition of a symmetric positive
325  * definite matrix m
326  *
327  * this is the constructor to uses with an SMatrix (and objects
328  * that behave like an SMatrix in terms of using
329  * operator()(int i, int j) for access to elements)
330  */
331  template<class M> CholeskyDecompGenDim(unsigned N, const M& m) :
332  fN(N), fL(new F[(fN * (fN + 1)) / 2]), fOk(false)
333  {
335  fOk = _decomposerGenDim<F, M>()(fL, m, fN);
336  }
337 
338  /// perform a Cholesky decomposition
339  /** perfrom a Cholesky decomposition of a symmetric positive
340  * definite matrix m
341  *
342  * this is the constructor to use in special applications where
343  * plain arrays are used
344  *
345  * NOTE: the matrix is given in packed representation, matrix
346  * element m(i,j) (j <= i) is supposed to be in array element
347  * (i * (i + 1)) / 2 + j
348  */
349  template<typename G> CholeskyDecompGenDim(unsigned N, G* m) :
350  fN(N), fL(new F[(fN * (fN + 1)) / 2]), fOk(false)
351  {
354  fOk = _decomposerGenDim<F, PackedArrayAdapter<G> >()(
355  fL, PackedArrayAdapter<G>(m), fN);
356  }
357 
358  /// destructor
359  ~CholeskyDecompGenDim() { delete[] fL; }
360 
361  /// returns true if decomposition was successful
362  /** @returns true if decomposition was successful */
363  bool ok() const { return fOk; }
364  /// returns true if decomposition was successful
365  /** @returns true if decomposition was successful */
366  operator bool() const { return fOk; }
367 
368  /** @brief solves a linear system for the given right hand side
369  *
370  * Note that you can use both SVector classes and plain arrays for
371  * rhs. (Make sure that the sizes match!). It will work with any vector
372  * implementing the operator [i]
373  *
374  * @returns if the decomposition was successful
375  */
376  template<class V> bool Solve(V& rhs) const
377  {
379  if (fOk) _solverGenDim<F,V>()(rhs, fL, fN);
380  return fOk;
381  }
382 
383  /** @brief place the inverse into m
384  *
385  * This is the method to use with an SMatrix.
386  *
387  * @returns if the decomposition was successful
388  */
389  template<class M> bool Invert(M& m) const
390  {
392  if (fOk) _inverterGenDim<F,M>()(m, fL, fN);
393  return fOk;
394  }
395 
396  /** @brief place the inverse into m
397  *
398  * This is the method to use with a plain array.
399  *
400  * @returns if the decomposition was successful
401  *
402  * NOTE: the matrix is given in packed representation, matrix
403  * element m(i,j) (j <= i) is supposed to be in array element
404  * (i * (i + 1)) / 2 + j
405  */
406  template<typename G> bool Invert(G* m) const
407  {
410  if (fOk) {
411  PackedArrayAdapter<G> adapted(m);
412  _inverterGenDim<F,PackedArrayAdapter<G> >()(adapted, fL, fN);
413  }
414  return fOk;
415  }
416 
417  /** @brief obtain the decomposed matrix L
418  *
419  * This is the method to use with a plain array.
420  *
421  * @returns if the decomposition was successful
422  */
423  template<class M> bool getL(M& m) const
424  {
425  if (!fOk) return false;
426  for (unsigned i = 0; i < fN; ++i) {
427  // zero upper half of matrix
428  for (unsigned j = i + 1; j < fN; ++j)
429  m(i, j) = F(0);
430  // copy the rest
431  for (unsigned j = 0; j <= i; ++j)
432  m(i, j) = fL[i * (i + 1) / 2 + j];
433  // adjust the diagonal - we save 1/L(i, i) in that position, so
434  // convert to what caller expects
435  m(i, i) = F(1) / m(i, i);
436  }
437  return true;
438  }
439 
440  /** @brief obtain the decomposed matrix L
441  *
442  * @returns if the decomposition was successful
443  *
444  * NOTE: the matrix is given in packed representation, matrix
445  * element m(i,j) (j <= i) is supposed to be in array element
446  * (i * (i + 1)) / 2 + j
447  */
448  template<typename G> bool getL(G* m) const
449  {
450  if (!fOk) return false;
451  // copy L
452  for (unsigned i = 0; i < (fN * (fN + 1)) / 2; ++i)
453  m[i] = fL[i];
454  // adjust diagonal - we save 1/L(i, i) in that position, so convert to
455  // what caller expects
456  for (unsigned i = 0; i < fN; ++i)
457  m[(i * (i + 1)) / 2 + i] = F(1) / fL[(i * (i + 1)) / 2 + i];
458  return true;
459  }
460 
461  /** @brief obtain the inverse of the decomposed matrix L
462  *
463  * This is the method to use with a plain array.
464  *
465  * @returns if the decomposition was successful
466  */
467  template<class M> bool getLi(M& m) const
468  {
469  if (!fOk) return false;
470  for (unsigned i = 0; i < fN; ++i) {
471  // zero lower half of matrix
472  for (unsigned j = i + 1; j < fN; ++j)
473  m(j, i) = F(0);
474  // copy the rest
475  for (unsigned j = 0; j <= i; ++j)
476  m(j, i) = fL[i * (i + 1) / 2 + j];
477  }
478  // invert the off-diagonal part of what we just copied
479  for (unsigned i = 1; i < fN; ++i) {
480  for (unsigned j = 0; j < i; ++j) {
481  typename M::value_type tmp = F(0);
482  for (unsigned k = i; k-- > j;)
483  tmp -= m(k, i) * m(j, k);
484  m(j, i) = tmp * m(i, i);
485  }
486  }
487  return true;
488  }
489 
490  /** @brief obtain the inverse of the decomposed matrix L
491  *
492  * @returns if the decomposition was successful
493  *
494  * NOTE: the matrix is given in packed representation, matrix
495  * element m(j,i) (j <= i) is supposed to be in array element
496  * (i * (i + 1)) / 2 + j
497  */
498  template<typename G> bool getLi(G* m) const
499  {
500  if (!fOk) return false;
501  // copy L
502  for (unsigned i = 0; i < (fN * (fN + 1)) / 2; ++i)
503  m[i] = fL[i];
504  // invert the off-diagonal part of what we just copied
505  G* base1 = &m[1];
506  for (unsigned i = 1; i < fN; base1 += ++i) {
507  for (unsigned j = 0; j < i; ++j) {
508  G tmp = F(0);
509  const G *base2 = &m[(i * (i - 1)) / 2];
510  for (unsigned k = i; k-- > j; base2 -= k)
511  tmp -= base1[k] * base2[j];
512  base1[j] = tmp * base1[i];
513  }
514  }
515  return true;
516  }
517 };
518 
519 namespace CholeskyDecompHelpers {
520  /// adapter for packed arrays (to SMatrix indexing conventions)
521  template<typename G> class PackedArrayAdapter
522  {
523  private:
524  G* fArr; ///< pointer to first array element
525  public:
526  /// constructor
527  PackedArrayAdapter(G* arr) : fArr(arr) {}
528  /// read access to elements (make sure that j <= i)
529  const G operator()(unsigned i, unsigned j) const
530  { return fArr[((i * (i + 1)) / 2) + j]; }
531  /// write access to elements (make sure that j <= i)
532  G& operator()(unsigned i, unsigned j)
533  { return fArr[((i * (i + 1)) / 2) + j]; }
534  };
535  /// struct to do a Cholesky decomposition (general dimensionality)
536  template<class F, class M> struct _decomposerGenDim
537  {
538  /// method to do the decomposition
539  /** @returns if the decomposition was successful */
540  bool operator()(F* dst, const M& src, unsigned N) const
541  {
542  // perform Cholesky decomposition of matrix: M = L L^T
543  // only thing that can go wrong: trying to take square
544  // root of negative number or zero (matrix is
545  // ill-conditioned or singular in these cases)
546 
547  // element L(i,j) is at array position (i * (i+1)) / 2 + j
548 
549  // quirk: we may need to invert L later anyway, so we can
550  // invert elements on diagonale straight away (we only
551  // ever need their reciprocals!)
552 
553  // cache starting address of rows of L for speed reasons
554  F *base1 = &dst[0];
555  for (unsigned i = 0; i < N; base1 += ++i) {
556  F tmpdiag = F(0.0); // for element on diagonale
557  // calculate off-diagonal elements
558  F *base2 = &dst[0];
559  for (unsigned j = 0; j < i; base2 += ++j) {
560  F tmp = src(i, j);
561  for (unsigned k = j; k--; )
562  tmp -= base1[k] * base2[k];
563  base1[j] = tmp *= base2[j];
564  // keep track of contribution to element on diagonale
565  tmpdiag += tmp * tmp;
566  }
567  // keep truncation error small
568  tmpdiag = src(i, i) - tmpdiag;
569  // check if positive definite
570  if (tmpdiag <= F(0.0)) return false;
571  else base1[i] = std::sqrt(F(1.0) / tmpdiag);
572  }
573  return true;
574  }
575  };
576 
577  /// struct to do a Cholesky decomposition
578  template<class F, unsigned N, class M> struct _decomposer
579  {
580  /// method to do the decomposition
581  /** @returns if the decomposition was successful */
582  bool operator()(F* dst, const M& src) const
583  { return _decomposerGenDim<F, M>()(dst, src, N); }
584  };
585 
586  /// struct to obtain the inverse from a Cholesky decomposition (general dimensionality)
587  template<class F, class M> struct _inverterGenDim
588  {
589  /// method to do the inversion
590  void operator()(M& dst, const F* src, unsigned N) const
591  {
592  // make working copy
593  F * l = new F[N * (N + 1) / 2];
594  std::copy(src, src + ((N * (N + 1)) / 2), l);
595  // ok, next step: invert off-diagonal part of matrix
596  F* base1 = &l[1];
597  for (unsigned i = 1; i < N; base1 += ++i) {
598  for (unsigned j = 0; j < i; ++j) {
599  F tmp = F(0.0);
600  const F *base2 = &l[(i * (i - 1)) / 2];
601  for (unsigned k = i; k-- > j; base2 -= k)
602  tmp -= base1[k] * base2[j];
603  base1[j] = tmp * base1[i];
604  }
605  }
606 
607  // Li = L^(-1) formed, now calculate M^(-1) = Li^T Li
608  for (unsigned i = N; i--; ) {
609  for (unsigned j = i + 1; j--; ) {
610  F tmp = F(0.0);
611  base1 = &l[(N * (N - 1)) / 2];
612  for (unsigned k = N; k-- > i; base1 -= k)
613  tmp += base1[i] * base1[j];
614  dst(i, j) = tmp;
615  }
616  }
617  delete [] l;
618  }
619  };
620 
621  /// struct to obtain the inverse from a Cholesky decomposition
622  template<class F, unsigned N, class M> struct _inverter
623  {
624  /// method to do the inversion
625  void operator()(M& dst, const F* src) const
626  { return _inverterGenDim<F, M>()(dst, src, N); }
627  };
628 
629  /// struct to solve a linear system using its Cholesky decomposition (generalised dimensionality)
630  template<class F, class V> struct _solverGenDim
631  {
632  /// method to solve the linear system
633  void operator()(V& rhs, const F* l, unsigned N) const
634  {
635  // solve Ly = rhs
636  for (unsigned k = 0; k < N; ++k) {
637  const unsigned base = (k * (k + 1)) / 2;
638  F sum = F(0.0);
639  for (unsigned i = k; i--; )
640  sum += rhs[i] * l[base + i];
641  // elements on diagonale are pre-inverted!
642  rhs[k] = (rhs[k] - sum) * l[base + k];
643  }
644  // solve L^Tx = y
645  for (unsigned k = N; k--; ) {
646  F sum = F(0.0);
647  for (unsigned i = N; --i > k; )
648  sum += rhs[i] * l[(i * (i + 1)) / 2 + k];
649  // elements on diagonale are pre-inverted!
650  rhs[k] = (rhs[k] - sum) * l[(k * (k + 1)) / 2 + k];
651  }
652  }
653  };
654 
655  /// struct to solve a linear system using its Cholesky decomposition
656  template<class F, unsigned N, class V> struct _solver
657  {
658  /// method to solve the linear system
659  void operator()(V& rhs, const F* l) const
660  { _solverGenDim<F, V>()(rhs, l, N); }
661  };
662 
663  /// struct to do a Cholesky decomposition (specialized, N = 6)
664  template<class F, class M> struct _decomposer<F, 6, M>
665  {
666  /// method to do the decomposition
667  bool operator()(F* dst, const M& src) const
668  {
669  if (src(0,0) <= F(0.0)) return false;
670  dst[0] = std::sqrt(F(1.0) / src(0,0));
671  dst[1] = src(1,0) * dst[0];
672  dst[2] = src(1,1) - dst[1] * dst[1];
673  if (dst[2] <= F(0.0)) return false;
674  else dst[2] = std::sqrt(F(1.0) / dst[2]);
675  dst[3] = src(2,0) * dst[0];
676  dst[4] = (src(2,1) - dst[1] * dst[3]) * dst[2];
677  dst[5] = src(2,2) - (dst[3] * dst[3] + dst[4] * dst[4]);
678  if (dst[5] <= F(0.0)) return false;
679  else dst[5] = std::sqrt(F(1.0) / dst[5]);
680  dst[6] = src(3,0) * dst[0];
681  dst[7] = (src(3,1) - dst[1] * dst[6]) * dst[2];
682  dst[8] = (src(3,2) - dst[3] * dst[6] - dst[4] * dst[7]) * dst[5];
683  dst[9] = src(3,3) - (dst[6] * dst[6] + dst[7] * dst[7] + dst[8] * dst[8]);
684  if (dst[9] <= F(0.0)) return false;
685  else dst[9] = std::sqrt(F(1.0) / dst[9]);
686  dst[10] = src(4,0) * dst[0];
687  dst[11] = (src(4,1) - dst[1] * dst[10]) * dst[2];
688  dst[12] = (src(4,2) - dst[3] * dst[10] - dst[4] * dst[11]) * dst[5];
689  dst[13] = (src(4,3) - dst[6] * dst[10] - dst[7] * dst[11] - dst[8] * dst[12]) * dst[9];
690  dst[14] = src(4,4) - (dst[10]*dst[10]+dst[11]*dst[11]+dst[12]*dst[12]+dst[13]*dst[13]);
691  if (dst[14] <= F(0.0)) return false;
692  else dst[14] = std::sqrt(F(1.0) / dst[14]);
693  dst[15] = src(5,0) * dst[0];
694  dst[16] = (src(5,1) - dst[1] * dst[15]) * dst[2];
695  dst[17] = (src(5,2) - dst[3] * dst[15] - dst[4] * dst[16]) * dst[5];
696  dst[18] = (src(5,3) - dst[6] * dst[15] - dst[7] * dst[16] - dst[8] * dst[17]) * dst[9];
697  dst[19] = (src(5,4) - dst[10] * dst[15] - dst[11] * dst[16] - dst[12] * dst[17] - dst[13] * dst[18]) * dst[14];
698  dst[20] = src(5,5) - (dst[15]*dst[15]+dst[16]*dst[16]+dst[17]*dst[17]+dst[18]*dst[18]+dst[19]*dst[19]);
699  if (dst[20] <= F(0.0)) return false;
700  else dst[20] = std::sqrt(F(1.0) / dst[20]);
701  return true;
702  }
703  };
704  /// struct to do a Cholesky decomposition (specialized, N = 5)
705  template<class F, class M> struct _decomposer<F, 5, M>
706  {
707  /// method to do the decomposition
708  bool operator()(F* dst, const M& src) const
709  {
710  if (src(0,0) <= F(0.0)) return false;
711  dst[0] = std::sqrt(F(1.0) / src(0,0));
712  dst[1] = src(1,0) * dst[0];
713  dst[2] = src(1,1) - dst[1] * dst[1];
714  if (dst[2] <= F(0.0)) return false;
715  else dst[2] = std::sqrt(F(1.0) / dst[2]);
716  dst[3] = src(2,0) * dst[0];
717  dst[4] = (src(2,1) - dst[1] * dst[3]) * dst[2];
718  dst[5] = src(2,2) - (dst[3] * dst[3] + dst[4] * dst[4]);
719  if (dst[5] <= F(0.0)) return false;
720  else dst[5] = std::sqrt(F(1.0) / dst[5]);
721  dst[6] = src(3,0) * dst[0];
722  dst[7] = (src(3,1) - dst[1] * dst[6]) * dst[2];
723  dst[8] = (src(3,2) - dst[3] * dst[6] - dst[4] * dst[7]) * dst[5];
724  dst[9] = src(3,3) - (dst[6] * dst[6] + dst[7] * dst[7] + dst[8] * dst[8]);
725  if (dst[9] <= F(0.0)) return false;
726  else dst[9] = std::sqrt(F(1.0) / dst[9]);
727  dst[10] = src(4,0) * dst[0];
728  dst[11] = (src(4,1) - dst[1] * dst[10]) * dst[2];
729  dst[12] = (src(4,2) - dst[3] * dst[10] - dst[4] * dst[11]) * dst[5];
730  dst[13] = (src(4,3) - dst[6] * dst[10] - dst[7] * dst[11] - dst[8] * dst[12]) * dst[9];
731  dst[14] = src(4,4) - (dst[10]*dst[10]+dst[11]*dst[11]+dst[12]*dst[12]+dst[13]*dst[13]);
732  if (dst[14] <= F(0.0)) return false;
733  else dst[14] = std::sqrt(F(1.0) / dst[14]);
734  return true;
735  }
736  };
737  /// struct to do a Cholesky decomposition (specialized, N = 4)
738  template<class F, class M> struct _decomposer<F, 4, M>
739  {
740  /// method to do the decomposition
741  bool operator()(F* dst, const M& src) const
742  {
743  if (src(0,0) <= F(0.0)) return false;
744  dst[0] = std::sqrt(F(1.0) / src(0,0));
745  dst[1] = src(1,0) * dst[0];
746  dst[2] = src(1,1) - dst[1] * dst[1];
747  if (dst[2] <= F(0.0)) return false;
748  else dst[2] = std::sqrt(F(1.0) / dst[2]);
749  dst[3] = src(2,0) * dst[0];
750  dst[4] = (src(2,1) - dst[1] * dst[3]) * dst[2];
751  dst[5] = src(2,2) - (dst[3] * dst[3] + dst[4] * dst[4]);
752  if (dst[5] <= F(0.0)) return false;
753  else dst[5] = std::sqrt(F(1.0) / dst[5]);
754  dst[6] = src(3,0) * dst[0];
755  dst[7] = (src(3,1) - dst[1] * dst[6]) * dst[2];
756  dst[8] = (src(3,2) - dst[3] * dst[6] - dst[4] * dst[7]) * dst[5];
757  dst[9] = src(3,3) - (dst[6] * dst[6] + dst[7] * dst[7] + dst[8] * dst[8]);
758  if (dst[9] <= F(0.0)) return false;
759  else dst[9] = std::sqrt(F(1.0) / dst[9]);
760  return true;
761  }
762  };
763  /// struct to do a Cholesky decomposition (specialized, N = 3)
764  template<class F, class M> struct _decomposer<F, 3, M>
765  {
766  /// method to do the decomposition
767  bool operator()(F* dst, const M& src) const
768  {
769  if (src(0,0) <= F(0.0)) return false;
770  dst[0] = std::sqrt(F(1.0) / src(0,0));
771  dst[1] = src(1,0) * dst[0];
772  dst[2] = src(1,1) - dst[1] * dst[1];
773  if (dst[2] <= F(0.0)) return false;
774  else dst[2] = std::sqrt(F(1.0) / dst[2]);
775  dst[3] = src(2,0) * dst[0];
776  dst[4] = (src(2,1) - dst[1] * dst[3]) * dst[2];
777  dst[5] = src(2,2) - (dst[3] * dst[3] + dst[4] * dst[4]);
778  if (dst[5] <= F(0.0)) return false;
779  else dst[5] = std::sqrt(F(1.0) / dst[5]);
780  return true;
781  }
782  };
783  /// struct to do a Cholesky decomposition (specialized, N = 2)
784  template<class F, class M> struct _decomposer<F, 2, M>
785  {
786  /// method to do the decomposition
787  bool operator()(F* dst, const M& src) const
788  {
789  if (src(0,0) <= F(0.0)) return false;
790  dst[0] = std::sqrt(F(1.0) / src(0,0));
791  dst[1] = src(1,0) * dst[0];
792  dst[2] = src(1,1) - dst[1] * dst[1];
793  if (dst[2] <= F(0.0)) return false;
794  else dst[2] = std::sqrt(F(1.0) / dst[2]);
795  return true;
796  }
797  };
798  /// struct to do a Cholesky decomposition (specialized, N = 1)
799  template<class F, class M> struct _decomposer<F, 1, M>
800  {
801  /// method to do the decomposition
802  bool operator()(F* dst, const M& src) const
803  {
804  if (src(0,0) <= F(0.0)) return false;
805  dst[0] = std::sqrt(F(1.0) / src(0,0));
806  return true;
807  }
808  };
809  /// struct to do a Cholesky decomposition (specialized, N = 0)
810  template<class F, class M> struct _decomposer<F, 0, M>
811  {
812  private:
813  _decomposer() { };
814  bool operator()(F* dst, const M& src) const;
815  };
816 
817  /// struct to obtain the inverse from a Cholesky decomposition (N = 6)
818  template<class F, class M> struct _inverter<F,6,M>
819  {
820  /// method to do the inversion
821  inline void operator()(M& dst, const F* src) const
822  {
823  const F li21 = -src[1] * src[0] * src[2];
824  const F li32 = -src[4] * src[2] * src[5];
825  const F li31 = (src[1] * src[4] * src[2] - src[3]) * src[0] * src[5];
826  const F li43 = -src[8] * src[9] * src[5];
827  const F li42 = (src[4] * src[8] * src[5] - src[7]) * src[2] * src[9];
828  const F li41 = (-src[1] * src[4] * src[8] * src[2] * src[5] +
829  src[1] * src[7] * src[2] + src[3] * src[8] * src[5] - src[6]) * src[0] * src[9];
830  const F li54 = -src[13] * src[14] * src[9];
831  const F li53 = (src[13] * src[8] * src[9] - src[12]) * src[5] * src[14];
832  const F li52 = (-src[4] * src[8] * src[13] * src[5] * src[9] +
833  src[4] * src[12] * src[5] + src[7] * src[13] * src[9] - src[11]) * src[2] * src[14];
834  const F li51 = (src[1]*src[4]*src[8]*src[13]*src[2]*src[5]*src[9] -
835  src[13]*src[8]*src[3]*src[9]*src[5] - src[12]*src[4]*src[1]*src[2]*src[5] - src[13]*src[7]*src[1]*src[9]*src[2] +
836  src[11]*src[1]*src[2] + src[12]*src[3]*src[5] + src[13]*src[6]*src[9] -src[10]) * src[0] * src[14];
837  const F li65 = -src[19] * src[20] * src[14];
838  const F li64 = (src[19] * src[13] * src[14] - src[18]) * src[9] * src[20];
839  const F li63 = (-src[8] * src[13] * src[19] * src[9] * src[14] +
840  src[8] * src[18] * src[9] + src[12] * src[19] * src[14] - src[17]) * src[5] * src[20];
841  const F li62 = (src[4]*src[8]*src[13]*src[19]*src[5]*src[9]*src[14] -
842  src[18]*src[8]*src[4]*src[9]*src[5] - src[19]*src[12]*src[4]*src[14]*src[5] -src[19]*src[13]*src[7]*src[14]*src[9] +
843  src[17]*src[4]*src[5] + src[18]*src[7]*src[9] + src[19]*src[11]*src[14] - src[16]) * src[2] * src[20];
844  const F li61 = (-src[19]*src[13]*src[8]*src[4]*src[1]*src[2]*src[5]*src[9]*src[14] +
845  src[18]*src[8]*src[4]*src[1]*src[2]*src[5]*src[9] + src[19]*src[12]*src[4]*src[1]*src[2]*src[5]*src[14] +
846  src[19]*src[13]*src[7]*src[1]*src[2]*src[9]*src[14] + src[19]*src[13]*src[8]*src[3]*src[5]*src[9]*src[14] -
847  src[17]*src[4]*src[1]*src[2]*src[5] - src[18]*src[7]*src[1]*src[2]*src[9] - src[19]*src[11]*src[1]*src[2]*src[14] -
848  src[18]*src[8]*src[3]*src[5]*src[9] - src[19]*src[12]*src[3]*src[5]*src[14] - src[19]*src[13]*src[6]*src[9]*src[14] +
849  src[16]*src[1]*src[2] + src[17]*src[3]*src[5] + src[18]*src[6]*src[9] + src[19]*src[10]*src[14] - src[15]) *
850  src[0] * src[20];
851 
852  dst(0,0) = li61*li61 + li51*li51 + li41*li41 + li31*li31 + li21*li21 + src[0]*src[0];
853  dst(1,0) = li61*li62 + li51*li52 + li41*li42 + li31*li32 + li21*src[2];
854  dst(1,1) = li62*li62 + li52*li52 + li42*li42 + li32*li32 + src[2]*src[2];
855  dst(2,0) = li61*li63 + li51*li53 + li41*li43 + li31*src[5];
856  dst(2,1) = li62*li63 + li52*li53 + li42*li43 + li32*src[5];
857  dst(2,2) = li63*li63 + li53*li53 + li43*li43 + src[5]*src[5];
858  dst(3,0) = li61*li64 + li51*li54 + li41*src[9];
859  dst(3,1) = li62*li64 + li52*li54 + li42*src[9];
860  dst(3,2) = li63*li64 + li53*li54 + li43*src[9];
861  dst(3,3) = li64*li64 + li54*li54 + src[9]*src[9];
862  dst(4,0) = li61*li65 + li51*src[14];
863  dst(4,1) = li62*li65 + li52*src[14];
864  dst(4,2) = li63*li65 + li53*src[14];
865  dst(4,3) = li64*li65 + li54*src[14];
866  dst(4,4) = li65*li65 + src[14]*src[14];
867  dst(5,0) = li61*src[20];
868  dst(5,1) = li62*src[20];
869  dst(5,2) = li63*src[20];
870  dst(5,3) = li64*src[20];
871  dst(5,4) = li65*src[20];
872  dst(5,5) = src[20]*src[20];
873  }
874  };
875  /// struct to obtain the inverse from a Cholesky decomposition (N = 5)
876  template<class F, class M> struct _inverter<F,5,M>
877  {
878  /// method to do the inversion
879  inline void operator()(M& dst, const F* src) const
880  {
881  const F li21 = -src[1] * src[0] * src[2];
882  const F li32 = -src[4] * src[2] * src[5];
883  const F li31 = (src[1] * src[4] * src[2] - src[3]) * src[0] * src[5];
884  const F li43 = -src[8] * src[9] * src[5];
885  const F li42 = (src[4] * src[8] * src[5] - src[7]) * src[2] * src[9];
886  const F li41 = (-src[1] * src[4] * src[8] * src[2] * src[5] +
887  src[1] * src[7] * src[2] + src[3] * src[8] * src[5] - src[6]) * src[0] * src[9];
888  const F li54 = -src[13] * src[14] * src[9];
889  const F li53 = (src[13] * src[8] * src[9] - src[12]) * src[5] * src[14];
890  const F li52 = (-src[4] * src[8] * src[13] * src[5] * src[9] +
891  src[4] * src[12] * src[5] + src[7] * src[13] * src[9] - src[11]) * src[2] * src[14];
892  const F li51 = (src[1]*src[4]*src[8]*src[13]*src[2]*src[5]*src[9] -
893  src[13]*src[8]*src[3]*src[9]*src[5] - src[12]*src[4]*src[1]*src[2]*src[5] - src[13]*src[7]*src[1]*src[9]*src[2] +
894  src[11]*src[1]*src[2] + src[12]*src[3]*src[5] + src[13]*src[6]*src[9] -src[10]) * src[0] * src[14];
895 
896  dst(0,0) = li51*li51 + li41*li41 + li31*li31 + li21*li21 + src[0]*src[0];
897  dst(1,0) = li51*li52 + li41*li42 + li31*li32 + li21*src[2];
898  dst(1,1) = li52*li52 + li42*li42 + li32*li32 + src[2]*src[2];
899  dst(2,0) = li51*li53 + li41*li43 + li31*src[5];
900  dst(2,1) = li52*li53 + li42*li43 + li32*src[5];
901  dst(2,2) = li53*li53 + li43*li43 + src[5]*src[5];
902  dst(3,0) = li51*li54 + li41*src[9];
903  dst(3,1) = li52*li54 + li42*src[9];
904  dst(3,2) = li53*li54 + li43*src[9];
905  dst(3,3) = li54*li54 + src[9]*src[9];
906  dst(4,0) = li51*src[14];
907  dst(4,1) = li52*src[14];
908  dst(4,2) = li53*src[14];
909  dst(4,3) = li54*src[14];
910  dst(4,4) = src[14]*src[14];
911  }
912  };
913  /// struct to obtain the inverse from a Cholesky decomposition (N = 4)
914  template<class F, class M> struct _inverter<F,4,M>
915  {
916  /// method to do the inversion
917  inline void operator()(M& dst, const F* src) const
918  {
919  const F li21 = -src[1] * src[0] * src[2];
920  const F li32 = -src[4] * src[2] * src[5];
921  const F li31 = (src[1] * src[4] * src[2] - src[3]) * src[0] * src[5];
922  const F li43 = -src[8] * src[9] * src[5];
923  const F li42 = (src[4] * src[8] * src[5] - src[7]) * src[2] * src[9];
924  const F li41 = (-src[1] * src[4] * src[8] * src[2] * src[5] +
925  src[1] * src[7] * src[2] + src[3] * src[8] * src[5] - src[6]) * src[0] * src[9];
926 
927  dst(0,0) = li41*li41 + li31*li31 + li21*li21 + src[0]*src[0];
928  dst(1,0) = li41*li42 + li31*li32 + li21*src[2];
929  dst(1,1) = li42*li42 + li32*li32 + src[2]*src[2];
930  dst(2,0) = li41*li43 + li31*src[5];
931  dst(2,1) = li42*li43 + li32*src[5];
932  dst(2,2) = li43*li43 + src[5]*src[5];
933  dst(3,0) = li41*src[9];
934  dst(3,1) = li42*src[9];
935  dst(3,2) = li43*src[9];
936  dst(3,3) = src[9]*src[9];
937  }
938  };
939  /// struct to obtain the inverse from a Cholesky decomposition (N = 3)
940  template<class F, class M> struct _inverter<F,3,M>
941  {
942  /// method to do the inversion
943  inline void operator()(M& dst, const F* src) const
944  {
945  const F li21 = -src[1] * src[0] * src[2];
946  const F li32 = -src[4] * src[2] * src[5];
947  const F li31 = (src[1] * src[4] * src[2] - src[3]) * src[0] * src[5];
948 
949  dst(0,0) = li31*li31 + li21*li21 + src[0]*src[0];
950  dst(1,0) = li31*li32 + li21*src[2];
951  dst(1,1) = li32*li32 + src[2]*src[2];
952  dst(2,0) = li31*src[5];
953  dst(2,1) = li32*src[5];
954  dst(2,2) = src[5]*src[5];
955  }
956  };
957  /// struct to obtain the inverse from a Cholesky decomposition (N = 2)
958  template<class F, class M> struct _inverter<F,2,M>
959  {
960  /// method to do the inversion
961  inline void operator()(M& dst, const F* src) const
962  {
963  const F li21 = -src[1] * src[0] * src[2];
964 
965  dst(0,0) = li21*li21 + src[0]*src[0];
966  dst(1,0) = li21*src[2];
967  dst(1,1) = src[2]*src[2];
968  }
969  };
970  /// struct to obtain the inverse from a Cholesky decomposition (N = 1)
971  template<class F, class M> struct _inverter<F,1,M>
972  {
973  /// method to do the inversion
974  inline void operator()(M& dst, const F* src) const
975  {
976  dst(0,0) = src[0]*src[0];
977  }
978  };
979  /// struct to obtain the inverse from a Cholesky decomposition (N = 0)
980  template<class F, class M> struct _inverter<F,0,M>
981  {
982  private:
983  _inverter();
984  void operator()(M& dst, const F* src) const;
985  };
986 
987  /// struct to solve a linear system using its Cholesky decomposition (N=6)
988  template<class F, class V> struct _solver<F,6,V>
989  {
990  /// method to solve the linear system
991  void operator()(V& rhs, const F* l) const
992  {
993  // solve Ly = rhs
994  const F y0 = rhs[0] * l[0];
995  const F y1 = (rhs[1]-l[1]*y0)*l[2];
996  const F y2 = (rhs[2]-(l[3]*y0+l[4]*y1))*l[5];
997  const F y3 = (rhs[3]-(l[6]*y0+l[7]*y1+l[8]*y2))*l[9];
998  const F y4 = (rhs[4]-(l[10]*y0+l[11]*y1+l[12]*y2+l[13]*y3))*l[14];
999  const F y5 = (rhs[5]-(l[15]*y0+l[16]*y1+l[17]*y2+l[18]*y3+l[19]*y4))*l[20];
1000  // solve L^Tx = y, and put x into rhs
1001  rhs[5] = y5 * l[20];
1002  rhs[4] = (y4-l[19]*rhs[5])*l[14];
1003  rhs[3] = (y3-(l[18]*rhs[5]+l[13]*rhs[4]))*l[9];
1004  rhs[2] = (y2-(l[17]*rhs[5]+l[12]*rhs[4]+l[8]*rhs[3]))*l[5];
1005  rhs[1] = (y1-(l[16]*rhs[5]+l[11]*rhs[4]+l[7]*rhs[3]+l[4]*rhs[2]))*l[2];
1006  rhs[0] = (y0-(l[15]*rhs[5]+l[10]*rhs[4]+l[6]*rhs[3]+l[3]*rhs[2]+l[1]*rhs[1]))*l[0];
1007  }
1008  };
1009  /// struct to solve a linear system using its Cholesky decomposition (N=5)
1010  template<class F, class V> struct _solver<F,5,V>
1011  {
1012  /// method to solve the linear system
1013  void operator()(V& rhs, const F* l) const
1014  {
1015  // solve Ly = rhs
1016  const F y0 = rhs[0] * l[0];
1017  const F y1 = (rhs[1]-l[1]*y0)*l[2];
1018  const F y2 = (rhs[2]-(l[3]*y0+l[4]*y1))*l[5];
1019  const F y3 = (rhs[3]-(l[6]*y0+l[7]*y1+l[8]*y2))*l[9];
1020  const F y4 = (rhs[4]-(l[10]*y0+l[11]*y1+l[12]*y2+l[13]*y3))*l[14];
1021  // solve L^Tx = y, and put x into rhs
1022  rhs[4] = (y4)*l[14];
1023  rhs[3] = (y3-(l[13]*rhs[4]))*l[9];
1024  rhs[2] = (y2-(l[12]*rhs[4]+l[8]*rhs[3]))*l[5];
1025  rhs[1] = (y1-(l[11]*rhs[4]+l[7]*rhs[3]+l[4]*rhs[2]))*l[2];
1026  rhs[0] = (y0-(l[10]*rhs[4]+l[6]*rhs[3]+l[3]*rhs[2]+l[1]*rhs[1]))*l[0];
1027  }
1028  };
1029  /// struct to solve a linear system using its Cholesky decomposition (N=4)
1030  template<class F, class V> struct _solver<F,4,V>
1031  {
1032  /// method to solve the linear system
1033  void operator()(V& rhs, const F* l) const
1034  {
1035  // solve Ly = rhs
1036  const F y0 = rhs[0] * l[0];
1037  const F y1 = (rhs[1]-l[1]*y0)*l[2];
1038  const F y2 = (rhs[2]-(l[3]*y0+l[4]*y1))*l[5];
1039  const F y3 = (rhs[3]-(l[6]*y0+l[7]*y1+l[8]*y2))*l[9];
1040  // solve L^Tx = y, and put x into rhs
1041  rhs[3] = (y3)*l[9];
1042  rhs[2] = (y2-(l[8]*rhs[3]))*l[5];
1043  rhs[1] = (y1-(l[7]*rhs[3]+l[4]*rhs[2]))*l[2];
1044  rhs[0] = (y0-(l[6]*rhs[3]+l[3]*rhs[2]+l[1]*rhs[1]))*l[0];
1045  }
1046  };
1047  /// struct to solve a linear system using its Cholesky decomposition (N=3)
1048  template<class F, class V> struct _solver<F,3,V>
1049  {
1050  /// method to solve the linear system
1051  void operator()(V& rhs, const F* l) const
1052  {
1053  // solve Ly = rhs
1054  const F y0 = rhs[0] * l[0];
1055  const F y1 = (rhs[1]-l[1]*y0)*l[2];
1056  const F y2 = (rhs[2]-(l[3]*y0+l[4]*y1))*l[5];
1057  // solve L^Tx = y, and put x into rhs
1058  rhs[2] = (y2)*l[5];
1059  rhs[1] = (y1-(l[4]*rhs[2]))*l[2];
1060  rhs[0] = (y0-(l[3]*rhs[2]+l[1]*rhs[1]))*l[0];
1061  }
1062  };
1063  /// struct to solve a linear system using its Cholesky decomposition (N=2)
1064  template<class F, class V> struct _solver<F,2,V>
1065  {
1066  /// method to solve the linear system
1067  void operator()(V& rhs, const F* l) const
1068  {
1069  // solve Ly = rhs
1070  const F y0 = rhs[0] * l[0];
1071  const F y1 = (rhs[1]-l[1]*y0)*l[2];
1072  // solve L^Tx = y, and put x into rhs
1073  rhs[1] = (y1)*l[2];
1074  rhs[0] = (y0-(l[1]*rhs[1]))*l[0];
1075  }
1076  };
1077  /// struct to solve a linear system using its Cholesky decomposition (N=1)
1078  template<class F, class V> struct _solver<F,1,V>
1079  {
1080  /// method to solve the linear system
1081  void operator()(V& rhs, const F* l) const
1082  {
1083  // solve LL^T x = rhs, put y into rhs
1084  rhs[0] *= l[0] * l[0];
1085  }
1086  };
1087  /// struct to solve a linear system using its Cholesky decomposition (N=0)
1088  template<class F, class V> struct _solver<F,0,V>
1089  {
1090  private:
1091  _solver();
1092  void operator()(V& rhs, const F* l) const;
1093  };
1094 }
1095 
1096 
1097  } // namespace Math
1098 
1099 } // namespace ROOT
1100 
1101 #endif // ROOT_Math_CHOLESKYDECOMP
1102 
bool ok() const
returns true if decomposition was successful
bool getL(G *m) const
obtain the decomposed matrix L
const G operator()(unsigned i, unsigned j) const
read access to elements (make sure that j <= i)
void operator()(V &rhs, const F *l) const
method to solve the linear system
Namespace for new ROOT classes and functions.
Definition: TFoamSampler.h:19
bool Invert(G *m) const
place the inverse into m
CholeskyDecompGenDim(unsigned N, G *m)
perform a Cholesky decomposition
bool Invert(M &m) const
place the inverse into m
bool fOk
flag indicating a successful decomposition
bool getL(M &m) const
obtain the decomposed matrix L
unsigned fN
dimensionality dimensionality of the problem
bool operator()(F *dst, const M &src, unsigned N) const
method to do the decomposition
void operator()(M &dst, const F *src) const
method to do the inversion
CholeskyDecomp(G *m)
perform a Cholesky decomposition
bool operator()(F *dst, const M &src) const
method to do the decomposition
adapter for packed arrays (to SMatrix indexing conventions)
void operator()(M &dst, const F *src) const
method to do the inversion
void operator()(M &dst, const F *src) const
method to do the inversion
bool getLi(M &m) const
obtain the inverse of the decomposed matrix L
void operator()(M &dst, const F *src, unsigned N) const
method to do the inversion
void operator()(V &rhs, const F *l) const
method to solve the linear system
void operator()(M &dst, const F *src) const
method to do the inversion
void operator()(V &rhs, const F *l) const
method to solve the linear system
bool fOk
flag indicating a successful decomposition
void operator()(M &dst, const F *src) const
method to do the inversion
class to compute the Cholesky decomposition of a matrix
double sqrt(double)
bool operator()(F *dst, const M &src) const
method to do the decomposition
bool getL(G *m) const
obtain the decomposed matrix L
bool getLi(G *m) const
obtain the inverse of the decomposed matrix L
void operator()(V &rhs, const F *l) const
method to solve the linear system
#define N
void operator()(M &dst, const F *src) const
method to do the inversion
bool getLi(G *m) const
obtain the inverse of the decomposed matrix L
void operator()(V &rhs, const F *l) const
method to solve the linear system
struct to solve a linear system using its Cholesky decomposition (generalised dimensionality) ...
bool operator()(F *dst, const M &src) const
method to do the decomposition
struct to do a Cholesky decomposition (general dimensionality)
struct to obtain the inverse from a Cholesky decomposition
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t Int_t m
Definition: TRolke.cxx:637
bool getLi(M &m) const
obtain the inverse of the decomposed matrix L
G & operator()(unsigned i, unsigned j)
write access to elements (make sure that j <= i)
constexpr Double_t G()
Gravitational constant in: .
Definition: TMath.h:139
bool Invert(M &m) const
place the inverse into m
struct to do a Cholesky decomposition
bool operator()(F *dst, const M &src) const
method to do the decomposition
bool operator()(F *dst, const M &src) const
method to do the decomposition
bool Solve(V &rhs) const
solves a linear system for the given right hand side
void operator()(M &dst, const F *src) const
method to do the inversion
class to compute the Cholesky decomposition of a matrix
CholeskyDecompGenDim(unsigned N, const M &m)
perform a Cholesky decomposition
bool getL(M &m) const
obtain the decomposed matrix L
CholeskyDecomp(const M &m)
perform a Cholesky decomposition
bool operator()(F *dst, const M &src) const
method to do the decomposition
F * fL
lower triangular matrix L
void operator()(V &rhs, const F *l) const
method to solve the linear system
void operator()(V &rhs, const F *l) const
method to solve the linear system
bool Invert(G *m) const
place the inverse into m
Namespace for new Math classes and functions.
bool Solve(V &rhs) const
solves a linear system for the given right hand side
bool operator()(F *dst, const M &src) const
method to do the decomposition
bool ok() const
returns true if decomposition was successful
struct to solve a linear system using its Cholesky decomposition
void operator()(V &rhs, const F *l, unsigned N) const
method to solve the linear system
struct to obtain the inverse from a Cholesky decomposition (general dimensionality) ...