I'm writing a class for Hermitian matrices. This is a complex matrix that has only n*(n+1)/2
independent complex numbers (ignoring details about the diagonal being exactly real).
My plan is to write only the upper triangular elements, where row number compared to column number satisfy the condition satisfy the rule: row >= column
. However, this requires something like a proxy? I'm not sure how to implement this. Here's the problem:
Say I implement the member function at(int row, int column)
to access an element.
template<typename T>
std::complex<T>& HermitianMatrix<T>::at(long row, long column)
{
if(row >= column)
return this->_matrix[ElementIndex(row,column)];
else
return std::conj(this->_matrix[ElementIndex(column,row)]);
}
where ElementIndex
converts the row
and column
input to the the position in the array std::complex<T>* _matrix = new std::complex<T>(...)
. Of course, this method returns a reference. The code you see above doesn't work for the lower triangular part of the matrix because the reference is gone after returning.
What is the right and most efficient way to implement this, such that I have some kind of "pipe" for the lower triangular matrix part always goes through std::conj
for both set and get?
Please ask for more information if required. Thank you.
Following the Franck's example, I propose to return a wrapper class (or struct) that wrap the reference to the element and memorize a boolean flag to remember if it's neccessary to coniugate the number.
Something like [caution: not tested]
template <typename T>
struct cWrapper
{
bool c;
std::complex<T> & r;
cWrapper (bool c0, std::complex<T> & r0) : c{c0}, r{r0}
{ }
operator std::complex<T>() const
{ return c ? std::conj(r) : r; }
cWrapper & operator= (const std::complex<T> & r0)
{
r = ( c ? std::conj(r0) : r0 );
return *this;
}
};
and your function could become [edit: modified after the corresponding edit in the question (row/column inversion for else case)]
template<typename T>
cWrapper<T> HermitianMatrix<T>::at(long row, long column)
{
if(row >= column)
return cWrapper<T>(false, this->_matrix[ElementIndex(row,column)]);
else
return cWrapper<T>(true, this->_matrix[ElementIndex(column,row)]);
}