add 32bit support for GreenSolvers

fix omega coversion

cleanup

sanitize_eigen also on vectors; allow non-BlasComplex

test new deterministic dual phi algorithm

tighten version for 32-bit test skip

skip ci cache for the moment

ci fix?

remove ubuntu-latest x86 from ci matrix

exclude x86 altogether (broken on nightly due to new Memory)

ComplexF32 fixes to rest of GreenSolvers

type-stability fixes

.github/workflows/ci.yml

@@ -5,6 +5,9 @@ on:
 concurrency:
   group: ${{ github.workflow }}-${{ github.ref }}
   cancel-in-progress: true
+permissions:
+  actions: write
+  contents: read
 jobs:
   test:
     if: github.event_name == 'push' && contains(toJson(github.event.commits), '[skip test]') == false && contains(toJson(github.event.commits), '[skip tests]') == false
@@ -14,23 +17,24 @@ jobs:
       matrix:
         version:
           - '1'
-          # - 'nightly'
+          - 'nightly'
         os:
           - ubuntu-latest
           - macOS-latest
           - windows-latest
         arch:
           - x64
-          - x86
-        exclude:
-          - os: macOS-latest
-            arch: x86
+          # - x86   ## Currently broken on nightly due to OutOfMemoryError()
+        # exclude:
+        #   - os: macOS-latest
+        #     arch: x86
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
+      - uses: julia-actions/cache@v1
       - uses: julia-actions/julia-buildpkg@latest
       - uses: julia-actions/julia-runtest@latest
       - uses: julia-actions/julia-uploadcodecov@latest

src/bands.jl

@@ -61,8 +61,9 @@ end
 
 bands(rng, rngs...; kw...) = h -> bands(h, rng, rngs...; kw...)
 
-bands(h::Union{Function,AbstractHamiltonian}, rng, rngs...; kw...) =
-    bands(h, mesh(rng, rngs...); kw...)
+bands(h::Function, rng, rngs...; kw...) = bands(h, mesh(rng, rngs...); kw...)
+bands(h::AbstractHamiltonian{T}, rng, rngs::Vararg{Any,L´}; kw...) where {T,L´} =
+    bands(h, convert(Mesh{SVector{L´+1,T}}, mesh(rng, rngs...)); kw...)
 
 bands(h::AbstractHamiltonian{<:Any,<:Any,L}; kw...) where {L} =
     bands(h, default_band_ticks(Val(L))...; kw...)
@@ -660,11 +661,14 @@ end
 function simplex_projector(hkav, verts, vind, εav, mindeg)
     φ = states(verts[vind])
     hproj = φ' * hkav * φ
-    _, P = eigen!(Hermitian(hproj), sortby = ε -> abs(ε - εav))
+    _, P = maybe_eigen!(Hermitian(hproj), sortby = ε -> abs(ε - εav))
     Pthin = view(P, :, 1:mindeg)
     return φ * Pthin
 end
 
+maybe_eigen!(A::AbstractMatrix{<:LinearAlgebra.BlasComplex}; kw...) = eigen!(A; kw...)
+maybe_eigen!(A; kw...) = eigen(A; kw...)    # generic fallback for e.g. Complex16
+
 #endregion
 
 ############################################################################################

src/convert.jl

@@ -15,6 +15,9 @@ Base.convert(::Type{T}, l::CellSites) where T<:CellSites = T(l)
 Base.convert(::Type{T}, l::T) where T<:AbstractHamiltonian = l
 Base.convert(::Type{T}, l::AbstractHamiltonian) where T<:AbstractHamiltonian = T(l)
 
+Base.convert(::Type{T}, l::T) where T<:Mesh = l
+Base.convert(::Type{T}, l::Mesh) where T<:Mesh = T(l)
+
 # Constructors for conversion
 Sublat{T,E}(s::Sublat, name = s.name) where {T<:AbstractFloat,E} =
     Sublat{T,E}([sanitize_SVector(SVector{E,T}, site) for site in sites(s)], name)
@@ -40,6 +43,9 @@ function ParametricHamiltonian{E}(ph::ParametricHamiltonian) where {E}
     return ParametricHamiltonian(hparent, h, ms, ptrs, pars)
 end
 
+Mesh{S}(m::Mesh) where {S} =
+    Mesh(convert(Vector{S}, vertices(m)), neighbors(m), simplices(m))
+
 # We need this to promote different sublats into common dimensionality and type to combine
 # into a lattice
 Base.promote(ss::Sublat{T,E}...) where {T,E} = ss

src/greenfunction.jl

@@ -31,23 +31,30 @@ default_green_solver(::AbstractHamiltonian) = GS.Bands()
 (g::GreenSlice)(; params...) = minimal_callsafe_copy(call!(g; params...))
 (g::GreenSlice)(ω; params...) = copy(call!(g, ω; params...))
 
-call!(g::GreenFunction, ω::Real; params...) = call!(g, retarded_omega(ω, solver(g)); params...)
+function call!(g::GreenFunction{T}, ω::Real; params...) where {T}
+    ω´ = retarded_omega(real_or_complex_typed(T, ω), solver(g))
+    return call!(g, ω´; params...)
+end
 
-function call!(g::GreenFunction, ω::Complex; params...)
+function call!(g::GreenFunction{T}, ω::Complex; params...) where {T}
+    ω´ = real_or_complex_typed(T, ω)
     h = parent(g)
     contacts´ = contacts(g)
     call!(h; params...)
-    Σblocks = call!(contacts´, ω; params...)
+    Σblocks = call!(contacts´, ω´; params...)
     corbs = contactorbitals(contacts´)
-    slicer = solver(g)(ω, Σblocks, corbs)
+    slicer = solver(g)(ω´, Σblocks, corbs)
     return GreenSolution(g, slicer, Σblocks, corbs)
 end
 
 call!(g::GreenSlice; params...) =
     GreenSlice(call!(greenfunction(g); params...), slicerows(g), slicecols(g))
 
-call!(g::GreenSlice, ω; params...) =
-    call!(greenfunction(g), ω; params...)[slicerows(g), slicecols(g)]
+call!(g::GreenSlice{T}, ω; params...) where {T} =
+    call!(greenfunction(g), real_or_complex_typed(T, ω); params...)[slicerows(g), slicecols(g)]
+
+real_or_complex_typed(::Type{T}, ω::Real) where {T<:Real} = convert(T, ω)
+real_or_complex_typed(::Type{T}, ω::Complex) where {T<:Real} = convert(Complex{T}, ω)
 
 retarded_omega(ω::T, s::AppliedGreenSolver) where {T<:Real} =
     ω + im * sqrt(eps(float(T))) * needs_omega_shift(s)
@@ -434,7 +441,7 @@ Base.view(s::TMatrixSlicer, ::Colon, ::Colon) = view(s.gcontacts, :, :)
 
 function Base.getindex(s::TMatrixSlicer, i::CellOrbitals, j::CellOrbitals)
     g0 = s.g0slicer
-    g0ij = g0[i, j]
+    g0ij = ensure_mutable_matrix(g0[i, j])
     tkk´ = s.tmatrix
     isempty(tkk´) && return g0ij
     k = s.contactorbs
@@ -444,6 +451,9 @@ function Base.getindex(s::TMatrixSlicer, i::CellOrbitals, j::CellOrbitals)
     return gij
 end
 
+ensure_mutable_matrix(m::SMatrix) = Matrix(m)
+ensure_mutable_matrix(m::AbstractMatrix) = m
+
 minimal_callsafe_copy(s::TMatrixSlicer) = TMatrixSlicer(minimal_callsafe_copy(s.g0slicer),
     s.tmatrix, s.gcontacts, s.contactorbs)
 

src/sanitizers.jl

@@ -82,17 +82,16 @@ end
 
 #endregion
 
-# ############################################################################################
-# # Block sanitizers
-# #region
+############################################################################################
+# CellIndices sanitizers
+#region
 
-# sanitize_block(S::Type{<:Number}, s, _) = convert(S, first(s))
-# sanitize_block(S::Type{<:SMatrix}, s::SMatrix, size) = sanitize_SMatrix(S, s, size)
-# sanitize_block(::Type{S}, s::Number, size) where {S<:SMatrix} = sanitize_SMatrix(S, S(s*I), size)
-# sanitize_block(::Type{S}, s::UniformScaling, size) where {S<:SMatrix} =
-#     sanitize_SMatrix(S, S(s), size)
+# an inds::Tuple fails some tests because convert(Tuple, ::UnitRange) doesnt exist, but
+# convert(SVector, ::UnitRange) does. Used e.g. in compute_or_retrieve_green @ sparselu.jl
+sanitize_cellindices(inds::Tuple) = SVector(inds)
+sanitize_cellindices(inds) = inds
 
-# #endregion
+#endregion
 
 ############################################################################################
 # Supercell sanitizers
@@ -115,11 +114,11 @@ sanitize_supercell(::Val{L}, v) where {L} =
 # Eigen sanitizers
 #region
 
-sanitize_eigen(ε, Ψ) = Eigen(sorteigs!(ε, Ψ)...)
 sanitize_eigen(ε::AbstractVector, Ψs::AbstractVector{<:AbstractVector}) =
     sanitize_eigen(ε, hcat(Ψs...))
-sanitize_eigen(ε::AbstractVector{<:Real}, Ψ::AbstractMatrix) =
-    sanitize_eigen(complex.(ε), Ψ)
+sanitize_eigen(ε, Ψ) = Eigen(sorteigs!(sanitize_eigen(ε), sanitize_eigen(Ψ))...)
+sanitize_eigen(x::AbstractArray{<:Real}) = complex.(x)
+sanitize_eigen(x::AbstractArray{<:Complex}) = x
 
 function sorteigs!(ϵ::AbstractVector, ψ::AbstractMatrix)
     p = Vector{Int}(undef, length(ϵ))

src/slices.jl

@@ -157,7 +157,7 @@ combine_subcells(c::C, cs::C...) where {C<:CellOrbitals} =
 function combine_subcells(c::C, cs::C...) where {C<:CellOrbitalsGrouped}
     groups´ = merge(orbgroups(c), orbgroups.(cs)...)
     indices´ = union(orbindices(c), orbindices.(cs)...)
-    return CellIndices(cell(c), indices´, OrbitalLikeGrouped(groups´))
+    return CellOrbitalsGrouped(cell(c), indices´, groups´)
 end
 
 #endregion
@@ -299,7 +299,7 @@ sites_to_orbs(c::AnyCellOrbitalsDict, _) = c
 sites_to_orbs(c::AnyCellOrbitals, _) = c
 
 ## convert SiteSlice -> OrbitalSliceGrouped/OrbitalSlice
-Contacts
+
 sites_to_orbs(s::SiteSelector, g) = sites_to_orbs(lattice(g)[s], g)
 sites_to_orbs(kw::NamedTuple, g) = sites_to_orbs(getindex(lattice(g); kw...), g)
 sites_to_orbs(i::Integer, g) = orbslice(selfenergies(contacts(g), i))
@@ -335,13 +335,13 @@ function sites_to_orbs(cs::CellSites, os::OrbitalBlockStructure)
     sites = siteindices(cs)
     groups = _groups(sites, os) # sites, orbranges
     orbinds = _orbinds(sites, groups, os)
-    return CellIndices(cell(cs), orbinds, OrbitalLikeGrouped(Dictionary(groups...)))
+    return CellOrbitalsGrouped(cell(cs), orbinds, Dictionary(groups...))
 end
 
 function sites_to_orbs_flat(cs::CellSites, os::OrbitalBlockStructure)
     sites = siteindices(cs)
     orbinds = _orbinds(sites, os)
-    return CellIndices(cell(cs), orbinds, OrbitalLike())
+    return CellOrbitals(cell(cs), orbinds)
 end
 
 _groups(i::Integer, os) = [i], [flatrange(os, i)]

src/solvers/green/bands.jl

@@ -143,7 +143,7 @@ struct BandSimplex{D,T,S1<:SVector{<:Any,<:Real},S2<:SMatrix{<:Any,D,<:Real},S3<
     ei::S1        # eᵢ::SVector{D´,T} = energy of vertex i
     kij::S2       # kᵢ[j]::SMatrix{D´,D,T,DD´} = coordinate j of momentum for vertex i
     eij::S3       # ϵᵢʲ::SMatrix{D´,D´,T,D´D´} = e_j - e_i
-    dual::S1      # dual::SVector{D´,T}, first hyperdual coefficient
+    dualphi::S1   # dual::SVector{D´,T}, first hyperdual coefficient
     VD::T         # D!V = |det(kᵢʲ - kᵢ⁰)|
     refex::R      # Ref to Series expansions for necessary functions
 end
@@ -178,9 +178,10 @@ function BandSimplex(ei::SVector{D´,T}, kij::SMatrix{D´,D,T}, refex = Ref(Expa
     ei, eij = snap_and_diff(ei)
     k0 = kij[1, :]
     U = kij[SVector{D}(2:D´),:]' .- k0          # edges as columns
-    VD = abs(det(U)) / (2π)^D
-    dual = generate_dual(eij)
-    return BandSimplex(ei, kij, eij, dual, VD, refex)
+    VD = T(abs(det(U)) / (2π)^D)
+    ones = [one(T) for _ in Combinations(D´, 3)]
+    dualphi = generate_dual_phi(eij, ones)
+    return BandSimplex(ei, kij, eij, dualphi, VD, refex)
 end
 
 # make similar ei[i] exactly the same, and compute the pairwise difference
@@ -196,15 +197,32 @@ function snap_and_diff(es)
     return ess, chop(ess' .- ess)
 end
 
-function generate_dual(eij::SMatrix{D´,D´,T}) where {D´,T}
-    dual = rand(SVector{D´,T})
-    iszero(eij) && return dual
-    while !is_valid_dual(dual, eij)
-        dual = rand(SVector{D´,T})
-    end
+# e_j such that e^j_k are all nonzero
+generate_dual_e(::Type{SVector{D´,T}}) where {D´,T} = SVector(ntuple(i -> T(i^2), Val(D´)))
+
+# dϕ such that d_ijk = e^j_k dϕ^j_l - e^j_l dϕ^j_k are all as close to 1 as possible
+function generate_dual_phi(eij::SMatrix{D´,D´,T}, ones) where {D´,T}
+    As = (SVector(ntuple(i -> dual_face(eij, face, i), Val(D´))) for face in Combinations(D´, 3))
+    As´ = Iterators.filter(!iszero, As)
+    isempty(As´) && return generate_dual_e(SVector{D´,T})
+    A = transpose(stack(As´))
+    ones´ = size(A, 1) == length(ones) ? ones : ones[1:size(A,1)]
+    dual = SVector{D´, T}(A \ ones´)
     return dual
 end
 
+function dual_face(emat, (j, k, l), i)
+    if i == j
+        return emat[k,l]
+    elseif i  == k
+        return emat[l,j]
+    elseif i == l
+        return emat[j,k]
+    else
+        return zero(eltype(emat))
+    end
+end
+
 # check whether iszero(eʲₖφʲₗ - φʲₖeʲₗ) for nonzero e's
 function is_valid_dual(phi, es)
     phis = phi' .- phi
@@ -244,13 +262,13 @@ NaN_to_Inf(x::T) where {T} = ifelse(isnan(x), T(Inf), x)
 # g_integrals_local: zero-dn g₀(ω) and gⱼ(ω) with normal or hyperdual numbers for φ
 #region
 
-function g_integrals_local(s::BandSimplex{D,T}, ω, ::Val{N} = Val(0)) where {D,T,N}
+function g_integrals_local(s::BandSimplex{<:Any,<:Any,SVector{D´,T}}, ω, ::Val{N} = Val(0)) where {D´,T,N}
     eⱼ = s.ei
     eₖʲ = s.eij
     g₀, gⱼ = begin
         if N > 0 || is_degenerate(eₖʲ)
-            eⱼ´ = s.dual
-            order = ifelse(N > 0, N, D+1)
+            eⱼ´ = generate_dual_e(SVector{D´,T})
+            order = ifelse(N > 0, N, D´)
             eⱼseries = Series{order}.(eⱼ, eⱼ´)
             g_integrals_local_e(s, ω, eⱼseries)
         else
@@ -314,7 +332,7 @@ end
 
 # imaginary log with branchcut in the lower plane
 logim(x::Complex) = iszero(imag(x)) ? logim(real(x)) : log(-im * x)
-logim(x::Real) =  log(abs(x)) - im * 0.5π * sign(x)
+logim(x::T) where {T<:Real} = log(abs(x)) - im * T(0.5π) * sign(x)
 logim(x, ex) = logim(x)
 
 # required for local degenerate case (expansion of logim(Δ::Series, ex))
@@ -344,7 +362,7 @@ function g_integrals_nonlocal(s::BandSimplex{D,T}, ω, dn, ::Val{N} = Val(0)) wh
         ## This dynamical estimate of the order is not type-stable. Not worth it
         # order = N == 0 ? simplex_degeneracy(ϕₖʲ, eₖʲ) + 1 : N
         # if order > 1
-            ϕⱼ´ = s.dual
+            ϕⱼ´ = s.dualphi
             ϕⱼseries = Series{order}.(ϕⱼ, ϕⱼ´)
             g_integrals_nonlocal_ϕ(s, ω, ϕⱼseries)
         else
@@ -416,7 +434,7 @@ function g_integrals_nonlocal_ϕ(s::BandSimplex{D,T}, ω::Number, ϕⱼ) where {
         Δ0 = chop(first(Δⱼ))
         if iszero(Δ0)
             g₀ = zero(complex(T))
-            gⱼ = SVector(ntuple(Returns(g₀), Val(D)))
+            gⱼ = ntuple(Returns(g₀), Val(D))
         else
             Δ0⁻¹ = inv(Δ0)
             γⱼ = αₖʲγⱼ[1,:]                         # if eₖʲ == 0, then αₖʲ == 1
@@ -426,7 +444,7 @@ function g_integrals_nonlocal_ϕ(s::BandSimplex{D,T}, ω::Number, ϕⱼ) where {
             g₀ = q * sum(scalar.(λⱼ))
             gⱼ = ntuple(Val(D)) do j
                 q * scalar(λⱼ[j+1] + im * sum(λₖʲ[:,j+1] - transpose(λₖʲ)[:,j+1]))
-            end |> SVector
+            end
         end
     else
         αₖʲγⱼeϕⱼ = αₖʲγⱼ .* transpose(eϕⱼ)          # αₖʲγⱼeϕⱼ :: SMatrix{D´,D´}
@@ -438,9 +456,9 @@ function g_integrals_nonlocal_ϕ(s::BandSimplex{D,T}, ω::Number, ϕⱼ) where {
         g₀ = q´ * sum(scalar.(Λⱼ))
         gⱼ = ntuple(Val(D)) do j
             q´ * scalar(Λⱼ[j+1] + im * sum(Λₖʲ[:,j+1] - transpose(Λₖʲ)[:,j+1]))
-        end |> SVector
+        end
     end
-    return g₀, gⱼ
+    return Complex{T}(g₀), SVector{D,Complex{T}}(gⱼ)
 end
 
 # Series of cis(ϕ)
@@ -531,13 +549,13 @@ function J_scalar(t::T, e, Δ, ex) where {T<:Real}
     iszero(e) && return zero(complex(T))
     tΔ = t * Δ
     J = iszero(tΔ) ? logim(Δ) : cis(tΔ) * J_integral(tΔ, t, Δ)
-    return J
+    return Complex{T}(J)
 end
 
 function J_scalar(t::Series{N,T}, e, Δ, ex) where {N,T<:Real}
     iszero(e) && return zero(Series{N,Complex{T}})
     N´ = N - 1
-    C = complex(T)
+    C = Complex{T}
     iszero(Δ) && return Series{N}(C(Inf))
     t₀ = t[0]
     tΔ = t₀ * Δ
@@ -567,11 +585,11 @@ function J_scalar(t::Series{N,T}, e, Δ, ex) where {N,T<:Real}
     return rescale(EJ, t[1] * Δ)
 end
 
-function J_integral(tΔ, t, Δ)
+function J_integral(tΔ, t::T, Δ) where {T}
     J = iszero(imag(tΔ)) ?
-        cosint(abs(tΔ)) - im*sinint(real(tΔ)) - im*0.5π*sign(Δ) :
-        -gamma(0, im*tΔ) - im*0.5π*(sign(real(Δ))+sign(real(tΔ)))
-    return J
+        cosint(abs(tΔ)) - im*sinint(real(tΔ)) - im*T(0.5π)*sign(Δ) :
+        -gamma(0, im*tΔ) - im*T(0.5π)*(sign(real(Δ))+sign(real(tΔ)))
+    return Complex{T}(J)
 end
 
 #endregion
@@ -841,7 +859,7 @@ function muladd_ψPψ⁺!(gmat, α, ψ, ψPdict, pind, orbs)
 end
 
 function muladd_ψPψ⁺!(gmat, α, ψ, (rows, cols)::Tuple)
-    if size(ψ,1) == length(rows) == length(cols)
+    if size(ψ, 1) == length(rows) == length(cols)
         mul!(gmat, ψ, ψ', α, 1)
     else
         ψrows = view_or_copy(ψ, rows, :)