In this work, a derivation and implementation of the relativistic time-dependent configuration interaction singles (RTDCIS) method is presented. Various observables for krypton and xenon atoms obtained by RTDCIS are compared with experimental data and alternative relativistic calculations. This includes energies of occupied orbitals in the Dirac-Fock ground state, Rydberg state energies, Fano resonances and photoionization cross sections. Diagrammatic many-body perturbation theory, based on the relativistic random phase approximation, is used as a benchmark with excellent agreement between RTDCIS reported at the Tamm-Dancoff level. Results from RTDCIS are computed in the length gauage, where the negative energy states can be omitted with acceptable loss of accuracy. A complex absorbing potential, that is used to remove photoelectrons far from the ion, is implemented as a scalar potential and validated for RTDCIS. The RTDCIS methodology presented here opens for future studies of strong-field processes, such as attosecond transient absorption and high-order harmonic generation, with electron and hole spin dynamics and other relativistic effects described by first principle via the Dirac equation.