On Monday 10 December 2007, Andreas Harnack wrote:
Dave Steffen schrieb:
On Thursday 06 December 2007, Andreas Harnack wrote:
Hi @,
shouldn't there be at least an elementary matrix class available in the standard? I notice there have been some discussions in the past, but there doesn't seem to be anything in the pipeline. So, is there any interest in a (really!) light weight statically typed matrix template class?
No, there shouldn't be. :-)
Well, I'm looking forward to a passionate discussion ;-)
Having been out for a few days with stomach flu, I don't know how passionate I'm gonna be. :-)
The C++ landscape is littered with such things, from simpler than your example to fantastically elaborate. See <http://www.oonumerics.org/oon/> for a sample.
I'd take that as a proof that there is a need for standardization, however the result might look like.
No, I meant precicely the opposite. Many many people have started out doing exactly what you're doing, for exactly the same reasons. No "one true way" has emerged, because it's a *hard* problem. We, as an industry, don't yet have a solution that's acceptable as a possible standard. Don't take my word for it. Go look. Look at Blit++, and it's offshoot tvmet (which I quite like, and sounds like what you want); look at MTL, and Boost's own UBlas, and TNT, and... The closest thing we have to a standard linear algebra library (that I know of) is BLAS + LAPACK, which has been around for decades (and was originally written in Fortran). And even that's not appropriate for everyone's needs. Again: doing this is much harder than it looks at first, and should only be attempted by people who know what they're doing. I'll go out on a limb and state that nobody who hasn't looked at some of the attempts at OO Linear Algebra out there knows what they're doing. :-)
One of the big problems is that everyone wants something slightly
different. Example:
#include <algorithm> #include <functional>
template<typename T, unsigned int M, unsigned int N> class matrix { protected: T m[M][N];
OK, this makes the matrix size a compile-time decision. This is ideal for some cases and anathema for others.
Isn't that argument a bit like asking, why C++ does have arrays at all, instead of just pointers to dynamically allocated memory? The huge advantage I see here is that the size is part of the type and hence subject to compile-time checks. (That doesn't exclude the option to have alternative implementations as well.)
It isn't even remotely like asking that question. Of course the advantage here is compile time checking, also no dynamic memory allocation, and (most importantly) huge performance gains from using template metaprogramming techniques (a la Blitz++ and tvmet). ... except that this doesn't support copy-on-write semantics (which you *can* do if the internal storage is dynamically allocated and stored e.g. in a smart pointer). And this doesn't scale well to large matrices (e.g. thousands of 400x400 matrices) on some machines. And this doesn't allow for compile-time sized matrices. And *all* of these issues can be deal breakers. In fact, they *are* all deal breakers for applications I work on.
T& operator()(unsigned int i, unsigned int j) { return m[i][j]; } T const& operator()(unsigned int i, unsigned int j) const { return m[i][j]; }
OK, is that one-based or zero-based? Some people demand zero-based, others demand one-based.
That's zero based, like any other C++ array or container.
Right, but sometimes you really want 1-based, like when you're also working with Fortran or Matlab code, or with Fortran or Matlab programmers, or even just Mathematics professors who write all their papers counting from 1 like all humans except C programmers. :-)
There are still a few essential but fairly basic things missing, like matrix multiplication and transposing a matrix. Elementary row operations might be useful as well for those, who want to provide complex algorithms, but that's it, basically. There shouldn't be anything in there that's more complicated than that. (Matrices are such a general concept that it can be very tricky, if not impossible, to provide any guaranties for anything above the basic stuff.)
OK, what does operator* mean? Inner product? Outer product? Direct product?
As far as I'm aware there's only one product to matrices; we're not talking about vectors here ;-). Of course, you can use row- or column-matrices to represent vectors, in which case the result (i.e. scalar or matrix) depends on the order of operands.
You're kidding, right? There are many ways to multiply two matrices together. I listed three above. Cross products are a fourth. The correct interface depends *heavily* on the intended uses.
Any interest? Comments most welcome!
I'm not saying your code is wrong, it's not; it's just fine for what you want. But it's not at all fine for what the general matrix-using population wants.
Am I really the only one who needs this?
No! Many people need what you're building (and either build it themselves, or use one of the available packages). But many more people need something *like* what you're proposing, except for this, that, and that other bit. Which is why you don't standardize this.
As an example of a library that *does* try to fill all those needs, look at Boost's very own UBlas. Flexible, has all the options covered... and IMHO rather difficult to use, and also rather slow.
Here, I think, we have a huge misunderstanding. I never tried, don't pretend I could and don't even want to provide a solution that suits anybodies need.
Then it's not suitable for standardization.
I think that's impossible. I totally agree that there are plenty of good libraries out there, some of them even excellent.
Not impossible, just difficult, and nobody in the community has a solution that they feel is good enough for standardization.
But I refuse to believe, that the C++ community consists only of people who do high performance scientific computing. For them is reasonably well cared for. But what about the rest of us?
That's nonsensical. The high performance scientific people are, alas, such a minority in the C++ community that the only attempt at any high-performance stuff in the C++ standard is "valarray", which is effectively a failure, IIRC because the only guy on the standardization committee who really understood the issues involved couldn't continue participating. But what about the rest of you? ...
An an example of a library that probably does exactly what you want, but way way fast, look at 'tvmet'.
Agreed, it does what I want. In fact, any library I looked and does what I want because all I want is really simple stuff. I just don't like to pay the price (i.e. the dependencies it introduces, the complexity that gets involved, the compile time overhead etc.).
And elsethread you wrote
I do plan indeed to add a more template parameters, first and foremost to allow to specify special matrix forms, like diagonal matrices, which have some special properties, that might be interesting to the user. Different storage models might interesting too.
... which means you're going to reinvent tvmet.
My point is, I guess, that general-purpose one-size-fits-all libraries have historically been hard to use and slow; libraries that are easy to use and fast are not applicable to all uses. I personally think this issue is rather fundamental to the "computer linear algebra" problem, and nobody has yet come up with The Good Solution (not for lack of work).
I totally agree, but as I said, that wasn't my intention. And by the way, talking about matrices does by no means imply to talk about linear algebra. It just happens that linear algebra peoples are the ones that talk most about matrices. :-)
That distinction isn't nearly as clear-cut as you think. But that aside, my fundamental statement stands: there's good reason not to standardize such a library (yet). I'll add another strong statement: don't roll your own unless you really know what you're doing. Use something already out there, and even then, test the hell out of your numerical results. (There are IIRC bugs in tvmet that can produce wrong answers. And poorly conditioned matrices can wreck havoc under any circumstances.)
(BTW, I think that UBlas + template typedefs [in whatever form we get them] will grow up to be a big, strong, fast, and very useful library some day. It doesn't fit my needs at all yet, so I'm using something else).
Nicely phrased ;-). What do you use then?
Interesting question. For years, we've been using dynamically-sized matrices with BLAS and LAPACK underneath, mainly because that was the only real option when our code base was born. It's not ideally suited for our application (LAPACK is particularly good for big matrices, but we almost always deal with small ones). We're looking at writing our own library right now. We may use tvmet, but probably not. The only reason we would consider doing our own from scratch is that we have the right kind of people to do so. (PhDs in math, computer science, and physics, with loads of scientific computing experience thrown in.) And even then, we're limiting the project scope severely, and building large testing frameworks around it. :-) -- Dave Steffen - Software Engineer 4 Numerica Corporation (www.numerica.us <http://www.numerica.us/> ) Email: dgsteffen @ numerica.us