From af6f61474f379fd75f2e164b12dc3c852d8d4042 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Thu, 12 Dec 2013 19:10:00 +0000
Subject: [PATCH] A direct dense solver constructing an energy-orthonormal
basis.
Nice, because it works for rank-deficient matrices.
And still fast enough for the systems I am considering here.
[[Imported from SVN: r9580]]
---
dune/gfe/gramschmidtsolver.hh | 94 +++++++++++++++++++++++++++++++++++
1 file changed, 94 insertions(+)
create mode 100644 dune/gfe/gramschmidtsolver.hh
diff --git a/dune/gfe/gramschmidtsolver.hh b/dune/gfe/gramschmidtsolver.hh
new file mode 100644
index 00000000..7e9be4d6
--- /dev/null
+++ b/dune/gfe/gramschmidtsolver.hh
@@ -0,0 +1,94 @@
+#ifndef DUNE_GFE_GRAMSCHMIDTSOLVER_HH
+#define DUNE_GFE_GRAMSCHMIDTSOLVER_HH
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+
+/** \brief Direct solver for a dense symmetric linear system, using an orthonormal basis
+ *
+ * This solver computes an A-orthonormal basis, and uses that to compute the solution
+ * of a linear system. The advantage of this is that it works even if the matrix is
+ * known to have a non-trivial kernel. This happens in GFE applications when the Hessian
+ * of a functional on a manifold is given in coordinates of the surrounding space.
+ * The method if efficient for the small systems that we consider in GFE applications.
+ *
+ * \tparam field_type Type used to store scalars
+ * \tparam embeddedDim Number of rows and columns of the linear system
+ * \tparam dim The rank of the matrix
+ */
+template <class field_type, int dim, int embeddedDim>
+class GramSchmidtSolver
+{
+ /** \brief Normalize a vector to unit length, measured in a matrix norm
+ * \param matrix The matrix inducing the matrix norm
+ * \param[in,out] v The vector to normalize
+ */
+ static void normalize(const Dune::FieldMatrix<field_type,embeddedDim,embeddedDim>& matrix,
+ Dune::FieldVector<field_type,embeddedDim>& v)
+ {
+ field_type energyNormSquared = 0;
+
+ for (int i=0; i<v.size(); i++)
+ for (int j=0; j<v.size(); j++)
+ energyNormSquared += matrix[i][j]*v[i]*v[j];
+
+ v /= std::sqrt(energyNormSquared);
+ }
+
+
+ /** \brief Project vj on vi, and subtract the result from vj
+ *
+ * \param matrix The matrix the defines the scalar product
+ */
+ static void project(const Dune::FieldMatrix<field_type,embeddedDim,embeddedDim>& matrix,
+ const Dune::FieldVector<field_type,embeddedDim>& vi,
+ Dune::FieldVector<field_type,embeddedDim>& vj)
+ {
+
+ field_type energyScalarProduct = 0;
+
+ for (int i=0; i<vi.size(); i++)
+ for (int j=0; j<vj.size(); j++)
+ energyScalarProduct += matrix[i][j]*vi[i]*vj[j];
+
+ for (int i=0; i<vj.size(); i++)
+ vj[i] -= energyScalarProduct * vi[i];
+
+ }
+
+public:
+ /** Solve linear system by constructing an energy-orthonormal basis */
+ static void solve(const Dune::FieldMatrix<field_type,embeddedDim,embeddedDim>& matrix,
+ Dune::FieldVector<field_type,embeddedDim>& x,
+ const Dune::FieldVector<field_type,embeddedDim>& rhs,
+ const Dune::FieldMatrix<field_type,dim,embeddedDim>& basis)
+ {
+ // 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);
+
+ for (int i=0; i<dim; i++) {
+
+ normalize(matrix, orthonormalBasis[i]);
+
+ for (int j=0; j<i; j++)
+ project(matrix, orthonormalBasis[i], orthonormalBasis[j]);
+
+ }
+
+ // Solve the system in the orthonormal basis
+ Dune::FieldVector<field_type,dim> orthoCoefficient;
+ for (int i=0; i<dim; i++)
+ orthoCoefficient[i] = rhs*orthonormalBasis[i];
+
+ // Solution in canonical basis
+ x = 0;
+ for (int i=0; i<dim; i++)
+ x.axpy(orthoCoefficient[i], orthonormalBasis[i]);
+
+ }
+
+};
+
+#endif
\ No newline at end of file
--
GitLab