Gentoo Forums
Gentoo Forums
Gentoo Forums
Quick Search: in
C++: 'constexpr' initialize a matrix
View unanswered posts
View posts from last 24 hours

 
Reply to topic    Gentoo Forums Forum Index Portage & Programming
View previous topic :: View next topic  
Author Message
Akkara
Administrator
Administrator


Joined: 28 Mar 2006
Posts: 6605
Location: &akkara

PostPosted: Tue Feb 26, 2019 10:22 pm    Post subject: C++: 'constexpr' initialize a matrix Reply with quote

How does one go about initializing a matrix such that it creates a array in memory, filled with the specified numbers? I have tried searching the net, and find numerous variations of the following, using std::initializer_list:

g++ -std=c++14 -O -Wall -S -o - try-matrix-1.c++:
#include <initializer_list>

/* A fixed-size matrix with R rows and C columns of elements having type T.
 * Each row is stored in consecutive addresses.  Stride S advances from one
 * row to the element in the same column the next row down.
 */
template <typename T, unsigned R, unsigned C, unsigned S = C>
struct  Matrix
{
    T           _elem[R][S];

    constexpr   Matrix(std::initializer_list< std::initializer_list<T> > x)
    {
        for(unsigned r = 0;  r < x.size();  ++r) {
            for(unsigned c = 0;  c < x.begin()[r].size();  ++c) {
                _elem[r][c] = x.begin()[r].begin()[c];
            }
        }
    }
};

using   Mat2x2  = Matrix<int, 2, 2>;

Mat2x2  m       {{1,2},{3,4}};
This has some problems:
  • In C++-11, constexpr constructors can only initialize variables using the '':' construct outside the body. The body must be empty.

  • This compiles in C++14 and onwards, but produces horrible code. It doesn't look at all like a 'constexpr'. I was expecting to see four storage locations initialized to their values. Instead it generates a long procedure that constructs the matrix from a private copy of the numbers stored in an array of pointers to arrays.

I'm now trying various template solutions. I have a row initialized but having trouble extending that to the full matrix. Here's a helper class handling one row:
g++ -std=c++11 -O -Wall -S -o - try-matrix-2.c++:
template <typename T, unsigned N>
struct  MatrixRow
{
    T           _elem[N];

    // Version 1: doesn't compile: invalid conversion from ‘T*’ to ‘T’
//  constexpr   MatrixRow(T (&&x)[N])   : _elem{x}      {}
   
    // Version 2: compiles but only initializes the first element
//  constexpr   MatrixRow(T (&&x)[N])   : _elem{*x}     {}
   
    // can integer_sequence be used somehow?
//  template <int... seq>
//  constexpr   MatrixRow(T (&&x)[N])   : _elem{x[seq]...}      {}

    // This one seems to work
    template <typename... Ta>
    constexpr   MatrixRow(Ta... x)      : _elem{x...}   {}
};

// Version 1,2 uses this constructor
// MatrixRow<int, 3>    try_row{{5,6,7}};

MatrixRow<int, 3>       try_row{5,6,7};


And here's what I have so far of a matrix done as a series of rows:
g++ -std=c++11 -O -Wall -S -o - try-matrix-3.c++:
template <typename T, unsigned N>
struct  MatrixRow
{
    T           _elem[N];

    template <typename... Ta>
    constexpr   MatrixRow(Ta... x)
    : _elem{x...}
    {}
};

template <typename T, unsigned Nrows, unsigned Ncolumns, unsigned Stride = Ncolumns>
struct  Matrix
{
    using       Trow    = MatrixRow<T,Stride>;

    Trow        _row[Nrows];

    constexpr   Matrix(Trow row0)
    : _row{row0}
    {}

    constexpr   Matrix(Trow row0, Trow row1) // Ugly to list all possibilities
    : _row{row0, row1}
    {}

    // Something like this would be nice but can't get it working yet
    template <typename... Tr>
    constexpr   Matrix(Tr... rows)
    : _row{rows...}
    {
    }
};

using   Mat3x3  = Matrix<int, 3, 3>;

Mat3x3  m1      {{1,2,3}};
Mat3x3  m2      {{1,2,3},{4,5,6}};

//Mat3x3        m3      {{1,2,3},{4,5,6},{7,8,9}};

Thanks!
_________________
Isn't it odd that no job openings ever list being potty-trained as one of the requirements?
Back to top
View user's profile Send private message
Hu
Moderator
Moderator


Joined: 06 Mar 2007
Posts: 13512

PostPosted: Wed Feb 27, 2019 3:10 am    Post subject: Reply with quote

You're working too hard. Skip the constructors and use aggregate initialization.
Code:
#include <array>

template <typename T, unsigned R, unsigned C, unsigned S = C>
struct  Matrix
{
   using row_type = std::array<T, S>;
   std::array<row_type, R> _elem;
};

using   Mat2x2  = Matrix<int, 2, 2>;
Mat2x2  m      ={Mat2x2::row_type{1,2},Mat2x2::row_type{3,4}};
As a bonus, this works in C++11 mode.
Back to top
View user's profile Send private message
Akkara
Administrator
Administrator


Joined: 28 Mar 2006
Posts: 6605
Location: &akkara

PostPosted: Thu Feb 28, 2019 11:33 am    Post subject: Reply with quote

Looks like I forgot one could do that :) Thanks! (again!)

What had me thinking about constructors, is this is meant as the start of a "constexpr"-capable matrix class, with the usual matrix operators + - and *. Constexpr inverse would be nice, if possible (ha!), but I'll be happy with just the basics for now. I thought I would need some sort of constructors while working up to those operators - and if I recall - the presence of any constructors in a class makes direct initialization of that class impossible. (Just tested this: defining a constructor gives errors in direct assignment.)

Proceeding with direct initialization anyway, my next thought was to write a "submatrix" operation. The idea being, to try defining +, -, * on 2x2 matrices then write the larger matrix ops recursively in terms of matrices of half the size. (Deferring for later how to handle edge cases when presented with an odd-sized matrix.) To make a long story short, that didn't seem to work. Apparently, can't reinterpret_cast a constexpr:
Quote:
error: accessing value of ‘m1.Matrix<int, 4u, 4u>::_elem.std::array<std::array<int, 4u>, 4u>::_M_elems[1u].std::array<int, 4u>::_M_elems[1u]’ through a ‘const Matrix<int, 2u, 2u, 4u>’ glvalue in a constant expression

