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.
Anyway, it may sound like Iâm rejecting your feedback, but Iâm not : - ) . I just mean to show that weâve thought about this issue. I do welcome a discussion about this.
SG14 list run by email@example.com
Older Archives on Google Groups