Improving the Performance of sparse LU Decomposition in GEMPACK

In GEMPACK the computational general equilibrium (CGE) model is solved in its linearized form.In this form, the solution to the CGE model is computed by solving a sparse non symmetric linear system of equations(LSE). In most computer programs, the solution to a LSE is computed either by using matrix decomposition methods (direct), iterative methods or subspace methods. GEMPACK employs a sparse LU decomposition with iterative refinement to solve the LSE. Until GEMPACK 12 these solvers purely relied on algorithms from the Harwell Subroutine Library (HSL) with MA48 as the default solver. While computationally efficient, computing the LU decomposition using these off-the-shelf routines remained the time limiting step for many CGE models. In this paper, we will present a new algorithm for the analysis phase and other improvements to the original HSL routines which greatly improve the solution time for most CGE problems. The preparation stage for the analysis phase involves the reorganisation of the matrix into strong components. By replacing the original modified Hopcroft-Karp algorithm with a tree grafting algorithm to determine a suitable bipartite matching the relative time spent can be reduced by a large factor. The most time consuming part in the analysis phase is the determination of a suitable pivot structure. We developed series of buffered sorted data structures, lower bound estimates and an on-the-fly partitioning of the matrix into sparse and full processed parts which efficiently handles this problem. The new algorithm is able to find a suitable pivot structure with lower number of fill-ins at a significantly reduced cost than the original MA48 solver. For the most common CGE problems we see reductions in solution time by a factor 5 to 10.