All geek questions in one place

Template function overloading and SFINAE implementation

I spend some time learning how to use templates in C ++. I have never used them before, and I'm not always sure what can be or what cannot be achieved in another situation.

As an exercise, I wrap some of the Blas and Lapack functions that I use for my actions, and I'm currently working on the ?GELS wrapper (which evaluates the solution to a linear system of equations).

  A x + b = 0

?GELS function (for real values ​​only) exists with two names: SGELS , for vectors with one precision and DGELS for double precision.

My understanding of the interface is the solve function as follows:

  const std::size_t rows = /* number of rows for A */; const std::size_t cols = /* number of cols for A */; std::array< double, rows * cols > A = { /* values */ }; std::array< double, ??? > b = { /* values */ }; // ??? it can be either // rows or cols. It depends on user // problem, in general // max( dim(x), dim(b) ) = // max( cols, rows ) solve< double, rows, cols >(A, b); // the solution x is stored in b, thus b // must be "large" enough to accomodate x

Depending on the user's requirements, the problem may be overridden or indefinite, which means:

  • if it is redefined dim(b) > dim(x) (the solution is pseudo-inverse)
  • if it is not defined dim(b) < dim(x) (the solution is a minimization of LSQ)
  • or the normal case in which dim(b) = dim(x) (the solution is the inverse of A )

(excluding singular cases).

Since ?GELS stores the result in the input vector b , std::array shouold has enough space to accommodate the solution, as described in the code comments ( max(rows, cols) ).

I want (compilation time) to determine which decision to make (is this changing the paraming in the call ?GELS ). I have two functions (I simplify for the sake of the question) that handle precision and already know what size b and number of rows / cols :

 namespace wrap { template <std::size_t rows, std::size_t cols, std::size_t dimb> void solve(std::array<float, rows * cols> & A, std::array<float, dimb> & b) { SGELS(/* Called in the right way */); } template <std::size_t rows, std::size_t cols, std::size_t dimb> void solve(std::array<double, rows * cols> & A, std::array<double, dimb> & b) { DGELS(/* Called in the right way */); } }; /* namespace wrap */

which are part of the inner wrapper. User function, specify the required size in vector b through the patterns:

 #include <type_traits> /** This struct makes the max between rows and cols */ template < std::size_t rows, std::size_t cols > struct biggest_dim { static std::size_t const value = std::conditional< rows >= cols, std::integral_constant< std::size_t, rows >, std::integral_constant< std::size_t, cols > >::type::value; }; /** A type for the b array is selected using "biggest_dim" */ template < typename REAL_T, std::size_t rows, std::size_t cols > using b_array_t = std::array< REAL_T, biggest_dim< rows, cols >::value >; /** Here we have the function that allows only the call with b of * the correct size to continue */ template < typename REAL_T, std::size_t rows, std::size_t cols > void solve(std::array< REAL_T, cols * rows > & A, b_array_t< REAL_T, cols, rows > & b) { static_assert(std::is_floating_point< REAL_T >::value, "Only float/double accepted"); wrap::solve< rows, cols, biggest_dim< rows, cols >::value >(A, b); }

So it actually works . But I want to take one more step, and I really do not know how to do it. If the user tries to call solve with b too small, an extremely difficult to read error occurs by the compiler.

I am trying to insert a static_assert which helps the user understand his error. But any direction that comes to my mind requires the use of two functions with the same signature (does this look like template overloading?), For which I cannot find the SFINAE strategy (and they do not compile at all).

Do you consider it possible to make a static statement for the case of incorrect b measurement without changing the user interface at compile time ? Hope this question is clear enough.

@Caninonos . For me, the user interface is how the user calls the solver, that is:

  solve< type, number of rows, number of cols > (matrix A, vector b)

This is a limitation that I used for my exercises to improve my skills. This means that I do not know if a solution can really be reached. Type b should match the function call, and it is easy if I add another template parameter and I change the user interface, violating my restriction.

Minimal complete and working example

This is a minimal complete and working example. Upon request, I removed any reference to the concepts of linear algebra. This is a number problem. Cases:

  • N1 = 2, N2 =2 . Since N3 = max(N1, N2) = 2 everything works
  • N1 = 2, N2 =1 . Since N3 = max(N1, N2) = N1 = 2 everything works
  • N1 = 1, N2 =2 . Since N3 = max(N1, N2) = N2 = 2 everything works
  • N1 = 1, N2 =2 . Since N3 = N1 = 1 < N2 it correctly causes a compilation error. I want to catch a compilation error with a static statement explaining the fact that the N3 dimension is incorrect. At the moment, the error is difficult to read and understand.

You can view and test it online here.

+5
source share
3 answers

First, a few improvements that simplify the design a bit and help increase its scalability:

  • no biggest_dim . std::max is constexpr sice C ++ 14. You should use it instead.

  • no need for b_array_t . You can simply write std::array< REAL_T, std::max(N1, N2)>

