# Re: [SG14] [EXTERNAL] Re: Linear algebra library proposal

From: Hoemmen, Mark <mhoemme_at_[hidden]>
Date: Tue, 18 Jun 2019 18:59:32 +0000
On 6/18/19, 10:19 AM, "p groarke" <philippe.groarke_at_[hidden].com<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.