This is the dead-end code:
g++ -std=c++11 -O2 -Wall -S -o - matrix.c++:
#include <array>

template <typename T, unsigned R, unsigned C, unsigned S = C>
struct  Matrix
{   
    using row_type = std::array<T, S>;
    std::array<row_type, R> _elem;

    constexpr T const &operator()(int r, int c) const
    {
        return _elem[r][c];   
    }

    template<unsigned Rsub, unsigned Csub>
    constexpr Matrix<T, Rsub, Csub, S> const &submatrix(int r0, int c0)
    {
        return reinterpret_cast<Matrix<T, Rsub, Csub, S> const &>((*this)(r0,c0));
    }
};
     
using   Mat4x4  = Matrix<int, 4, 4>;

constexpr Mat4x4 m1 = {
    Mat4x4::row_type{1,2,3,4},
    Mat4x4::row_type{5,6,7,8},
    Mat4x4::row_type{9,10,11,12},
    Mat4x4::row_type{13,14,15,16}
};

static_assert(  m1(1,1) == 6,   "m1(1,1) should be 6"); // OK

//constexpr auto m2 = m1.submatrix<2,2>(1,1);   // Error

auto    m2 = m1.submatrix<2,2>(1,1); // Works, but isn't constexpr, and uses initialization procedure

//static_assert(  m2(1,1) == 11,   "m2(1,1) should be 11"); // m2 isn't constexpr


