Simplifying and unifying the design of user-defined mathematical containers

[VinceRev]

Hello.

I have rewritten the beginning of my proposal avoiding the abusive word "vectorization" (https://groups.google.com/a/isocpp.org/forum/#!topic/std-proposals/pnZyN07uAak) and I tried to improve it a little. Before writing the technical specifications, I would like to have some feedback about the motivation and the small example at the beginning of the proposal. Do you think that a tool that simplifies the design of user-defined mathematical containers could fit in the standard library?

Thanks,

Vincent

[stackm...@hotmail.com]

I would prefer the names vector, static_vector and dynamic_vector., but since vector is already taken... What about introducing a namespace? std::math::vector comes to my mind.

[Jens Maurer]

I'm mostly looking at Figure 1 in your paper.
Here are a few concerns:

- Please show a few lines of code that actually make use of
MyMatrix, highlighting the facilities your <mathematizer>
provides "for free".

- Today's trend is towards SIMD vectorization. The interface
you need for your static_mathematizer<...> seems be called
"access", and does not convey enough information so that an
implementation of static_mathematizer<...> could actually
employ SIMD vector instructions, simply because "access"
could use non-contiguous memory for storage of the elements.

- Standardizing a facility that is actually adverse to
SIMD vectorization, yet tries to supply mathematical vectors
and matrices that ought to be perfectly suited for SIMD, seems
to be a non-starter for me.

- If there is a need in the marketplace for standardization of
2D, 3D, 4D vectors and 2x2, 3x3, 4x4 matrices as well as
dynamically-sized mathematical vectors and matrices, we should
rather invest our energy in providing standard types for these
instead of helping to produce even more user-specific variants.

- On the other hand, if people are looking for something
like the boost.iterator library that helps to create a "complete"
iterator type from just a few basic operations, I think
the choice of "access" as the only basic operation is not
sufficiently performance-friendly.

- Is there a concise summary why std::valarray doesn't work?

- In any case, this is complicated enough that it should go
to boost (or some other public place) first before coming before
the committee as a proposal.

Jens

[Vincent R]

Thanks for your comments!

---------------------------------------
1) I will update a more detailed example later.
But for example, you can type:

MyMatrix<double, 3, 3> a(1); // A matrix of one

MyMatrix<double, 3, 3> b(2); // A matrix of two

MyMatrix<double, 3, 3> c = a+4*b; // The calculation is done thanks to the operators provided by static_mathematizer
---------------------------------------

2) I do not master SIMD optimizations enough, but in a very naive way, I thought that a loop like that (where Size is a constant, and where access is inlined) could be vectorized by a compiler.

template<class Kind, std::size_t Size, template<class, Parameters...> class Crtp, class Type, Kind... Parameters>
struct static_mathematizer
{

    template<class OtherType>
    Crtp<Type, Parameters...>& operator+=(const static_mathematizer<Kind, Size, Crtp, OtherType, Parameters...>& rhs)
    {
        for (std::size_t i = 0; i < Size; ++i) {
            static_cast<Crtp<Type, Parameters...>&>(*this).access(i) += rhs.access(i);
        }
        return static_cast<Crtp<Type, Paramters...>&>(*this);
    }
};

I am wrong ?

---------------------------------------

3) I totally agree (so the current design is only ok if you can do SIMD-ization through CRTP)
---------------------------------------

4) On this point, I do not agree, and in fact that was the starting point of this proposal. If we provide some Vector2D, Vector3D etc... in a standardized way and that I want an order 3 tensor with 4*4*4 components with some tensor_product function, I will have to also reimplement all the operator+,-,*... from scratch (and that is the whole problem of the std::valarray) because I cannot derive from the provided containers safely (absence of virtual destructor/protected destructor) and if I use composition, I have to provide an acessor function to get the underlying array, and the user will have to use this accessor for some operations (like mytensor1.data() += mytensor2.data()) and not for other ones (tensor_product(mytensor1, mytensor2)).

---------------------------------------

5) I will answer later

---------------------------------------

6) See 4)

---------------------------------------

7) You are probably right. After having finished my library (that use a mechanism close to that, but with too many functions for a standardization), I will provide a basic implementation of the classes showed in the proposal.

---------------------------------------

2013/1/1 Jens Maurer <Jens....@gmx.net>

Jens

--

--
Vincent Reverdy
Phd Student @ Laboratory Universe and Theories
Cosmology and General Relativity Group
Observatory of Paris-Meudon, France

[Róbert Dávid]

Hello Vincent,

I got a few thoughts after reading your PDF:

- I don't see a clear goal for the library. Is it getting a bunch of "mathy" functions "for free"? Is it (SIMD) vectorization? Is it both?

