RE: sparce matrices

[Mikael Adlers <mikael at mathcore dot com>]

Hi all,
I have to add some remarks to this discussion. First of all I would say that
a matrix of the size 1e8 x 1e8 seem quite big to me. Even if only the
diagonal 
of this matrix is stored it would require about 3 Gbytes of memory. However,
on a large cluster of computers and if the matrix has a lot of structure
(band matrix, block triangular or other nive structure) who knows.

I do not agree that a matrix of the order 10.000 "is on shaky numerical
ground".
A sparse matrix of this size is easily solvable for a sparse solver (even a
dense
solver can handle a matrix of this size). There is no shaky numerical
behaviour
for matrices of this size and if it is sparse it is even possible to sharpen

the usual error estimates that are quite pessimistic in most cases (there is

a growth factor in the error estimate that are 2^(n-1), however numerical
analasists
agree that this grows is very unlikely (but attainable) and the growth
factor is in most 
practical cases around 10). I have worked with a sparse factorization
algorithm, 
the QR factorization (computed the least squares solution to a system) and
my 
test problems was about 30.000 x 10.000 with 130.000 nonzero elements and
the 
factorization took a couple of seconds. I considered my test problems as
small ones.
In practice the density of a sparse matrix decreases as the dimension
increases
(like meshes, the number of neighbours do not increase as the mesh
increase).

The really largescale FEM problems can generate sparse matrices of the order
1e6 - 1e7
and the solution of there system is possible on a supercomputer or cluster
of ordinary pc
using iterative methods with sufficient accuracy. I know of an example with
an 3D grid
if the space shuttle consisting of 16.000.000 grid points, each point
generation at least 
one equation.

It is to simplify to say that a goal of sparse computation is to minimize
the movement
of zeros to and from the memory. The goal is more to structure the
computation to avoid 
operations on zeros elements and to minimize the creation of new nonzero
elements
(happens in all direct methods, i.e. methods based on factorization).

The reason why arithmetic matrix libraries like gsl or LAPACK (I have no
knowledge about VSIPL) 
do not support any sparse matrix data type is the following. To be able to
use the sparsity
of the problem it is important that all intermediate results also are
sparse. For example
when solving a linear system with LU decomposition both L and U has to be
sparse
otherwise all is for nothing. To do this, a graph representing the structure
is
created. Then some graph algorithms are added to analyse the graph and try
to
minimize any destruction of sparsity.

However, in iterative libraries like IML++ only matrix - vector products is
needed and 
therefore it is possible to use it with sparse data structures.

The different storage formats that exists is to maximize the performance of
the
specific operation that you would like to do to the matrix. Solving a
tridiagonal
matrix you would like to have it stored compressed rowwize format and
computing the 
QR factorization it is more useful to have it column wise. 

If you would like to solve a problem that you cannot load into memory you
have
to use iterative methods like IML++ otherwise I would recommend the SuperLU
package
(http://www.nersc.gov/~xiaoye/SuperLU/).


Regards,
/Mikael Adlers


------------------------------------------------------------------ 
 Mikael Adlers, Ph.D.          email: mikael@mathcore.com 
 MathCore AB                   phone: +4613 32 85 07 
 Wallenbergsgata 4             fax:         21 27 01
 SE-583 35 Linköping, Sweden   http://www.mathcore.com



> -----Original Message-----
> From: Edwin Robert Tisdale [mailto:E.Robert.Tisdale@jpl.nasa.gov] 
> Sent: den 3 oktober 2001 20:31
> To: gsl-discuss@sources.redhat.com
> Subject: Re: sparce matrices
> 
> 
> Matthew J. Doller wrote:
> 
> > I am about to begin a project at my school using gsl
> > and my advising professor seems to be very interested in gsl.
> > one question he did have was,
> > "Are there any constructs to store sparces matrices?"
> > It might be beneficial to have such a way to store them
> > especially when dealing with matrices that are 10e8 x 10e8 in size.
> 
> First, let me point out that a matrix -- even a sparse matrix --
> of that size has no physical relevance.
> Even a matrix of size 10,000 by 10,000 is on shaky numerical ground.
> More typically, large sparse matrices are on the order of 1,000 by
> 1,000.
> If they are much larger than that,
> the user is probably doing something very, very wrong.
> 
> There are two distinct considerations in sparse matrix libraries.
> The primary consideration is to save time by avoiding algorithms
> that move zeros in and out of memory.
> The secondary consideration is to save memory
> by storing only the non zero entries of the sparse matrix.
> Modern computer architectures have fast floating-point pipelines
> and large, fast and cheap memory so neither of these considerations
> is worth pursuing unless the matrix is very large and very sparse
> (i.e., less than 10% of the entries are non zero).
> It is especially important to keep this in mind
> when considering triangular and symmetric matrix representations.
> 
> There doesn't appear to be any suitable abstraction for 
> sparse matrices.
> Sparse matrix libraries are all about the specific data 
> representation.
> Take a look at Survey of Sparse Matrix Storage Formats
> 
>http://www.netlib.org/linalg/old_html_templates/subsection2.8.3.1.html
>
>Low level vector and matrix arithmetic libraries like the GSL and
>the Vector, Signal and Image Processing Library (VSIPL)
>
>	http://www.vsipl.org/
>
>don't support sparse matrix data types directly but the VSIPL API
>specifies a gather operation which can be used to collect
>the non zero elements of a sparse matrix into a dense vector
>and a scatter operation which can be used to insert
>the elements of a dense vector into a sparse matrix.
>Otherwise, the only other truly useful operation
>is sparse matrix - dense vector multiplication
>which can be used to solve a system of equations iteratively --
>see IML++ (Iterative Methods Library) v. 1.2a
>
>	http://math.nist.gov/iml++/
>
>for example.