This paper shows that an appropriate form of Gaussian elimination is applicable to H-matrices. If A is an arbitrary real matrix it is called an M-matrix if \(a_{ij}\leq 0\) for \(i\neq j\) and if all its eigenvalues have non-negative real part. The comparison matrix \({\mathcal M}(A)\equiv m_{ij}\) of any square complex matrix A is defined by \(m_{ii}=| a_{ii}|\), \(m_{ij}=-| a_{ij}|\), \(i\neq j\). Such an A is an H-matrix if \({\mathcal M}(A)\) is an M-matrix. Any square matrix A is column diagonally dominant (cdd) if \(\exists\) scalars \(d_ i>0\) such that \(d_ j| a_{jj}| \geq \sum_{i\neq j}d_ i| a_{ij}|,\) \(1\leq j\leq n\) I. Kuo [ibid. 16, 217-220 (1977; Zbl 0353.15021)] has shown that \(\exists\) a permutation matrix P such that \(PAP^ T\) admits an LU factorization into M-matrices and a numerically stable algorithm for the computation is given by the first and the third author [SIAM J. Algebraic Discrete Methods 7, 368-378 (1986; Zbl 0613.65027)]. The current authors extend this result to H-matrices and also show that the computed factorization is numerically stable. Further, it is shown that the growth factor \(\gamma\) resulting from application of Gaussian elimination with cdd pivoting is bounded above by n, the size of A.