  We introduce a software package called HIFIR for preconditioning sparse,
unsymmetric, ill-conditioned, and potentially singular systems. HIFIR computes
a hybrid incomplete factorization, which combines multilevel incomplete LU
factorization with a truncated, rank-revealing QR factorization on the final
Schur complement. This novel hybridization is based on the new theory of
approximate generalized inverse and $\epsilon$-accuracy. It enables
near-optimal preconditioners for consistent systems and enables flexible GMRES
to solve inconsistent systems when coupled with iterative refinement. In this
paper, we focus on some practical algorithmic and software issues of HIFIR. In
particular, we introduce a new inverse-based rook pivoting into ILU, which
improves the robustness and the overall efficiency for some ill-conditioned
systems by significantly reducing the size of the final Schur complement for
some systems. We also describe the software design of HIFIR in terms of its
efficient data structures for supporting rook pivoting in a multilevel setting,
its template-based generic programming interfaces for mixed-precision real and
complex values in C++, and its user-friendly high-level interfaces in MATLAB
and Python. We demonstrate the effectiveness of HIFIR for ill-conditioned or
singular systems arising from several applications, including the Helmholtz
equation, linear elasticity, stationary incompressible Navier--Stokes
equations, and time-dependent advection-diffusion equation.
