  We present an efficient numerical method for computing Hamiltonian matrix
elements between non-orthogonal Slater determinants, focusing on the most
time-consuming component of the calculation that involves a sparse array. In
the usual case where many matrix elements should be calculated, this
computation can be transformed into a multiplication of dense matrices. It is
demonstrated that the present method based on the matrix-matrix multiplication
attains $\sim$80% of the theoretical peak performance measured on systems
equipped with modern microprocessors, a factor of 5-10 better than the normal
method using indirectly indexed arrays to treat a sparse array. The reason for
such different performances is discussed from the viewpoint of memory access.
