Re: [boost] [rfc] A light weight statically typed matrix template class (again?)

From: Dave Steffen (dgsteffen_at_[hidden])
Date: 2007-12-12 12:08:51

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