# SG14

Subject: Re: [EXTERNAL] Re: Linear algebra library proposal
From: Hoemmen, Mark (mhoemme_at_[hidden])
Date: 2019-06-18 13:59:32

On 6/18/19, 10:19 AM, "p groarke" <philippe.groarke_at_[hidden]<mailto:philippe.groarke_at_[hidden]>> wrote:
> I guess I shouldâ€™ve prefaced my question with, great work! Certainly a good step forward for C++.

Thanks so much! : - D Donâ€™t worry, weâ€™re ready for critiques!

> 3. Yes ðŸ˜Š Iâ€™d be happy to be that person.

Thank you! Weâ€™ll be happy to accept your feedback.

> The reason I bring up quaternions is, they should be treated quite differently than lets say a vec4 (mdarray<float, 1, 4> I believe). Affine matrices can also be treated differently. So the point isnâ€™t so much about features (though quaternions should be considered a requirement for 3D), but how the APIs might be too general in accepting any rank-1 types, for example.

This point distills to the following:

> There needs to be a way to specify extra mathematical properties that have no bearing on storage.

For example: â€œthis mdarray<float, extents<4>> is actually a quaternionâ€ or â€œthis matrix is actually Hermitian.â€ We had a lot of discussions about the latter. In the BLAS, â€œHermitian matrix typeâ€ conflates a mathematical property (A(j,i) == conj(A(i,j))) with a promise by the algorithm (will only access the upper / lower triangle of the matrix). You canâ€™t represent either of these ideas with mdspan or mdarray. Thatâ€™s intentional; mdspan and mdarray are meant to be low level and easy for compilers and libraries to optimize.

Another example: the BLAS has a â€œTriangular matrix type.â€ â€œTriangularâ€ again combines a mathematical property with a promise by the algorithm. A data structure that ignores the promise and tries to represent a â€œtriangular matrixâ€ would behave like a sparse matrix with a fixed structure. (You can read from any entry, but you can only write to entries in the triangle). The implementation would ignore all this anyway â€“ it could just grab the underlying mdspan and just remember the promise and mathematical properties.

If we wanted to represent â€œHermitian matrixâ€ or â€œTriangular matrixâ€ in our design, without writing a complicated data structure that the algorithms would never use anyway, we would have needed some way to â€œtagâ€ a matrix with mathematical properties. For example:

> Have you considered tag dispatching for algorithm overloads?

We considered this idea. You can see a brief discussion at the end of P1673 (â€œOptions and votesâ€). We actually want to bring it up for a vote.

The â€œtag dispatchingâ€ approach calls for a way for users to add their own tags, since library canâ€™t possibly cover all mathematical properties of interest. Thus, many (perhaps all) of the functions would need to be customization points. That would add a LOT of customization points to the Standard Library. LEWG and LWG probably wonâ€™t accept that without a strong argument. Hence, the vote: If people have a strong argument for customization points, we want to hear it.

We made the decision to exclude tags and custom mathematical properties because our proposal intends to be a C++ BLAS binding. It does what the BLAS does. Itâ€™s a low-level building block, not a â€œMatlab in C++.â€ (Cleve Moler originally wrote Matlab as a higher-level interface to the BLAS, LINPACK, and EISPACK. We think that both low-level and higher-level libraries have their place and can cooperate.)

We left out tensors for the same reason. People writing dense tensor libraries often just write loops around optimized BLAS or BLAS-like functions. Tensors come with a mathematical structure that we didnâ€™t want to try to represent, but people writing tensor libraries could call into our library to get good performance.

> Just as we implement iterator contracts, the same can be enforced here with traits. Food for thought.

There are relatively few kinds of iterators, and they all describe variations on one-dimensional iteration. That makes it easier to write generic algorithms based on the iterator categories. I think itâ€™s a lot harder to write generic algorithms on matrices with arbitrary properties, that conflate mathematical and access assertions. If weâ€™re just providing hooks for usersâ€™ implementations, then weâ€™re writing the interface for a higher-level library.

The more we worked on the design, the more we realized this. User customization of mathematical properties properly belongs to a higher-level library, that dispatches to a lower-level library as needed for performance.