Axel Ruhe September 15, 2006 |
2D1252,
Computational Algebra, Fall 2006, |

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 an L-membrane. This is the one
that you see as logo for Matlab. 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?

Note: The Minimum Degree
ordering of G does not look as I have described.

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.

You can use the routine

to get a matrix with the vector x spread over the grid G and plot it by the command

- Factorization: Compute the LU factors of the reordered matrix
- Solution: Solve a system for a new right hand side b

Start with a moderate grid, say 20*20 points. When your code works increase the grid size until you either run out of memory or the experiment takes more than say 3 minutes. Record for which n that happens, and report which machine you use! It may be different on different workstations or PC:s.

Compare to a full matrix code

and you get a full matrix. Do the same computations as in the previous task. Record the time and space needed, and report for which matrix sizes the full matrix code is faster than the sparse one. Be careful, the full matrix may need too much memory for a much smaller n than a sparse matrix! On my account, I got

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 twodimensional finite
element stiffness
and mass matrix generated by FEMLAB. It will probably behave like the
finite
difference matrices of previous task.

- 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 15, 2006 by Axel Ruhe