[Dune-devel] Fieldmatrices of different shapes

[Elias Pipping]

On Mon, 26 Oct 2015 at 17:22 Elias Pipping <elias.pipping at fu-berlin.de>
wrote:

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

I've now turned this into

  https://gitlab.dune-project.org/core/dune-common/issues/5


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