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.