From 408c0310f0b645a595d3245b8d773c9b536d7ace Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Thu, 12 Dec 2013 21:09:11 +0000
Subject: [PATCH] Start using a dedicated implementation of a symmetric matrix
[[Imported from SVN: r9582]]
---
dune/gfe/gramschmidtsolver.hh | 22 +++--------
dune/gfe/symmetricmatrix.hh | 68 ++++++++++++++++++++++++++++++++
dune/gfe/targetspacertrsolver.cc | 8 +++-
3 files changed, 81 insertions(+), 17 deletions(-)
create mode 100644 dune/gfe/symmetricmatrix.hh
diff --git a/dune/gfe/gramschmidtsolver.hh b/dune/gfe/gramschmidtsolver.hh
index 7e9be4d6..3623ff6b 100644
--- a/dune/gfe/gramschmidtsolver.hh
+++ b/dune/gfe/gramschmidtsolver.hh
@@ -2,8 +2,8 @@
#define DUNE_GFE_GRAMSCHMIDTSOLVER_HH
#include <dune/common/fvector.hh>
-#include <dune/common/fmatrix.hh>
+#include <dune/gfe/symmetricmatrix.hh>
/** \brief Direct solver for a dense symmetric linear system, using an orthonormal basis
*
@@ -24,16 +24,10 @@ class GramSchmidtSolver
* \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,
+ static void normalize(const Dune::SymmetricMatrix<field_type,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);
+ v /= std::sqrt(matrix.energyScalarProduct(v,v));
}
@@ -41,16 +35,12 @@ class GramSchmidtSolver
*
* \param matrix The matrix the defines the scalar product
*/
- static void project(const Dune::FieldMatrix<field_type,embeddedDim,embeddedDim>& matrix,
+ static void project(const Dune::SymmetricMatrix<field_type,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];
+ field_type energyScalarProduct = matrix.energyScalarProduct(vi,vj);
for (int i=0; i<vj.size(); i++)
vj[i] -= energyScalarProduct * vi[i];
@@ -59,7 +49,7 @@ class GramSchmidtSolver
public:
/** Solve linear system by constructing an energy-orthonormal basis */
- static void solve(const Dune::FieldMatrix<field_type,embeddedDim,embeddedDim>& matrix,
+ static void solve(const Dune::SymmetricMatrix<field_type,embeddedDim>& matrix,
Dune::FieldVector<field_type,embeddedDim>& x,
const Dune::FieldVector<field_type,embeddedDim>& rhs,
const Dune::FieldMatrix<field_type,dim,embeddedDim>& basis)
diff --git a/dune/gfe/symmetricmatrix.hh b/dune/gfe/symmetricmatrix.hh
new file mode 100644
index 00000000..c4611231
--- /dev/null
+++ b/dune/gfe/symmetricmatrix.hh
@@ -0,0 +1,68 @@
+#ifndef DUNE_GFE_SYMMETRICMATRIX_HH
+#define DUNE_GFE_SYMMETRICMATRIX_HH
+
+#include <dune/common/fvector.hh>
+
+namespace Dune {
+
+/** \brief A class implementing a symmetric matrix with compile-time size
+ *
+ * A \f$ dim\times dim \f$ matrix is stored internally as a <tt>Dune::FieldVector<double, dim*(dim+1)/2></tt>
+ * The components are assumed to be \f$ [ E(1,1),\ E(2,1),\ E(2,2),\ E(3,1),\ E(3,2),\ E(3,3) ]\f$
+ * and analogous for other dimensions
+ * \tparam dim number of lines/columns of the matrix
+ */
+template <class T, int N>
+class SymmetricMatrix
+{
+
+public:
+ /** \brief Default constructor
+ *
+ * Tensor is initialized containing zeros if no argument is given.
+ * \param eye if true tensor is initialized as identity
+ */
+ SymmetricMatrix()
+ {}
+
+ /** \brief Matrix style random read/write access to components
+ * \param i line index
+ * \param j column index
+ * \note You need to know what you are doing: You can only access the lower
+ * left triangular submatrix using this methods. It requires i<=j.
+ */
+ T& operator() (int i, int j)
+ {
+ assert(i>=j);
+ return data_[i*(i+1)/2+j];
+ }
+
+ /** \brief Matrix style random read access to components
+ * \param i line index
+ * \param j column index
+ * \note You need to know what you are doing: You can only access the lower
+ * left triangular submatrix using this methods. It requires i<=j.
+ */
+ const T& operator() (int i, int j) const
+ {
+ assert(i>=j);
+ return data_[i*(i+1)/2+j];
+ }
+
+ T energyScalarProduct (const FieldVector<T,N>& v1, const FieldVector<T,N>& v2) const
+ {
+ T result = 0;
+ for (size_t i=0; i<N; i++)
+ for (size_t j=0; j<=i; j++)
+ result += (1+(i!=j)) * operator()(i,j) * v1[i] * v2[j];
+
+ return result;
+ }
+
+
+private:
+ Dune::FieldVector<T,N*(N+1)/2> data_;
+};
+
+}
+#endif
diff --git a/dune/gfe/targetspacertrsolver.cc b/dune/gfe/targetspacertrsolver.cc
index 25909b44..ca3d9e05 100644
--- a/dune/gfe/targetspacertrsolver.cc
+++ b/dune/gfe/targetspacertrsolver.cc
@@ -107,8 +107,14 @@ void TargetSpaceRiemannianTRSolver<TargetSpace>::solve()
#ifdef USE_GAUSS_SEIDEL_SOLVER
innerSolver_->solve();
#else
+ Dune::SymmetricMatrix<field_type, embeddedBlocksize> symmetricHessian;
+ for (size_t j=0; j<embeddedBlocksize; j++)
+ for (size_t k=0; k<=j; k++)
+ symmetricHessian(j,k) = hesseMatrix[0][0][j][k];
+
+
Dune::FieldMatrix<field_type,blocksize,embeddedBlocksize> basis = x_.orthonormalFrame();
- GramSchmidtSolver<field_type, blocksize, embeddedBlocksize>::solve(hesseMatrix[0][0], corr[0], rhs[0], basis);
+ GramSchmidtSolver<field_type, blocksize, embeddedBlocksize>::solve(symmetricHessian, corr[0], rhs[0], basis);
#endif
#ifdef USE_GAUSS_SEIDEL_SOLVER
--
GitLab