I'm dealing with some variable size square sparse matrices resulting from a FEM analysis, and my next step is optimizing the system solving in terms of speed. This is a visualization of some aspects of a typical matrix using MATLAB's spy:

MATLAB visualization

Please note I'm not using MATLAB (or this question would be useless), I just used it for generating the above graph.

The top left graph show my matrix as-is (non zero elements are blue), and bottom left shows the LU decomposition of said matrix. On the top right I have applied the COLAMD algorithm (through MATLAB function colamd) to generate a permutation to my matrix that would result in sparser LU factorized matrix, which indeed happens as the bottom right graph suggests (less than a third of the non zero elements I had before).


Now my problem is understanding how to implement the COLAMD function by hand in the language I'm using. I've been browsing through the COLAMD project page, read the paper on it, and tried to decypher the C source code to no avail.

What I'm looking for is an explanation of the algorithm in the likes of the ones provided by Rosetta Code (example of an algorithm), or some pseudo-code, basically anything that clearly shows what steps must be implemented.

I'm not sure if this belongs in Stack Overflow or here, please advise if it's in the wrong place. Thanks.


1 Answer 1


In the following paper

the authors provided COLAMD preordering algorithm that relies on the same strategy of Matlab's COLMMD preordering algorithm, but with a better heuristic ordering. Their algorithm is faster and produces better orderings, with less nonzeros in the factors of the matrix.

They used COLAMD for square unsymmetric matrices to provide a column preordering for sparse partial pivoting. The aim is to pick a column preordering relies only on the nonzero pattern $A$ such that the factorization stays as sparse as possible.

At step, say $l$, they pick a pivot column $v$ to minimize some metric on all the pivot columns as trial to decrease the fill. Then columns $v$ and $l$ are exchanged.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .