In This Package:

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

inverts a symmetric matrix (dimension n, general, implementation) More...

#include <SymMatrixInverter.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::SymMatrixInverter::inverterForAllN< T >

inverts a symmetric matrix (dimension n, general, implementation)

Definition at line 21 of file SymMatrixInverter.h.

---

Member Typedef Documentation

template<class T>

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

Definition at line 23 of file SymMatrixInverter.h.

---

Member Function Documentation

template<class T>

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

Definition at line 24 of file SymMatrixInverter.h.

        {
          // perform variant cholesky decomposition of matrix: M = L O L^T
          // L is lower triangular matrix with positive eigenvalues, thus same
          // use as in proper cholesky decomposition
          // O is diagonal, elements are either +1 or -1, to account for potential
          // negative eigenvalues of M (name is o because it contains "ones")
          // only thing that can go wrong: trying to divide by zero on diagonale
          // (matrix is singular in this case)

          // 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];
          // sign on diagonale
          F o[N];

          // 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] * o[k];
              base1[j] = tmp *= base2[j] * o[j];
              // keep track of contribution to element on diagonale
              tmpdiag -= tmp * tmp * o[j];
            }
            // keep truncation error small
            tmpdiag = m(i, i) + tmpdiag;
            // save sign to have L positive definite
            if (tmpdiag < F(0)) {
              tmpdiag = -tmpdiag;
              o[i] = F(-1);
            } else {
              // tmpdiag has right sign
              o[i] = F(1);
            }
            if (tmpdiag == F(0)) return false;
            else base1[i] = std::sqrt(F(1) / tmpdiag);
          }

          // ok, next step: invert off-diagonal part of matrix L
          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 O 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] * o[k];
              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/SymMatrixInverter.h

---