Recent SIMD changes break tests in dune-fem.

dune-common/dune/common/simd/interface.hh:247:72: error: invalid use of incomplete type ‘struct Dune::Simd::Overloads::IndexType<Dune::FieldMatrix<double, 1, 1>*, void>’
     using Index = typename Overloads::IndexType<std::decay_t<V> >::type;

What would be a proper fix for that?

I'm assuming you're referring to something like the error when compiling dune/fem/operator/linear/test/CMakeFiles/test_istlmatrix.dir/test-matrices.cc.o in e.g. this log. The error isn't exactly the same, but I could not find an occurrence of invalid use of incomplete struct in that log, so I'll go with the one I found.

The complete error I found from the log is below. First thing I noticed is that the DynamicMatrix appears to be instantiated with a pointer-to-FieldMatrix field type, which strikes me as odd. Is that something that dune-fem does, and what is it used for?

I can see this being the the cause, though I'll have to investigate closer. The old interface just treated any unknown type as a scalar, which might have been sufficient to make whatever dune-fem does with pointers work. The new interface has an explicit whitelist for types to treat as scalar types in an attempt to report errors earlier. Naturally, pointers are not in that whitelist, since they aren't arithmetic types and don't emulate them either.

[ 59%] Building CXX object dune/fem/operator/linear/test/CMakeFiles/test_istlmatrix.dir/test-matrices.cc.o
In file included from /duneci/modules/dune-common/dune/common/simd/simd.hh:11:0,
                 from /duneci/modules/dune-common/dune/common/densematrix.hh:19,
                 from /duneci/modules/dune-common/dune/common/fmatrix.hh:15,
                 from /duneci/modules/dune-grid/dune/grid/common/geometry.hh:12,
                 from /duneci/modules/dune-grid/dune/grid/common/grid.hh:1181,
                 from /duneci/modules/dune-grid/dune/grid/yaspgrid.hh:20,
                 from /builds/dune-fem/dune-fem/build-cmake/config.h:621,
                 from /builds/dune-fem/dune-fem/dune/fem/operator/linear/test/test-matrices.cc:2:
