[Dune-devel] Fieldmatrices of different shapes

[Martin Nolte]

Hi Elias,

except for the infinite loop, this behavior is to be expected. FieldMatrix 
defines a method named rightmultiplyany without "using" the implementation in 
the base class. Therefore, the method on the base class is never called.

To me, it seems you have run into two separate errors:
(a) the method on the DenseMatrix should be accessible on the FieldMatrix
(b) FieldMatrix< 2, 3 > should not cast into FieldMatrix< 3, 3 > (which you 
already wrote a bug report on).

While (a) is trivial to fix, (b) looks like a chicken-egg problem to me and 
might take some time to fix.

Best,

Martin

On 12/07/2015 11:45 PM, Elias Pipping wrote:
> On Mon, 26 Oct 2015 at 17:22 Elias Pipping <elias.pipping at fu-berlin.de
> <mailto: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
>
>
> _______________________________________________
> Dune-devel mailing list
> Dune-devel at dune-project.org
> http://lists.dune-project.org/mailman/listinfo/dune-devel
>

-- 
Dr. Martin Nolte <nolte at mathematik.uni-freiburg.de>

Universität Freiburg                                   phone: +49-761-203-5630
Abteilung für angewandte Mathematik                    fax:   +49-761-203-5632
Hermann-Herder-Straße 10
79104 Freiburg, Germany