Mismatches in matrix size cause infinite loop
In the reduction of an issue I was facing, I arrived at the following ill-formed piece of code which surprisingly compiles for me:
using M23 = Dune::FieldMatrix<double, 2, 3>;
using M32 = Dune::FieldMatrix<double, 3, 2>;
M23 A = {{1, 2, 3}, {10, 20, 30}};
M32 const B = {{1, 2}, {10, 20}, {100, 200}};
M23 M(B); // should this really compile?
At runtime, I do get a segmentation fault, but the cause is all too curious: A backtrace reveals the following loop:
#0 0x000000000050e7ea in Dune::(anonymous namespace)::DenseMatrixAssignerImplementation<Dune::FieldMatrix<double, 2, 3>, Dune::FieldMatrix<double, 3, 2>, false>::DefaultImplementation<Dune::FieldMatrix<double, 2, 3>, Dune::FieldMatrix<double, 3, 2>, false>::apply(Dune::FieldMatrix<double, 2, 3>&, Dune::FieldMatrix<double, 3, 2> const&) ()
#1 0x000000000050e805 in Dune::(anonymous namespace)::DenseMatrixAssignerImplementation<Dune::FieldMatrix<double, 2, 3>, Dune::FieldMatrix<double, 3, 2>, false>::DefaultImplementation<Dune::FieldMatrix<double, 2, 3>, Dune::FieldMatrix<double, 3, 2>, false>::apply(Dune::FieldMatrix<double, 2, 3>&, Dune::FieldMatrix<double, 3, 2> const&) ()
It seems that the piece of code from densematrix.hh
// static_cast
template< class M, class T >
struct DefaultImplementation< M, T, false >
{
static void apply ( M &m, const T &t )
{
static_assert( (Conversion< const T, const M >::exists), "No template specialization of \
DenseMatrixAssigner found" );
m = static_cast< const M & >( t );
}
};
calls itself over and over again because for some reason, the Conversion trait determines that there is a conversion between Dune::FieldMatrix<double, 2, 3> and Dune::FieldMatrix<double, 3, 2>.
The original issue
Similarly, there seems to be a conversion from Dune::FieldMatrix<double, 3, 3> to Dune::FieldMatrix<double, 3, 2>: I originally wanted to test bounds checking and wrote the intentionally ill-formed piece of code
Dune::FieldMatrix<double, 2, 3> A = {{1, 2, 3}, {10, 20, 30}};
Dune::FieldMatrix<double, 3, 2> const B = {{1, 2}, {10, 20}, {100, 200}};
A.rightmultiply(B);
Since rightmultiply (in contrast to rightmultiplyany) makes sense only if the matrix B is square, with a dimension that matches the number of columns of A, the above example is an incorrect use of this function: The result would be a 2x2 matrix yet rightmultiply() expects the resulting matrix to have the same shape as A (i.e. 2x3).
There are in fact two rightmultiply() functions: One in fmatrix.hh:
FieldMatrix& rightmultiply (const FieldMatrix<K,cols,cols>& M)
This one requires B to be a FieldMatrix of the right size. There’s another function in densematrix.hh
template<typename M2>
MAT& rightmultiply (const DenseMatrix<M2>& M)
that does not have such a constraint, hence I expected the call A.rightmultiply(B) in my example to use this function. That’s not what happens, though. The DenseMatrix implementation is only used if I comment out the FieldMatrix implementation. I found this surprising: The FieldMatrix implementation obviously cannot be used. What happens instead (with the FieldMatrix implementation in place): Neither one is called. Here, too, I get a segmentation fault, for which the backtrace reveals the all too familiar infinite loop (albeit with different sizes):
#0 0x0000000000407074 in Dune::(anonymous namespace)::DenseMatrixAssignerImplementation<Dune::FieldMatrix<double, 3, 3>, Dune::FieldMatrix<double, 3, 2>, false>::DefaultImplementation<Dune::FieldMatrix<double, 3, 3>, Dune::FieldMatrix<double, 3, 2>, false>::apply(Dune::FieldMatrix<double, 3, 3>&, Dune::FieldMatrix<double, 3, 2> const&) ()
#1 0x0000000000407094 in Dune::(anonymous namespace)::DenseMatrixAssignerImplementation<Dune::FieldMatrix<double, 3, 3>, Dune::FieldMatrix<double, 3, 2>, false>::DefaultImplementation<Dune::FieldMatrix<double, 3, 3>, Dune::FieldMatrix<double, 3, 2>, false>::apply(Dune::FieldMatrix<double, 3, 3>&, Dune::FieldMatrix<double, 3, 2> const&) ()