/duneci/modules/dune-common/dune/common/simd/interface.hh: In substitution of 'template<class V> using Index = typename Dune::Simd::Overloads::IndexType<typename std::decay<_Tp>::type>::type [with V = Dune::DenseMatVecTraits<Dune::DynamicMatrix<Dune::FieldMatrix<double, 1, 1>*> >::value_type]':
/duneci/modules/dune-common/dune/common/densematrix.hh:206:52:   required from 'class Dune::DenseMatrix<Dune::DynamicMatrix<Dune::FieldMatrix<double, 1, 1>*> >'
/duneci/modules/dune-common/dune/common/dynmatrix.hh:58:9:   required from 'class Dune::DynamicMatrix<Dune::FieldMatrix<double, 1, 1>*>'
/builds/dune-fem/dune-fem/dune/fem/operator/matrix/istlmatrix.hh:599:47:   required from 'class Dune::Fem::ISTLLocalMatrix<Dune::Fem::ISTLMatrixObject<Dune::Fem::LagrangeDiscreteFunctionSpace<Dune::Fem::FunctionSpace<double, double, 2, 1>, Dune::Fem::LeafGridPart<Dune::YaspGrid<2> >, 1>, Dune::Fem::LagrangeDiscreteFunctionSpace<Dune::Fem::FunctionSpace<double, double, 2, 1>, Dune::Fem::LeafGridPart<Dune::YaspGrid<2> >, 1> > >'
/builds/dune-fem/dune-fem/dune/fem/storage/objectstack.hh:72:9:   required from 'Dune::Fem::ObjectStackEntry<ObjectFactoryImp>::~ObjectStackEntry() [with ObjectFactoryImp = Dune::Fem::ISTLMatrixObject<Dune::Fem::LagrangeDiscreteFunctionSpace<Dune::Fem::FunctionSpace<double, double, 2, 1>, Dune::Fem::LeafGridPart<Dune::YaspGrid<2> >, 1>, Dune::Fem::LagrangeDiscreteFunctionSpace<Dune::Fem::FunctionSpace<double, double, 2, 1>, Dune::Fem::LeafGridPart<Dune::YaspGrid<2> >, 1> >]'
/builds/dune-fem/dune-fem/dune/fem/storage/objectstack.hh:161:13:   required from 'Dune::Fem::ObjectStack<ObjectFactoryImp>::~ObjectStack() [with ObjectFactoryImp = Dune::Fem::ISTLMatrixObject<Dune::Fem::LagrangeDiscreteFunctionSpace<Dune::Fem::FunctionSpace<double, double, 2, 1>, Dune::Fem::LeafGridPart<Dune::YaspGrid<2> >, 1>, Dune::Fem::LagrangeDiscreteFunctionSpace<Dune::Fem::FunctionSpace<double, double, 2, 1>, Dune::Fem::LeafGridPart<Dune::YaspGrid<2> >, 1> >]'
/builds/dune-fem/dune-fem/dune/fem/operator/matrix/istlmatrix.hh:619:11:   required from 'Dune::Fem::ISTLLinearOperator<DomainFunction, RangeFunction>::ISTLLinearOperator(const string&, const DomainSpaceType&, const RangeSpaceType&, const Dune::Fem::MatrixParameter&) [with DomainFunction = Dune::Fem::ISTLBlockVectorDiscreteFunction<Dune::Fem::LagrangeDiscreteFunctionSpace<Dune::Fem::FunctionSpace<double, double, 2, 1>, Dune::Fem::LeafGridPart<Dune::YaspGrid<2> >, 1>, Dune::FieldVector<double, 1> >; RangeFunction = Dune::Fem::ISTLBlockVectorDiscreteFunction<Dune::Fem::LagrangeDiscreteFunctionSpace<Dune::Fem::FunctionSpace<double, double, 2, 1>, Dune::Fem::LeafGridPart<Dune::YaspGrid<2> >, 1>, Dune::FieldVector<double, 1> >; std::__cxx11::string = std::__cxx11::basic_string<char>; Dune::Fem::ISTLLinearOperator<DomainFunction, RangeFunction>::DomainSpaceType = Dune::Fem::LagrangeDiscreteFunctionSpace<Dune::Fem::FunctionSpace<double, double, 2, 1>, Dune::Fem::LeafGridPart<Dune::YaspGrid<2> >, 1>; Dune::Fem::ISTLLinearOperator<DomainFunction, RangeFunction>::RangeSpaceType = Dune::Fem::LagrangeDiscreteFunctionSpace<Dune::Fem::FunctionSpace<double, double, 2, 1>, Dune::Fem::LeafGridPart<Dune::YaspGrid<2> >, 1>]'
/builds/dune-fem/dune-fem/dune/fem/operator/linear/test/test-matrices.cc:128:34:   required from here
/duneci/modules/dune-common/dune/common/simd/interface.hh:247:72: error: invalid use of incomplete type 'struct Dune::Simd::Overloads::IndexType<Dune::FieldMatrix<double, 1, 1>*, void>'
     using Index = typename Overloads::IndexType<std::decay_t<V> >::type;
                                                                        ^

changed the description

Oh, after protecting the error message you posted from markdown mangling, it seems to be the one I found.

OK, DynamicMatrix does not like a pointer field type, because it is derived from DenseMatrix. DenseMatrix has a private typedef simd_index_type (size_t in the case of scalars), but that typedef cannot be determined because there is no SIMD::Index for pointers, so instantiation fails.

The pointer field type is due to ISTLLocalMatrix in istlmatrix.hh, line 599:

      Dune::DynamicMatrix< LittleBlockType* > matrices_;

Apparently, DynamicMatrix is used here to store a bunch of pointers to FieldMatrixes stored elsewhere, which (I guess) are themselves logically arranged in a matrix-like way.

