Incorporate ldlt from StaticMatrix

dune/matrix-vector/CMakeLists.txt

@@ -4,6 +4,7 @@ install(FILES
   axpy.hh
   componentwisematrixmap.hh
   genericvectortools.hh
+  ldlt.hh
   matrixtraits.hh
   scalartraits.hh
   singlenonzerocolumnmatrix.hh

dune/matrix-vector/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