diff --git a/dune/fufem/symmetricmatrix.hh b/dune/fufem/symmetricmatrix.hh index c4913f1f61883d916aa70a5db4940135a07c635c..dbc351f18b330a74f18b96823fe13fd7e0bb32cf 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 */