Axel Ruhe
September  13, 2004
2D1252, Computational Algebra, Fall 2004, 

Programming Assignment 2

Hand in by September 28 at lecture!

Sparse Matrices, direct algorithms

A short introduction given in lecture notes

The Matlab command
help/ MATLAB/Using Matlab/Mathematics/Sparse matrices
gives an instructive description. Some details are given in  Extra tips.

If you are really interested read the  original paper:
Gilbert, John R., Cleve Moler, and Robert Schreiber, "Sparse Matrices in MATLAB: Design and Implementation," SIAM J. Matrix Anal. Appl., Vol. 13, No. 1, January 1992, pp. 333-356.  gilbert92sparse.ps gilbert92sparse.pdf

Test matrices are found in Matrix Market  at National Institute of Standards and Technology (NIST).

1. Simple example: Compare reorderings!

Start with the matrix you get from a finite difference approximation of the Laplace equation. You get a grid by the command  G=numgrid and a matrix by  A=delsq(G)! Look at the matrix with spy(A).

Take an appropriate size, start with a square. Look at the grid matrix G and see how the points are ordered by columns.

You will get  permutation vectors for band width reduction (RCM) and minimal degree (MMD) by the calls pr=symrcm(A) and pm=symmmd(A). Substructuring, nested dissection, is not implemented, but there is an option in numgrid that gives that order for a square. Look at the reordered matrices with spy  and compare to what I showed in the lecture!

Look at the grid G reordered by Reversed Cuthill Mc Kee. You can use the routine  vector2grid as described in Extra tips. Note that it cannot be used directly on the permutation vector pr, it must have the inverse permutation. It is simple to get the inverse permutation by doing the call
iv=1:n ; rp(pr)=iv ;
which gives the inverse permutation in the vector rp.

Does your result look like what I showed in the lecture?

The Minimum Degree ordering of G does not look as I have described. Take a look at it and compare to Nested Dissection, which is precisely as described. Look also at the matrices with spy, here there is a qualitative similarity between MMD and ND.


2. When is sparse factorization of advantage?

The sparse factorization will use much less storage space and fewer arithmetic operations than a full matrix code. On the other hand it needs quite a lot of book keeping to track at which places fill in occurs. In industrial sparse matrix codes, one switches to handle a full matrix when a certain percentage of the elements are nonzero. Look at the LU factors of a matrix ordered with minimum degree, they are rather dense in the lower right corner!

Now we want to see when it is of advantage to use a sparse matrix code.

Matrix:  Takethe matrix A as finite difference matrices over a square as in previous task.
Right hand side b: Take the vector b as a 1 (one) in one or a set of contigous positions in the center of the grid G.
Solution x: Plotted over the square it will be a tent like shape with poles in those positions where b is one. You can use the routine
X=v2g(x,G);
to get a matrix with the vector x spread over the grid G.

Experiment: Let the number of grid points vary. Make up a table of the number of nonzeros in the original matrix A and the LU factors for the different reorderings RCM, MMD and ND. Measure the time needed for
  1. Analyse: Finding the permutation
  2. Factor:    Compute the LU factors of the reordered matrix
  3. Solve:     Solve a system for a new right hand side b
  4. Check solution: Compute r=Ax-b and list ||r||
Compare to the space and time requirements for a full matrix code. Just call
     AF=full(A); 
and you get a full matrix. Be careful, the full matrix may need too much memory!

A note on timing:
The Matlab clock ticks very slowly. You might need to run the shorter operations several times and divide the total time by the number of repeats.

3. General Sparse Matrices

I havc stored some quite large matrices as Matlab files. Look at them and see how much fill in there is in Gaussian eleimination with different reorderings. Note that these matrices are too large to store as full matrices! Do not print the result of the spy command in the report, it can be very space consuming.

Note that:

 Take one or more of the matrices:

  1. fleigAB.mat A finite element stiffness and mass matrix generated by FEMLAB. Will probably behave like the finite difference matrices of previous task.
  2. Try a 3 dimensional Laplace generated by delsq3d.m i. e. for  20*20*40 points.
  3. A power network for western  USA bcspwr
  4. A VLSI circuit of a clock clock.mat Only one of the matrices, G, has many elements enough to be interesting, I think each element stands for a resistor. A VLSI circuit has a tree like structure with few cycles. This makes the matrices very sparse and algorithms that are built on level structures like RCM will give very little fill in.
  5. You may find another interesting matrix on Matrix Market. Get it and study it!
Good Tour!


Finished September 13, 2004 by Axel Ruhe