From 917fe5b9b8ab770afec83e2fcf6dc78e53635253 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Thu, 10 Oct 2019 21:57:14 +0200
Subject: [PATCH] Fix GramSchmidtSolver
The computation of the A-orthogonal matrices was wrong,
and for some reason nobody ever really noticed this.
One reason could be that in GFE applications A is frequently
very close to the identity, and then you wouldn't notice
the bug that much.
Anyway, localgeodesicfefunctiontest.cc passes now.
---
dune/gfe/gramschmidtsolver.hh | 15 ++++++++-------
1 file changed, 8 insertions(+), 7 deletions(-)
diff --git a/dune/gfe/gramschmidtsolver.hh b/dune/gfe/gramschmidtsolver.hh
index 5b70d665..0577e358 100644
--- a/dune/gfe/gramschmidtsolver.hh
+++ b/dune/gfe/gramschmidtsolver.hh
@@ -63,13 +63,14 @@ public:
{
// Use the Gram-Schmidt algorithm to compute a basis that is orthonormal
// with respect to the given matrix.
- for (int i=0; i<dim; i++) {
+ normalize(matrix, orthonormalBasis_[0]);
- normalize(matrix, orthonormalBasis_[i]);
+ for (int i=1; i<dim; i++) {
for (int j=0; j<i; j++)
- project(matrix, orthonormalBasis_[i], orthonormalBasis_[j]);
+ project(matrix, orthonormalBasis_[j], orthonormalBasis_[i]);
+ normalize(matrix, orthonormalBasis_[i]);
}
}
@@ -100,14 +101,14 @@ public:
// Use the Gram-Schmidt algorithm to compute a basis that is orthonormal
// with respect to the given matrix.
Dune::FieldMatrix<field_type,dim,embeddedDim> orthonormalBasis(basis);
+ normalize(matrix, orthonormalBasis[0]);
- for (int i=0; i<dim; i++) {
-
- normalize(matrix, orthonormalBasis[i]);
+ for (int i=1; i<dim; i++) {
for (int j=0; j<i; j++)
- project(matrix, orthonormalBasis[i], orthonormalBasis[j]);
+ project(matrix, orthonormalBasis[j], orthonormalBasis[i]);
+ normalize(matrix, orthonormalBasis[i]);
}
// Solve the system in the orthonormal basis
--
GitLab