I'm not going to fight the whole "should only contain entities that represent physical fields"-cause. Though not using DynamicMatrix for pointer storage should be one way to solve this.

If you're not willing to do that, here is what I'll probably do instead: simd_index_type is used within DenseMatrix only within luDecomposition(), or in policy classes meant for use with luDecomposition(). If your field type is a pointer, you're not going to use luDecomposition() directly or indirectly anyway. It should be possible to remove simd_index_type from the main DenseMatrix class, and have it as a template parameter or similar for the policy objects. Then, it will only be instantiated by anything that uses luDecomposition(), e.g. determinant(), solve() or invert().

I should probably phrase the question explicitly:

@robert.kloefkorn Are you going to replace the DynamicMatrix used for pointer storage in ISTLLocalMatrix by something better suited for that task, or do you want me to change DenseMatrix so it won't croak on pointers?

@joe: Thanks for the investigation. I will fix it in dune-fem but I also like your suggestions to move the simd_index_type to the methods. Some code in OPM also breaks with this changes and I'll keep the issue open for now.

dune-common/dune/common/simd/interface.hh:247:72: error: invalid use of incomplete type ‘struct Dune::Simd::Overloads::IndexType<Opm::DenseAd::Evaluation<__float128, 3>, void>’
     using Index = typename Overloads::IndexType<std::decay_t<V> >::type;

@joe: That is the problem in OPM, when a class called Evaluation is used instead of a POD. Evaluation is supposed to overload all math operations. That code was running before, and I'm sure it was also using some of the methods mentioned above. How should that be fixed?

Some additional information on the issues encountered by the OPM simulators: they break because they implement automatic differentiation using quasi-scalar objects and these objects need to be used as the scalar type for FieldMatrix in some places. Note that using SIMD in the linear algebra code for these AD object does not really make sense because they already use SIMD instructions internally. To solve this issue, Dune's dense linear algebra classes could keep the behavior which they exhibited since the beginning of time, except for scalar types that have been explicitly marked to be SIMD compatible. (i.e., use whitelisting for SIMD?)

That sounds like a good solution to me.

As long as they behave like fancy scalar types (possibly use vectorization internally but don't expose that), it should be enough to specialize Dune::Simd::IsStandard for them.

Though I'll grant that you are probably right, and I am considering treating every unknown type like a standard scalar type. The downside would be delayed error detection, i.e. if something goes wrong it'll be detected only deep inside the SIMD machinery, making error messages even harder to understand. The upside is also delayed error detection, i.e. if the parts of the machinery that will fail with the strange types are never actually instantiated, there will be no error.

Oh, and if I treat unknown types like scalars, Dune::Simd::IsStandard will go away, so you would not be able to specialize it in that case...

why not just change the default of that trait or rename it to Dune::Simd::IsSimd?

Yeah, as a quick fix, I'll default the trait to always true. That should mean any type will be accepted. I'll throw the traits class out completely in a later MR.

(IsSimd would be a misnomer -- IsStandard says whether or not a type should be accepted by the adaption layer indented for types like double, float, std::complex. If a type is a true simd type it should have it's own adaption layer, which will take precedence over the adaption layer for standard types no matter what value IsStandard has for that type.)

mentioned in commit 4e48f7c5

mentioned in merge request !509 (merged)

mentioned in issue #122 (closed)

closed via merge request !509 (merged)

mentioned in commit 0daf88cd

simd_index_type is used within DenseMatrix only within luDecomposition(), or in policy classes meant for use with luDecomposition(). If your field type is a pointer, you're not going to use luDecomposition() directly or indirectly anyway.

On a side note: Actually, I think it would be a useful extension to allow users to compute an luDecomposition once and reuse it for solving.

On a side note: Actually, I think it would be a useful extension to allow users to compute an luDecomposition once and reuse it for solving.

Could you plaese open a feature-request-issue so that we can keep track of this. In this issues your side note will surely get lost...