Matrix(H, slaters) is slow for multi-configurational problems

[jagot]

Finding the NBodyMatrixElement between two (one-electron) Slater determinants costs on average 5 μs:

BenchmarkTools.Trial:
  memory estimate:  5.67 KiB
  allocs estimate:  92
  --------------
  minimum time:     2.958 μs (0.00% GC)
  median time:      4.366 μs (0.00% GC)
  mean time:        6.285 μs (26.06% GC)
  maximum time:     1.100 ms (98.67% GC)
  --------------
  samples:          10000
  evals/sample:     8

When generating equations of motion, one needs to find the energy expression matrix, the complexity of which is quadratic with the number of Slater determinants. In single-active electron calculations, the number of Slater determinants grows quadratically with ℓ (~2ℓ^2), which means that for ℓmax=30, it takes about 15 seconds to just generate the energy expression matrix (which only has ~6ℓ^2 non-zero matrix elements).

Of course, for such cases it is wiser to implement an optimized version, but it would be useful if the general case could be made faster too. However, that would require making NBodyMatrixElement(a, op, b, overlaps) at least two orders of magnitude faster.

[jagot]

Quick profiling shows that a lot of time is spent in cofactor:

I have the feeling that in the same fashion nonzero_minors can efficiently find which minors are nonzero, by looking at the sparsity pattern, cofactor should be able to efficiently calculate their values, without having to recurse to recursion. An obvious optimization is of course the very common diagonal case.