The discretization scheme results in an algebraic system with
4*n*_{c} + *n*_{b}/2 unknowns, where *n*_{c} is the number of cells in the problems and
*n*_{b} is the number of boundary faces.^{2} This is roughly four times the number
of unknowns that a method like the standard seven-point orthogonal operator
produces, but the new method cannot be expressed in terms of only the cell
center variables without obtaining a dense diffusion matrix. There are no
methods which use only cell center unknowns that are second-order accurate
on skewed meshes. The matrix system for this method is unsymmetric, which
necessitates an unsymmetric solver. There is a maximum of 11 non-zero elements
in any row, due to the cell face equation stencil.

When a Krylov subspace solver is used, a specialized preconditioning system has been developed to speed convergence. The preconditioner consists of a low-order version of the new method itself, which is derived by setting the minor direction terms for each face to zero (see section entitled ``Flux Terms''). This results in a matrix which can be reduced to a system involving the cell center unknowns only. In addition to having four times fewer unknowns, this system has a maximum of seven non-zero elements per row and is symmetric. After solution with a symmetric solver (which tends to be faster than an unsymmetric solver), the face unknowns are obtained via back-substitution.

The Krylov subspace solvers used for the solution of the main system are GMRES and TFQMR, among others. The low-order preconditioning system is solved using the Conjugate Gradients method, which is in turn preconditioned using SSOR. Alternately, an incomplete direct method called the unstructured multi-frontal method may be used. This solver gives very exact answers rapidly for small systems, but ultimately loses with the respect to the Krylov solvers when large time-dependent problems with loose tolerances are solved.