So then I tried a different approach, using integer_sequence to write matrix addition directly. It took some doing to find a workable way and get all the syntax right, but this seems to work. It also uses make_integer_sequence, which is a C++14 thing. But it's not hard to write one so I'm not worried about it. (If I do need to write one, is it considered bad practice to put it in the std:: namespace if it duplicates functionality that would otherwise be there?) (Also, how is it possible that integer_sequence made it into C++11, yet they somehow forgot make_integer_sequence? And that's not the only one, witness also unique_ptr and make_unique.)

g++ -std=c++14 -O2 -Wall -S -o - matrix.c++:
#include <array>
#include <utility>

template <typename T, unsigned R, unsigned C, unsigned S = C,
        typename rowseq = std::make_integer_sequence<int, R>,
        typename colseq = std::make_integer_sequence<int, C>
    >
struct  Matrix;

template <typename T, unsigned R, unsigned C, unsigned S, int... rowseq, int... colseq>
struct  Matrix<T, R, C, S,
        std::integer_sequence<int, rowseq...>,
        std::integer_sequence<int, colseq...>
    >
{   
    using row_type = std::array<T, S>;
    std::array<row_type, R> _elem;
   
    constexpr T const &operator()(int r, int c) const
    {
        return _elem[r][c];
    }


    static constexpr row_type addrow(row_type const &a, row_type const &b)
    {
        return row_type{(a[colseq] + b[colseq])...};
    }

    constexpr Matrix    operator+(Matrix const &b) const
    {
        return Matrix{addrow(_elem[rowseq], b._elem[rowseq])...};
    }
};
     
using   Mat4x4  = Matrix<int, 4, 4>;

constexpr Mat4x4 m1 = {
    Mat4x4::row_type{1,2,3,4},
    Mat4x4::row_type{5,6,7,8},
    Mat4x4::row_type{9,10,11,12},
    Mat4x4::row_type{13,14,15,16}
};

static_assert(  m1(1,1) == 6,   "m1(1,1) should be 6"); // OK

constexpr Mat4x4 m2 = m1 + m1;
static_assert(  m2(2,1) == 20,  "m2(2,1) should be 20"); // OK

extern const Mat4x4     m3 = m1 + m1;   // OK: allocated in .rodata, with no procedures

Multiplication is expected to turn out to be a lot harder to massage into constexpr-friendly terms. Hopefully the rowseq... and colseq... theme can be made to work for it as well. I'll need a way of doing dot-products of two vectors. The '...' thing can multiply corresponding elements of two arrays but it is less clear how to sum them all into a single result.

[to be continued]
_________________
Isn't it odd that no job openings ever list being potty-trained as one of the requirements?
Back to top
View user's profile Send private message
Hu
Moderator
Moderator


Joined: 06 Mar 2007
Posts: 13512

PostPosted: Fri Mar 01, 2019 3:31 am    Post subject: Reply with quote

Akkara wrote:
I thought I would need some sort of constructors while working up to those operators - and if I recall - the presence of any constructors in a class makes direct initialization of that class impossible.
You may or may not need the constructors for those operators. You could have them return their results using aggregate syntax too, which should avoid the need for a constructor. Nonetheless, if you want them, and you are willing to tolerate slightly worse initialization syntax, you can get them easily enough:
Code:
#include <array>

template <typename T, unsigned R, unsigned C, unsigned S = C>
struct  Matrix
{
   using row_type = std::array<T, S>;
   using array_type = std::array<row_type, R>;
   array_type _elem;
   constexpr Matrix(const array_type &e) :
      _elem(e)
   {
   }
   constexpr Matrix(std::nullptr_t) : _elem() {}
};

using   Mat2x2  = Matrix<int, 2, 2>;
extern const Mat2x2 m, n;
constexpr Mat2x2  m = {Mat2x2::array_type{Mat2x2::row_type{1,2},Mat2x2::row_type{3,4}}};
constexpr Mat2x2 n = nullptr;
You can add other constructors, including non-constexpr ones, as needed. I included a dummy that explicitly zero-fills the array in response to the otherwise meaningless use of nullptr.

Changes from the prior version:
  • Add convenience typedef array_type.
  • Add constexpr constructor taking a reference to a ready-to-use data array.
  • Change the initialization to construct an anonymous array_type and pass that to the constructor. For me, gcc still manages to evaluate all that at compile time and emit a pre-initialized object.

Akkara wrote:
Apparently, can't reinterpret_cast a constexpr:
Right. constexpr is intentionally a constrained subset of the language because people expect it to all be compile-time solved.
Akkara wrote:
It also uses make_integer_sequence, which is a C++14 thing. But it's not hard to write one so I'm not worried about it. (If I do need to write one, is it considered bad practice to put it in the std:: namespace if it duplicates functionality that would otherwise be there?)
Yes, adding your own symbols to namespace std is discouraged. I would write it in a separate namespace (possibly the top level one). Then, if C++14 is available, use a using declaration to alias the std::make_integer_sequence implementation instead of using yours.
Code:
#ifdef HAVE_CXX14_MAKE_INTEGER_SEQUENCE
using std::make_integer_sequence;
#else
template <typename T, unsigned S>
class make_integer_sequence
{
   /* your implementation and helpers here */
};
#endif
Akkara wrote:
(Also, how is it possible that integer_sequence made it into C++11, yet they somehow forgot make_integer_sequence? And that's not the only one, witness also unique_ptr and make_unique.)
I don't know. Perhaps they underestimated how popular those features would be outside the implementation of the standard library.
Akkara wrote:
g++ -std=c++14 -O2 -Wall -S -o - matrix.c++:
template <typename T, unsigned R, unsigned C, unsigned S = C,
        typename rowseq = std::make_integer_sequence<int, R>,
        typename colseq = std::make_integer_sequence<int, C>
    >
Personally, I would not make rowseq and colseq part of the main signature. Instead, I would make them available only in the specific functions where you need them. With a bit of care, you could even hoist inner helpers out far enough that they can be reused by different matrix types that have the same T and C.
Code:
#include <array>
#include <utility>

template <typename T, unsigned R, unsigned C, unsigned S = C>
struct  Matrix
{
   using row_type = std::array<T, S>;
   using array_type = std::array<row_type, R>;
   array_type _elem;

   constexpr T const &operator()(int r, int c) const
   {
      return _elem[r][c];
   }

   template <std::size_t... colseq>
   static constexpr row_type addrow(row_type const &a, row_type const &b, std::index_sequence<colseq...>)
   {
      return row_type{(a[colseq] + b[colseq])...};
   }

   template <std::size_t... rowseq>
   static constexpr Matrix addmatrix(Matrix const &a, Matrix const &b, std::index_sequence<rowseq...>)
   {
      return Matrix{addrow(a._elem[rowseq], b._elem[rowseq], std::make_index_sequence<C>())...};
   }
   constexpr Matrix    operator+(Matrix const &b) const
   {
      return addmatrix(*this, b, std::make_index_sequence<R>());
   }
};

using   Mat2x2  = Matrix<int, 2, 2>;
extern const Mat2x2 m, n;
constexpr Mat2x2  m = {Mat2x2::array_type{Mat2x2::row_type{1,2},Mat2x2::row_type{3,4}}};
constexpr Mat2x2 n = m + m;
Akkara wrote:
Multiplication is expected to turn out to be a lot harder to massage into constexpr-friendly terms. Hopefully the rowseq... and colseq... theme can be made to work for it as well. I'll need a way of doing dot-products of two vectors. The '...' thing can multiply corresponding elements of two arrays but it is less clear how to sum them all into a single result.
You could use a C++17 fold operator, or recursive descent peeling off and adding values at each layer.
Back to top
View user's profile Send private message
Akkara
Administrator
Administrator


Joined: 28 Mar 2006
Posts: 6605
Location: &akkara

PostPosted: Wed Mar 13, 2019 5:22 am    Post subject: Reply with quote

Found some time to work on this some more. I have what appears to be a working matrix multiply.

I took your advice, Hu, about not making rowseq and colseq part of the main class signature. (Thanks!) I also liked the idea of passing a unused parameter as a way to get the integer_sequence template filled out.

I also removed the stride parameter. I had thought it would be helpful when partitioning off sub-matrices in a recursive divide-and-conquer approach, but it turned out not to be useful since one can't reinterpret-cast a constexpr. And I accidentally used the wrong version of the code when starting to prototype multiplication and the std:array changes were lost. I was too far along to stop and revert so I'm just going to use this code to discuss for now.

There's what I have so far. This does matrix addition and multiplication.
matrix.h:
#include <utility>

template <typename T, int Nrow, int Ncol>
struct  Matrix
{
    template <int... seq>
    using       sequence= std::integer_sequence<int, seq...>;

    using       rowseq  = std::make_integer_sequence<int, Nrow>;
    using       colseq  = std::make_integer_sequence<int, Ncol>;

    struct      Row
    {
        T       _elem[Ncol];

        template <typename... Ta>
        constexpr       Row(Ta... x)    : _elem{x...}   {}

        template <int... seq>   
        static constexpr Row    add(Row const &a, Row const &b, sequence<seq...>)
        {   
            return Row{(a._elem[seq] + b._elem[seq])...};
        }
    };

    Row _row[Nrow];
           
    // Not ideal, will have to do for now
    constexpr   Matrix(Row r)                           : _row{r}               {}   
    constexpr   Matrix(Row r0, Row r1)                  : _row{r0,r1}           {}
    constexpr   Matrix(Row r0, Row r1, Row r2)          : _row{r0,r1,r2}        {}   
    constexpr   Matrix(Row r0, Row r1, Row r2, Row r3)  : _row{r0,r1,r2,r3}     {}

    constexpr   Matrix(Row r0, Row r1, Row r2, Row r3, Row r4) 
    :   _row{r0,r1,r2,r3,r4}
    {}
    constexpr   Matrix(Row r0, Row r1, Row r2, Row r3, Row r4, Row r5) 
    :   _row{r0,r1,r2,r3,r4,r5}
    {}
    constexpr   Matrix(Row r0, Row r1, Row r2, Row r3, Row r4, Row r5, Row r6) 
    :   _row{r0,r1,r2,r3,r4,r5,r6}
    {}
    constexpr   Matrix(Row r0, Row r1, Row r2, Row r3, Row r4, Row r5, Row r6, Row r7) 
    :   _row{r0,r1,r2,r3,r4,r5,r6,r7}
    {}

//      Doesn't work
//    template <typename... Tr>
//    constexpr Matrix(Tr... rows)      : _row{Row(rows)...}    {}

//      Also doesn't work
//    template <typename... Tr>
//    constexpr Matrix(Row r0, Tr... rows)      : _row{r0, Row(rows)...}        {}

    // Fetch element at r,c
    constexpr T const   &operator()(int r, int c)       const
    {
        return _row[r]._elem[c];
    }

    // Matrix add: operates element-by-element
    template <int... seq>
    static constexpr Matrix     add_rows(Matrix const &a, Matrix const &b, sequence<seq...>)
    {
        return Matrix{Row::add(a._row[seq], b._row[seq], colseq())...};
    }

    constexpr Matrix    operator+(Matrix const &rhs)    const
    {
        return add_rows(*this, rhs, rowseq());
    }

    // Sum a list of arguments
    // g++-5.4 missing fold-expressions even with -std=c++17 flag (g++6 and above is fine)
    // clang-4.0 doesn't accept -std=c++17 flag
    // Testing this macro seems to catch all the cases
#if __cpp_fold_expressions
    template <typename... Ts>
    static constexpr T  sum(Ts... args)
    {
        return (0 + ... + args);
    }
#else
    template <typename... Ts>
    static constexpr T  sum(T first, Ts... rest)
    {
        return first + sum(rest...);
    }

    static constexpr T  sum()
    {
        return 0;
    }
#endif

    // Scalar sum of element-by-element products of row of a times column of b
    template <int Ndot, int... seq>
    static constexpr T          row_dot_col(Matrix<T,Nrow,Ndot> const &a, Matrix<T,Ndot,Ncol> const &b, int row, int col, sequence<seq...>)
    {
        return sum(a._row[row]._elem[seq] * b._row[seq]._elem[col]...);
    }

    // row_dot_col above, iterated over all columns of b giving a row result
    template <int Ndot, int... seq>
    static constexpr Row        matmul_do_row(Matrix<T,Nrow,Ndot> const &a, Matrix<T,Ndot,Ncol> const &b, int row, sequence<seq...>)
    {
        return Row(row_dot_col(a,b,row,seq,std::make_integer_sequence<int,Ndot>())...);
    }

    // matmul_do_row above, iterated over all rows of a giving matrix result
    template <int Ndot, int... seq>
    static constexpr Matrix     matmul(Matrix<T,Nrow,Ndot> const &a, Matrix<T,Ndot,Ncol> const &b, sequence<seq...>)
    {
        return Matrix{matmul_do_row(a,b,seq,colseq())...};
    }

    // Product of two matrices giving this class template instance as a result
    template <int Ndot>
    static constexpr Matrix     matmul(Matrix<T,Nrow,Ndot> const &a, Matrix<T,Ndot,Ncol> const &b)
    {
        return matmul(a,b,rowseq());
    }
};

template <typename T, int Nrow, int Ncol, int Ndot>
constexpr Matrix<T,Nrow,Ncol>   operator*(Matrix<T,Nrow,Ndot> const &a, Matrix<T,Ndot,Ncol> const &b)
{
    return Matrix<T,Nrow,Ncol>::matmul(a, b);
}


matrix-test.c++:
#include "matrix.h"

using   Matrix1x1       = Matrix<int, 1, 1>;
using   Matrix1x2       = Matrix<int, 1, 2>;
using   Matrix1x3       = Matrix<int, 1, 3>;
using   Matrix1x4       = Matrix<int, 1, 4>;

using   Matrix2x1       = Matrix<int, 2, 1>;
using   Matrix2x2       = Matrix<int, 2, 2>;
using   Matrix2x3       = Matrix<int, 2, 3>;
using   Matrix2x4       = Matrix<int, 2, 4>;

using   Matrix3x1       = Matrix<int, 3, 1>;
using   Matrix3x2       = Matrix<int, 3, 2>;
using   Matrix3x3       = Matrix<int, 3, 3>;
using   Matrix3x4       = Matrix<int, 3, 4>;

using   Matrix4x1       = Matrix<int, 4, 1>;
using   Matrix4x2       = Matrix<int, 4, 2>;
using   Matrix4x3       = Matrix<int, 4, 3>;
using   Matrix4x4       = Matrix<int, 4, 4>;

#define CHECK_ENTRY(mat,row,col, xpect) static_assert(  \
        mat(row,col) == xpect,                          \
        #mat "(" #row "," #col ") is not " #xpect)

#define CHECK_ROW2(mat,row, c0,c1)      \
        CHECK_ENTRY(mat,row,0,c0);      \
        CHECK_ENTRY(mat,row,1,c1)       \

#define CHECK_ROW3(mat,row, c0,c1,c2)   \
        CHECK_ROW2(mat,row, c0,c1);     \
        CHECK_ENTRY(mat,row,2,c2)

#define CHECK_ROW4(mat,row, c0,c1,c2,c3)\
        CHECK_ROW3(mat,row, c0,c1,c2);  \
        CHECK_ENTRY(mat,row,3,c3)

constexpr Matrix4x4     x {
    { 1, 2, 3, 4},
    { 5, 6, 7, 8},
    { 9,10,11,12},
    {13,14,15,16}
};

CHECK_ROW4(x,0,          1, 2, 3, 4);
CHECK_ROW4(x,1,          5, 6, 7, 8);
CHECK_ROW4(x,2,          9,10,11,12);
CHECK_ROW4(x,3,         13,14,15,16);

constexpr Matrix4x4     y {
    { 1,-1, 2,-2},
    { 1, 2,-1,-2},
    { 2,-1,-2, 1},
    {-1,-2, 2,-1}
};

constexpr Matrix4x4     z {
    {15, -6, 3,-15},
    { 5, -2,-9,-15},
    {15,-12,-9,-15},
    { 5,-14,-3,-15}
};

// Check matrix addition
constexpr Matrix4x4     xplusy  = x + y;

CHECK_ROW4(xplusy,0,     2, 1, 5, 2);
CHECK_ROW4(xplusy,1,     6, 8, 6, 6);
CHECK_ROW4(xplusy,2,    11, 9, 9,13);
CHECK_ROW4(xplusy,3,    12,12,17,15);

// Check matrix multiplication
constexpr Matrix4x4     xy      = x * y;

CHECK_ROW4(xy,0,          5, -8,  2, -7);
CHECK_ROW4(xy,1,         17,-16,  6,-23);
CHECK_ROW4(xy,2,         29,-24, 10,-39);
CHECK_ROW4(xy,3,         41,-32, 14,-55);


// Matrix multiplication is not commutative
constexpr Matrix4x4     yx      = y * x;

CHECK_ROW4(yx,0,        -12,-12,-12,-12);
CHECK_ROW4(yx,1,        -24,-24,-24,-24);
CHECK_ROW4(yx,2,         -8, -8, -8, -8);
CHECK_ROW4(yx,3,         -6, -8,-10,-12);


// z is 30 * the inverse of y
constexpr Matrix4x4     yz      = y * z;

CHECK_ROW4(yz,0,        30, 0, 0, 0);
CHECK_ROW4(yz,1,         0,30, 0, 0);
CHECK_ROW4(yz,2,         0, 0,30, 0);
CHECK_ROW4(yz,3,         0, 0, 0,30);

constexpr Matrix4x4     zy      = z * y;

CHECK_ROW4(zy,0,        30, 0, 0, 0);
CHECK_ROW4(zy,1,         0,30, 0, 0);
CHECK_ROW4(zy,2,         0, 0,30, 0);
CHECK_ROW4(zy,3,         0, 0, 0,30);


// Check products of non-square arguments
constexpr Matrix4x3     p {{1,2,3},{4,5,6},{7,8,9},{10,11,12}};
constexpr Matrix3x4     q {{1,2,3,4},{5,6,7,8},{9,10,11,12}};
constexpr Matrix3x2     r {{3,2},{1,6},{5,4}};


constexpr Matrix4x4     pq      = p * q;

CHECK_ROW4(pq,0,         38,  44,  50,  56);
CHECK_ROW4(pq,1,         83,  98, 113, 128);
CHECK_ROW4(pq,2,        128, 152, 176, 200);
CHECK_ROW4(pq,3,        173, 206, 239, 272);

constexpr Matrix3x3     qp      = q * p;

CHECK_ROW3(qp,0,         70,  80,  90);
CHECK_ROW3(qp,1,        158, 184, 210);
CHECK_ROW3(qp,2,        246, 288, 330);


constexpr Matrix4x2     pr      = p * r;

CHECK_ROW2(pr,0,         20, 26);
CHECK_ROW2(pr,1,         47, 62);
CHECK_ROW2(pr,2,         74, 98);
CHECK_ROW2(pr,3,        101,134);

// Invalid matrix multiplication: columns of r doesn't match rows of p
// clang gives nice error message; g++ (at least 6.3) not so much
// constexpr auto       rp      = r * p;

Compile like one of these or pick something else of your choice:
Code:
g++ -std=c++14 -O -S -o - matrix-test.c++
g++ -std=c++17 -O -S -o - matrix-test.c++
clang -std=c++14 -O -S -o - matrix-test.c++


Questions currently on mind:
  • Is it advisable to overload '*' for matrix multiplication, in light that it is not commutative? Many equations write better using '*' but I don't know if the compiler is free to re-order things.

  • The addition operator was defined inside the class. Multiplication was defined outside the class, mostly because of the complicated matrix-size relation between arguments and result it seemed easier that way. Is there a style reason to prefer doing it one way or the other if there isn't a compelling technical reason forcing it?

  • I have to admit, I prefer the look of the initialization using the Row subclass. It just suffers from severe restrictivitis. Is there a way of getting the std::array version to write so cleanly?
      Sub-question: why does
      Code:
      template <typename... Ta>
      constexpr Row(Ta... x)    : _elem{x...}   {}
      work in the Row class, but a similar idiom doesn't work in the main Matrix class?

  • Will the std::array direct-initialization version still work once this is turned into a "proper" class with private members for the portions that aren't part of the API?

  • What's the declaration to use to cause a constexpr matrix expression to also produce the result as a const matrix in memory, for use at runtime (for example, to further multiply by a runtime-variable matrix)?

  • A big one: This results in a huge amount of code being generated, of very poor quality (load/store from memory greatly dominates the actual arithmetic)
    Code:
    g++ -std=c++17 -x c++ -O2 -S -o - - <<END
    #include "matrix.h"

    Matrix<int,8,8> f(Matrix<int,8,8> const &a, Matrix<int,8,8> const &b)
    {
        return a * b;
    }
    END
    How can I define an overload so that when the operands aren't constexpr, it calls out to a version that uses the ordinary looped version of matrix multiply found in every textbook?

Finally, thanks in advance for whatever assistance may arrive, as well as apologies offered for another huge post and one which probably comes off as slightly rude. I don't mean it that way. But it has been difficult context switching between writing and getting the code to work after many false starts, and now trying to put it all into English. I can only hope if someone else ever needs a constexpr-friendly matrix class, perhaps part of the way might be paved for them. (And I probably need to change the thread title since this is now much more than just about initialization.)
_________________
Isn't it odd that no job openings ever list being potty-trained as one of the requirements?
Back to top
View user's profile Send private message
Hu
Moderator
Moderator


Joined: 06 Mar 2007
Posts: 13512

PostPosted: Thu Mar 14, 2019 4:33 am    Post subject: Reply with quote

Akkara wrote:
matrix.h:
    template <int... seq>
    using       sequence= std::integer_sequence<int, seq...>;
s/integer/index and you can eliminate the type. std::index_sequence is a sequence of std::size_t, but std::integer_sequence is a sequence of whatever type you specify.
Akkara wrote:
matrix.h:
    // Fetch element at r,c
    constexpr T const   &operator()(int r, int c)       const
Personally, I'd use unsigned for indexes instead of int. It makes range-checking easier, if you decide to add that.
Akkara wrote:
matrix.h:
    template <typename... Ts>
    static constexpr T  sum(Ts... args)
    {
        return (0 + ... + args);
    }
If you have a type where there is no well-defined identity value for addition, you can avoid the need for one by requiring at least one argument and using that as your initial term. For here, where 0 can be assumed to be the correct identity, this form is fine.
Akkara wrote:
  • Is it advisable to overload '*' for matrix multiplication, in light that it is not commutative? Many equations write better using '*' but I don't know if the compiler is free to re-order things.
I think overloaded operators are called as written. The compiler cannot swap the parameters.
Akkara wrote:
  • The addition operator was defined inside the class. Multiplication was defined outside the class, mostly because of the complicated matrix-size relation between arguments and result it seemed easier that way. Is there a style reason to prefer doing it one way or the other if there isn't a compelling technical reason forcing it?
Some projects have coding style rules, but I don't know of a good technical reason to pick one when not forced.
Akkara wrote:
  • I have to admit, I prefer the look of the initialization using the Row subclass. It just suffers from severe restrictivitis. Is there a way of getting the std::array version to write so cleanly?
Yes. I think the style I showed in prior posts will work here. Either define no Row constructors and rely on aggregate initialization or define one which takes a const std::array<T, Ncol> & and use copy-initialization.
Akkara wrote:
      Sub-question: why does
      Code:
      template <typename... Ta>
      constexpr Row(Ta... x)    : _elem{x...}   {}
      work in the Row class, but a similar idiom doesn't work in the main Matrix class?
If I were to guess without seeing error messages, it's because the working version does not reiterate the type.
Akkara wrote:
  • Will the std::array direct-initialization version still work once this is turned into a "proper" class with private members for the portions that aren't part of the API?
Private member variables or private methods? If you add private member variables, then the class is not an aggregate.
Akkara wrote:
  • What's the declaration to use to cause a constexpr matrix expression to also produce the result as a const matrix in memory, for use at runtime (for example, to further multiply by a runtime-variable matrix)?
Code:
constexpr Matrix4x4 x = constexpr_initialization_of_Matrix4x4();
Akkara wrote:
  • A big one: ...
    How can I define an overload so that when the operands aren't constexpr, it calls out to a version that uses the ordinary looped version of matrix multiply found in every textbook?
As far as I know, there is no fully-conformant way to test whether an expression is constexpr, although there are some interesting limited purpose tricks for it discussed on Stack Overflow. Generally, you would want to write it as:
Code:
constexpr T compute_result_as_constexpr() {
  ...
}

constexpr std::enable_if<is_constexpr<...>, T>::type compute_result_wrapper() {
   return compute_result_as_constexpr();
}

constexpr std::enable_if<!is_constexpr<...>, T>::type compute_result_wrapper() {
   return compute_result_as_library();
}
Caveat: is_constexpr is not part of the standard library and, as above, I didn't find any fully robust (no false positives and no false negatives, for an arbitrary constexpr expression) implementation. If you can find or make one that is accurate for your use cases, the skeleton above should do the rest.
Back to top
View user's profile Send private message
Akkara
Administrator
Administrator


Joined: 28 Mar 2006
Posts: 6605
Location: &akkara

PostPosted: Fri Mar 15, 2019 2:04 am    Post subject: Reply with quote

Thanks as always for your insightful advice.

Quote:
Personally, I'd use unsigned for indexes instead of int.

I do tend to use unsigned for the most part. But some of the lines can get mighty long with a constexpr up front, a templated-something return, the function name and its template arguments, several unsigned arguments, maybe some const references to additional templated arguments, and often a const and noexcept in the back. It can get hard to see the function name from among all those adjectives :). So I've been trying int to see if that would help any. (You might have noticed, I also tend to line up the function names, to help see at an easier glance what names are being defined.) Also, I recently I read in Google's style guide that they require using 'int' for indexes, and that put some additional doubt in my mind (I don't work there, but I try to look around for ideas and best practices and adopt the ones I find help the most.)

[Edit/note: I figured out some of what follows; see my next post. The friend question (toward the end of the code) is still open.]

Getting back to initialization for a moment, there's something I'm just not understanding. Here's a trivial example class that uses a varargs constructor:
Code:
template <typename T, unsigned N>
class   Vector
{
public:
    Vector()    {}

    constexpr   T       operator()(unsigned i)  const
    {
        return _elem[i];
    }

    template <typename... Ta>   // Varargs constructor
    constexpr   Vector(Ta... args)
    :    _elem  {args...}
    {}

private:
    std::array<T,N>     _elem[N];
};

constexpr Vector<int,4> v4      {1,2,3,4};

// Force allocation so we can see it
Vector<int,4>           x4      = v4;

But when I try the same varargs constructor idea in Matrix, I can't get it to compile.
Code:
#include <utility>
#include <array>

template <typename T, unsigned Nrow, unsigned Ncol>
class   Matrix
{
    template <unsigned... seq>
    using       sequence        = std::integer_sequence<unsigned, seq...>;

    template <unsigned N>
    using       make_sequence   = std::make_integer_sequence<unsigned, N>;

    using       sub_matrix      = Matrix<T, Nrow-1, Ncol>;

    using       Trow            = std::array<T,Ncol>;

public:
    Matrix()    {}

    constexpr   T       operator()(unsigned r, unsigned c)      const
    {
        return _elem[r][c];
    }

    constexpr   Matrix(Trow r0)         // Constructor #1a
    :   _elem   {r0}
    {}

#if 1
    constexpr   Matrix(Trow r0, Trow r1)// Constructor #1b: fine, but limited to 2 rows
    :    Matrix(r0, sub_matrix(r1), make_sequence<Nrow-1>())
    {}
#elif 1
    template <typename... Ta>           // Constructor #1c: try to generalize 1b to more rows (doesn't work)
    constexpr   Matrix(Trow r0, Ta... rows)
    :    Matrix(r0, sub_matrix(rows...), make_sequence<Nrow-1>())
    {}
#else   // I continue to be unpleasantly surprised that something like this was never made standard
    constexpr   Matrix(Trow r0, Trow rows...)
    :    Matrix(r0, sub_matrix(rows...), make_sequence<Nrow-1>())
    {}
#endif

    template <unsigned... seq>          // Constructor #3
    constexpr   Matrix(Trow r0, Matrix<T,Nrow-1,Ncol> mat)
    :   Matrix(r0, mat, make_sequence<Nrow-1>())
    {}

private:
    template <unsigned... seq>          // Constructor #4
    constexpr   Matrix(Trow r0, sub_matrix mat, sequence<seq...>)
    :   _elem   {r0, mat._elem[seq]...}
    {}

// Constructor #4 complains Matrix<int, 1u, 2u>::_elem [1] is private unless _elem is made public.
//      friend  Matrix; // Doesn't work.  "warning: class ‘Matrix<T, Nrow, Ncol>’ is implicitly friends with itself"
//      friend  Matrix<T, Nrow-1, Ncol>;        // Doesn't work either
public:
    Trow        _elem[Nrow];
};

constexpr Matrix<int, 1,2>      m12     {{3,4}};                // OK: uses construtor #1a
constexpr Matrix<int, 2,2>      m22a    {{1, 2}, m12};          // OK: uses construtor #3, calls #4
constexpr Matrix<int, 2,2>      m22b    {{5, 6}, {7, 8}};       // OK: uses construtor #1b, calls #4
                                                                // Fails if try to use constructor #1c
Matrix<int,2,2> x22a    = m22a; // Force allocation so we can see it
Matrix<int,2,2> x22b    = m22b;

g++ -std=c++14 -O2 -Wall -S -o - matrix.c++ as usual.

Comments show most of the questions. I have to run so can't go into more detail now, perhaps add/edit later... Thanks in advance regardless.
_________________
Isn't it odd that no job openings ever list being potty-trained as one of the requirements?
Back to top
View user's profile Send private message
Akkara
Administrator
Administrator


Joined: 28 Mar 2006
Posts: 6605
Location: &akkara

PostPosted: Fri Mar 15, 2019 9:22 am    Post subject: Reply with quote

I think I figured it out. The answer was implicit in one of Hu's earlier comments I just didn't understand it fully at the time. It appears that {...} sequences don't get auto-deduced types in templates (or something like that - I'm probably not being precise here). This change lets it work with Constructor #1c:
Code:
constexpr Matrix<int, 2,2>      m22b    {{5, 6}, std::array<int,2>{{7, 8}}};
I just wish I didn't have to uglyfy the code so much. Typedefs can help but I'd still be repeating the type on every row. There'd be more boilerplate code than actual data.

Tried "hinting" at the type by specifying the first two argument types in the template:
Code:
template <typename T, unsigned Nrow, unsigned Ncol>
class   Matrix
{
    /* ... same as before ... */

    template <typename... Ta>           // Constructor #1c: try to generalize 1b to more rows
    constexpr   Matrix(Trow r0, Trow r1, Ta... rows)
    :    Matrix(r0, sub_matrix(r1, rows...), make_sequence<Nrow-1>())
    {}

    /* ... same as before ... */
};
This now works for the 2x2 case, but still fails for larger number of rows.

Oh, well. :( I think I'll have to stick with direct initialization and no private members. That probably means I can't have matrix constants directly in an expression, is that right? I'll have to declare 'constexpr' variables to hold each different matrix and use those names in the expression instead. Even a compromise solution such as using a series of specific-length constructors for the small-matrix cases (like I did in the multiplication code above) doesn't work: defining any constructor loses the ability to use direct initialization for larger matrices.

It still seems I must be missing something. It can't be that writing an initialized instance is that hard to express cleanly... is it?
_________________
Isn't it odd that no job openings ever list being potty-trained as one of the requirements?
Back to top
View user's profile Send private message
Hu
Moderator
Moderator


Joined: 06 Mar 2007
Posts: 13512

PostPosted: Sat Mar 16, 2019 2:08 am    Post subject: Reply with quote

Akkara wrote:
But some of the lines can get mighty long
I enable virtual line wrapping (Vim: set wrap) and live with the long lines.
Akkara wrote:
Also, I recently I read in Google's style guide that they require using 'int' for indexes, and that put some additional doubt in my mind (I don't work there, but I try to look around for ideas and best practices and adopt the ones I find help the most.)
I read that too. That stood out to me as a bad practice. They make vague claims about unidentified standards body members disliking use of unsigned math, add in some vague claims that the less-defined nature of signed integer overflow improves error checking and optimization relative to the more-defined nature of unsigned integer overflow, and conclude that indices should be signed. They fail to address that range checking is easier with unsigned values. I don't want every range check to need to test for a negative index. In unsigned math, if an integer underflows, it will become an extremely large index and the test for exceeding container size will handle it automatically. Compare:
Code:
if (si < 0 || si >= max_size) { error; }
Code:
if (ui >= max_size) { error; }
Unsigned range checks are shorter and encourage the compiler to generate a single test. In environments where you cannot count on a code review to catch a missing si < 0, I prefer a convention that indices are unsigned.
Akkara wrote:

Getting back to initialization for a moment, there's something I'm just not understanding. Here's a trivial example class that uses a varargs constructor:
But when I try the same varargs constructor idea in Matrix, I can't get it to compile.

#1c fails because you cannot deduce a matching type from a brace-enclosed initializer. In the case of a simple integer (for Vector), you aren't using brace-enclosed initializers, so it deduces a type of int and works. Generally, mixing templates with constructs that require deduction (like brace initializers) will end in frustration.
Akkara wrote:
This change lets it work with Constructor #1c:
Code:
constexpr Matrix<int, 2,2>      m22b    {{5, 6}, std::array<int,2>{{7, 8}}};
I just wish I didn't have to uglyfy the code so much.
So write it the way I showed in my second post:
Code:
#include <utility>
#include <array>

template <typename T, unsigned Nrow, unsigned Ncol>
class   Matrix
{
    using       Trow            = std::array<T,Ncol>;
public:
   using array_type = std::array<Trow, Nrow>;
   constexpr Matrix(const array_type &a) : x(200), _elem(a) {}
private:
   int x;
   array_type _elem;
};

using M22 = Matrix<int, 2,2>;
constexpr M22      m22b    (M22::array_type{{{5, 6}, {7, 8}}});
Matrix<int,2,2> x22b    = m22b;
One extra type::array_type at the outermost scope of the constructor call, instead of repeating the row type. If that's your only constructor, you can even omit the type and let brace-deduction solve it, because the constructor is not a template.
Akkara wrote:
I think I'll have to stick with direct initialization and no private members. That probably means I can't have matrix constants directly in an expression, is that right?
I don't understand how the limitation follows from the lack of private member variables. Whether or not you have constructors, you can declare anonymous immediately-initialized instances for single-shot usage. Could you show an example of the expression you think cannot work?
Back to top
View user's profile Send private message
Akkara
Administrator
Administrator


Joined: 28 Mar 2006
Posts: 6605
Location: &akkara

PostPosted: Sun Mar 17, 2019 5:22 am    Post subject: Reply with quote

Quote:
They fail to address that range checking is easier with unsigned values. I don't want every range check to need to test for a negative index. In unsigned math, if an integer underflows, it will become an extremely large index and the test for exceeding container size will handle it automatically.

I agree with your sentiment. I'd add that using unsigned - or better - size_t is more correct because it lets the size span the entire range of possibilities. Although arguably with 64-bit sizes this is much less important since I don't foresee memory size getting anywhere near (2**63 - 1) bytes, nevermind twice as much. But in 32-bit userland the difference between limiting a size to 2gig vs 4gig can definitely be noticeable.

Regarding the optimization, however, it seems the compiler already knows about it and applies it for us: Both
    Code:
    echo 'bool f(int a) { return a < 0 || a > 100; }'  |  g++ -O2 -x c++ -S -o - -
    and
    Code:
    echo 'bool f(unsigned a) { return a > 100; }'  |  g++ -O2 -x c++ -S -o - -
    gives
    arm::g++-6.3:
       cmp   r0, #100
       movls   r0, #0
       movhi   r0, #1
       bx   lr
It even knows to offset the argument to use the unsigned trick:
    Code:
    echo 'bool f(int a) { return a < 3 || a > 100; }'  |  g++ -O2 -x c++ -S -o - -
    gives
    arm::g++-6.3:
       sub   r0, r0, #3
       cmp   r0, #97
       movls   r0, #0
       movhi   r0, #1
       bx   lr
Any of -O1, -O2, and -Os gives the same result.

The optimization is also used with 'if' statements:
    Code:
    echo 'void g(), h(); void f(int a) { if(a < 3 || a > 100) g(); else h(); }'  |  g++ -O2 -x c++ -S -o - -
    gives
    arm::g++-6.3:
       sub   r0, r0, #3
       cmp   r0, #97
       bls   .L2
       b   _Z1gv
    .L2:
       b   _Z1hv
As before any of the optimization flags (other than -O0) applies the "unsigned compare" transformation, however g++ doesn't do tail-call elimination on -O1 so the code here isn't as good as with -O2 or -Os. Clang does a similar transformation, and it does tail-call elimination on all optimization levels other than 0.
_________________
Isn't it odd that no job openings ever list being potty-trained as one of the requirements?
Back to top
View user's profile Send private message
Akkara
Administrator
Administrator


Joined: 28 Mar 2006
Posts: 6605
Location: &akkara

PostPosted: Sun Mar 17, 2019 5:22 am    Post subject: Reply with quote

Quote:
So write it the way I showed in my second post: [...]

:oops: How did I not remember that one?? I was at a different computer, and instead of connecting and reviewing again, I thought I'd try to do it from memory, to see if it "stuck". Clearly it hadn't. :) Hopefully now it will. That is a fine solution, thanks again and sorry you had to repeat yourself.

