From 1888608dd11d672026a6bb11a520335e4c7675bf Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Mon, 25 Jul 2016 13:32:22 +0200
Subject: [PATCH] ldlt moved to dune-matrix-vector
---
dune/solvers/common/staticmatrixtools.hh | 33 ++-------------------
dune/solvers/iterationsteps/blockgssteps.hh | 7 +++--
2 files changed, 7 insertions(+), 33 deletions(-)
diff --git a/dune/solvers/common/staticmatrixtools.hh b/dune/solvers/common/staticmatrixtools.hh
index 733407ce..1d159a94 100644
--- a/dune/solvers/common/staticmatrixtools.hh
+++ b/dune/solvers/common/staticmatrixtools.hh
@@ -11,6 +11,7 @@
#include <dune/matrix-vector/addtodiagonal.hh>
#include <dune/matrix-vector/axpy.hh>
+#include <dune/matrix-vector/ldlt.hh>
#include <dune/matrix-vector/scalartraits.hh>
#include "arithmetic.hh"
@@ -458,19 +459,7 @@ namespace StaticMatrix
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];
- }
+ Dune::MatrixVector::ldlt(A,L,D);
}
/** \brief Solve linear system using a LDL^T decomposition.
@@ -492,24 +481,8 @@ namespace StaticMatrix
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];
- }
+ Dune::MatrixVector::ldltSolve(L,D,b,x);
}
-
-
-
-
}
#endif
diff --git a/dune/solvers/iterationsteps/blockgssteps.hh b/dune/solvers/iterationsteps/blockgssteps.hh
index 6835c835..01d17c64 100644
--- a/dune/solvers/iterationsteps/blockgssteps.hh
+++ b/dune/solvers/iterationsteps/blockgssteps.hh
@@ -6,7 +6,8 @@
#include <dune/common/parametertree.hh>
-#include <dune/solvers/common/staticmatrixtools.hh>
+#include <dune/matrix-vector/ldlt.hh>
+
#include <dune/solvers/common/defaultbitvector.hh>
#include <dune/solvers/iterationsteps/lineariterationstep.hh>
@@ -92,8 +93,8 @@ VBlock direct(const MBlock& m, const VBlock& b) {
template <class MBlock, class VBlock>
VBlock ldlt(MBlock m, const VBlock& b) {
VBlock x;
- StaticMatrix::ldlt(m, m, m);
- StaticMatrix::ldltSolve(m, m, b, x);
+ Dune::MatrixVector::ldlt(m, m, m);
+ Dune::MatrixVector::ldltSolve(m, m, b, x);
return x;
}
--
GitLab