-1

I am relatively new to C++ and am building a matrix class for carrying out with it common linear algebra operations. Some snippet of the hpp:

template<const int rows, const int columns, typename T=float>
struct Matrix3D
{
    const int shape[2] = {rows, columns};
    
    float content[rows][columns]{};
    
    Matrix3D()
    {
        for (int i = 0; i < rows; i++)
        {
            for (int j = 0; j < columns; j++)
            {
                this->content[i][j] = 0.;
            }
        }
    }
    
    explicit Matrix3D(const T (&arr)[rows][columns])
    {
        for (int i = 0; i < rows; i++)
        {
            for (int j = 0; j < columns; j++)
            {
                this->content[i][j] = arr[i][j];
            }
        }
    }
    
    bool operator==( const Matrix3D<rows, columns> &mat ) const noexcept
    {
        int equal = 0;
        
        for (int i = 0; i < rows; i++)
        {
            for (int j = 0; j < columns; j++)
            {
                equal += this->content[i][j] == mat.content[i][j];
            }
        }
        return  equal == rows * columns;
    }
    
    bool operator!=( const Matrix3D<rows, columns> &mat ) const noexcept
    {
        int equal = 0;
        
        for (int i = 0; i < rows; i++)
        {
            for (int j = 0; j < columns; j++)
            {
                equal += this->content[i][j] == mat.content[i][j];
            }
        }
        return  equal != rows * columns;
    }

    ...

};

You get the idea: for pretty much all the operations that I will need to manipulate that matrix object I am going to need these double loops in every method.

Is there a way to abstract that, for example, to a static method, and pass the core looping expression as argument? It is worth noting that my solution must be compatible with CUDA, meaning I cannot use std::vector, or Eigen, or similar.

Thanks in advance.

mosegui
  • 664
  • 6
  • 15
  • 3
    The issue more or less solves itself once you solve another issue: using an array of arrays is for several reasons not a suitable representation of matrices or higher-dimensional arrays. Use a one-dimensional array and manually calculate the index into that array. Once you have that, all your functions with nested loops can instead use a single loop over the array, and they can use the standard library algorithms (such as `std::fill`, `std::transform`, `std::equal_range`). And if you use `std::array` instead of C arrays you get even more benefits for free (e.g. copy construction). – Konrad Rudolph Mar 19 '21 at 17:01
  • 1
    It's usually better to use a 1D array and emulate the 2D part of it, that means fewer loops. It also means you can do a simple linear scan for differences or use `std::copy` for cloning. – tadman Mar 19 '21 at 17:02
  • Since all of the loop indexes are constant odds are good the compiler will do whatever unrolling it sees as beneficial for you. You could make a function that does the looping and accepts a function that does the work as an as an argument. Don't know if you\ll gain much other than less typing. – user4581301 Mar 19 '21 at 17:03
  • You could implement a two-indexed range with iterators. With c++17 it'll look like this `for(auto [i,j] : index_range(cols,rows)){...}`. Don't know if compilers can completely optimize this structure away. – ALX23z Mar 19 '21 at 17:03
  • @KonradRudolph and tadman as I stated in my original question I cannot use STL as the code has to be CUDA-compatible – mosegui Mar 20 '21 at 10:22
  • @mosegui You can replace the algorithms by their Thrust equivalent then. In general it’s still a good idea to write your `Matrix3D` code as if you were writing for the C++ standard library, and then supplement the missing parts one by one: the algorithms you need are all trivial to implement, if need be. – Konrad Rudolph Mar 20 '21 at 13:12

1 Answers1

2

You could write a for_each_element member template:

template <typename F>
void for_each_element(F f) {
    for (int i = 0; i < rows; i++)
    {
        for (int j = 0; j < columns; j++)
        {
            f(i,j);
        }
    }
}

Such that for example the operator!= can be written as:

bool operator!=( const Matrix3D<rows, columns> &mat ) const noexcept
    {
        int equal = 0;
        for_each_element( [&](int i,int j) {                                  
                equal += this->content[i][j] == mat.content[i][j];
        });
        return  equal != rows * columns;
    }

However, eg operator== is a good example for why this isnt the best idea: Instead of looping through all elements you should stop as soon as you encounter one that isn't equal. for_each_element does not allow you to do that (at least not in its current form). As others suggested in comments, for a matrix you should rather "abstract away" the second dimension, by using a 1d array plus an index transformation:

T& get(size_t i,size_t j) { return this->content[i*columns + j]; }

where content is a properly sized 1d array.

You can then write most loops as

for (auto& element : content) {
    //...
}

which mitigates the need to make it any shorter.

463035818_is_not_an_ai
  • 109,796
  • 11
  • 89
  • 185