In This Package:

Gaudi::Math::SymPosDefMatrixInverter::inverterForAllN< T > Struct Template Reference

invert symmetric positive definite matrix (dimension n, general, implementation) More...

#include <SymPosDefMatrixInverter.h>

List of all members.

Public Types
typedef T::value_type	F
Public Member Functions
bool	operator() (T &m, int N) const

---

Detailed Description

template<class T>
struct Gaudi::Math::SymPosDefMatrixInverter::inverterForAllN< T >

invert symmetric positive definite matrix (dimension n, general, implementation)

Definition at line 21 of file SymPosDefMatrixInverter.h.

---

Member Typedef Documentation

template<class T>

typedef T::value_type Gaudi::Math::SymPosDefMatrixInverter::inverterForAllN< T >::F

Definition at line 23 of file SymPosDefMatrixInverter.h.

---

Member Function Documentation

template<class T>

bool Gaudi::Math::SymPosDefMatrixInverter::inverterForAllN< T >::operator()	(	T &	m,
		int	N
	)			const [inline]

Definition at line 24 of file SymPosDefMatrixInverter.h.

        {
          // 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)

          // we only need N * (N + 1) / 2 elements to store L
          // element L(i,j) is at array position (i * (i + 1)) / 2 + j
          F l[(N * (N + 1)) / 2];

          // quirk: we need to invert L later anyway, so we can invert elements
          // on diagonale straight away (we only ever need their reciprocals!)

          // cache starting address of rows of L for speed reasons
          F *base1 = &l[0];
          for (int i = 0; i < N; base1 += ++i) {
            F tmpdiag = F(0);   // for element on diagonale
            // calculate off-diagonal elements
            F *base2 = &l[0];
            for (int j = 0; j < i; base2 += ++j) {
              F tmp = m(i, j);
              for (int k = j; k--; )
                tmp -= base1[k] * base2[k];
              base1[j] = tmp *= base2[j];
              // keep track of contribution to element on diagonale
              tmpdiag += tmp * tmp;
            }
            // keep truncation error small
            tmpdiag = m(i, i) - tmpdiag;
            // check if positive definite
            if (tmpdiag <= F(0)) return false;
            else base1[i] = std::sqrt(F(1) / tmpdiag);
          }

          // ok, next step: invert off-diagonal part of matrix
          base1 = &l[1];
          for (int i = 1; i < N; base1 += ++i) {
            for (int j = 0; j < i; ++j) {
              F tmp = F(0);
              F *base2 = &l[(i * (i - 1)) / 2];
              for (int 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 (int i = N; i--; ) {
            for (int j = i + 1; j--; ) {
              F tmp = F(0);
              base1 = &l[(N * (N - 1)) / 2];
              for (int k = N; k-- > i; base1 -= k)
                tmp += base1[i] * base1[j];
              m(i,j) = tmp;
            }
          }

          return true;
        }

---

The documentation for this struct was generated from the following file:

• /home/dayabay/installs/trunk_2011_04_11_opt/NuWa-trunk/lhcb/Kernel/LHCbMath/LHCbMath/SymPosDefMatrixInverter.h

---