C++ linear algebra phone meeting minutes 01 May 2019 Present: Cem Bassoy (CEB) Charles Bay (CB) mbutler Guy Davidson (GD) Matt Harrington (MH) (4154764628) Mark Hoemmen (MFH) Jens Maurer (JM) Bob Steagall (BS) BS: Latest release of paper: https://wg21.link/p1385 Latest version of the code: https://github.com/BobSteagall/wg21/tree/la_r2 BS: Giving a linear algebra talk at C++Now next week. Update on progress on P1385 since our last call in April. Implemented most if not all points in implementation guidance at end of R1 version of P1385. Added several things over the last month: - const and mutable iterator types for vectors - const and mutable begin for vector - row and column view types, for rows resp. columns of a matrix: read-only view type - work in strawman implementation: filled in incomplete member functions just to make everything work; Guy working on impl for arithmetic traits - applied constexpr to almost everything, except for dynamically resizable engines and arith traits (not gotten to it yet) - noexcept'd most of the things - added more debugging helpers and stubs mdspan, executors: still on the to-do list. We wanted basic implementation done first. One major development in how traits work: if interest, I can pop up a screen and show how it works. We've modified operation traits in such a way that it's possible to create new operation traits type, and can override traits that are part of default operation traits. Example in code repo: Overwriting the element promotion traits so that float + float gives a double. Also examples of overriding other things, like arithmetic operator. New traits system makes it easier for user to override just what they want. It doesn't change interface much. Implementation more complex, but this is hidden from user. GD: I'm nearly done with the arithmetic traits; just need to do the complicated ones. Started looking at specializations for geometry, but that's beyond scope of this meeting. BS: Kona feedback: Merge column and row vector types. This is all in std::math namespace. Matrix * vector means matrix * (column vector); vector * matrix means (row vector) * matrix. vector * vector is inner product. Outer product will be named function. JM: Are you just working on implementation, or any progress on updating the paper? The paper does not mention std::math, for example. BS: R2 version of paper will propose std::math. Not part of R1. JM: I would encourage that the R2 update be more verbose in the changelog than what I'm currently seeing in the revision history of the paper. This is where std::math, made compatible with mdspan, etc. changes would go. BS: OK. JM: Changelog would help when rereading paper. Also mention things like zero-based indexing in its part of the paper, not in a remote place. BS: OK, we'll include that in the main text for the R2 paper. BS: You mentioned earlier, JM, that you had some questions. JM: I started to read through the paper for this call, and made a few notes. I would be happy to talk about them with you. BS: Sure. JM: Concern with entire way of approaching this: Non-concept-constrained template parameters. For example: fs_matrix_engine just takes a template class T. We're in a world now where we want to have concepts. Presumably T is an element type that must satisfy the element type constraints. BS: Yes, adding concepts is on the to-do list. We haven't attempted to do it yet. JM: It will change some of the presentation in here. Would be good to have an open issues list in the paper to see that it's on the to-do list, so people won't repeat tha5t concern. GD: You're quite right, JM. In an ideal world, we would have started on this sooner, but concepts are fresh. JM: Fine, you have a lot of concepts in the text of the paper, just not in the C++. BS: A lot of this is just my unfamiliarity with concepts. GD: Me too. BS: For example: matrix element. Could even be noncommutative ring for some matrix arithmetic, but determinant requires more refinement. How do we write C++ concepts that embody these mathematical concepts? I don't know how to do that yet. JM: My view: Ignore the question. We have historically punted on defining what we mean by "usual floating-point type." For example, with std::complex, we only allow float, double, long double. Ditto for random number generation. I understand, you don't want to constrain your library similarly. People have continually failed to describe "a type that acts like double" accurately. For linear algebra, maybe enough that it's a ring. Limited operations. Don't need sine, cosine, log. Maybe possible to describe set of required operations. Again, std::complex, people have tried and failed to make that reasonable. MFH: ... it's a function of the operations you want to do, not of the matrix data structure, so constraining "matrix element type" is not appropriate. JM: Right... BS: ... it will take us a while. JM: The C++ concepts we have mostly embody syntax, not semantics. I need type with the following operators, etc. You're not saying "this is a ring" or "this is a field." BS: Right. A potential fly in the ointment is that of computing the result type. For example, John Mcfarlane has his elastic integer proposal, where the class is parameterized by the number of bits. The result type of addition is determined at compile time by the interaction of the template parameters of the inputs. One thing I don't understand very well is how one could cover that in a C++ concept. JM: You don't necessarily want to ... [Cem Bassoy enters] BS,JM: ... BS: Add two vectors of float16, get vector of float32. JM: Two concept things, must distinguish them. First, two element types, may differ, the two argument types you have, they must have syntactically valid operator+. Second: btw, the result of this operation needs to be this type or convertible to this type. Maybe you should punt on the latter part. The concept requirement for you would be that you can put operator+ between those two types. For the result, you can only hope that the result type makes sense. "Matrix element" in 6.2 -- not sure whether ... essentially an explicit opt-in. If I write a new type, I must explicitly opt in to make it a matrix element type. That's not usually something we do. We don't usually have a special switch to enable this. Any technical reason why you need this? BS: R2 version has changed this. Instead of explicitly specializing a traits class, there's a traits like allocator_traits or pointer_types, in which one describes the capabilities of the numeric type -- it's called number_traits. You specialize it for your new numeric type. MFH: explains history JM: Gneeral nuemric props class, tries ot embody math properties -- I am hesitant to find that a good idea. If you go there, SG6 has an opinion. Second, I undersatnd why you need for adnvanced operations a division, for others you don't. You can check with constexpr... Currently, I am amiss which other props of T which are not checkable by syntactical checks in a concepts... would be needed to help you to implement a matrix library. If there are such properties, I owuld like to have them listed in the paper. BS: OK. That's a good question. I admit that one motivation for number_traits is to constrain set of types of which matrices can be instantiation. JM: Why the additional burden on the user? CEB: Because you don't want to allow e.g., bitwise operations -- you might like to exclude some types, to help back-end developers. JM: The T we're talking about here, element type, is supposed to be airhmetic in shape and form, not std::pair. It's totally OK if I have a back-end developer who decides that his engine is unsuitable for Ts that are not trivially copyable or larger than 8 bytes or not double. That's stuff you can formulate with concepts. [Michael Wong] BS: Without constraint, one could add two matrices of std::string. JM: What's wrong with that? You can add two strings and it does something. CEB: IF you have a matrix, you would expect as a mathematician that you would have a nuemrical type. GD: Not necessarily. You can have matrices of functions. BS: Philosophical point. We want linear algebra to represent math. MFH: What about matrices of polynomials? BS: That's math; strings are not. JM: You could further constrain by requiring plus, minus, and multiply -- that would exclude string -- if you're worried about it. Anyone who defines these operations on a type without actually implementing a ring, they may do so. GD: Yes, they have dug their own grave. BS: OK. Your points on concepts are well taken; I am writing them down. JM: Next: 6.3: talks about element promotion traits. Talks about a config mechanism to find out what the compound element type would be, for input element types. It talks about partial specializations. I would drop partial specialziations -- all impl -- and show what kind of type unification you want here. For plain types T1 and T2, you pick maybe the larger one, and for std::complex and complex you pick complex. That would be actual rules, not partial specs without rules that don't give me clue what's going on . BS. OK. Diff in R2. This was a hallway point I got in Kona. Element promotion begins with decltype of addition operator -- let the types say it -- JM: If that actually gets rid of element propmotion thing ... fine with me. BS: We can't actulaly get rid of it, since it defeats flexibility of the system -- ability to override the dfeault behavior -- e.g., element prom traits were float + flaot -> double. JM: I hope, as written here, it seems the matrix elemeent promotion is a global struct. BS: In R1, it is a global struct. R2 removes that and lets one create partrial specs in your own namespace and insert them into the operand types so they overrid3e the default traints in namespace std. JM: OK, good. My concern here is if we have global setting, two diff parts of program get different outcome when I add two half-floats, then taost. BS: ... [R2 fixes that] ... JM: Yes, that's nice. I believe the only sensible behavior for namespace std is decltype(a+b). BS: Right, that's the default. JM: That is the only sensible thing we would have in namespace std, so we don't have to call it out. ... fine. Not user customizable. Good. JM: Fixed-size engine in 6.4.1. First, there is no showing on the paper of the requirements of a basic engine type. It "provides a basic indexing ..." and then you go on to fixed-size engine.... which part of class declaration is part of basic engine requirements, and which are part of fixed-size engine. For purpose of exposition: would help to have class italics basic_engine italics, that show the required names. JM: 6.4.1: First concern: Two template parameters rows and cols do not seem to be reflected as public static const data members. Template should expose them as public things. BS: In the R2 version, the cols and rows -- all the member functions of fsengine are constexpr -- rows and cols member functions return the template parameter values. JM: that doesn't address the problem -- if you have a run-time fsmatrixengine value, having a constexpr member function just means you can call that function in a constexpr context, but I don't think it will work ... if it's a static member function, that will work, but those functions are nonstatic here, for no good reason actually. That means I need a this pointer, and this pointer is nonconstexpr in general. That means I can't get cols and rows as compile-time constants. BS: OK. The problem goes away if the functions are static? JM: Yes. Also goes away for a different reason if you use mdspan extent. JM: is_dense, is_rectangular as true_type / false_type -- instead of static constexpr bool. What's the benefit of the type as opposed to a simple compile-time boolean? I understand why traits -- must derive from true_type etc. -- but not really necessary here. BS: Reason is precedence set by things in containers ... JM: I'm curious to see precedence having true_type and false_type as member typedef in some larger structure. BS: in allocator_traits: propagate on container copy assignment ... is always equal ... all four of those nested type aliases are assigned to either true_type or false_type. That precdence is what I was attempting to follow, though you're right, having static constexpr bool as variable would be cleaner and easier. GD: [inaudible] JM: I would expressly ask LEWG and LWG whether that is their design guidance in this area, instead of deducing anything from allocator_traits. BS: I was [just] trying to follow an existing model. JM: numeric_limits is closer to your domain and follows a different model. Again, I am not a library expert in this area, just that this strikes me as odd, and we should ask for library guidance. GD: Absolutely; your input is welcome in this regard. JM: when you have single-argument constructors, e.g., drmatrixengine, that takes size utpule ... BS: That should be epxlicit; that's an oversight. JM: Yes. JM: There is a curious amount ... you have this transpose engine here, which claims to be just a view on something else, where operator() exchanges i and j. BS: right. JM: Does that cuase ownership concerns? No, because you don't own the other engine type, you have a pointer to it, not a by-value copy. BS: Correct. JM: Do you get mmeory management concerns here? Can your original matrix go out of scope while you hold the transpose matrix? BS: Yes, it's possible for that to occur, just like a string to which a stirntg_view refers can go out of scope. JM: It feels odd to call this "engine," because I understood engine to own the memory of a matrix. BS: We distinguish between owning and nonowning engines -- transpose engine is nonowning. JM: I get a bit of stomachache ... until you introduced the transpose engine, everything was value based, e.g., return ... by value. Having a view has all the scope and dangling pointers trouble. HOwever, I'm being told that maybe this will change in R2 so I won't press the point. BS: This design is meant to let one create expression template enginees. The enginees would be the expression templates themselves, all of which themslevs would be nonowning view types. Possible to contrive a situaiton where one instantiates a matrix whose engine is an epxxresion, and the matrices to which those expressions refer all go out of scope. The problem as I see it is no worse than string_view and string. The ability to have expression templates and create your own ... is a critical part of the flexigblity we're trying to build in. Expression templates would improve performance ... JM: You eventually need values, but I understand that there are efficiency gains, e.g., if multiplication by a constant is folded into matrix-vector multiplication. Maybe my main thing here is that the 6.4 intro to engine types does not even mention that there are owning and nonowning [engines]. it says, the underarching purpose is to perform "resource management" -- that suggests storage management. Maybe it should be an "engine view" to highlight that it's not a storage management ... BS: OK. We allude to that in Definition 4, of what an engine is. We haven't clearly spelled out that storage management could be owning or nonowning. JM: Frankly, with this wording, I'm not getting that this is where expression templates fit in. If engines are, then this needs to be much more explicit. Expression templates are not just about conveying access to the math elements. The expr templates , in their core meaning, just say, you're performing all computations at the operator=. There's no direct element access at that point, it's just a bulk operation that you do at one point, rather than several operations with temporaries in between. I don't read this in section ... BS: OK, that's fine. This is good feedback. JM: One thing I have not understood in here, is in 6.7.1, we have e.g., negation traits. When I call unary operator- on a matrix, this tells me how to actually compute that, and what the result is, ..., right? BS: Yes. JM: I understand that there's a class operand1, which is the matrix or mathematical vector type. Where does optraits R come from? BS: A better way to look at the problem: Skip to Section 7.1. This explains how the various traits work together. It makes more sense from the top down. JM: You're computing new operation traits from the two arguments ... and those traits are the ones you pass into ... BS: Exactly [explains relation to element type promotion]. JM: OK, the thing I was missing is that the op traits were doing these promotions -- unary negation, it's a noop. BS: It's almost trivial, but there just in case. JM: Fine. JM: What I was wondering is, you don't necessarily have to show this in the paper as code, but I undersatnd there are substantial algorithm gains if you special-case diagonal and band matrices? BS: Yes. There are special cases for upper and lower triangular and also for banded matrices. JM: Those are not yet shown here, in terms of code. BS: that's a question of scope for P1385, whether it should include additional engines. JM: I think P1385 will already be large enough once you're done with all the details. I'm fine with leaving those out. What I'm concerned about, is whether the infrastructure works out when those get added later. BS: We haven't implemented those yet. That's a question of, for a given engine type that represents one of those specialized layouts, that its corresponding arithmetic operation is written correctly. There's no conceptual difference between that and say, the transpose engine. [explains] That transpose engine works with our test code, for example. MFH: Historically, this was possible, and done in different ways. JM: Yes, I'm not worried that it can be done. I just want to know whether the infrastructure works with it. I have the same concern... what's conspicuously absenst here is std::valarray. it was targeted to solve some of these problems, but was a dead end. I don't want a dead end like that, that does not do what it promised to do. Same for the expression templates. Do you have experimental code that does expression templates and shows that the engine mechanics work? BS: I don't have that yet. ... binary addition. If I call m.t and return a view of the transpose, that is, in effect, an expression. JM: That's fine. When do you trigger materializing that expression when you see the operator=? ONe of the unsolved problems with expression templates, is "auto x =" is not a value type, it's just a view. However, "specific_matrix_type x = " would hopefully make a real copy of the matrix. What's the mechanic to make that happen, and is that another operation... that needs to be in the infrastructure -- the materialization or copy operation. BS: That's a good question: I don't have an answer yet. We will think about that. I don't want to speak for GD. GD: We'll just have to keep writing experimental code. CEB: In Boost, we have a rudimentary expression template system. I think we have done it with CRTP. You have a base expression structure where you derive all your types -- matrix and vector become expressions -- as JM said, you need an operator= that evaluates that expression. This is how we did it in Boost uBLAS. I have never seen that like what you have done here, with an engine ... this could be reevaluated ... BS: A long time ago, when I was working in medical imaging, I wrote a medical imaging library for my company. some of my experiences went into this. I was looking at whether to derive matrices and vectors from base expression type or use engines. I tried the engine approach and it seems to work well. Engines avoid conceptual baggage of inheritance. JM: I have no doubt that not using inheritance in this area is beneficial. All I'm saying is that until I've seen it [expression templates etc.] or have credible hearsay that it works, I wouldn't count on it. For example, operator= materialization. there is an operator= that takes a different engine and a different optraits and does something -- maybe that's enough to deal with the expression templates business, but I don't know. CEB: What do you mean by "materialization"? JM: Evaluation of the accumulated instructions -- expression templates work by accumulating mathematical operations into a complicated C++ data structure, and materialization says "go get the values." BS: Your points are well taken, about documenting the desire to have expression templates, and proving that it works. GD: That's good high-quality feedback, JM; thank you very much. BS: Yes, thank you. This is the first paper I've worked on. BS: MFH: blas lapack wrapper paper suggestions? BS: mdarray? BS would like paper proposal BS: strawman mostly not designed for optimal perf, designed to test how pieces of interface work together, to get the compile-time behavior we expect. We had talked about mdspan as common language for element layout, perhaps within and certainly beyond the lin alge librayry ith the low-level stuff. I like it just haven't had time to think about it. could use mdarray accessor for expression templates --- CB: Lurking. I'm learning allt he things about writing a paper ... MW: marks paper is good use of expression tempaltes in library has been borught up do you think it could be problem with library evolution not in the C++ standard BS: I don't know how poeple will react to that. the operators would be value based at this point. Our original thinking was to leave it up to implementer, whether operators work by expression templatse or values Values are easier to program. We could provide impl of expression template engines in std::experimental to get more experience with them MFH: valarray: optional whether it uses expression tempaltes. MW: don't get me wrong, i am supportive of using expression templates, but you're not making it a requirement. BS: We are not. BS: Thank you JM! MW: Yes, thank you JM! next call is beofre paper deadline. Deadline is 17 June-ish. BS: Wed 05 June.