Commit efd45e10 authored by Patrick Jaap's avatar Patrick Jaap

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
parent 409103ac
......@@ -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 */
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment