Commit 7e1e1b3f authored by graeser's avatar graeser

Merge branch 'feature/isometric-symmetric-matrix' into 'master'

Make symmetrixMatrix isometric

See merge request !47
parents 409103ac 09c8b4ac
Pipeline #16288 passed with stage
in 69 minutes and 51 seconds
......@@ -38,6 +38,27 @@ public:
data_[i] = *it;
}
/** \brief Construct identityMatrix */
static constexpr SymmetricMatrix<T,N> identityMatrix()
{
SymmetricMatrix<T,N> id;
for( size_t i=0; i<N*(N+1)/2; ++i)
id.setEntry(i,i,1);
return id;
}
/** \brief return the trace map in vector representation
* the returned vector is such that
* <data,a> = trace(A) where A.data = a
*/
static constexpr Data traceMap()
{
// trace(A) = <Id,A>, therefore return data of Id
// use isometry!
return identityMatrix().data();
}
SymmetricMatrix<T,N>& operator=(const T& s)
{
data_ = s;
......@@ -50,32 +71,36 @@ 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
* \note: You need to know what you are doing: Modifying an off-diagonal
* entry modifies two entries in the matrix!
*/
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
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 +117,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