  Two numerical algorithms are proposed for computing an interval matrix
containing the matrix gamma function. In 2014, the author presented algorithms
for enclosing all the eigenvalues and basis of invariant subspaces of $A \in
\mathbb{C}^{n \times n}$. As byproducts of these algorithms, we can obtain
interval matrices containing small matrices whose spectrums are included in
that of $A$. In this paper, we interpret the interval matrices containing the
basis and small matrices as a result of verified block diagonalization (VBD),
and establish a new framework for enclosing matrix functions using the VBD. To
achieve enclosure for the gamma function of the small matrices, we derive
computable perturbation bounds. We can apply these bounds if input matrices
satisfy conditions. We incorporate matrix argument reductions (ARs) to force
the input matrices to satisfy the conditions, and develop theories for
accelerating the ARs. The first algorithm uses the VBD based on a numerical
spectral decomposition, and involves only cubic complexity under an assumption.
The second algorithm adopts the VBD based on a numerical Jordan decomposition,
and is applicable even for defective matrices. Numerical results show
efficiency and robustness of the algorithms.
