Skip to content

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