And now to your problem. One good way in C ++ 17:

 template < typename REAL_T, std::size_t N1, std::size_t N2, std::size_t N3> void solve(std::array< REAL_T, N1 * N2 > & A, std::array< REAL_T, N3> & b) { if constexpr (N3 == std::max(N1, N2)) wrap::internal< N1, N2, N3 >(A, b); else static_assert(N3 == std::max(N1, N2), "invalid 3rd dimmension"); // don't write static_assert(false) // this would make the program ill-formed (*) }

Or as indicated by @ max66

 template < typename REAL_T, std::size_t N1, std::size_t N2, std::size_t N3> void solve(std::array< REAL_T, N1 * N2 > & A, std::array< REAL_T, N3> & b) { static_assert(N3 == std::max(N1, N2), "invalid 3rd dimmension"); if constexpr (N3 == std::max(N1, N2)) wrap::internal< N1, N2, N3 >(A, b); }

Tadaa !! Simple, elegant, nice error message.

The difference between the constexpr if version and just static_assert Ie:

 void solve(...) { static_assert(...); wrap::internal(...); }

it is that only with static_assert compiler will try to create an instance of wrap::internal even if static_assert fails, polluting the error output. With constexpr, if the wrap::internal call is not part of the body under the condition of failure, therefore, the error output is clean.

(*) The reason why I did not just write static_asert(false, "error msg) is because it would make the program poorly formed, no diagnostics were required. See constexpr if and static_assert

You can also franchise float / double if you wish by moving the template argument after non-deductible:

 template < std::size_t N1, std::size_t N2, std::size_t N3, typename REAL_T> void solve(std::array< REAL_T, N1 * N2 > & A, std::array< REAL_T, N3> & b) {

Thus, the call becomes:

 solve< n1_3, n2_3>(A_3, b_3);
+2
source

Why don't you try combining tag sending with some static_assert s? Hope this is one way to achieve what you want to solve. I mean that all three correct cases are correctly transferred to the correct blas calls, different types and mismatches of measurements are processed, and violation of the order of float and double also processed, all in a way convenient for the user, thanks to static_assert .

EDIT. I'm not sure about your requirement for a C++ version, but below C++11 .

 #include <algorithm> #include <iostream> #include <type_traits> template <class value_t, int nrows, int ncols> struct Matrix {}; template <class value_t, int rows> struct Vector {}; template <class value_t> struct blas; template <> struct blas<float> { static void overdet(...) { std::cout << __PRETTY_FUNCTION__ << std::endl; } static void underdet(...) { std::cout << __PRETTY_FUNCTION__ << std::endl; } static void normal(...) { std::cout << __PRETTY_FUNCTION__ << std::endl; } }; template <> struct blas<double> { static void overdet(...) { std::cout << __PRETTY_FUNCTION__ << std::endl; } static void underdet(...) { std::cout << __PRETTY_FUNCTION__ << std::endl; } static void normal(...) { std::cout << __PRETTY_FUNCTION__ << std::endl; } }; class overdet {}; class underdet {}; class normal {}; template <class T1, class T2, int nrows, int ncols, int dim> void solve(const Matrix<T1, nrows, ncols> &lhs, Vector<T2, dim> &rhs) { static_assert(std::is_same<T1, T2>::value, "lhs and rhs must have the same value types"); static_assert(dim >= nrows && dim >= ncols, "rhs does not have enough space"); static_assert(std::is_same<T1, float>::value || std::is_same<T1, double>::value, "Only float or double are accepted"); solve_impl(lhs, rhs, typename std::conditional<(nrows < ncols), underdet, typename std::conditional<(nrows > ncols), overdet, normal>::type>::type{}); } template <class value_t, int nrows, int ncols, int dim> void solve_impl(const Matrix<value_t, nrows, ncols> &lhs, Vector<value_t, dim> &rhs, underdet) { /* get the pointers and dimension information from lhs and rhs */ blas<value_t>::underdet( /* trans, m, n, nrhs, A, lda, B, ldb, work, lwork, info */); } template <class value_t, int nrows, int ncols, int dim> void solve_impl(const Matrix<value_t, nrows, ncols> &lhs, Vector<value_t, dim> &rhs, overdet) { /* get the pointers and dimension information from lhs and rhs */ blas<value_t>::overdet( /* trans, m, n, nrhs, A, lda, B, ldb, work, lwork, info */); } template <class value_t, int nrows, int ncols, int dim> void solve_impl(const Matrix<value_t, nrows, ncols> &lhs, Vector<value_t, dim> &rhs, normal) { /* get the pointers and dimension information from lhs and rhs */ blas<value_t>::normal( /* trans, m, n, nrhs, A, lda, B, ldb, work, lwork, info */); } int main() { /* valid types */ Matrix<float, 2, 4> A1; Matrix<float, 4, 4> A2; Matrix<float, 5, 4> A3; Vector<float, 4> b1; Vector<float, 5> b2; solve(A1, b1); solve(A2, b1); solve(A3, b2); Matrix<int, 4, 4> A4; Vector<int, 4> b3; // solve(A4, b3); // static_assert for float & double Matrix<float, 4, 4> A5; Vector<int, 4> b4; // solve(A5, b4); // static_assert for different types // solve(A3, b1); // static_assert for dimension problem return 0; }
+1
source

You must consider why the interface offers this (confusing) mess of parameters. The author had several things. First of all, you can solve problems of the form A x + b == 0 and A^T x + b == 0 in one function. Secondly, data A and b can actually indicate memory in matrices larger than those needed by alg. This can be seen using the LDA and LDB .

This is subaddress that complicates the situation. If you need a simple but possibly quite useful API, you can ignore this part:

 using ::std::size_t; using ::std::array; template<typename T, size_t rows, size_t cols> using matrix = array<T, rows * cols>; enum class TransposeMode : bool { None = false, Transposed = true }; // See https://stackoverflow.com/questions/14637356/ template<typename T> struct always_false_t : std::false_type {}; template<typename T> constexpr bool always_false_v = always_false_t<T>::value; template < typename T, size_t rowsA, size_t colsA, size_t rowsB, size_t colsB , TransposeMode mode = TransposeMode::None > void solve(matrix<T, rowsA, colsA>& A, matrix<T, rowsB, colsB>& B) { // Since the algorithm works in place, b needs to be able to store // both input and output static_assert(rowsB >= rowsA && rowsB >= colsA, "b is too small"); // LDA = rowsA, LDB = rowsB if constexpr (::std::is_same_v<T, float>) { // SGELS(mode == TransposeMode::None ? 'N' : 'T', ....); } else if constexpr (::std::is_same_v<T, double>) { // DGELS(mode == TransposeMode::None ? 'N' : 'T', ....); } else { static_assert(always_false_v<T>, "Unknown type"); } }

Now, turning to subaddressing, perhaps with the help of LDA and LDB . I suggest you make this part of your data type, and not directly part of the template signature. You want to have your own type of matrix that can reference storage in the matrix. Maybe something like this:

 // Since we store elements in a column-major order, we can always // pretend that our matrix has less columns than it actually has // less rows than allocated. We can not equally pretend less rows // otherwise the addressing into the array is off. // Thus, we'd only four total parameters: // offset = columnSkipped * actualRows + rowSkipped), actualRows, rows, cols // We store the offset implicitly by adjusting our begin pointer template<typename T, size_t rows, size_t cols, size_t actualRows> class matrix_view { // Name derived from string_view :) static_assert(actualRows >= rows); T* start; matrix_view(T* start) : start(start) {} template<typename U, size_t r, size_t c, size_t ac> friend class matrix_view; public: template<typename U> matrix_view(matrix<U, rows, cols>& ref) : start(ref.data()) { } template<size_t rowSkipped, size_t colSkipped, size_t newRows, size_t newCols> auto submat() { static_assert(colSkipped + newCols <= cols, "can only shrink"); static_assert(rowSkipped + newRows <= rows, "can only shrink"); auto newStart = start + colSkipped * actualRows + rowSkipped; using newType = matrix_view<T, newRows, newCols, actualRows> return newType{ newStart }; } T* data() { return start; } };

Now you need to adapt your interface to this new data type, which basically just introduces a few new parameters. Checks remain basically the same.

 // Using this instead of just type-defing allows us to use deducation guides // Replaces: using matrix = std::array from above template<typename T, size_t rows, size_t cols> class matrix { public: std::array<T, rows * cols> storage; auto data() { return storage.data(); } auto data() const { return storage.data(); } }; extern void dgels(char TRANS , integer M, integer N , integer NRHS , double* A, integer LDA , double* B, integer LDB); // Mock, missing a few parameters at the end // Replaces the solve method from above template < typename T, size_t rowsA, size_t colsA, size_t actualRowsA , size_t rowsB, size_t colsB, size_t actualRowsB , TransposeMode mode = TransposeMode::None > void solve(matrix_view<T, rowsA, colsA, actualRowsA> A, matrix_view<T, rowsB, colsB, actualRowsB> B) { static_assert(rowsB >= rowsA && rowsB >= colsA, "b is too small"); char transMode = mode == TransposeMode::None ? 'N' : 'T'; // LDA = rowsA, LDB = rowsB if constexpr (::std::is_same_v<T, float>) { fgels(transMode, rowsA, colsA, colsB, A.data(), actualRowsA, B.data(), actualRowsB); } else if constexpr (::std::is_same_v<T, double>) { dgels(transMode, rowsA, colsA, colsB, A.data(), actualRowsA, B.data(), actualRowsB); // DGELS(, ....); } else { static_assert(always_false_v<T>, "Unknown type"); } }

Usage example:

 int main() { matrix<float, 5, 5> A; matrix<float, 4, 1> b; auto viewA = matrix_view{A}.submat<1, 1, 4, 4>(); auto viewb = matrix_view{b}; solve(viewA, viewb); // solve(viewA, viewb.submat<1, 0, 2, 1>()); // Error: b is too small // solve(matrix_view{A}, viewb.submat<0, 0, 5, 1>()); // Error: can only shrink (b is 4x1 and can not be viewed as 5x1) }
+1
source

Source: https://habr.com/ru/post/1275781/

More articles:


All Articles