Axel Ruhe
September 18, 2006
|
2D1252,
Computational Algebra, Fall 2006, |
 |
Programming Assignment 2
Hand in by Ocober 5 at lecture!
Sparse Matrices, direct algorithms
A short introduction given in the lecture notes Chapter 2!
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 an L-membrane of size n=15 . 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
approximate minimal degree (AMD) by the calls pr=symrcm(A) and
pm=symamd(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 v2g 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?
Note: The Minimum Degree
ordering of G does not look as I have described.
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.
Now we want to see when it is of advantage to use a sparse matrix code.
Matrix: Take the matrix A as finite difference matrices
over a L-membrane 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 membrane, it will look like a tent
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 and plot it by
the command mesh(X).
Experiment with sparse code: Let the number of grid points vary.
Make up a table of the number of nonzeros in the original matrix A and
the Cholesky factors for the original ordering and the reorderings RCM
and MMD. Determine the permutations and measure the time
needed for
- Factorization:
Compute the LU factors of the reordered matrix
- Solution:
Solve a system for a new right hand side b
Check the solution: Compute the residual r=Ax-b and list ||r||
Start with a moderate grid, say 15*15 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: Just call
AF=full(A);
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 Out of memory for n just above 4000.
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 have stored some quite large matrices as Matlab files. Look at them
and see how much fill in there is in Gaussian elimination 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:
- .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 bcspwr10.mtx
- 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!
Good Tour!
Finished September 18, 2006 by Axel Ruhe