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 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 2

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 3

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.