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

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).

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

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.

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!

to get a matrix with the vector x spread over the grid G.

- Analyse: Finding the permutation
- Factor: Compute the LU factors of the reordered matrix
- Solve: Solve a system for a new right hand side b
- Check solution: Compute r=Ax-b and list ||r||

and you get a full matrix. Be careful, the full matrix may need too much memory!

A note on timing:

Note that:

- .m files are executed in Matlab. Just write the file name!
- .mat files are read in by the command load
- .mtx files are in Matrix Market format and is read with the
command A=mmread('file name');

If you do not have the routine mmread, you may fetch it from Matrix Market.

Take one or more of the matrices:

- fleigAB.mat
A finite element stiffness and mass matrix generated by FEMLAB. Will probably
behave like the finite difference matrices of previous task.

- Try a 3 dimensional Laplace generated by delsq3d.m i. e. for 20*20*40 points.
- A power network for western USA bcspwr
- 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. - You may find another interesting matrix on Matrix Market. Get it and study it!

Finished September 13, 2004 by Axel Ruhe