Quote:
Akkara wrote:
I think I'll have to stick with direct initialization and no private members. That probably means I can't have matrix constants directly in an expression, is that right?
I don't understand how the limitation follows from the lack of private member variables. Whether or not you have constructors, you can declare anonymous immediately-initialized instances for single-shot usage. Could you show an example of the expression you think cannot work?

I think I mostly worked it out. It was confusion spurred by use of () or {} in making an instance, combined with a suspicion that if there's any constructors present, direct initialization can no longer be used. This code shows what I was exploring:
g++ -O2 -Wall -S -o - c.c++:
struct  Test1 {
    int x, y, z;

    constexpr Test1     operator+(Test1 const &b)       const
    {
        return Test1{x + b.x, y + b.y, z + b.z};
    }
};

struct  Test2 {
    constexpr Test2(int p, int q)       : x(p), y(q), z(0)      {}

    constexpr Test2     operator+(Test2 const &b)       const
    {
        return Test2(x + b.x, y + b.y); // (ignoring z for now)
    }

    int x, y, z;
};

constexpr Test1 a1      = {3,4};        // Requires {}
constexpr Test2 a2      = {3,4};        // Requires {}

constexpr Test1 b1      {5,6};          // Requires {}
constexpr Test2 b2      (5,6);          // Either {} or () work

