Further Dirac speedup

[certik]

TODO:

• Faster eigensolver for the benchmark #18
• The S matrix depends on kappa (it's an integral which contains r^(2*alpha)), but it doesn't change between SCF iterations; So we should precompute S(kappa) and factorize it (LU or similar); then in each SCF iteration only compute regular eigenproblem
• Once previous point is done, investigate if Gauss-Lobatto can be used to obtain a diagonal S, which would speedup the application of factorized S, similarly as the Schroedinger solver