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&) ()