constexpr Test1 c1      = Test1{3,4} + Test1{5,6};      // Requires {}
constexpr Test2 c2      = Test2{3,4} + Test2(5,6);      // Either {} or () work

constexpr Test1 d1      = a1 + b1;     
constexpr Test2 d2      = a2 + b2;     

constexpr Test1 e1      = {3,4,5};      // Direct initialize the "extra" element ("=" sign is optional)
//constexpr Test2 e2    = {3,4,5};      // Constructor precludes direct initialization

// Force allocation to see them all
Test1   A1      = a1;
Test1   B1      = b1;
Test1   C1      = c1;
Test1   D1      = d1;
Test1   E1      = e1;

Test2   A2      = a2;
Test2   B2      = b2;
Test2   C2      = c2;
Test2   D2      = d2;
//Test2 E2      = e2;

But, what is the difference between using an =, and not, when creating a variable? (As with constants a1, a2, b1, and b2 in the code above.) This is my current understanding of the situation (constexpr omitted for clarity):
Code:
Test1 a1 = {3,4};       // direct member initialization with list of initializers
Test2 a2 = {3,4};       // calls constructor making a temp, which is assigned to a2
Test1 b1 {5,6};         // ??? works anyway (this is what I thought wouldn't work)
Test2 b2 (3,4);         // calls constructor making it directly into b2, no temp involved
Test1{3,4} + Test1{5,6} // two temporary instances of Test1 created using direct member initialization, then added (didn't think that would work, but it does)

