From ba36bed073921c8b18f3786251ce4d3cb54b0d60 Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Fri, 26 Jan 2024 17:30:43 +0100 Subject: [PATCH] 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 fix graphene model fix graphene (redo) --- .github/workflows/ci.yml | 16 +++++--- src/bands.jl | 10 +++-- src/convert.jl | 6 +++ src/greenfunction.jl | 24 ++++++++---- src/presets/hamiltonians.jl | 8 ++-- src/sanitizers.jl | 23 ++++++------ src/slices.jl | 8 ++-- src/solvers/green/bands.jl | 70 +++++++++++++++++++++------------- src/solvers/green/kpm.jl | 4 +- src/solvers/green/schur.jl | 38 +++++++++++-------- src/solvers/green/sparselu.jl | 71 +++++++++++++++++++---------------- src/tools.jl | 2 + src/types.jl | 7 +++- test/test_greenfunction.jl | 32 +++++++++++++++- 14 files changed, 203 insertions(+), 116 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 841fcfc8..280411c2 100644 --- a/.github/workflows/ci.yml +++ b/.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 diff --git a/src/bands.jl b/src/bands.jl index 2461c0ea..7945dbec 100644 --- a/src/bands.jl +++ b/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 ############################################################################################ diff --git a/src/convert.jl b/src/convert.jl index 24b0928f..48a2570a 100644 --- a/src/convert.jl +++ b/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 diff --git a/src/greenfunction.jl b/src/greenfunction.jl index 509b0eb9..47ea78aa 100644 --- a/src/greenfunction.jl +++ b/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) diff --git a/src/presets/hamiltonians.jl b/src/presets/hamiltonians.jl index 2273922d..f1a26a59 100644 --- a/src/presets/hamiltonians.jl +++ b/src/presets/hamiltonians.jl @@ -6,16 +6,16 @@ module HamiltonianPresets using Quantica, LinearAlgebra -function graphene(; a0 = 0.246, range = a0/sqrt(3), t0 = 2.7, β = 3, dim = 2, type = Float64, names = (:A, :B), kw...) +function graphene(; a0 = 0.246, range = neighbors(1), t0 = 2.7, β = 3, dim = 2, type = Float64, names = (:A, :B), kw...) lat = LatticePresets.honeycomb(; a0, dim, type, names) h = hamiltonian(lat, - hopping((r, dr) -> t0 * exp(-β*(sqrt(3) * norm(dr)/a0 - 1)) * I,range = range); kw...) + hopping((r, dr) -> t0 * exp(-β*(sqrt(3) * norm(dr)/a0 - 1)) * I, range = range); kw...) return h end function twisted_bilayer_graphene(; twistindex = 1, twistindices = (twistindex, 1), a0 = 0.246, - interlayerdistance = 1.36a0, rangeintralayer = a0/sqrt(3), rangeinterlayer = 4a0/sqrt(3), + interlayerdistance = 1.36a0, rangeintralayer = neighbors(1), rangeinterlayer = 4a0/sqrt(3), hopintra = 2.70 * I, hopinter = 0.48, modelintra = hopping(hopintra, range = rangeintralayer), type = Float64, names = (:Ab, :Bb, :At, :Bt), kw...) @@ -59,4 +59,4 @@ end # module const HP = HamiltonianPresets -#endregion \ No newline at end of file +#endregion diff --git a/src/sanitizers.jl b/src/sanitizers.jl index a3cb6d79..28b87408 100644 --- a/src/sanitizers.jl +++ b/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(ϵ)) diff --git a/src/slices.jl b/src/slices.jl index 0f964f70..137a0591 100644 --- a/src/slices.jl +++ b/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)] diff --git a/src/solvers/green/bands.jl b/src/solvers/green/bands.jl index 63caf2d5..8310e6a0 100644 --- a/src/solvers/green/bands.jl +++ b/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, :) diff --git a/src/solvers/green/kpm.jl b/src/solvers/green/kpm.jl index f6d75ed7..76af50aa 100644 --- a/src/solvers/green/kpm.jl +++ b/src/solvers/green/kpm.jl @@ -174,10 +174,10 @@ needs_omega_shift(s::AppliedKPMGreenSolver) = false #region ## apply ## -function apply(s::GS.KPM, h::Hamiltonian{<:Any,<:Any,0}, cs::Contacts) +function apply(s::GS.KPM, h::Hamiltonian{T,<:Any,0}, cs::Contacts) where {T} isempty(cs) && argerror("The KPM solver requires at least one contact to be added that defiens where the Green function will be computed. A dummy contact can be created with `attach(nothing; sites...)`.") hmat = h(()) - bandCH = band_ceter_halfwidth(hmat, s.bandrange, s.padfactor) + bandCH = T.(band_ceter_halfwidth(hmat, s.bandrange, s.padfactor)) ket = contact_basis(h, cs) momenta = momentaKPM(hmat, ket, bandCH; order = s.order, kernel = s.kernel) return AppliedKPMGreenSolver(momenta, bandCH) diff --git a/src/solvers/green/schur.jl b/src/solvers/green/schur.jl index 61c8fa0e..66bfe0cf 100644 --- a/src/solvers/green/schur.jl +++ b/src/solvers/green/schur.jl @@ -7,8 +7,8 @@ #region struct SchurWorkspace{C} - GL::Matrix{C} - GR::Matrix{C} + GL::Matrix{ComplexF64} + GR::Matrix{ComplexF64} LG::Matrix{C} RG::Matrix{C} A::Matrix{C} @@ -34,10 +34,10 @@ struct SchurFactorsSolver{T,B} linds::Vector{Int} # orbital indices on left surface rinds::Vector{Int} # orbital indices on right surface sinds::Vector{Int} # orbital indices on the smallest surface (left for l<=r, right for l>r) - L::Matrix{Complex{T}} # l<=r ? PL : PL*H' === hp PR (n × min(l,r)) - R::Matrix{Complex{T}} # l<=r ? PR*H === hm PL : PR (n × min(l,r)) - R´L´::Matrix{Complex{T}} # [R'; -L']. L and R must be dense for iG \ (L,R) - tmp::SchurWorkspace{Complex{T}} + L::Matrix{ComplexF64} # l<=r ? PL : PL*H' === hp PR (n × min(l,r)) + R::Matrix{ComplexF64} # l<=r ? PR*H === hm PL : PR (n × min(l,r)) + R´L´::Matrix{ComplexF64} # [R'; -L']. L and R must be dense for iG \ (L,R) + tmp::SchurWorkspace{Complex{T}} # L, R, R´L´ need 64bit end #region ## Constructors ## @@ -58,8 +58,8 @@ function SchurFactorsSolver(h::Hamiltonian{T,<:Any,1}, shift = one(Complex{T})) end function SchurWorkspace{C}((n, d), l, r) where {C} - GL = Matrix{C}(undef, n, d) - GR = Matrix{C}(undef, n, d) + GL = Matrix{ComplexF64}(undef, n, d) + GR = Matrix{ComplexF64}(undef, n, d) LG = Matrix{C}(undef, d, n) RG = Matrix{C}(undef, d, n) A = Matrix{C}(undef, 2d, 2d) @@ -94,19 +94,19 @@ function left_right_projectors(hm::SparseMatrixCSC, hp::SparseMatrixCSC) linds = stored_cols(hm) rinds = stored_cols(hp) # dense projectors - o = one(eltype(hp)) * I + o = one(ComplexF64) * I allrows = 1:size(hp,1) l_leq_r = length(linds) <= length(rinds) PR = o[allrows, rinds] PL = o[allrows, linds] if l_leq_r sinds = linds - R = Matrix(hm[:, linds]) # R = PR H = hm PL + R = Matrix{ComplexF64}(hm[:, linds]) # R = PR H = hm PL L = PL else sinds = rinds R = PR - L = Matrix(hp[:, rinds]) # L = PL H' = hp PR + L = Matrix{ComplexF64}(hp[:, rinds]) # L = PL H' = hp PR end return linds, rinds, L, R, sinds, l_leq_r end @@ -481,10 +481,10 @@ function Base.getproperty(s::SchurGreenSlicer, f::Symbol) return getfield(s, f) end -# note that g.source is taller than L, R, due to extended sites, but of >= witdth +# note that g.sourceC is taller than L, R, due to extended sites, but of >= witdth # size(L, 2) = size(R, 2) = min(l, r) = d (deflated surface) function extended_ldiv!(gL::Matrix{C}, g::SparseLUSlicer, L) where {C} - Lext = view(g.source, :, axes(L, 2)) + Lext = view(g.source64, :, axes(L, 2)) fill!(Lext, zero(C)) copyto!(Lext, CartesianIndices(L), L, CartesianIndices(L)) copy!(gL, view(ldiv!(g.fact, Lext), axes(L)...)) @@ -492,7 +492,7 @@ function extended_ldiv!(gL::Matrix{C}, g::SparseLUSlicer, L) where {C} end function extended_rdiv!(L´g::Matrix{C}, L, g::SparseLUSlicer) where {C} - Lext = view(g.source, :, axes(L, 2)) + Lext = view(g.source64, :, axes(L, 2)) fill!(Lext, zero(C)) copyto!(Lext, CartesianIndices(L), L, CartesianIndices(L)) copy!(L´g, view(ldiv!(g.fact', Lext), axes(L)...)') @@ -503,8 +503,11 @@ end #region ## API ## -Base.getindex(s::SchurGreenSlicer, i::CellOrbitals, j::CellOrbitals) = - isinf(s.boundary) ? inf_schur_slice(s, i, j) : semi_schur_slice(s, i, j) +function Base.getindex(s::SchurGreenSlicer, i::CellOrbitals, j::CellOrbitals) + G = isinf(s.boundary) ? inf_schur_slice(s, i, j) : semi_schur_slice(s, i, j) + # for type-stability with SVector indices + return maybe_SMatrix(G, orbindices(i), orbindices(j)) +end function inf_schur_slice(s::SchurGreenSlicer, i::CellOrbitals, j::CellOrbitals) rows, cols = orbindices(i), orbindices(j) @@ -570,6 +573,9 @@ function semi_schur_slice(s::SchurGreenSlicer{C}, i, j) where {C} end end +maybe_SMatrix(G::Matrix, rows::SVector{L}, cols::SVector{L´}) where {L,L´} = SMatrix{L,L´}(G) +maybe_SMatrix(G, rows, cols) = G + # TODO: Perhaps too conservative function minimal_callsafe_copy(s::SchurGreenSlicer) s´ = SchurGreenSlicer(s.ω, minimal_callsafe_copy(s.solver)) diff --git a/src/solvers/green/sparselu.jl b/src/solvers/green/sparselu.jl index 2c71db1a..91433656 100644 --- a/src/solvers/green/sparselu.jl +++ b/src/solvers/green/sparselu.jl @@ -7,26 +7,25 @@ struct AppliedSparseLUGreenSolver{C} <:AppliedGreenSolver end mutable struct SparseLUSlicer{C} <:GreenSlicer{C} - fact::SparseArrays.UMFPACK.UmfpackLU{C,Int} # of full system plus extended orbs - nonextrng::UnitRange{Int} # range of non-extended orbital indices - unitcinds::Vector{Vector{Int}} # non-extended fact indices per contact - unitcindsall::Vector{Int} # merged and uniqued unitcinds - source::Matrix{C} # preallocation for ldiv! source @ contacts - unitg::Matrix{C} # lazy storage of a full ldiv! solve of all (nonextrng) sites - function SparseLUSlicer{C}(fact, nonextrng, unitcinds, unitcindsall, source) where {C} + fact::SparseArrays.UMFPACK.UmfpackLU{ComplexF64,Int} # of full system plus extended orbs + nonextrng::UnitRange{Int} # range of non-extended orbital indices + unitcinds::Vector{Vector{Int}} # non-extended fact indices per contact + unitcindsall::Vector{Int} # merged and uniqued unitcinds + source64::Matrix{ComplexF64} # preallocation for ldiv! source @ contacts + sourceC::Matrix{C} # alias of source64 or C conversion + unitg::Matrix{C} # lazy storage of a full ldiv! solve of all (nonextrng) sites + function SparseLUSlicer{C}(fact, nonextrng, unitcinds, unitcindsall, source64) where {C} s = new() - s.fact = fact + s.fact = fact # Note that there is no SparseArrays.UMFPACK.UmfpackLU{ComplexF32} s.nonextrng = nonextrng s.unitcinds = unitcinds s.unitcindsall = unitcindsall - s.source = source + s.source64 = source64 + s.sourceC = convert(Matrix{C}, source64) return s end end -SparseLUSlicer(fact::Factorization{C}, nonextrng, unitcinds, unitcindsall, source) where {C} = - SparseLUSlicer{C}(fact, nonextrng, unitcinds, unitcindsall, source) - #region ## API ## inverse_green(s::AppliedSparseLUGreenSolver) = s.invgreen @@ -40,8 +39,8 @@ unitcellinds_contacts_merged(s::SparseLUSlicer) = s.unitcindsall minimal_callsafe_copy(s::AppliedSparseLUGreenSolver) = AppliedSparseLUGreenSolver(minimal_callsafe_copy(s.invgreen)) -minimal_callsafe_copy(s::SparseLUSlicer) = - SparseLUSlicer(s.fact, s.nonextrng, s.unitcinds, s.unitcindsall, copy(s.source)) +minimal_callsafe_copy(s::SparseLUSlicer{C}) where {C} = + SparseLUSlicer{C}(s.fact, s.nonextrng, s.unitcinds, s.unitcindsall, copy(s.source64)) #endregion @@ -65,7 +64,7 @@ function (s::AppliedSparseLUGreenSolver{C})(ω, Σblocks, contactorbitals) where nonextrng = orbrange(invgreen) unitcinds = invgreen.unitcinds unitcindsall = invgreen.unitcindsall - source = s.invgreen.source + source64 = convert(Matrix{ComplexF64}, s.invgreen.source) # the H0 and Σs inside invgreen have already been updated by the parent call!(g, ω; ...) update!(invgreen, ω) igmat = matrix(invgreen) @@ -76,7 +75,7 @@ function (s::AppliedSparseLUGreenSolver{C})(ω, Σblocks, contactorbitals) where argerror("Encountered a singular G⁻¹(ω) at ω = $ω, cannot factorize") end - so = SparseLUSlicer{C}(fact, nonextrng, unitcinds, unitcindsall, source) + so = SparseLUSlicer{C}(fact, nonextrng, unitcinds, unitcindsall, source64) return so end @@ -91,42 +90,48 @@ end function Base.view(s::SparseLUSlicer, i::Integer, j::Integer) dstinds = unitcellinds_contacts(s, i) srcinds = unitcellinds_contacts(s, j) - source = view(s.source, :, 1:length(srcinds)) - return compute_or_retrieve_green(s, dstinds, srcinds, source) + source64 = view(s.source64, :, 1:length(srcinds)) + sourceC = view(s.sourceC, :, 1:length(srcinds)) + return compute_or_retrieve_green(s, dstinds, srcinds, source64, sourceC) end Base.view(s::SparseLUSlicer, ::Colon, ::Colon) = - compute_or_retrieve_green(s, s.unitcindsall, s.unitcindsall, s.source) + compute_or_retrieve_green(s, s.unitcindsall, s.unitcindsall, s.source64, s.sourceC) -function Base.view(s::SparseLUSlicer, i::CellOrbitals, j::CellOrbitals) - # cannot use s.source, because it has only ncols = number of orbitals in contacts - source = similar_source(s, j) - v = compute_or_retrieve_green(s, orbindices(i), orbindices(j), source) +function Base.view(s::SparseLUSlicer{C}, i::CellOrbitals, j::CellOrbitals) where {C} + # cannot use s.source64, because it has only ncols = number of orbitals in contacts + source64 = similar_source64(s, j) + sourceC = convert(Matrix{C}, source64) # will alias if C == ComplexF64 + v = compute_or_retrieve_green(s, orbindices(i), orbindices(j), source64, sourceC) return v end # Implements cache for full ldiv! solve (unitg) -function compute_or_retrieve_green(s::SparseLUSlicer{C}, dstinds, srcinds, source) where {C} +# source64 and sourceC are of size (total_orbs_including_extended, length(srcinds)) +function compute_or_retrieve_green(s::SparseLUSlicer{C}, dstinds, srcinds, source64, sourceC) where {C} if isdefined(s, :unitg) g = view(s.unitg, dstinds, srcinds) else fact = s.fact - allinds = 1:size(fact, 1) - one!(source, srcinds) - gext = ldiv!(fact, source) + allinds = 1:size(fact, 1) # axes not defined on SparseArrays.UMFPACK.UmfpackLU + one!(source64, srcinds) + gext = ldiv!(fact, source64) + sourceC === gext || copy!(sourceC, gext) # required when C != ComplexF64 + # allinds may include extended orbs -> exclude them from the view dstinds´ = ifelse(dstinds === allinds, s.nonextrng, dstinds) - g = view(gext, dstinds´, :) + srcinds´ = convert(typeof(srcinds), 1:length(srcinds)) + g = view(sourceC, dstinds´, srcinds´) if srcinds === allinds - s.unitg = copy(view(gext, s.nonextrng, s.nonextrng)) + s.unitg = copy(view(gext, s.nonextrng, s.nonextrng)) # exclude extended orbs end end return g end -similar_source(s::SparseLUSlicer, ::CellOrbitals{<:Any,Colon}) = - similar(s.source, size(s.source, 1), maximum(s.nonextrng)) -similar_source(s::SparseLUSlicer, j::CellOrbitals) = - similar(s.source, size(s.source, 1), norbitals(j)) +similar_source64(s::SparseLUSlicer, ::CellOrbitals{<:Any,Colon}) = + similar(s.source64, size(s.source64, 1), maximum(s.nonextrng)) +similar_source64(s::SparseLUSlicer, j::CellOrbitals) = + similar(s.source64, size(s.source64, 1), norbitals(j)) # getindex must return a Matrix Base.getindex(s::SparseLUSlicer, i::CellOrbitals, j::CellOrbitals) = copy(view(s, i, j)) diff --git a/src/tools.jl b/src/tools.jl index bcdc6877..f243eec1 100644 --- a/src/tools.jl +++ b/src/tools.jl @@ -15,6 +15,8 @@ rdr((r1, r2)::Pair) = (0.5 * (r1 + r2), r2 - r1) @inline tupleflatten(x0::Tuple, x1, xs...) = tupleflatten((x0..., x1), xs...) @inline tupleflatten(x0::Tuple, x1::Tuple, xs...) = tupleflatten((x0..., x1...), xs...) +tuplefill(x, ::Val{N}) where {N} = ntuple(Returns(x), Val(N)) + padtuple(t, x, N) = ntuple(i -> i <= length(t) ? t[i] : x, N) @inline tupletake(x, ::Val{N}) where {N} = ntuple(i -> x[i], Val(N)) diff --git a/src/types.jl b/src/types.jl index 27c0194b..cfb743fb 100644 --- a/src/types.jl +++ b/src/types.jl @@ -316,15 +316,18 @@ const AnyOrbitalSlice = Union{OrbitalSlice,OrbitalSliceGrouped} #region ## Constructors ## CellSite(cell, ind::Int) = CellIndices(sanitize_SVector(Int, cell), ind, SiteLike()) -CellSites(cell, inds = Int[]) = CellIndices(sanitize_SVector(Int, cell), inds, SiteLike()) +CellSites(cell, inds = Int[]) = CellIndices(sanitize_SVector(Int, cell), sanitize_cellindices(inds), SiteLike()) # exported lowercase constructor for general inds cellsites(cell, inds) = CellSites(cell, inds) CellOrbitals(cell, inds = Int[]) = - CellIndices(sanitize_SVector(Int, cell), inds, OrbitalLike()) + CellIndices(sanitize_SVector(Int, cell), sanitize_cellindices(inds), OrbitalLike()) CellOrbital(cell, ind::Int) = CellIndices(sanitize_SVector(Int, cell), ind, OrbitalLike()) +CellOrbitalsGrouped(cell, inds, groups::Dictionary) = + CellIndices(sanitize_SVector(Int, cell), sanitize_cellindices(inds), OrbitalLikeGrouped(groups)) + # LatticeSlice from an AbstractVector of CellIndices LatticeSlice(lat::Lattice, cs::AbstractVector{<:CellIndices}) = LatticeSlice(lat, CellIndicesDict(cs)) diff --git a/test/test_greenfunction.jl b/test/test_greenfunction.jl index d2b7c28d..85ac8228 100644 --- a/test/test_greenfunction.jl +++ b/test/test_greenfunction.jl @@ -10,7 +10,7 @@ function testgreen(h, s; kw...) L = Quantica.latdim(lattice(h)) z = zero(SVector{L,Int}) o = Quantica.unitvector(1, SVector{L,Int}) - locs = (cellsites(z, :), cellsites(z, 2:3), cellsites(z, 2), cellsites(o, :), CellOrbitals(o, 1), CellOrbitals(z, 2:3)) + locs = (cellsites(z, :), cellsites(z, 2), cellsites(z, 1:2), cellsites(z, (1,2)), cellsites(o, :), CellOrbitals(o, 1), CellOrbitals(z, 1:2), CellOrbitals(z, SA[2,1])) for loc in locs, loc´ in locs gs = g[loc, loc´] @test gs isa GreenSlice @@ -166,6 +166,36 @@ end end end +@testset "greenfunction 32bit" begin + # GS.Bands + # skip for older versions due to https://github.com/JuliaLang/julia/issues/53054 + # which makes solve not deterministic when adding sufficiently many band simplices + if VERSION >= v"1.11.0-DEV.632" + h = HP.graphene(type = Float32) + s = GS.Bands() + testgreen(h, s) + end + h = HP.graphene(type = Float32, a0 = 1) |> supercell(10) |> supercell + s = GS.SparseLU() + testgreen(h, s) + + # GS.Schur + h = HP.graphene(type = Float32, a0 = 1) |> supercell((1,-1), region = r -> 0<=r[2]<=3) + s = GS.Schur() + testgreen(h, s) + s = GS.Schur(boundary = -1) + testgreen(h, s) + + # GS.KPM + g = HP.graphene(a0 = 1, t0 = 1, orbitals = (2,1), type = Float32) |> + supercell(region = RP.circle(20)) |> + attach(nothing, region = RP.circle(1)) |> greenfunction(GS.KPM(order = 300, bandrange = (-3.1, 3.1))) + ρs = ldos(g[1], kernel = missing) + for ω in -3:0.1:3 + @test all(>=(0), ρs(ω)) + end +end + @testset "diagonal slicing" begin g = HP.graphene(a0 = 1, t0 = 1, orbitals = (2,1)) |> supercell((2,-2), region = r -> 0<=r[2]<=5) |> attach(nothing, cells = 2, region = r -> 0<=r[2]<=2) |> attach(nothing, cells = 3) |>