Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
30namespace ROOT {
31
32 namespace Math {
33
34/// helpers for CholeskyDecomp
35namespace 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 */
76template<class F, unsigned N> class CholeskyDecomp
77{
78private:
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;
85public:
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 */
310template<class F> class CholeskyDecompGenDim
311{
312private:
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 */
320 /// flag indicating a successful decomposition
321 bool fOk;
322public:
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
519namespace 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:
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:
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:
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
#define N
class to compute the Cholesky decomposition of a matrix
bool fOk
flag indicating a successful decomposition
bool getL(G *m) const
obtain the decomposed matrix L
bool getL(M &m) const
obtain the decomposed matrix L
CholeskyDecompGenDim(unsigned N, G *m)
perform a Cholesky decomposition
bool Invert(G *m) const
place the inverse into m
bool Invert(M &m) const
place the inverse into m
bool Solve(V &rhs) const
solves a linear system for the given right hand side
bool getLi(M &m) const
obtain the inverse of the decomposed matrix L
F * fL
lower triangular matrix L
unsigned fN
dimensionality dimensionality of the problem
CholeskyDecompGenDim(unsigned N, const M &m)
perform a Cholesky decomposition
bool getLi(G *m) const
obtain the inverse of the decomposed matrix L
bool ok() const
returns true if decomposition was successful
adapter for packed arrays (to SMatrix indexing conventions)
G & operator()(unsigned i, unsigned j)
write access to elements (make sure that j <= i)
const G operator()(unsigned i, unsigned j) const
read access to elements (make sure that j <= i)
class to compute the Cholesky decomposition of a matrix
bool fOk
flag indicating a successful decomposition
bool ok() const
returns true if decomposition was successful
bool getL(G *m) const
obtain the decomposed matrix L
bool Invert(G *m) const
place the inverse into m
CholeskyDecomp(const M &m)
perform a Cholesky decomposition
bool getLi(G *m) const
obtain the inverse of the decomposed matrix L
bool getL(M &m) const
obtain the decomposed matrix L
bool Solve(V &rhs) const
solves a linear system for the given right hand side
F fL[N *(N+1)/2]
lower triangular matrix L
CholeskyDecomp(G *m)
perform a Cholesky decomposition
bool getLi(M &m) const
obtain the inverse of the decomposed matrix L
bool Invert(M &m) const
place the inverse into m
#define F(x, y, z)
#define G(x, y, z)
Namespace for new Math classes and functions.
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
struct to do a Cholesky decomposition (general dimensionality)
bool operator()(F *dst, const M &src, unsigned N) const
method to do the 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 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 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
struct to do a Cholesky decomposition
bool operator()(F *dst, const M &src) const
method to do the decomposition
struct to obtain the inverse from a Cholesky decomposition (general dimensionality)
void operator()(M &dst, const F *src, unsigned N) const
method to do the inversion
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
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
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
struct to obtain the inverse from a Cholesky decomposition
void operator()(M &dst, const F *src) const
method to do the inversion
struct to solve a linear system using its Cholesky decomposition (generalised dimensionality)
void operator()(V &rhs, const F *l, unsigned N) const
method to solve the linear system
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
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
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
struct to solve a linear system using its Cholesky decomposition
void operator()(V &rhs, const F *l) const
method to solve the linear system
auto * m
Definition textangle.C:8
auto * l
Definition textangle.C:4
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345