_________________
Isn't it odd that no job openings ever list being potty-trained as one of the requirements?
Back to top
View user's profile Send private message
Hu
Moderator
Moderator


Joined: 06 Mar 2007
Posts: 13512

PostPosted: Sun Mar 17, 2019 4:49 pm    Post subject: Reply with quote

Akkara wrote:
Regarding the optimization, however, it seems the compiler already knows about it and applies it for us:
Yes, but the signed form also has the downside that the programmer must write both checks whenever you want a range validation. The compiler won't warn if you forget the < 0 test, and if your project lacks picky code reviewers, the omission may go unnoticed.
Akkara wrote:
That is a fine solution, thanks again and sorry you had to repeat yourself.
It's ok. It wasn't that much to repeat.
Akkara wrote:
It was confusion spurred by use of () or {} in making an instance, combined with a suspicion that if there's any constructors present, direct initialization can no longer be used.
If constructors are present, aggregate initialization is not allowed. However, a syntax that would otherwise be aggregate initialization will instead try to become a constructor call. If that finds a constructor, you get a constructor call. Otherwise, you get a compile error about the lack of an acceptable constructor. This is very convenient, but hardly an application of the principle of least surprise.
Akkara wrote:
But, what is the difference between using an =, and not, when creating a variable?
Without =, you request list initialization. Aggregate initialization is one form of list initialization, but is only used when the left hand side is actually an aggregate type. The full rules are rather lengthy. See https://en.cppreference.com/w/cpp/language/list_initialization#Explanation for the full logic tree of what happens when you request list initialization. I think this covers your examples (not quoted), but if you want further elaboration, post back and I will step through them.
Back to top
View user's profile Send private message
Display posts from previous:   
Reply to topic    Gentoo Forums Forum Index Portage & Programming All times are GMT
Page 1 of 1

 
Jump to:  
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum