[Dune-devel] Fieldmatrices of different shapes

[Elias Pipping]

Dear list,

I’m very confused and my impression is that at the bottom of my confusion
lies something that’s not working the way it should. I’d appreciate it if
you followed my journey and helped me to undo the knot in my head :)

To test bounds checking in dune-common I wrote a piece of code that I
expected to compile but not run:

    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. I
get a segmentation fault instead (this is with clang 3.7, I have not tested
any other compiler). I don’t understand why but apparently a cast of A to a
3x3 matrix is attempted, which fails, apparently in an infinite loop. A
backtrace looks like this:

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

(over and over again)

A similar issue can be triggered with a simpler piece of code:

  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);

This, too, gives me a segmentation fault, and the backtrace reveals a
similar infinite loop:

0x0000000000406daa 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> (as well as Dune::FieldMatrix<double, 3, 3>
and Dune::FieldMatrix<double, 3, 2> in my original example). That doesn't
seem correct to me... (is it?)

Now, of course in both cases I started out with an illegal piece of code.
But somehow I expected the compiler to handle this more gracefully than
with an infinite loop and a segfault. Is this a dune bug?


Elias
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.dune-project.org/pipermail/dune-devel/attachments/20151026/f834ecb6/attachment.htm>