Numerically stable LDLT-factorization of F-type saddle point matrices

We present a new algorithm that constructs a fill-reducing ordering for a special class of saddle point matrices: the F-matrices. This class contains the matrix occurring after discretization of the Stokes equation on a C-grid. The commonly used approach is to construct a fill-reducing ordering for the whole matrix followed by an adaptation of the ordering such that it becomes feasible. We propose to compute first a fill-reducing ordering for an extension of the definite submatrix. This ordering can be easily extended to an ordering for the whole matrix. In this manner, the construction of the ordering is straightforward and it can be computed efficiently. We show that much of the structure of the matrix is preserved during Gaussian elimination. For an F-matrix, the preserved structure allows us to prove that any feasible ordering obtained in this way is numerically stable. The growth factor of this factorization is much smaller than the one for general indefinite matrices and is bounded by a number that depends linearly on the number of indefinite nodes. The algorithm allows for generalization to saddle point problems that are not of F-type and are nonsymmetric, e.g. the incompressible Navier-Stokes equations (with Coriolis force) on a C-grid. Numerical results for F-matrices show that the algorithm is able to produce a factorization with low fill.