A fast method to compute higher-order approximate inverses is constructed for multidiagonal band matrices with strong pivots. It is based on concepts such as the order of an element in a strong-pivot matrix, the order matrix, and the zone of influence in an elimination step. This method may be applied to nonsymmetric and indefinite matrices, where there are strong pivots. Together with the block-preconditioned conjugate gradient method, it may be used for the numerical solution of elliptic partial differential equations and related problems. Through numerical experiments it is shown that the method is more efficient than other methods.