Re: matrix library ?

From: lums_at_[hidden]
Date: 2001-03-14 08:59:52

--- In boost_at_y..., walter_at_g... wrote:

> If I understand TNT correctly, it supports some sparse storage 
> schemes. But I think another goal to achieve is also some 
performance 
> enhancement using the right iterator concept. So if you require a 
> matrix library to deal with dense and sparse matrices supporting 
> expression templates too, things seem to get very complicated. I 
> don't know whether any existing library implements all these 
concepts 
> together.

Again at the risk of being somewhat tacky, MTL does in fact support 
sparse and dense in a unified fashion along with unified iterators 
and very careful attention to performance.  The present release does 
not have expression templates.  However, Jeremy has created a high-
performance expression template layer for MTL.  So, there is an 
available library that implements most of these things, and one in 
development that implements them all.  (I am trying to get more 
person-power to devote to this so hopefully things will accelerate 
soon.)

> I have got another problem with sparse matrices also, which I would 
> like to mention here: IMHO dense matrices and classical linear 
> equation solvers together form a pair, for example. Sparse matrices 
> are normally used in the context of iterative solvers, which are 
> quite subtle to handle AFAIK.

Sparse versus dense is solely an issue of underlying representation.  
Both types of representation represent the same thing -- a matrix 
(or, more precisely a finite dimensional linear operator).  It is 
quite an easy matter to write a linear solver in a few lines of code 
for a dense matrix (the classic solver you mention), and one can even 
write a stable solver in a few more lines of code.  Doing LU 
factorization, then back and forward solves are together called 
a 'direct' solution process.

However, one can still perform direct solution on a sparse system of 
equations.  In many cases, such an approach will give you performance 
comparable to, or better than, using an iterative solver.  The 
drawback is that sparse solvers are extraordinarily complicated.  One 
has to represent only the non-zero elements of the matrix and then, 
when factorizing, add fill-in elements where required by the 
factorization operations.  To make this most efficient, one has to 
first order the sparse matrix for low fill (finding the order for 
absolute minimum fill is an NP complete problem).  FWIW, there are a 
few heuristic ordering algorithms in the BGL.

There are some domains where iterative solvers just don't work (or 
haven't been made to work) because of conditioning of the problem -- 
circuit simulation is one notable example.

When you start talking about iterative solvers you also have to start 
talking about preconditioners (incomplete factorizations), a black 
art unto itself.  The subtlety in iterative solvers all comes down to 
preconditioning -- the iterative solvers are not all that 
complicated.  The Iterative Template Library has almost all iterative 
solvers of consequence (and a number of preconditioners).  The 
problem is that no one has ever found a one-size fits all solver (and 
in fact, there is a very famous result by Faber and Manteuffel that 
shows that one solver can't fit all problems) and no one has ever 
found a one-size fits-all preconditioner.  The best preconditioners 
are those that exploit some part of the problem domain (think of the 
preconditioner as an approximate solver).  Usually a good 
preconditioner will be almost as complicated as a direct solver.

A broad over-generalization about performance.  The performance of 
direct vs. iterative methods is quite controversial.  However, for 
certain model problems (discretized elliptic PDEs), the rule of thumb 
is that direct is a win for 1-D problems, it is a wash for 2-D, and 
iterative solver will be a win for 3-D.  Iterative here means an 
optimal solver (CG or GMRES) with an incomplete factorization and 
natrual ordering.

In short.  Sparse matrices are used for efficiency of 
representation.  Direct methods can be used with sparse matrices 
although ordering and factorization are complicated and there will be 
some growth in memory as a result of fill-in.  Iterative solvers are 
easy to write, but subtle and quick to anger -- use the right 
preconditioner.  

For the curious, the link for ITL is 
http://www.lsc.nd.edu/research/itl