From efd45e108bb8541e7af584459f890d1217deeec3 Mon Sep 17 00:00:00 2001
From: Patrick Jaap <patrick.jaap@tu-dresden.de>
Date: Fri, 29 Mar 2019 11:41:17 +0100
Subject: [PATCH] Make symmetrixMatrix isometric
- We have to add the factor 1./sqrt(2) in read access for off-diagonal entries
- The random write access is removed
- As replacement a setEntry method is provided that ensures isometry
---
dune/fufem/symmetricmatrix.hh | 41 +++++++++++++++++------------------
1 file changed, 20 insertions(+), 21 deletions(-)
diff --git a/dune/fufem/symmetricmatrix.hh b/dune/fufem/symmetricmatrix.hh
index c4913f1f..dbc351f1 100644
--- a/dune/fufem/symmetricmatrix.hh
+++ b/dune/fufem/symmetricmatrix.hh
@@ -50,32 +50,34 @@ public:
return *this;
}
- /** \brief Matrix style random read/write access to components
+ /** \brief Random 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.
+ * The off-diagonal entries are scaled by sqrt(2) to ensure isometry
*/
- T& operator() (int i, int j)
+ void setEntry(int i, int j, const T& entry)
{
- if( i >= j )
- return data_[i*(i+1)/2+j];
- else
- return data_[j*(j+1)/2+i];
+ if( i > j )
+ data_[i*(i+1)/2+j] = std::sqrt(2.) * entry;
+ else if ( i < j )
+ data_[j*(j+1)/2+i] = std::sqrt(2.) * entry;
+ else // ( i == j )
+ data_[i*(i+1)/2+i] = entry;
}
/** \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.
+ * The off-diagonal entries are scaled by 1/sqrt(2) to ensure isometry
*/
- const T& operator() (int i, int j) const
+ const T operator() (int i, int j) const
{
- if( i >= j )
- return data_[i*(i+1)/2+j];
- else
- return data_[j*(j+1)/2+i];
+ if( i > j )
+ return 1./std::sqrt(2.) * data_[i*(i+1)/2+j];
+ else if ( i < j )
+ return 1./std::sqrt(2.) * data_[j*(j+1)/2+i];
+ else // ( i == j )
+ return data_[i*(i+1)/2+i];
}
T energyScalarProduct (const FieldVector<T,N>& v1, const FieldVector<T,N>& v2) const
@@ -92,14 +94,11 @@ public:
data_.axpy(a,other.data_);
}
- /** \brief Compute the Frobenius norm of the matrix */
+ /** \brief Compute the Frobenius norm of the matrix
+ * Make use the isometry */
T frobenius_norm2() const
{
- T result = 0;
- for (size_t i=0; i<N; ++i)
- for (size_t j=0; j<N; ++j)
- result += operator()(i,j)*operator()(i,j);
- return result;
+ data_.two_norm2();
}
/** \brief Compute the Frobenius norm of the matrix */
--
GitLab