We address persistence and stability of three-dimensional vortex configurations in the discrete nonlinear Schrödinger (NLS) equation and develop a symbolic package based on Wolfram's MATHEMATICA for computations of the Lyapunov--Schmidt reduction method. The Lyapunov--Schmidt reduction method is a theoretical tool which enables us to study continuations and terminations of the discrete vortices for small coupling between lattice nodes as well as the spectral stability of the persistent configurations. The method was developed earlier in the context of the two-dimensional NLS lattice and applied to the on-site and off-site configurations (called the vortex cross and the vortex cell) by using semi-analytical computations. The present treatment develops a full symbolic computational package which takes a desired waveform at the anti-continuum limit of uncoupled sites, performs a required number of Lyapunov--Schmidt reductions and outputs the predictions on whether the configuration persists, for finite coupling, in the three-dimensional lattice and whether it is stable or unstable. It also provides approximations for the eigenvalues of the linearized stability problem. We report a number of applications of the algorithm to important multi-site configurations, such as the simple cube, the double cross, and the diamond. For each three-dimensional configuration, we identify exactly one solution, which is stable for small coupling between lattice nodes.