ldlt.hh

#ifndef DUNE_MATRIX_VECTOR_LDLT_HH
#define DUNE_MATRIX_VECTOR_LDLT_HH

namespace Dune {
  namespace MatrixVector {

    /** \brief Compute an LDL^T decomposition
             *
             * The methods computes a decomposition A=LDL^T of a given dense
             * symmetric matrix A such that L is lower triangular with all
             * diagonal elements equal to 1 and D is diagonal. If A is positive
             * definite then A=(LD^0.5)(LD^0.5)^T is the Cholesky decomposition.
             * However, the LDL^T decomposition does also work for indefinite
             * symmetric matrices and is more stable than the Cholesky
     * decomposition
             * since no square roots are required.
             *
             * The method does only change the nontrivial entries of the given
     * matrix
             * L and D, i.e., it does not set the trivial 0 and 1 entries.
             * Thus one can store both in a single matrix M and use
             * M as argument for L and D.
             *
             * The method can furthermore work in-place, i.e., it is safe to
             * use A as argument for L and D. In this case the entries of A
             * below and on the diagonal are changed to those to those of
             * L and D, respectively.
             *
             * \param A Matrix to be decomposed. Only the lower triangle is
     * used.
             * \param L Matrix to store the lower triangle. Only entries below
     * the diagonal are written.
             * \param D Matrix to store the diagonal. Only diagonal entries are
     * written.
             */
    template <class SymmetricMatrix, class LowerTriangularMatrix,
              class DiagonalMatrix>
    static void ldlt(const SymmetricMatrix& A, LowerTriangularMatrix& L,
                     DiagonalMatrix& D) {
      for (unsigned int i = 0; i < A.N(); ++i) {
        D[i][i] = A[i][i];
        for (unsigned int j = 0; j < i; ++j) {
          L[i][j] = A[i][j];
          for (unsigned int k = 0; k < j; ++k)
            L[i][j] -= L[i][k] * L[j][k] * D[k][k];
          L[i][j] /= D[j][j];
        }
        for (unsigned int k = 0; k < i; ++k)
          D[i][i] -= L[i][k] * L[i][k] * D[k][k];
      }
    }

    /** \brief Solve linear system using a LDL^T decomposition.
     *
     * The methods solves a linear system Mx=b where A is given
     * by a decomposition A=LDL^T. The method does only use
     * the values of L and D below and on the diagonal, respectively.
     * The 1 entries on the diagonal of L are not required.
     * If L and D are stored in a single matrix it is safe
     * the use this matrix as argument for both.
     *
     * Note that the solution vector must already have the correct size.
     *
     * \param L Matrix containing the lower triangle of the decomposition
     * \param D Matrix containing the diagonal of the decomposition
     * \param b Right hand side on the linear system
     * \param x Vector to store the solution of the linear system
     */
    template <class LowerTriangularMatrix, class DiagonalMatrix,
              class RhsVector, class SolVector>
    static void ldltSolve(const LowerTriangularMatrix& L,
                          const DiagonalMatrix& D, const RhsVector& b,
                          SolVector& x) {
      for (unsigned int i = 0; i < x.size(); ++i) {
        x[i] = b[i];
        for (unsigned int j = 0; j < i; ++j)
          x[i] -= L[i][j] * x[j];
      }
      for (unsigned int i = 0; i < x.size(); ++i)
        x[i] /= D[i][i];
      for (int i = x.size() - 1; i >= 0; --i) {
        for (unsigned int j = i + 1; j < x.size(); ++j)
          x[i] -= L[j][i] * x[j];
      }
    }
  }
}
#endif