1# Get all eigensolutions.
 2function _diagonalize(ham::Hamiltonian)
 3    kpoints = ham.basis.kpoints
 4    results = Vector{Any}(undef, length(kpoints))
 5
 6    for ik in eachindex(kpoints)
 7        Afull = Hermitian(Array(ham[ik]))
 8        E = eigen(Afull)
 9        X = E.vectors
10        λ = E.values
11        results[ik] = (; λ, X)
12    end
13    (; λ=[real.(result.λ) for result in results],
14       X=[result.X for result in results])
15end
16
17function diagonalize(bases; n_bands=false)
18    data = map(bases) do basis
19        ham   = Hamiltonian(basis)
20        if n_bands isa Int
21            diagonalize_all_kblocks(diag_full, ham, n_bands)
22        else
23            _diagonalize(ham)
24        end
25    end
26    (; λ=getfield.(data, :λ), X=getfield.(data, :X))
27end