EIC code displayed by LXR master/​include/​root/​Math/​CholeskyDecomp.h	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2025-09-17 09:13:56

 // @(#)root/smatrix:$Id$
 // Author: M. Schiller    2009
 
 #ifndef ROOT_Math_CholeskyDecomp
 #define ROOT_Math_CholeskyDecomp
 
 /** @file
  * header file containing the templated implementation of matrix inversion
  * routines for use with ROOT's SMatrix classes (symmetric positive
  * definite case)
  *
  * @author Manuel Schiller
  * @date Aug 29 2008
  *    initial release inside LHCb
  * @date May 7 2009
  * factored code to provide a nice Cholesky decomposition class, along
  * with separate methods for solving a single linear system and to
  * obtain the inverse matrix from the decomposition
  * @date July 15th 2013
  *    provide a version of that class which works if the dimension of the
  *    problem is only known at run time
  * @date September 30th 2013
  *    provide routines to access the result of the decomposition L and its
  *    inverse
  * @date January 20th 2024
  *    CholeskyDecompGenDim allowed copying and moving, but did not provide copy
  *    and move constructors or assignment operators, resulting in incorrect
  *    code and potential memory corruption, if somebody tried to copy/move a
  *    CholeskyDecompGenDim instance. Since we're in a C++17 world now, use the
  *    excellent movable std::vector<F> to get rid of explicit memory management
  *    in CholeskyDecompGenDim, and thus provide correct code for copy/move
  *    construction and assignment. Moreover, the provided releaseWorkspace
  *    method allows reuse of the allocated vector inside a CholeskyDecompGenDim
  *    to reduce the cost of repeated matrix decomposition.
  */
 
 #include <cmath>
 #include <vector>
 #include <algorithm>
 
 namespace ROOT {
 
 namespace Math {
 
 /// helpers for CholeskyDecomp
 namespace CholeskyDecompHelpers {
 // forward decls
 template <class F, class M>
 struct _decomposerGenDim;
 template <class F, unsigned N, class M>
 struct _decomposer;
 template <class F, class M>
 struct _inverterGenDim;
 template <class F, unsigned N, class M>
 struct _inverter;
 template <class F, class V>
 struct _solverGenDim;
 template <class F, unsigned N, class V>
 struct _solver;
 template <typename G>
 class PackedArrayAdapter;
 } // namespace CholeskyDecompHelpers
 
 /// class to compute the Cholesky decomposition of a matrix
 /** class to compute the Cholesky decomposition of a symmetric
  * positive definite matrix
  *
  * provides routines to check if the decomposition succeeded (i.e. if
  * matrix is positive definite and non-singular), to solve a linear
  * system for the given matrix and to obtain its inverse
  *
  * the actual functionality is implemented in templated helper
  * classes which have specializations for dimensions N = 1 to 6
  * to achieve a gain in speed for common matrix sizes
  *
  * usage example:
  * @code
  * // let m be a symmetric positive definite SMatrix (use type float
  * // for internal computations, matrix size is 4x4)
  * CholeskyDecomp<float, 4> decomp(m);
  * // check if the decomposition succeeded
  * if (!decomp) {
  *   std::cerr << "decomposition failed!" << std::endl;
  * } else {
  *   // let rhs be a vector; we seek a vector x such that m * x = rhs
  *   decomp.Solve(rhs);
  *   // rhs now contains the solution we are looking for
  *
  *   // obtain the inverse of m, put it into m itself
  *   decomp.Invert(m);
  * }
  * @endcode
  */
 template <class F, unsigned N>
 class CholeskyDecomp {
 private:
    /// lower triangular matrix L
    /** lower triangular matrix L, packed storage, with diagonal
     * elements pre-inverted */
    F fL[N * (N + 1) / 2];
    /// flag indicating a successful decomposition
    bool fOk;
 
 public:
    /// perform a Cholesky decomposition
    /** perform a Cholesky decomposition of a symmetric positive
     * definite matrix m
     *
     * this is the constructor to uses with an SMatrix (and objects
     * that behave like an SMatrix in terms of using
     * operator()(int i, int j) for access to elements)
     */
    template <class M>
    CholeskyDecomp(const M &m) : fOk(false)
    {
       using CholeskyDecompHelpers::_decomposer;
       fOk = _decomposer<F, N, M>()(fL, m);
    }
 
    /// perform a Cholesky decomposition
    /** perform a Cholesky decomposition of a symmetric positive
     * definite matrix m
     *
     * this is the constructor to use in special applications where
     * plain arrays are used
     *
     * NOTE: the matrix is given in packed representation, matrix
     * element m(i,j) (j <= i) is supposed to be in array element
     * (i * (i + 1)) / 2 + j
     */
    template <typename G>
    CholeskyDecomp(G *m) : fOk(false)
    {
       using CholeskyDecompHelpers::_decomposer;
       using CholeskyDecompHelpers::PackedArrayAdapter;
       fOk = _decomposer<F, N, PackedArrayAdapter<G>>()(fL, PackedArrayAdapter<G>(m));
    }
 
    /// returns true if decomposition was successful
    /** @returns true if decomposition was successful */
    bool ok() const { return fOk; }
    /// returns true if decomposition was successful
    /** @returns true if decomposition was successful */
    operator bool() const { return fOk; }
 
    /** @brief solves a linear system for the given right hand side
     *
     * Note that you can use both SVector classes and plain arrays for
     * rhs. (Make sure that the sizes match!). It will work with any vector
     * implementing the operator [i]
     *
     * @returns if the decomposition was successful
     */
    template <class V>
    bool Solve(V &rhs) const
    {
       using CholeskyDecompHelpers::_solver;
       if (fOk)
          _solver<F, N, V>()(rhs, fL);
       return fOk;
    }
 
    /** @brief place the inverse into m
     *
     * This is the method to use with an SMatrix.
     *
     * @returns if the decomposition was successful
     */
    template <class M>
    bool Invert(M &m) const
    {
       using CholeskyDecompHelpers::_inverter;
       if (fOk)
          _inverter<F, N, M>()(m, fL);
       return fOk;
    }
 
    /** @brief place the inverse into m
     *
     * This is the method to use with a plain array.
     *
     * @returns if the decomposition was successful
     *
     * NOTE: the matrix is given in packed representation, matrix
     * element m(i,j) (j <= i) is supposed to be in array element
     * (i * (i + 1)) / 2 + j
     */
    template <typename G>
    bool Invert(G *m) const
    {
       using CholeskyDecompHelpers::_inverter;
       using CholeskyDecompHelpers::PackedArrayAdapter;
       if (fOk) {
          PackedArrayAdapter<G> adapted(m);
          _inverter<F, N, PackedArrayAdapter<G>>()(adapted, fL);
       }
       return fOk;
    }
 
    /** @brief obtain the decomposed matrix L
     *
     * This is the method to use with a plain array.
     *
     * @returns if the decomposition was successful
     */
    template <class M>
    bool getL(M &m) const
    {
       if (!fOk)
          return false;
       for (unsigned i = 0; i < N; ++i) {
          // zero upper half of matrix
          for (unsigned j = i + 1; j < N; ++j)
             m(i, j) = F(0);
          // copy the rest
          for (unsigned j = 0; j <= i; ++j)
             m(i, j) = fL[i * (i + 1) / 2 + j];
          // adjust the diagonal - we save 1/L(i, i) in that position, so
          // convert to what caller expects
          m(i, i) = F(1) / m(i, i);
       }
       return true;
    }
 
    /** @brief obtain the decomposed matrix L
     *
     * @returns if the decomposition was successful
     *
     * NOTE: the matrix is given in packed representation, matrix
     * element m(i,j) (j <= i) is supposed to be in array element
     * (i * (i + 1)) / 2 + j
     */
    template <typename G>
    bool getL(G *m) const
    {
       if (!fOk)
          return false;
       // copy L
       for (unsigned i = 0; i < (N * (N + 1)) / 2; ++i)
          m[i] = fL[i];
       // adjust diagonal - we save 1/L(i, i) in that position, so convert to
       // what caller expects
       for (unsigned i = 0; i < N; ++i)
          m[(i * (i + 1)) / 2 + i] = F(1) / fL[(i * (i + 1)) / 2 + i];
       return true;
    }
 
    /** @brief obtain the inverse of the decomposed matrix L
     *
     * This is the method to use with a plain array.
     *
     * @returns if the decomposition was successful
     */
    template <class M>
    bool getLi(M &m) const
    {
       if (!fOk)
          return false;
       for (unsigned i = 0; i < N; ++i) {
          // zero lower half of matrix
          for (unsigned j = i + 1; j < N; ++j)
             m(j, i) = F(0);
          // copy the rest
          for (unsigned j = 0; j <= i; ++j)
             m(j, i) = fL[i * (i + 1) / 2 + j];
       }
       // invert the off-diagonal part of what we just copied
       for (unsigned i = 1; i < N; ++i) {
          for (unsigned j = 0; j < i; ++j) {
             typename M::value_type tmp = F(0);
             for (unsigned k = i; k-- > j;)
                tmp -= m(k, i) * m(j, k);
             m(j, i) = tmp * m(i, i);
          }
       }
       return true;
    }
 
    /** @brief obtain the inverse of the decomposed matrix L
     *
     * @returns if the decomposition was successful
     *
     * NOTE: the matrix is given in packed representation, matrix
     * element m(j,i) (j <= i) is supposed to be in array element
     * (i * (i + 1)) / 2 + j
     */
    template <typename G>
    bool getLi(G *m) const
    {
       if (!fOk)
          return false;
       // copy L
       for (unsigned i = 0; i < (N * (N + 1)) / 2; ++i)
          m[i] = fL[i];
       // invert the off-diagonal part of what we just copied
       G *base1 = &m[1];
       for (unsigned i = 1; i < N; base1 += ++i) {
          for (unsigned j = 0; j < i; ++j) {
             G tmp = F(0);
             const G *base2 = &m[(i * (i - 1)) / 2];
             for (unsigned k = i; k-- > j; base2 -= k)
                tmp -= base1[k] * base2[j];
             base1[j] = tmp * base1[i];
          }
       }
       return true;
    }
 };
 
 /// class to compute the Cholesky decomposition of a matrix
 /** class to compute the Cholesky decomposition of a symmetric
  * positive definite matrix when the dimensionality of the problem is not known
  * at compile time
  *
  * provides routines to check if the decomposition succeeded (i.e. if
  * matrix is positive definite and non-singular), to solve a linear
  * system for the given matrix and to obtain its inverse
  *
  * the actual functionality is implemented in templated helper
  * classes which have specializations for dimensions N = 1 to 6
  * to achieve a gain in speed for common matrix sizes
  *
  * usage example:
  * @code
  * // let m be a symmetric positive definite SMatrix (use type float
  * // for internal computations, matrix size is 4x4)
  * CholeskyDecompGenDim<float> decomp(4, m);
  * // check if the decomposition succeeded
  * if (!decomp) {
  *   std::cerr << "decomposition failed!" << std::endl;
  * } else {
  *   // let rhs be a vector; we seek a vector x such that m * x = rhs
  *   decomp.Solve(rhs);
  *   // rhs now contains the solution we are looking for
  *
  *   // obtain the inverse of m, put it into m itself
  *   decomp.Invert(m);
  * }
  * @endcode
  */
 template <class F>
 class CholeskyDecompGenDim {
 private:
    /** @brief dimensionality
     * dimensionality of the problem */
    unsigned fN = 0;
    /// lower triangular matrix L
    /** lower triangular matrix L, packed storage, with diagonal
     * elements pre-inverted */
    std::vector<F> fL;
    /// flag indicating a successful decomposition
    bool fOk = false;
 
 public:
    /// perform a Cholesky decomposition
    /** perform a Cholesky decomposition of a symmetric positive
     * definite matrix m
     *
     * this is the constructor to uses with an SMatrix (and objects
     * that behave like an SMatrix in terms of using
     * operator()(int i, int j) for access to elements)
     *
     * The optional wksp parameter that can be used to avoid (re-)allocations
     * on repeated decompositions of the same size (by passing a workspace of
     * size N * (N + 1), or reusing one returned from releaseWorkspace by a
     * previous decomposition).
     */
    template <class M>
    CholeskyDecompGenDim(unsigned N, const M &m, std::vector<F> wksp = {}) : fN(N), fL(std::move(wksp))
    {
       using CholeskyDecompHelpers::_decomposerGenDim;
       if (fL.size() < fN * (fN + 1) / 2)
          fL.resize(fN * (fN + 1) / 2);
       fOk = _decomposerGenDim<F, M>()(fL.data(), m, fN);
    }
 
    /// perform a Cholesky decomposition
    /** perform a Cholesky decomposition of a symmetric positive
     * definite matrix m
     *
     * this is the constructor to use in special applications where
     * plain arrays are used
     *
     * NOTE: the matrix is given in packed representation, matrix
     * element m(i,j) (j <= i) is supposed to be in array element
     * (i * (i + 1)) / 2 + j
     *
     * The optional wksp parameter that can be used to avoid (re-)allocations
     * on repeated decompositions of the same size (by passing a workspace of
     * size N * (N + 1), or reusing one returned from releaseWorkspace by a
     * previous decomposition).
     */
    template <typename G>
    CholeskyDecompGenDim(unsigned N, G *m, std::vector<F> wksp = {}) : fN(N), fL(std::move(wksp))
    {
       using CholeskyDecompHelpers::_decomposerGenDim;
       using CholeskyDecompHelpers::PackedArrayAdapter;
       if (fL.size() < fN * (fN + 1) / 2)
          fL.resize(fN * (fN + 1) / 2);
       fOk = _decomposerGenDim<F, PackedArrayAdapter<G>>()(fL.data(), PackedArrayAdapter<G>(m), fN);
    }
 
    /// release the workspace of the decomposition for reuse
    /** release the workspace of the decomposition for reuse
     *
     * If you use CholeskyDecompGenDim repeatedly with the same size N, it is
     * possible to avoid repeatedly allocating the internal workspace. The
     * constructors take an optional wksp argument, and this routine can be
     * called to get the workspace out of a decomposition when you are done
     * with it.
     *
     * Please note that once you call releaseWorkspace, you cannot do further
     * operations on the decomposition object (and this is indicated by having
     * it return false when you call ok()).
     *
     * Here is an example snippet of (pseudo-)code which illustrates the use of
     * releaseWorkspace for inverting the covariance matrices of a number of
     * items (which must all be of the same size):
     *
     * @code
     * std::vector<float> wksp;
     * for (const auto& item: items) {
     *     CholeskyDecompGenDim decomp(item.cov().size(), item.cov(),
     *                                 std::move(wksp));
     *     if (decomp.ok()) {
     *         decomp.Invert(item.invcov());
     *     }
     *     wksp = std::move(decomp.releaseWorkspace());
     * }
     * @endcode
     *
     * In the above code snippet, only the first CholeskyDecompGenDim
     * constructor call will allocate memory. Subsequent calls in the loop will
     * reuse the buffer from the first iteration.
     */
    std::vector<F> releaseWorkspace()
    {
       fN = 0;
       fOk = false;
       return std::move(fL);
    }
 
    /// returns true if decomposition was successful
    /** @returns true if decomposition was successful */
    bool ok() const { return fOk; }
    /// returns true if decomposition was successful
    /** @returns true if decomposition was successful */
    operator bool() const { return fOk; }
 
    /** @brief solves a linear system for the given right hand side
     *
     * Note that you can use both SVector classes and plain arrays for
     * rhs. (Make sure that the sizes match!). It will work with any vector
     * implementing the operator [i]
     *
     * @returns if the decomposition was successful
     */
    template <class V>
    bool Solve(V &rhs) const
    {
       using CholeskyDecompHelpers::_solverGenDim;
       if (fOk)
          _solverGenDim<F, V>()(rhs, fL.data(), fN);
       return fOk;
    }
 
    /** @brief place the inverse into m
     *
     * This is the method to use with an SMatrix.
     *
     * @returns if the decomposition was successful
     */
    template <class M>
    bool Invert(M &m) const
    {
       using CholeskyDecompHelpers::_inverterGenDim;
       if (fOk) {
          auto wksp = fL;
          _inverterGenDim<F, M>()(m, fN, wksp.data());
       }
       return fOk;
    }
 
    /** @brief place the inverse into m
     *
     * This is the method to use with a plain array.
     *
     * @returns if the decomposition was successful
     *
     * NOTE: the matrix is given in packed representation, matrix
     * element m(i,j) (j <= i) is supposed to be in array element
     * (i * (i + 1)) / 2 + j
     */
    template <typename G>
    bool Invert(G *m) const
    {
       using CholeskyDecompHelpers::_inverterGenDim;
       using CholeskyDecompHelpers::PackedArrayAdapter;
       if (fOk) {
          auto wksp = fL;
          PackedArrayAdapter<G> adapted(m);
          _inverterGenDim<F, PackedArrayAdapter<G>>()(adapted, fN, wksp.data());
       }
       return fOk;
    }
 
    /** @brief obtain the decomposed matrix L
     *
     * This is the method to use with a plain array.
     *
     * @returns if the decomposition was successful
     */
    template <class M>
    bool getL(M &m) const
    {
       if (!fOk)
          return false;
       const F *L = fL.data();
       for (unsigned i = 0; i < fN; ++i) {
          // zero upper half of matrix
          for (unsigned j = i + 1; j < fN; ++j)
             m(i, j) = F(0);
          // copy the rest
          for (unsigned j = 0; j <= i; ++j)
             m(i, j) = L[i * (i + 1) / 2 + j];
          // adjust the diagonal - we save 1/L(i, i) in that position, so
          // convert to what caller expects
          m(i, i) = F(1) / m(i, i);
       }
       return true;
    }
 
    /** @brief obtain the decomposed matrix L
     *
     * @returns if the decomposition was successful
     *
     * NOTE: the matrix is given in packed representation, matrix
     * element m(i,j) (j <= i) is supposed to be in array element
     * (i * (i + 1)) / 2 + j
     */
    template <typename G>
    bool getL(G *m) const
    {
       if (!fOk)
          return false;
       const F *L = fL.data();
       // copy L
       for (unsigned i = 0; i < (fN * (fN + 1)) / 2; ++i)
          m[i] = L[i];
       // adjust diagonal - we save 1/L(i, i) in that position, so convert to
       // what caller expects
       for (unsigned i = 0; i < fN; ++i)
          m[(i * (i + 1)) / 2 + i] = F(1) / L[(i * (i + 1)) / 2 + i];
       return true;
    }
 
    /** @brief obtain the inverse of the decomposed matrix L
     *
     * This is the method to use with a plain array.
     *
     * @returns if the decomposition was successful
     */
    template <class M>
    bool getLi(M &m) const
    {
       if (!fOk)
          return false;
       const F *L = fL.data();
       for (unsigned i = 0; i < fN; ++i) {
          // zero lower half of matrix
          for (unsigned j = i + 1; j < fN; ++j)
             m(j, i) = F(0);
          // copy the rest
          for (unsigned j = 0; j <= i; ++j)
             m(j, i) = L[i * (i + 1) / 2 + j];
       }
       // invert the off-diagonal part of what we just copied
       for (unsigned i = 1; i < fN; ++i) {
          for (unsigned j = 0; j < i; ++j) {
             typename M::value_type tmp = F(0);
             for (unsigned k = i; k-- > j;)
                tmp -= m(k, i) * m(j, k);
             m(j, i) = tmp * m(i, i);
          }
       }
       return true;
    }
 
    /** @brief obtain the inverse of the decomposed matrix L
     *
     * @returns if the decomposition was successful
     *
     * NOTE: the matrix is given in packed representation, matrix
     * element m(j,i) (j <= i) is supposed to be in array element
     * (i * (i + 1)) / 2 + j
     */
    template <typename G>
    bool getLi(G *m) const
    {
       if (!fOk)
          return false;
       const F *L = fL.data();
       // copy L
       for (unsigned i = 0; i < (fN * (fN + 1)) / 2; ++i)
          m[i] = L[i];
       // invert the off-diagonal part of what we just copied
       G *base1 = &m[1];
       for (unsigned i = 1; i < fN; base1 += ++i) {
          for (unsigned j = 0; j < i; ++j) {
             G tmp = F(0);
             const G *base2 = &m[(i * (i - 1)) / 2];
             for (unsigned k = i; k-- > j; base2 -= k)
                tmp -= base1[k] * base2[j];
             base1[j] = tmp * base1[i];
          }
       }
       return true;
    }
 };
 
 namespace CholeskyDecompHelpers {
 /// adapter for packed arrays (to SMatrix indexing conventions)
 template <typename G>
 class PackedArrayAdapter {
 private:
    G *fArr; ///< pointer to first array element
 public:
    /// constructor
    PackedArrayAdapter(G *arr) : fArr(arr) {}
    /// read access to elements (make sure that j <= i)
    const G operator()(unsigned i, unsigned j) const { return fArr[((i * (i + 1)) / 2) + j]; }
    /// write access to elements (make sure that j <= i)
    G &operator()(unsigned i, unsigned j) { return fArr[((i * (i + 1)) / 2) + j]; }
 };
 /// struct to do a Cholesky decomposition (general dimensionality)
 template <class F, class M>
 struct _decomposerGenDim {
    /// method to do the decomposition
    /** @returns if the decomposition was successful */
    bool operator()(F *dst, const M &src, unsigned N) const
    {
       // perform Cholesky decomposition of matrix: M = L L^T
       // only thing that can go wrong: trying to take square
       // root of negative number or zero (matrix is
       // ill-conditioned or singular in these cases)
 
       // element L(i,j) is at array position (i * (i+1)) / 2 + j
 
       // quirk: we may need to invert L later anyway, so we can
       // invert elements on diagonal straight away (we only
       // ever need their reciprocals!)
 
       // cache starting address of rows of L for speed reasons
       F *base1 = &dst[0];
       for (unsigned i = 0; i < N; base1 += ++i) {
          F tmpdiag = F(0.0); // for element on diagonal
          // calculate off-diagonal elements
          F *base2 = &dst[0];
          for (unsigned j = 0; j < i; base2 += ++j) {
             F tmp = src(i, j);
             for (unsigned k = j; k--;)
                tmp -= base1[k] * base2[k];
             base1[j] = tmp *= base2[j];
             // keep track of contribution to element on diagonal
             tmpdiag += tmp * tmp;
          }
          // keep truncation error small
          tmpdiag = src(i, i) - tmpdiag;
          // check if positive definite
          if (tmpdiag <= F(0.0))
             return false;
          else
             base1[i] = std::sqrt(F(1.0) / tmpdiag);
       }
       return true;
    }
 };
 
 /// struct to do a Cholesky decomposition
 template <class F, unsigned N, class M>
 struct _decomposer {
    /// method to do the decomposition
    /** @returns if the decomposition was successful */
    bool operator()(F *dst, const M &src) const { return _decomposerGenDim<F, M>()(dst, src, N); }
 };
 
 /// struct to obtain the inverse from a Cholesky decomposition (general dimensionality)
 template <class F, class M>
 struct _inverterGenDim {
    /// method to do the inversion
    void operator()(M &dst, unsigned N, F *wksp) const
    {
       // invert off-diagonal part of matrix
       F *base1 = &wksp[1];
       for (unsigned i = 1; i < N; base1 += ++i) {
          for (unsigned j = 0; j < i; ++j) {
             F tmp = F(0.0);
             const F *base2 = &wksp[(i * (i - 1)) / 2];
             for (unsigned k = i; k-- > j; base2 -= k)
                tmp -= base1[k] * base2[j];
             base1[j] = tmp * base1[i];
          }
       }
 
       // Li = L^(-1) formed, now calculate M^(-1) = Li^T Li
       for (unsigned i = N; i--;) {
          for (unsigned j = i + 1; j--;) {
             F tmp = F(0.0);
             base1 = &wksp[(N * (N - 1)) / 2];
             for (unsigned k = N; k-- > i; base1 -= k)
                tmp += base1[i] * base1[j];
             dst(i, j) = tmp;
          }
       }
    }
 };
 
 /// struct to obtain the inverse from a Cholesky decomposition
 template <class F, unsigned N, class M>
 struct _inverter {
    /// method to do the inversion
    void operator()(M &dst, const F *src) const
    {
       // make working copy
       F wksp[N * (N + 1) / 2];
       std::copy(src, src + ((N * (N + 1)) / 2), wksp);
       return _inverterGenDim<F, M>()(dst, N, wksp);
    }
 };
 
 /// struct to solve a linear system using its Cholesky decomposition (generalised dimensionality)
 template <class F, class V>
 struct _solverGenDim {
    /// method to solve the linear system
    void operator()(V &rhs, const F *l, unsigned N) const
    {
       // solve Ly = rhs
       for (unsigned k = 0; k < N; ++k) {
          const unsigned base = (k * (k + 1)) / 2;
          F sum = F(0.0);
          for (unsigned i = k; i--;)
             sum += rhs[i] * l[base + i];
          // elements on diagonal are pre-inverted!
          rhs[k] = (rhs[k] - sum) * l[base + k];
       }
       // solve L^Tx = y
       for (unsigned k = N; k--;) {
          F sum = F(0.0);
          for (unsigned i = N; --i > k;)
             sum += rhs[i] * l[(i * (i + 1)) / 2 + k];
          // elements on diagonal are pre-inverted!
          rhs[k] = (rhs[k] - sum) * l[(k * (k + 1)) / 2 + k];
       }
    }
 };
 
 /// struct to solve a linear system using its Cholesky decomposition
 template <class F, unsigned N, class V>
 struct _solver {
    /// method to solve the linear system
    void operator()(V &rhs, const F *l) const { _solverGenDim<F, V>()(rhs, l, N); }
 };
 
 /// struct to do a Cholesky decomposition (specialized, N = 6)
 template <class F, class M>
 struct _decomposer<F, 6, M> {
    /// method to do the decomposition
    bool operator()(F *dst, const M &src) const
    {
       if (src(0, 0) <= F(0.0))
          return false;
       dst[0] = std::sqrt(F(1.0) / src(0, 0));
       dst[1] = src(1, 0) * dst[0];
       dst[2] = src(1, 1) - dst[1] * dst[1];
       if (dst[2] <= F(0.0))
          return false;
       else
          dst[2] = std::sqrt(F(1.0) / dst[2]);
       dst[3] = src(2, 0) * dst[0];
       dst[4] = (src(2, 1) - dst[1] * dst[3]) * dst[2];
       dst[5] = src(2, 2) - (dst[3] * dst[3] + dst[4] * dst[4]);
       if (dst[5] <= F(0.0))
          return false;
       else
          dst[5] = std::sqrt(F(1.0) / dst[5]);
       dst[6] = src(3, 0) * dst[0];
       dst[7] = (src(3, 1) - dst[1] * dst[6]) * dst[2];
       dst[8] = (src(3, 2) - dst[3] * dst[6] - dst[4] * dst[7]) * dst[5];
       dst[9] = src(3, 3) - (dst[6] * dst[6] + dst[7] * dst[7] + dst[8] * dst[8]);
       if (dst[9] <= F(0.0))
          return false;
       else
          dst[9] = std::sqrt(F(1.0) / dst[9]);
       dst[10] = src(4, 0) * dst[0];
       dst[11] = (src(4, 1) - dst[1] * dst[10]) * dst[2];
       dst[12] = (src(4, 2) - dst[3] * dst[10] - dst[4] * dst[11]) * dst[5];
       dst[13] = (src(4, 3) - dst[6] * dst[10] - dst[7] * dst[11] - dst[8] * dst[12]) * dst[9];
       dst[14] = src(4, 4) - (dst[10] * dst[10] + dst[11] * dst[11] + dst[12] * dst[12] + dst[13] * dst[13]);
       if (dst[14] <= F(0.0))
          return false;
       else
          dst[14] = std::sqrt(F(1.0) / dst[14]);
       dst[15] = src(5, 0) * dst[0];
       dst[16] = (src(5, 1) - dst[1] * dst[15]) * dst[2];
       dst[17] = (src(5, 2) - dst[3] * dst[15] - dst[4] * dst[16]) * dst[5];
       dst[18] = (src(5, 3) - dst[6] * dst[15] - dst[7] * dst[16] - dst[8] * dst[17]) * dst[9];
       dst[19] = (src(5, 4) - dst[10] * dst[15] - dst[11] * dst[16] - dst[12] * dst[17] - dst[13] * dst[18]) * dst[14];
       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]);
       if (dst[20] <= F(0.0))
          return false;
       else
          dst[20] = std::sqrt(F(1.0) / dst[20]);
       return true;
    }
 };
 /// struct to do a Cholesky decomposition (specialized, N = 5)
 template <class F, class M>
 struct _decomposer<F, 5, M> {
    /// method to do the decomposition
    bool operator()(F *dst, const M &src) const
    {
       if (src(0, 0) <= F(0.0))
          return false;
       dst[0] = std::sqrt(F(1.0) / src(0, 0));
       dst[1] = src(1, 0) * dst[0];
       dst[2] = src(1, 1) - dst[1] * dst[1];
       if (dst[2] <= F(0.0))
          return false;
       else
          dst[2] = std::sqrt(F(1.0) / dst[2]);
       dst[3] = src(2, 0) * dst[0];
       dst[4] = (src(2, 1) - dst[1] * dst[3]) * dst[2];
       dst[5] = src(2, 2) - (dst[3] * dst[3] + dst[4] * dst[4]);
       if (dst[5] <= F(0.0))
          return false;
       else
          dst[5] = std::sqrt(F(1.0) / dst[5]);
       dst[6] = src(3, 0) * dst[0];
       dst[7] = (src(3, 1) - dst[1] * dst[6]) * dst[2];
       dst[8] = (src(3, 2) - dst[3] * dst[6] - dst[4] * dst[7]) * dst[5];
       dst[9] = src(3, 3) - (dst[6] * dst[6] + dst[7] * dst[7] + dst[8] * dst[8]);
       if (dst[9] <= F(0.0))
          return false;
       else
          dst[9] = std::sqrt(F(1.0) / dst[9]);
       dst[10] = src(4, 0) * dst[0];
       dst[11] = (src(4, 1) - dst[1] * dst[10]) * dst[2];
       dst[12] = (src(4, 2) - dst[3] * dst[10] - dst[4] * dst[11]) * dst[5];
       dst[13] = (src(4, 3) - dst[6] * dst[10] - dst[7] * dst[11] - dst[8] * dst[12]) * dst[9];
       dst[14] = src(4, 4) - (dst[10] * dst[10] + dst[11] * dst[11] + dst[12] * dst[12] + dst[13] * dst[13]);
       if (dst[14] <= F(0.0))
          return false;
       else
          dst[14] = std::sqrt(F(1.0) / dst[14]);
       return true;
    }
 };
 /// struct to do a Cholesky decomposition (specialized, N = 4)
 template <class F, class M>
 struct _decomposer<F, 4, M> {
    /// method to do the decomposition
    bool operator()(F *dst, const M &src) const
    {
       if (src(0, 0) <= F(0.0))
          return false;
       dst[0] = std::sqrt(F(1.0) / src(0, 0));
       dst[1] = src(1, 0) * dst[0];
       dst[2] = src(1, 1) - dst[1] * dst[1];
       if (dst[2] <= F(0.0))
          return false;
       else
          dst[2] = std::sqrt(F(1.0) / dst[2]);
       dst[3] = src(2, 0) * dst[0];
       dst[4] = (src(2, 1) - dst[1] * dst[3]) * dst[2];
       dst[5] = src(2, 2) - (dst[3] * dst[3] + dst[4] * dst[4]);
       if (dst[5] <= F(0.0))
          return false;
       else
          dst[5] = std::sqrt(F(1.0) / dst[5]);
       dst[6] = src(3, 0) * dst[0];
       dst[7] = (src(3, 1) - dst[1] * dst[6]) * dst[2];
       dst[8] = (src(3, 2) - dst[3] * dst[6] - dst[4] * dst[7]) * dst[5];
       dst[9] = src(3, 3) - (dst[6] * dst[6] + dst[7] * dst[7] + dst[8] * dst[8]);
       if (dst[9] <= F(0.0))
          return false;
       else
          dst[9] = std::sqrt(F(1.0) / dst[9]);
       return true;
    }
 };
 /// struct to do a Cholesky decomposition (specialized, N = 3)
 template <class F, class M>
 struct _decomposer<F, 3, M> {
    /// method to do the decomposition
    bool operator()(F *dst, const M &src) const
    {
       if (src(0, 0) <= F(0.0))
          return false;
       dst[0] = std::sqrt(F(1.0) / src(0, 0));
       dst[1] = src(1, 0) * dst[0];
       dst[2] = src(1, 1) - dst[1] * dst[1];
       if (dst[2] <= F(0.0))
          return false;
       else
          dst[2] = std::sqrt(F(1.0) / dst[2]);
       dst[3] = src(2, 0) * dst[0];
       dst[4] = (src(2, 1) - dst[1] * dst[3]) * dst[2];
       dst[5] = src(2, 2) - (dst[3] * dst[3] + dst[4] * dst[4]);
       if (dst[5] <= F(0.0))
          return false;
       else
          dst[5] = std::sqrt(F(1.0) / dst[5]);
       return true;
    }
 };
 /// struct to do a Cholesky decomposition (specialized, N = 2)
 template <class F, class M>
 struct _decomposer<F, 2, M> {
    /// method to do the decomposition
    bool operator()(F *dst, const M &src) const
    {
       if (src(0, 0) <= F(0.0))
          return false;
       dst[0] = std::sqrt(F(1.0) / src(0, 0));
       dst[1] = src(1, 0) * dst[0];
       dst[2] = src(1, 1) - dst[1] * dst[1];
       if (dst[2] <= F(0.0))
          return false;
       else
          dst[2] = std::sqrt(F(1.0) / dst[2]);
       return true;
    }
 };
 /// struct to do a Cholesky decomposition (specialized, N = 1)
 template <class F, class M>
 struct _decomposer<F, 1, M> {
    /// method to do the decomposition
    bool operator()(F *dst, const M &src) const
    {
       if (src(0, 0) <= F(0.0))
          return false;
       dst[0] = std::sqrt(F(1.0) / src(0, 0));
       return true;
    }
 };
 /// struct to do a Cholesky decomposition (specialized, N = 0)
 template <class F, class M>
 struct _decomposer<F, 0, M> {
 private:
    _decomposer(){};
    bool operator()(F *dst, const M &src) const;
 };
 
 /// struct to obtain the inverse from a Cholesky decomposition (N = 6)
 template <class F, class M>
 struct _inverter<F, 6, M> {
    /// method to do the inversion
    inline void operator()(M &dst, const F *src) const
    {
       const F li21 = -src[1] * src[0] * src[2];
       const F li32 = -src[4] * src[2] * src[5];
       const F li31 = (src[1] * src[4] * src[2] - src[3]) * src[0] * src[5];
       const F li43 = -src[8] * src[9] * src[5];
       const F li42 = (src[4] * src[8] * src[5] - src[7]) * src[2] * src[9];
       const F li41 =
          (-src[1] * src[4] * src[8] * src[2] * src[5] + src[1] * src[7] * src[2] + src[3] * src[8] * src[5] - src[6]) *
          src[0] * src[9];
       const F li54 = -src[13] * src[14] * src[9];
       const F li53 = (src[13] * src[8] * src[9] - src[12]) * src[5] * src[14];
       const F li52 = (-src[4] * src[8] * src[13] * src[5] * src[9] + src[4] * src[12] * src[5] +
                       src[7] * src[13] * src[9] - src[11]) *
                      src[2] * src[14];
       const F li51 =
          (src[1] * src[4] * src[8] * src[13] * src[2] * src[5] * src[9] - 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] +
           src[11] * src[1] * src[2] + src[12] * src[3] * src[5] + src[13] * src[6] * src[9] - src[10]) *
          src[0] * src[14];
       const F li65 = -src[19] * src[20] * src[14];
       const F li64 = (src[19] * src[13] * src[14] - src[18]) * src[9] * src[20];
       const F li63 = (-src[8] * src[13] * src[19] * src[9] * src[14] + src[8] * src[18] * src[9] +
                       src[12] * src[19] * src[14] - src[17]) *
                      src[5] * src[20];
       const F li62 = (src[4] * src[8] * src[13] * src[19] * src[5] * src[9] * src[14] -
                       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] + src[17] * src[4] * src[5] +
                       src[18] * src[7] * src[9] + src[19] * src[11] * src[14] - src[16]) *
                      src[2] * src[20];
       const F li61 = (-src[19] * src[13] * src[8] * src[4] * src[1] * src[2] * src[5] * src[9] * src[14] +
                       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] +
                       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] -
                       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] - 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] +
                       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]) *
                      src[0] * src[20];
 
       dst(0, 0) = li61 * li61 + li51 * li51 + li41 * li41 + li31 * li31 + li21 * li21 + src[0] * src[0];
       dst(1, 0) = li61 * li62 + li51 * li52 + li41 * li42 + li31 * li32 + li21 * src[2];
       dst(1, 1) = li62 * li62 + li52 * li52 + li42 * li42 + li32 * li32 + src[2] * src[2];
       dst(2, 0) = li61 * li63 + li51 * li53 + li41 * li43 + li31 * src[5];
       dst(2, 1) = li62 * li63 + li52 * li53 + li42 * li43 + li32 * src[5];
       dst(2, 2) = li63 * li63 + li53 * li53 + li43 * li43 + src[5] * src[5];
       dst(3, 0) = li61 * li64 + li51 * li54 + li41 * src[9];
       dst(3, 1) = li62 * li64 + li52 * li54 + li42 * src[9];
       dst(3, 2) = li63 * li64 + li53 * li54 + li43 * src[9];
       dst(3, 3) = li64 * li64 + li54 * li54 + src[9] * src[9];
       dst(4, 0) = li61 * li65 + li51 * src[14];
       dst(4, 1) = li62 * li65 + li52 * src[14];
       dst(4, 2) = li63 * li65 + li53 * src[14];
       dst(4, 3) = li64 * li65 + li54 * src[14];
       dst(4, 4) = li65 * li65 + src[14] * src[14];
       dst(5, 0) = li61 * src[20];
       dst(5, 1) = li62 * src[20];
       dst(5, 2) = li63 * src[20];
       dst(5, 3) = li64 * src[20];
       dst(5, 4) = li65 * src[20];
       dst(5, 5) = src[20] * src[20];
    }
 };
 /// struct to obtain the inverse from a Cholesky decomposition (N = 5)
 template <class F, class M>
 struct _inverter<F, 5, M> {
    /// method to do the inversion
    inline void operator()(M &dst, const F *src) const
    {
       const F li21 = -src[1] * src[0] * src[2];
       const F li32 = -src[4] * src[2] * src[5];
       const F li31 = (src[1] * src[4] * src[2] - src[3]) * src[0] * src[5];
       const F li43 = -src[8] * src[9] * src[5];
       const F li42 = (src[4] * src[8] * src[5] - src[7]) * src[2] * src[9];
       const F li41 =
          (-src[1] * src[4] * src[8] * src[2] * src[5] + src[1] * src[7] * src[2] + src[3] * src[8] * src[5] - src[6]) *
          src[0] * src[9];
       const F li54 = -src[13] * src[14] * src[9];
       const F li53 = (src[13] * src[8] * src[9] - src[12]) * src[5] * src[14];
       const F li52 = (-src[4] * src[8] * src[13] * src[5] * src[9] + src[4] * src[12] * src[5] +
                       src[7] * src[13] * src[9] - src[11]) *
                      src[2] * src[14];
       const F li51 =
          (src[1] * src[4] * src[8] * src[13] * src[2] * src[5] * src[9] - 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] +
           src[11] * src[1] * src[2] + src[12] * src[3] * src[5] + src[13] * src[6] * src[9] - src[10]) *
          src[0] * src[14];
 
       dst(0, 0) = li51 * li51 + li41 * li41 + li31 * li31 + li21 * li21 + src[0] * src[0];
       dst(1, 0) = li51 * li52 + li41 * li42 + li31 * li32 + li21 * src[2];
       dst(1, 1) = li52 * li52 + li42 * li42 + li32 * li32 + src[2] * src[2];
       dst(2, 0) = li51 * li53 + li41 * li43 + li31 * src[5];
       dst(2, 1) = li52 * li53 + li42 * li43 + li32 * src[5];
       dst(2, 2) = li53 * li53 + li43 * li43 + src[5] * src[5];
       dst(3, 0) = li51 * li54 + li41 * src[9];
       dst(3, 1) = li52 * li54 + li42 * src[9];
       dst(3, 2) = li53 * li54 + li43 * src[9];
       dst(3, 3) = li54 * li54 + src[9] * src[9];
       dst(4, 0) = li51 * src[14];
       dst(4, 1) = li52 * src[14];
       dst(4, 2) = li53 * src[14];
       dst(4, 3) = li54 * src[14];
       dst(4, 4) = src[14] * src[14];
    }
 };
 /// struct to obtain the inverse from a Cholesky decomposition (N = 4)
 template <class F, class M>
 struct _inverter<F, 4, M> {
    /// method to do the inversion
    inline void operator()(M &dst, const F *src) const
    {
       const F li21 = -src[1] * src[0] * src[2];
       const F li32 = -src[4] * src[2] * src[5];
       const F li31 = (src[1] * src[4] * src[2] - src[3]) * src[0] * src[5];
       const F li43 = -src[8] * src[9] * src[5];
       const F li42 = (src[4] * src[8] * src[5] - src[7]) * src[2] * src[9];
       const F li41 =
          (-src[1] * src[4] * src[8] * src[2] * src[5] + src[1] * src[7] * src[2] + src[3] * src[8] * src[5] - src[6]) *
          src[0] * src[9];
 
       dst(0, 0) = li41 * li41 + li31 * li31 + li21 * li21 + src[0] * src[0];
       dst(1, 0) = li41 * li42 + li31 * li32 + li21 * src[2];
       dst(1, 1) = li42 * li42 + li32 * li32 + src[2] * src[2];
       dst(2, 0) = li41 * li43 + li31 * src[5];
       dst(2, 1) = li42 * li43 + li32 * src[5];
       dst(2, 2) = li43 * li43 + src[5] * src[5];
       dst(3, 0) = li41 * src[9];
       dst(3, 1) = li42 * src[9];
       dst(3, 2) = li43 * src[9];
       dst(3, 3) = src[9] * src[9];
    }
 };
 /// struct to obtain the inverse from a Cholesky decomposition (N = 3)
 template <class F, class M>
 struct _inverter<F, 3, M> {
    /// method to do the inversion
    inline void operator()(M &dst, const F *src) const
    {
       const F li21 = -src[1] * src[0] * src[2];
       const F li32 = -src[4] * src[2] * src[5];
       const F li31 = (src[1] * src[4] * src[2] - src[3]) * src[0] * src[5];
 
       dst(0, 0) = li31 * li31 + li21 * li21 + src[0] * src[0];
       dst(1, 0) = li31 * li32 + li21 * src[2];
       dst(1, 1) = li32 * li32 + src[2] * src[2];
       dst(2, 0) = li31 * src[5];
       dst(2, 1) = li32 * src[5];
       dst(2, 2) = src[5] * src[5];
    }
 };
 /// struct to obtain the inverse from a Cholesky decomposition (N = 2)
 template <class F, class M>
 struct _inverter<F, 2, M> {
    /// method to do the inversion
    inline void operator()(M &dst, const F *src) const
    {
       const F li21 = -src[1] * src[0] * src[2];
 
       dst(0, 0) = li21 * li21 + src[0] * src[0];
       dst(1, 0) = li21 * src[2];
       dst(1, 1) = src[2] * src[2];
    }
 };
 /// struct to obtain the inverse from a Cholesky decomposition (N = 1)
 template <class F, class M>
 struct _inverter<F, 1, M> {
    /// method to do the inversion
    inline void operator()(M &dst, const F *src) const { dst(0, 0) = src[0] * src[0]; }
 };
 /// struct to obtain the inverse from a Cholesky decomposition (N = 0)
 template <class F, class M>
 struct _inverter<F, 0, M> {
 private:
    _inverter();
    void operator()(M &dst, const F *src) const;
 };
 
 /// struct to solve a linear system using its Cholesky decomposition (N=6)
 template <class F, class V>
 struct _solver<F, 6, V> {
    /// method to solve the linear system
    void operator()(V &rhs, const F *l) const
    {
       // solve Ly = rhs
       const F y0 = rhs[0] * l[0];
       const F y1 = (rhs[1] - l[1] * y0) * l[2];
       const F y2 = (rhs[2] - (l[3] * y0 + l[4] * y1)) * l[5];
       const F y3 = (rhs[3] - (l[6] * y0 + l[7] * y1 + l[8] * y2)) * l[9];
       const F y4 = (rhs[4] - (l[10] * y0 + l[11] * y1 + l[12] * y2 + l[13] * y3)) * l[14];
       const F y5 = (rhs[5] - (l[15] * y0 + l[16] * y1 + l[17] * y2 + l[18] * y3 + l[19] * y4)) * l[20];
       // solve L^Tx = y, and put x into rhs
       rhs[5] = y5 * l[20];
       rhs[4] = (y4 - l[19] * rhs[5]) * l[14];
       rhs[3] = (y3 - (l[18] * rhs[5] + l[13] * rhs[4])) * l[9];
       rhs[2] = (y2 - (l[17] * rhs[5] + l[12] * rhs[4] + l[8] * rhs[3])) * l[5];
       rhs[1] = (y1 - (l[16] * rhs[5] + l[11] * rhs[4] + l[7] * rhs[3] + l[4] * rhs[2])) * l[2];
       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];
    }
 };
 /// struct to solve a linear system using its Cholesky decomposition (N=5)
 template <class F, class V>
 struct _solver<F, 5, V> {
    /// method to solve the linear system
    void operator()(V &rhs, const F *l) const
    {
       // solve Ly = rhs
       const F y0 = rhs[0] * l[0];
       const F y1 = (rhs[1] - l[1] * y0) * l[2];
       const F y2 = (rhs[2] - (l[3] * y0 + l[4] * y1)) * l[5];
       const F y3 = (rhs[3] - (l[6] * y0 + l[7] * y1 + l[8] * y2)) * l[9];
       const F y4 = (rhs[4] - (l[10] * y0 + l[11] * y1 + l[12] * y2 + l[13] * y3)) * l[14];
       // solve L^Tx = y, and put x into rhs
       rhs[4] = (y4)*l[14];
       rhs[3] = (y3 - (l[13] * rhs[4])) * l[9];
       rhs[2] = (y2 - (l[12] * rhs[4] + l[8] * rhs[3])) * l[5];
       rhs[1] = (y1 - (l[11] * rhs[4] + l[7] * rhs[3] + l[4] * rhs[2])) * l[2];
       rhs[0] = (y0 - (l[10] * rhs[4] + l[6] * rhs[3] + l[3] * rhs[2] + l[1] * rhs[1])) * l[0];
    }
 };
 /// struct to solve a linear system using its Cholesky decomposition (N=4)
 template <class F, class V>
 struct _solver<F, 4, V> {
    /// method to solve the linear system
    void operator()(V &rhs, const F *l) const
    {
       // solve Ly = rhs
       const F y0 = rhs[0] * l[0];
       const F y1 = (rhs[1] - l[1] * y0) * l[2];
       const F y2 = (rhs[2] - (l[3] * y0 + l[4] * y1)) * l[5];
       const F y3 = (rhs[3] - (l[6] * y0 + l[7] * y1 + l[8] * y2)) * l[9];
       // solve L^Tx = y, and put x into rhs
       rhs[3] = (y3)*l[9];
       rhs[2] = (y2 - (l[8] * rhs[3])) * l[5];
       rhs[1] = (y1 - (l[7] * rhs[3] + l[4] * rhs[2])) * l[2];
       rhs[0] = (y0 - (l[6] * rhs[3] + l[3] * rhs[2] + l[1] * rhs[1])) * l[0];
    }
 };
 /// struct to solve a linear system using its Cholesky decomposition (N=3)
 template <class F, class V>
 struct _solver<F, 3, V> {
    /// method to solve the linear system
    void operator()(V &rhs, const F *l) const
    {
       // solve Ly = rhs
       const F y0 = rhs[0] * l[0];
       const F y1 = (rhs[1] - l[1] * y0) * l[2];
       const F y2 = (rhs[2] - (l[3] * y0 + l[4] * y1)) * l[5];
       // solve L^Tx = y, and put x into rhs
       rhs[2] = (y2)*l[5];
       rhs[1] = (y1 - (l[4] * rhs[2])) * l[2];
       rhs[0] = (y0 - (l[3] * rhs[2] + l[1] * rhs[1])) * l[0];
    }
 };
 /// struct to solve a linear system using its Cholesky decomposition (N=2)
 template <class F, class V>
 struct _solver<F, 2, V> {
    /// method to solve the linear system
    void operator()(V &rhs, const F *l) const
    {
       // solve Ly = rhs
       const F y0 = rhs[0] * l[0];
       const F y1 = (rhs[1] - l[1] * y0) * l[2];
       // solve L^Tx = y, and put x into rhs
       rhs[1] = (y1)*l[2];
       rhs[0] = (y0 - (l[1] * rhs[1])) * l[0];
    }
 };
 /// struct to solve a linear system using its Cholesky decomposition (N=1)
 template <class F, class V>
 struct _solver<F, 1, V> {
    /// method to solve the linear system
    void operator()(V &rhs, const F *l) const
    {
       // solve LL^T x = rhs, put y into rhs
       rhs[0] *= l[0] * l[0];
    }
 };
 /// struct to solve a linear system using its Cholesky decomposition (N=0)
 template <class F, class V>
 struct _solver<F, 0, V> {
 private:
    _solver();
    void operator()(V &rhs, const F *l) const;
 };
 } // namespace CholeskyDecompHelpers
 
 } // namespace Math
 
 } // namespace ROOT
 
 #endif // ROOT_Math_CHOLESKYDECOMP

This page was automatically generated by the 2.3.7 LXR engine. The LXR team