| View previous topic :: View next topic |
| Author |
Message |
Akkara
Bodhisattva
Joined: 28 Mar 2006
Posts: 6702
Location: &akkara
|
 Posted: Tue Feb 26, 2019 10:22 pm Post subject: C++: 'constexpr' initialize a matrix |
 |
|
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!
_________________
Many think that Dilbert is a comic. Unfortunately it is a documentary. |
|
| Back to top |
|
 |
Hu
Moderator
Joined: 06 Mar 2007
Posts: 21603
|
| Posted: Wed Feb 27, 2019 3:10 am Post subject: |
|
|
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 |
|
|
Akkara
Bodhisattva
Joined: 28 Mar 2006
Posts: 6702
Location: &akkara
|
| Posted: Thu Feb 28, 2019 11:33 am Post subject: |
|
|
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]
_________________
Many think that Dilbert is a comic. Unfortunately it is a documentary. |
|
| Back to top |
|
|
Hu
Moderator
Joined: 06 Mar 2007
Posts: 21603
|
| Posted: Fri Mar 01, 2019 3:31 am Post subject: |
|
|
| 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 |
|
|
Akkara
Bodhisattva
Joined: 28 Mar 2006
Posts: 6702
Location: &akkara
|
| Posted: Wed Mar 13, 2019 5:22 am Post subject: |
|
|
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?
| 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.)
_________________
Many think that Dilbert is a comic. Unfortunately it is a documentary. |
|
| Back to top |
|
|
Hu
Moderator
Joined: 06 Mar 2007
Posts: 21603
|
| Posted: Thu Mar 14, 2019 4:33 am Post subject: |
|
|
| 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 |
|
|
Akkara
Bodhisattva
Joined: 28 Mar 2006
Posts: 6702
Location: &akkara
|
| Posted: Fri Mar 15, 2019 2:04 am Post subject: |
|
|
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.
_________________
Many think that Dilbert is a comic. Unfortunately it is a documentary. |
|
| Back to top |
|
|
Akkara
Bodhisattva
Joined: 28 Mar 2006
Posts: 6702
Location: &akkara
|
| Posted: Fri Mar 15, 2019 9:22 am Post subject: |
|
|
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?
_________________
Many think that Dilbert is a comic. Unfortunately it is a documentary. |
|
| Back to top |
|
|
Hu
Moderator
Joined: 06 Mar 2007
Posts: 21603
|
| Posted: Sat Mar 16, 2019 2:08 am Post subject: |
|
|
| 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 |
|
|
Akkara
Bodhisattva
Joined: 28 Mar 2006
Posts: 6702
Location: &akkara
|
| Posted: Sun Mar 17, 2019 5:22 am Post subject: |
|
|
| 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.
_________________
Many think that Dilbert is a comic. Unfortunately it is a documentary. |
|
| Back to top |
|
|
Akkara
Bodhisattva
Joined: 28 Mar 2006
Posts: 6702
Location: &akkara
|
| Posted: Sun Mar 17, 2019 5:22 am Post subject: |
|
|
| Quote: | | So write it the way I showed in my second post: [...] |
 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)
|
_________________
Many think that Dilbert is a comic. Unfortunately it is a documentary. |
|
| Back to top |
|
|
Hu
Moderator
Joined: 06 Mar 2007
Posts: 21603
|
| Posted: Sun Mar 17, 2019 4:49 pm Post subject: |
|
|
| 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 |
|
|
|
|
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
|
|
Portage & Programming
