Let's say that I have a struct named like this one:
struct OneBodyMatrixElements {
static const bool is_selfadjoint = 1;
static const bool property2 = 0;
// ........
static value_type val( const std::size_t i, const std::size_t j ) {
value_type ret = 0.;
if( i == j ) ret = 1.;
if( i == j + 1 || i == j - 1 ) ret = -1.;
if( i == 0 && j == 6 ) ret = -1.;
if( i == 6 && j == 0 ) ret = -1.;
return static_cast<value_type>(ret);
}
};
The static variable define different properties of this class and the static function val return the value depending on (i,j). This actually represents a matrix of size 7 whose scalar type is value_type (either float, double, std::complex).
In another part of my code, I need to loop over (i,j) for i,j=0,...,N (here for example N=7):
for( std::size_t i = 0; i < N; i++ ) {
for( std::size_t j = 0; j < N; j++ ) {
if( "OneBodyMatrixElements::val( i, j ) is non zero" ) {
//...
OneBodyMatrixElements::val( i, j ) * func( i, j );
//...
}
}
I have an if close to check if the OneBodyMatrixElements is zero because the multiplication just below is expensive and does not need to be computed if the prefactor is zero.
My question is the following: Since I know the properties of OneBodyMatrixElements, i.e. selfadjoint or diagonal, I also know the location of zero elements, how can I optimise the two for loops to iterate only over non zero elements and also take into account the selfadjoint property?
The obvious answer to this question would be to use Eigen SparseMatrix but I is not what I want because
I did have a look to the visitor pattern but I'm not sure that this is what is need. Any comments on the code is also welcome and a solution which actually store the values in the memory in a sparse form is acceptable. A solution involving only iterations over non-zero elements is also acceptable. Thank you.
EDIT: Answer:
Thank you very much for your answers. I did implement a custom iterator as explained below. For the sake of completeness, here is my version:
class tridiag_iterator {
std::size_t i, j;
public:
struct value_type {
std::size_t i, j;
};
tridiag_iterator( std::size_t i_, std::size_t j_ ) : i( i_ ), j( j_ ) { }
value_type operator*() const {
return { i, j };
}
bool operator!=( const tridiag_iterator &other ) {
if( i >= other.i && j > other.j ) return 0;
else return 1;
}
tridiag_iterator &operator++() {
if( j == i ) j = i + 1;
else if( j == i - 1 ) j = i;
else if( j == i + 1 ) {
i = i + 1;
j = i - 1;
}
}
};
and I added two members to OneBodyMatrixElements:
static tridiag_iterator begin() {
return tridiag_iterator( 0, 0 );
}
static tridiag_iterator end() {
return tridiag_iterator( N - 1, N - 1 );
}
with N defined elsewhere. Usage:
for( auto it = OneBodyMatrixElements::begin(), end = OneBodyMatrixElements::end(); it != end; ++it ) {
OneBodyMatrixElements::val( (*it).i, (*it).j ) * func( (*it).i, (*it).j );
}
Any comments are still welcome.
You could have OneBodyMatrixElements provide an iterator over the non-zero elements. Something like this:
struct Iterator
{
struct value_type
{
std::size_t i, j;
OneBodyMatrixElements::value_type val;
};
Iterator() : i(7), j(7) {}
Iterator(std::size_t i, std::size_t j) : i(i), j(j) {}
value_type operator* () const
{
return {i, j, OneBodyMatrixElements::val(i, j)};
}
Iterator& operator++ ()
{
if (j == i || j == i - 1) j = i + 1;
else if (j == i + 1) { i = i +1; j = i - 1; }
else :::
}
private:
std::size_t i, j;
};
struct OneBodyMatrixElements
{
static Iterator begin()
{
return Iterator(0, 0);
}
static Iterator end()
{
return Iterator();
}
};
// Usage:
for (auto it = OneBodyMatrixElements::begin(), end = OneBodyMatrixElements::end(); it != end; ++it)
{
it->val * func(it->i, it->j);
}
The code above is intended to illustrate the idea, not provide a 100% ready solution. Other operator definitions, typedefs for the iterator, etc. have to be added.
Collected from the Internet
Please contact [email protected] to delete if infringement.
Comments