- What would be the list of functions that the library provides? Please consider the following questions while filling up the list: What functions? Why those? Why only those? And what do they mean? (For example, what would operator* do? Element-wise multiplication? I certainly don't want that happening for matrices..)

- Standard C++ should/must/have to work on all platforms, so a quite critical question is of implementability. What would this library do for a CPU that does not have any SMID functionality? What would it do on VLIW machine? What would it do on a 8 bit embedded controller without memory unit? What would it do on x86, where there are optional extensions for 128 and 256 bit "machine vectors" (SSEx and AVX) that either present or not? Would the library's implementation change if I turn on the --use-avx-instructions switch at the compiler? Even, 'best' implementation is different for Intel and AMD processors, even if they would have the exact same instruction set - both vendors currently provide a library implementation (called MKL and ACML), those provide implementation to (a subset of) these functions - what are your thoughts on such libraries in the light of this proposal?

- I too don't see why is std::valarray insufficient. The whole problem is written valarray all over it. The statement "it (std::valarray) is not well designed for composition or inheritance. Consequently, it cannot be used as a common base on which one can add new operators or functions." is just nonsense - sorry. Almost all C++ class' are perfectly suited for composition (only 'broken' classes like auto_ptr are not), you just need to provide a façade-like interface for all the operators you want to use (and only for those - think of operator* and matrices again), and provide additional operators/functions based on what the class means (like tensor product), and not fizzle with a .data() accessor function:

class MyDoubleMatrix {
    std::valarray<double> data;
public:
    MyDoubleMatrix& operator+=(const MyDoubleMatrix& right) {
        this->data += right.data;
        return *this;
    }
};

int main() {
    MyDoubleMatrix a,b;
    a += b;
}

One place I can see valarray failing is the support for compile-time-size containers -> need of a std::valarray<double, 3>?

- This whole thing seems to be good for one thing only. Then why not providing a standard library that allows us to write

typedef std::math::vector<double, 4> Vector4d;
typedef std::math::matrix<float, 3, 3> Matrix33s;

typedef std::math::matrix<double, 3, std::math::dynamic_size> Matrix3Xd;

typedef std::math::multi_matrix<int, 2,3,4> IntegerMatrixOf2x3x4Size;

etc? Or even, have those typedefs already in. Missing a tensor and tensor_product is not really good argument, add them to the list, and all your worries are gone - don't write a proposal that is missing something :) I can agree with Jens Maurer, let's solve the problem, not the problem of implementing a solution to the original problem. I would love to have such standard library (so I can stop maintaining my own :) ).

The basic idea is interesting, but in this phase this is more like a highly incomplete wish-list than a proposal, and there is a lot work to be done here.

Regards,
Robert

[VinceRev]

Hello Robert.

Thanks for your comments.

- On the point "I too don't see why is std::valarray insufficient." you provide an example of how forwarding the operation to the valarray. Of course that is the good way to do it (better than .data()), but if there is like ~70 operators (counting value/valarray lhs/rhs mixes) and it takes 5 lines of code per operator you have 350 lines to write just to be able to add the function you want to the container.

- On your last point, I agree 100% with you: it would be far better to have standardized small vector/small matrices/small tensors. So if it is the ultimate goal, we will have to standardize all the multilinear algebra operations (and even for a basic tensor you need a metric to raise and lower you indexes and you have to be able to work with covariant/contravariant/lie derivatives): and that is a very huge task. If that is the problem, what would be the first step?

Regards,
Vincent

[VinceRev]

And on the point about forwarding operations to valarray, I think that you break expression template optimizations by doing that, isn't it?

[Róbert Dávid]

If you want proxy results for (insert operator here) two matrices, you also need to have proxy matrix classes (that wrap the proxy from valarray operators) as well.

To answer the question in the previous mail, if we want a standard linear algebra library, first, we should list what we want from it. A look on the current implementations (boost::ublas, Eigen, etc.) is a must.

Regards, Robert

[Jens Maurer]

On 01/02/2013 12:52 AM, Vincent R wrote:
> Thanks for your comments!
>
> ---------------------------------------
> 1) I will update a more detailed example later.
> But for example, you can type:
> MyMatrix<double, 3, 3> a(1); // A matrix of one
> MyMatrix<double, 3, 3> b(2); // A matrix of two
> MyMatrix<double, 3, 3> c = a+4*b; // The calculation is done thanks to the operators provided by static_mathematizer

Ah. It seems you're discussing the mathematical concept of vector spaces here.

In particular, you seem to be trying to offer an easy way to define
the operations on a vector space, i.e.
- addition (subtraction) of two items of a vector space is element-wise addition (subtraction)
- scalar multiplication with an item of a vector space is element-wise multiplication

Are there any other operations you'd like to support in your "static_mathematizer"?
If not, please name it "vector_space_operations" and people will actually know
what you intend to mean.

I'm always looking for generalities, so are there any other mathematical
structures that would benefit from a treatment similar to what you attempt to
offer for vector spaces?

> ---------------------------------------
> 2) I do not master SIMD optimizations enough, but in a very naive way, I thought that a loop like that (where Size is a constant, and where access is inlined) could be vectorized by a compiler.
>
> template<class Kind, std::size_t Size, template<class, Parameters...> class Crtp, class Type, Kind... Parameters>
> struct static_mathematizer
> {
> template<class OtherType>
> Crtp<Type, Parameters...>& operator+=(const static_mathematizer<Kind, Size, Crtp, OtherType, Parameters...>& rhs)
> {
> for (std::size_t i = 0; i < Size; ++i) {
> static_cast<Crtp<Type, Parameters...>&>(*this).access(i) += rhs.access(i);
> }
> return static_cast<Crtp<Type, Paramters...>&>(*this);
> }
> };
>
> I am wrong ?

I don't know. It's at least two abstraction levels away from the C-style array.
Here's another concern:

- Let's assume we're talking about a 3x3 "float" matrix. The SIMD on the machine can deal
with 4 floats in one operation. Let's assume the derived class actually implements
a 4x3 matrix, with the last column ignored, so that there's a nice match to the SIMD
data width. I don't think we can expect the compiler to "know" that it can use
SIMD instructions right away, because those would modify the fourth column, which would
(in general) be visible to other parts of the code.

Jens