From 136507be3d831d0747c8ccbcebc63ebb3db3c484 Mon Sep 17 00:00:00 2001
From: Max Kahnt <max.kahnt@fu-berlin.de>
Date: Tue, 21 Jul 2015 16:00:44 +0200
Subject: [PATCH] Add LDLt decomposition to dune-solvers (synchronize with
 dune-fufem) staticmatrixtools.

---
 dune/solvers/common/staticmatrixtools.hh | 75 ++++++++++++++++++++++++
 1 file changed, 75 insertions(+)

diff --git a/dune/solvers/common/staticmatrixtools.hh b/dune/solvers/common/staticmatrixtools.hh
index 088bd901..576ea7cf 100644
--- a/dune/solvers/common/staticmatrixtools.hh
+++ b/dune/solvers/common/staticmatrixtools.hh
@@ -382,6 +382,81 @@ class StaticMatrix
             return *(reinterpret_cast< Dune::FieldMatrix<K,1,1>* > (&x));
         }
 
+        /** \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];
+          }
+        }
 
 };
 
-- 
GitLab