I'm spending some time in learning how to use templates in C++. I never used them before and I'm not always sure what can be or what cannot be achieved in different situation.
As an exercise I'm wrapping some of the Blas and Lapack functions that I use for my activities,
and I'm currently working on the wrapping of ?GELS
(that evaluates the solution of a linear set of equations).
A x + b = 0
?GELS
function (for real values only) exists with two names: SGELS
, for single precision vectors and
DGELS
for double precision.
My idea of interface is a function solve
in this way:
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 user requirements, the problem may be overdetermined or undetermined, that means:
- if it is overdetermined
dim(b) > dim(x)
(the solution is a pseudo-inverse) - if it is undetermined
dim(b) < dim(x)
(the solution is a LSQ minimization) - or the normal case in which
dim(b) = dim(x)
(the solution is the inverse ofA
)
(without considering singular cases).
Since ?GELS
stores the result in the input vector b
, the std::array
shouold
have enough space to accomodate the solution, as described in code comments (max(rows, cols)
).
I want to (compile time) determine wich kind of solution to adopt (it is a paramenter change
in ?GELS
call). I have two functions (I'm simplifying for the sake of the question),
that handle the precision and already know which is the dimension of b
and the 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 */
that are part of an internal wrapper. The user function, detemine the size required
in the b
vector through templates:
#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);
}
In this way it actually works. But I want to go one step further, and I really don't have a clue on how to do it.
If the user tries to call solve
with b
of a size that is too small an extremely difficult-to-read error is raised by the compiler.
I'm trying to insert
a static_assert
that helps the user to understand his error. But any direction that comes in my mind
requires the use of two function with the same signature (it is like a template overloading?) for which
I cannot find a SFINAE strategy (and they actually do not compile at all).
Do you think it is possible to raise a static assertion for the case of wrong b
dimension without changing the user interface at compile time?
I hope the 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 constraint that I put on my exercise, in order to improve my skills. That means, I don't know if it is actually possible to achieve the solution. The type of b
must match the function call, and it is easy if I add another template parameter and I change the user interface, violating my constraint.
Minimal complete and working example
This is a minimal complete and working example. As requested I removed any reference to linear algebra concepts. It is a problem of number. The cases are:
N1 = 2, N2 =2
. SinceN3 = max(N1, N2) = 2
everything worksN1 = 2, N2 =1
. SinceN3 = max(N1, N2) = N1 = 2
everything worksN1 = 1, N2 =2
. SinceN3 = max(N1, N2) = N2 = 2
everything worksN1 = 1, N2 =2
. SinceN3 = N1 = 1 < N2
it correctly raises a compilation error. Iwant to intercept the compilation error with a static assertion that explains the fact that the dimension ofN3
is wrong. As for now the error is difficult to read and understand.
You can view and test it online here
Why don't you try to combine tag dispatch together with some
static_assert
s? Below is one way of achieving what you want to solve, I hope. I mean, all the three correct cases are properly piped to the correctblas
calls, different types and dimension mismatches are handled, and the violation aboutfloat
anddouble
s is also handled, all in a user-friendly way, thanks tostatic_assert
.EDIT. I am not sure about your
C++
version requirement, but below isC++11
friendly.You have to consider why the interface offers this (convoluted) mess of parameters. The author had several things in mind. First of all, you can solve problems of the form
A x + b == 0
andA^T x + b == 0
in one function. Secondly, the givenA
andb
can actually point to memory in matrices larger than the ones needed by alg. This can be seen by theLDA
andLDB
parameters.It is the subaddressing that makes things complicated. If you want a simple but maybe useful enough API, you could chose to ignore that part:
Now, addressing the subaddressing possible with
LDA
andLDB
. I propose that you make that part of your data type, not directly part of the template signature. You want to have your own matrix type that can reference storage in a matrix. Perhaps something like this:Now, you'd need to adapt your interface to this new datatype, that's basically just introducing a few new parameters. The checks stay basically the same.
Example usage:
First some improvements that simplify the design a bit and help redability:
there is no need for
biggest_dim
.std::max
is constexpr sice C++14. You should use it instead.there is no need for
b_array_t
. You can just writestd::array< REAL_T, std::max(N1, N2)>
And now to your problem. One nice way in C++17 is:
Or, as pointed by @max66
Tadaa!! Simple, elegant, nice error message.
The difference between the constexpr if version and just a
static_assert
I.e.:is that with just the
static_assert
the compiler will try to instantiatewrap::internal
even onstatic_assert
fail, polluting the error output. With the constexpr if the call towrap::internal
is not part of the body on condition fail so the error output is clean.(*) The reason I didn't just write
static_asert(false, "error msg)
is because that would make the program ill-formed, no diagnostics required. See constexpr if and static_assertYou can also make the
float
/double
deductible if you want by moving the template argument after the non-deductible ones:So the call becomes: