Best base type to deal with linear algebra

Question

I'm writing a small and inadequate linear algebra library in C++ for a project (I'm sorry). I'm implementing matrices and operations using double precision numbers. I'm doing right? Should I implement a template class instead? Is there a more precise type around?

Answer 1

I've written a C++ linear algebra library using templates. My thought was that we might want to use complex numbers or extended precision numbers some day. That was maybe seven years ago, and we haven't done it yet. We almost always use doubles as the template type, and we have typedefs to make that easy.

A few times we've gone the other way, using types smaller than a double. For example, we've used float rather than double in a memory-bound application described here. But 99.9 percent of the time we use doubles.

If you do use a template argument, watch out for using an integer type but implicitly requiring a floating point type. For example, say you have a matrix whose entries are all integers and so you use a matrix<int> class. But then you pass that to a linear solver. Now your arithmetic is done using integer division, and your results are wrong. (I've done that!)

Answer 2

I would implement the class/struct using a template. In the beginning, you will most likely be satisfied with just double, but I have found that in every project where I didn't implement matrices as templates, I later regretted it.

Also, it gives you an opportunity to use more interesting element-algebras - interval arithmetic, probability distributions, complex math, fixed-point match, sub-matrices, simple math :-), etc.

Answer 3

There is no other type more precise than long double which also has hardware support. But you are free to make your own types if you feel the need for more precision. They will however be quite slower than the native double type, even with extensive optimization.

Answer 4

Final question answer: Yes, there is, it's called long double and is at least as precise as double. For whether to use templates or not, yes i would use templates. That's a great usecase for them and i think it will make porting to some other scalar number type easier. You can then also just typedef a float and/or a double matrix, depending on the system you are running on and which one works faster/better there.

Answer 5

Don't make any extra work for yourself. If you can get by with double (or long double) go with that.

It sounds like this is just a little project, in which case the template thing will just make work for you.

Another option that hasn't been discussed is using a template to define your element type. That doesn't cause much, if any, extra work, but allows some changes later on.

Answer 6

I'm writing a small and inadequate linear algebra library in C++ for a project (I'm sorry)

OUCH! Be careful, be very very careful... Check out JAMA/TNT -- it's got the NIST stamp-of-approval on it and they've already handled some of the "simpler" linear algebra math e.g. various factoring algorithms. Linear algebra involves lots of tricky issues with numerical precision (e.g. Hilbert matrices) and as much as I like doing my own thing this is one of those areas where you might want to use a good solid foundation that's already been well-tested.

You should be able to use long double with it (not quite sure of that) but the algorithms themselves are probably more critical than the precision of the matrices.