Fill-reducing is key to sparse matrix factorization. This course contains a very good introduction to this topic. This course is much earlier but still useful. And George & Liu’s book on sparse linear systems can be downloaded there. Minimum of fill-in is an NP-complete problem. There are two major heuristic method, minimum degree and nested dissection. AMD and Metis are two popular packages corresponding to these two methods. AMD and its variants are contained in Matlab. Metis can be used in Matlab with this mex interface.

What I concern is a special case, square grid(above). There is a Matlab package MESHND for nested dissection in this simple case. Its benchmark shows it does not perform well for 5-point stencil. I thought there should be some simple and efficient method for this simple symmetric case. I tried to use constrained AMD(CAMD) with some proper constraint to improve the result. I found that eliminating four triangle in Picture A below first will produce a better result. I contacted with Bora Uçar, and he suggested to use more diamonds to partition the grid, like Picture B.The method finds out the borders and constrain them to be ordered later in CAMD. This does further improve the result, though the best number of diamonds is not easy to determine. Borders of diamonds are circles under the distance defined by the connection of 5-point stencil.They separate the domain with the minimal size. Under the condition of 9-point stencil, circles turn out to be squares in Picture C. This is the way MESHND partitions the grid.

However, the size of separators just partially explain the result. This becomes obvious when this method is extended to 3D 7-point stencil case. In 3D, the best shape for the size of separators is truncated octahedron. But those truncated octahedra are not good for fill-reducing. The best method I found is simply extending 2D diamonds to 3D. They become octahedra. Octahedra can not fill the 3D domain like diamonds in 2D (ignoring boundaries). There leaves many tetrahedra. These separators are not optimized in size, but they are much more suitable for fill-reducing. The Matlab code to generate the constraint in this way for CAMD is here. 2D version can be easily reduced from that.

Here is a comparison among different methods on s×s regular 5-point stencil grid.

s | MESHND | AMD | Metis | Proposed Method |

64 | 8.81E+004 | 6.72E+004 | 7.06E+004 | 6.19E+004 |

128 | 4.62E+005 | 3.81E+005 | 3.51E+005 | 3.02E+005 |

256 | 2.32E+006 | 1.97E+006 | 1.57E+006 | 1.45E+006 |

512 | 1.12E+007 | 9.90E+006 | 7.83E+006 | 6.76E+006 |

1024 | 5.27E+007 | 4.75E+007 | 3.75E+007 | 3.11E+007 |

There are 4 diamonds per row like picture B above for s=64,128,256. And this number turns out to be 8 when s=512,1024. As a general method, Metis performs well, but it costs much more time than the others.

And for 3D s×s×s grid, the number of nonzeros by this method is as good as the result of Metis, while Metis is much slower and its memory consumption is much higher. Metis runs out of memory on my Win32 machine when s=128.

s | MESHND | AMD | Metis | Proposed Method |

16 | 3.86E+005 | 2.81E+005 | 2.21E+005 | 2.54E+005 |

24 | 2.32E+006 | 1.87E+006 | 1.52E+006 | 1.49E+006 |

32 | 8.11E+006 | 7.75E+006 | 5.52E+006 | 5.18E+006 |

48 | 4.63E+007 | 5.23E+007 | 2.71E+007 | 2.90E+007 |

64 | 1.57E+008 | 1.84E+008 | 9.55E+007 | 9.89E+007 |

128 | 2.83E+009 | 5.17E+009 | NA | 1.91E+009 |

If you have any better explanation or better method, please let me know.

There is a project applying this fill-reducing.