It looks to me to that many of the things one would do within the framework of density functional theory (DFT) ultimately boil down to solving a generalized eigenvalue problem $Av = \lambda Bv$ where $A$ and $B$ are matrices, $\lambda$ is the eigenvalue and $v$ is the corresponding generalized eigenvector.
What are the core numerical routines used in DFT, besides solving generalized eigenvalue problems? What about optimization routines?
DFT is a nonlinear eigenvalue problem. In the community of chemistry and material science, it is almost always solved via the self-consistent-field method (SCF).
If the DFT is solved on the level of GGA or LDA (referring to Jacob's Ladder), then, within each SCF iteration, the most expensive numerical routine is indeed (generalized) eigenvalue problem. In addition to the linear (generalized) eigenvalue problem, the density/orbital mixing also involves cubic scaling numerical routines. To my knowledge, they are still much cheaper than the eigensolver.
However, if you further move upon Jacob's ladder to hybrid-functionals or above, then you will encounter numerical routines of much higher order scaling than the eigensolver. In practice, some approximation and numerical tricks can be used to reduce the computational cost of hybrid-functional to the same level as eigensolver. Hence the discussion is more complicated in this case.
Jumping out of the chemistry and material science community, some mathematicians have tried to solve the nonlinear eigenvalue problem directly from an optimization point of view and not uses SCF+mixing. These methods have only been applied to LDA-level functionals. Hence the core numerical routine varies from method to method.
The answer to your question depends a little on the context you are talking about, most importantly the type of basis function you use. As I will describe especially this aspect crucially influences the mathematical structure of the problem you try to solve.
Being a bit sloppy on the mathematical details (for more see e.g. section 4.6 of my thesis Kohn-Sham DFT and related problems such as Hartree-Fock are concerned with solving the following minimisation problem
$$ \min_{\varphi_i}\left( E[\{\varphi_i\}_i] \right) \quad \text{under the constraint} \quad \int \varphi_i(r) \varphi_j(r) \text{d}r = \delta_{ij} $$
where one seeks to find the minimising set of orthonormal orbitals for the Kohn-Sham functional $E[\{\varphi_i\}_i]$. Its precise form depends, for example, on the type of exchange-correlation functional employed.
Now finding the minimum of this problem means as usual that we need to take the functional derivative of that equation and set it to zero. The resulting Euler-Lagrange equations are the famous Kohn-Sham equations
$$ \hat{F}[\{\varphi_i\}_i] \varphi_i = \epsilon_i \varphi_i $$
involving the orbital-dependent Kohn-Sham operator $\hat{F}[\{\varphi_i\}_i]$. Notice that we are still in the infinite-dimensional regime of functions and operators at this point and notice further that $\hat{F}$ is actually the derivative of the energy functional with respect to the orbitals.
As Yingzhou Li already mentioned there are two ways to tackle this problem. The most common approach is the self-consistent field (SCF) procedure, where our angle of attack is the second equation. Introducing a basis the spectral problem of the operator becomes a generalised matrix eigenvalue problem
$$ F[\{v_i\}_i] v_i = \epsilon_i S v_i. $$
Since $F$ depends on the eigenvectors $\{v_i\}_i$, the procedure is to first guess some $\{v_i\}_i$, build an initial $F$, diagonalise it to get new $\{v_i\}_i$ and repeat until convergence. Because the diagonalisation is done many times (order of 20 to 50) for a single ground-state DFT calculation this is where most time is spent. I should mention that I am glossing over details here. In practice you never follow this simple procedure but additional measures like "mixing", "damping" or various forms of preconditioning are employed, but they are not the time-dependent steps.
The structure of $F$ and $S$ depend on the basis you employ. If your basis is orthogonal (like properly normalised plane-waves) than $S$ is just the identity and can be dropped and you get a standard eigenvalue problem. If you use for example Gaussians, than $S$ needs to be kept. Also for plane waves you typically need lots of them, say like 1 million or more, whereas for Gaussians typically large bases are around 1000 functions. On the other hand the matrices $F$ and $S$ for Gaussians are pretty densely populated, whereas for plane-wave-discretised problems $F$ is very sparse and contains many zeros. (I am focusing on Gaussians versus plane waves here, but the picture is similar for other types of bases).
It is easy to imagine that storing a $10^6$ times $10^6$ matrix in memory is going to be a problem, whereas storing a $1000$ times $1000$ matrix is completely fine. As a result one usually does not use the same diagonalisation procedures for both type of basis sets: For Gaussians one assembles the matrix and uses a dense diagonalisation procedure from LAPACK, whereas for plane-waves one uses iterative methods like Davidson or LOBPCG, where the key step is to be able to apply $F$. That is form matrix-vector products $F x$ for trial eigenvectors $x$ that come up during the iterative procedure.
To summarise: For an SCF-based approach for Gaussians building $F$ is the key step, which primarily involves the computation of the two-electron integrals as the most expensive constituent of $F$. For plane waves, however, the application of $F$ is most crucial. Here the most costly steps is the computation of fast-fourier transforms, which are needed to transform between Fourier space (the space in which the elements of $v$ and $F$ live) and real space (the space in which the potential is usually stored).
The second approach to tackle Kohn-Sham is not to use a self-consistent field procedure, but instead to directly tackle the minimisation problem I mentioned first. For this one typically uses Newton or Quasi-Newton minimisation procedures to directly minimise the energy by varying the input orbitals (hence the name). These algorithms all require the derivative of the energy wrt. the variations, hence again $F$. But similar to iterative diagonalisations $F$ only needs to be applied during the Newton minimsation. Therefore these procedures can be used both for Gaussian as well as plane-wave basis sets and the most costly step is, as before, building $F$ or the application of it.
Regarding what Yingzhou Li said, there is no limit to using direct minimisation LDA I am aware of. In fact in DFTK we can easily do direct minimisation for both LDA and GGA functionals and I do not immediately see a reason why an extension to hybrids would not be possible, but we never tried it.
It should be said, however, that some restrictions do apply to direct minimisation. In the typical implementations it is assumed that only the orbitals need to be optimised and not the occupation numbers. For systems with zero or small gaps (such as metals) direct minimisation is therefore not suitable.
Well, arguably the core numerical routine in density functional theory is the evaluation of the exchange-correlation functional itself, and the evaluation by quadrature of the integrals where the exchange-correlation functional sits.
