From a5df211a79b6c56f2cbfa8179124995a96a3a47c Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Fri, 14 May 2021 18:31:06 +0200 Subject: [PATCH 1/4] green functions through deflation draft heavy wip debugging wip wip Jordan chain works new method for Jordan compromise: r0 = 1 wip one example down working! atol = 0 defaults to Wimmer refactor + fastpaths fixes improved rank-defect handling cleanup fix type instabilities from schur another fix in recursion consistency fix safety cutoff in jordan loop C1 deflate method remove debug print Deflating Schur hessenberg approach reducing allocations reduce more allocations elliminate more allocations cleanup fix memory overwrite + bugfix optimize 1=>n green by using dense invG simplify again, inplace lu! for g(1=>n) same in-place opt for local-infinite-fastpath in-place optimizations on top of the QR C1 deflate removed debug print tweak show greens remove redundant shiftw! remove SpecialFunctions dep cover the case with zero h_+- flatten/unflatten greens r x r Jordan reconstruction with tests unflatten! fix test fix --- Project.toml | 3 +- src/Quantica.jl | 4 +- src/bandstructure.jl | 9 +- src/diagonalizer.jl | 2 +- src/greens.jl | 783 ++++++++++++++++++++----------------- src/hamiltonian.jl | 176 +++++---- src/lattice.jl | 8 +- src/plot_vegalite.jl | 8 +- src/tools.jl | 134 ++++++- test/test_bandstructure.jl | 6 +- test/test_greens.jl | 67 ++-- 11 files changed, 711 insertions(+), 489 deletions(-) diff --git a/Project.toml b/Project.toml index 6615cc05..bcdf4683 100644 --- a/Project.toml +++ b/Project.toml @@ -14,9 +14,9 @@ ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Requires = "ae029012-a4dd-5104-9daa-d747884805df" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" -SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" [compat] ExprTools = "^0.1" @@ -26,7 +26,6 @@ NearestNeighbors = "0.4" OffsetArrays = "1.0" ProgressMeter = "1.2" Requires = "1.0" -SpecialFunctions = "0.10, 1.0" StaticArrays = "0.12.3, 0.13, 1.0" julia = "1.5" diff --git a/src/Quantica.jl b/src/Quantica.jl index 8812acae..509da604 100644 --- a/src/Quantica.jl +++ b/src/Quantica.jl @@ -11,7 +11,7 @@ function __init__() end using StaticArrays, NearestNeighbors, SparseArrays, LinearAlgebra, OffsetArrays, - ProgressMeter, LinearMaps, Random, SpecialFunctions + ProgressMeter, LinearMaps, Random, SuiteSparse using ExprTools @@ -29,7 +29,7 @@ export sublat, bravais, lattice, dims, supercell, unitcell, bands, vertices, energies, states, degeneracy, momentaKPM, dosKPM, averageKPM, densityKPM, bandrangeKPM, - greens, greensolver, SingleShot1D + greens, greensolver, Schur1D export RegionPresets, RP, LatticePresets, LP, HamiltonianPresets, HP diff --git a/src/bandstructure.jl b/src/bandstructure.jl index fcc478c8..32d73ff2 100644 --- a/src/bandstructure.jl +++ b/src/bandstructure.jl @@ -11,7 +11,7 @@ Subspace(h::Hamiltonian, args...) = Subspace(h.orbstruct, args...) Subspace(::Missing, energy, basis, basevert...) = Subspace(OrbitalStructure(eltype(basis), size(basis, 1)), energy, basis, basevert...) Subspace(o::OrbitalStructure, energy, basis, basevert...) = - Subspace(energy, unflatten_or_reinterpret(basis, o), o, basevert...) + Subspace(energy, unflatten_orbitals_or_reinterpret(basis, o), o, basevert...) Subspace(energy::T, basis, orbstruct) where {T} = Subspace(energy, basis, orbstruct, SVector{0,T}()) Subspace(energy::T, basis, orbstruct, basevert) where {T} = Subspace(energy, basis, orbstruct, SVector(T.(basevert))) @@ -38,9 +38,10 @@ degeneracy(m::AbstractMatrix) = isempty(m) ? 1 : size(m, 2) # To support sentin orbitalstructure(s::Subspace) = s.orbstruct -flatten(s::Subspace) = Subspace(s.energy, _flatten(s.basis, s.orbstruct), s.orbstruct, s.basevert) +flatten(s::Subspace) = + Subspace(s.energy, flatten(s.basis, s.orbstruct), s.orbstruct, s.basevert) -unflatten(s::Subspace, o::OrbitalStructure) = Subspace(s.energy, unflatten(s.basis, o), o, s.basevert) +unflatten(s::Subspace, o::OrbitalStructure) = Subspace(s.energy, unflatten_orbitals(s.basis, o), o, s.basevert) # destructuring Base.iterate(s::Subspace) = s.energy, Val(:basis) @@ -51,6 +52,8 @@ Base.first(s::Subspace) = s.energy Base.last(s::Subspace) = s.basis Base.length(s::Subspace) = 2 +Base.copy(s::Subspace) = deepcopy(s) + ####################################################################### # Spectrum ####################################################################### diff --git a/src/diagonalizer.jl b/src/diagonalizer.jl index 3f5d010f..02231449 100644 --- a/src/diagonalizer.jl +++ b/src/diagonalizer.jl @@ -77,7 +77,7 @@ end function (d::Diagonalizer)(φs) ϵ, ψ = d(φs, NoUnflatten()) - ψ´ = unflatten_or_reinterpret(ψ, d.orbstruct) + ψ´ = unflatten_orbitals_or_reinterpret(ψ, d.orbstruct) return ϵ, ψ´ end diff --git a/src/greens.jl b/src/greens.jl index 306fd642..2c717601 100644 --- a/src/greens.jl +++ b/src/greens.jl @@ -1,5 +1,5 @@ ####################################################################### -# Green's function +# Green function ####################################################################### abstract type AbstractGreensSolver end @@ -17,7 +17,7 @@ the provided `solveobject`. Currently valid `solveobject`s are - the `Bandstructure` of `h` (for an unbounded `h` or an `Hamiltonian{<:Superlattice}}`) - the `Spectrum` of `h` (for a bounded `h`) -- `SingleShot1D(; direct = false)` (single-shot generalized [or direct if `direct = true`] eigenvalue approach for 1D Hamiltonians) +- `Schur1D(; direct = false)` (single-shot generalized [or direct if `direct = true`] eigenvalue approach for 1D Hamiltonians) If a `boundaries = (n₁, n₂, ...)` is provided, a reflecting boundary is assumed for each non-missing `nᵢ` perpendicular to Bravais vector `i` at a cell distance `nᵢ` from the @@ -59,7 +59,7 @@ julia> m = similarmatrix(g); g(m, 0.2) ``` # See also - `greens!`, `SingleShot1D` + `greens!`, `Schur1D` """ greens(h::Hamiltonian{<:Any,L}, solverobject; boundaries = filltuple(missing, Val(L))) where {L} = GreensFunction(greensolver(solverobject, h), h, boundaries) @@ -68,29 +68,45 @@ greens(solver::Function, args...; kw...) = h -> greens(h, solver(h), args...; kw # solver fallback greensolver(s::AbstractGreensSolver) = s +# missing cells +(g::GreensFunction)(ω; kw...) = g(ω, default_cells(g); kw...) + # call API fallback -(g::GreensFunction)(ω, cells = default_cells(g)) = greens!(similarmatrix(g), g, ω, cells) +(g::GreensFunction)(ω, cells; kw...) = greens!(similarmatrix(g), g, ω, cells; kw...) similarmatrix(g::GreensFunction, type = Matrix{blocktype(g.h)}) = similarmatrix(g.h, type) -greens!(matrix, g, ω, cells) = greens!(matrix, g, ω, sanitize_cells(cells, g)) +greens!(matrix, g, ω, cells; kw...) = greens!(matrix, g, ω, sanitize_cells(cells, g); kw...) + +default_cells(g::GreensFunction) = _plusone.(g.boundaries) => _plusone.(g.boundaries) -default_cells(::GreensFunction{S,L}) where {S,L} = filltuple(1, Val(L)) => filltuple(1, Val(L)) +_plusone(::Missing) = 1 +_plusone(n) = n + 1 sanitize_cells((cell0, cell1)::Pair{<:Integer,<:Integer}, ::GreensFunction{S,1}) where {S} = SA[cell0] => SA[cell1] sanitize_cells((cell0, cell1)::Pair{<:NTuple{L,Integer},<:NTuple{L,Integer}}, ::GreensFunction{S,L}) where {S,L} = SVector(cell0) => SVector(cell1) sanitize_cells(cells, g::GreensFunction{S,L}) where {S,L} = - throw(ArgumentError("Cells should be of the form `cᵢ => cⱼ`, with each `c` an `NTuple{$L,Integer}`")) + throw(ArgumentError("Cells should be of the form `cᵢ => cⱼ`, with each `c` an `NTuple{$L,Integer}`, got $cells")) const SVectorPair{L} = Pair{SVector{L,Int},SVector{L,Int}} +Base.size(g::GreensFunction, args...) = size(g.h, args...) +Base.eltype(g::GreensFunction) = eltype(g.h) + +flatsize(g::GreensFunction, args...) = flatsize(g.h, args...) +blockeltype(g::GreensFunction) = blockeltype(g.h) +blocktype(g::GreensFunction) = blocktype(g.h) +orbitaltype(g::GreensFunction) = orbitaltype(g.h) +orbitalstructure(g::GreensFunction) = orbitalstructure(g.h) + ####################################################################### -# SingleShot1DGreensSolver +# Schur1DGreensSolver ####################################################################### + """ - SingleShot1D(; direct = false) + Schur1D(; deflation = default_tol(T)) Return a Greens function solver using the generalized eigenvalue approach, whereby given the energy `ω`, the eigenmodes of the infinite 1D Hamiltonian, and the corresponding infinite @@ -98,22 +114,20 @@ and semi-infinite Greens function can be computed by solving the generalized eig equation A⋅φχ = λ B⋅φχ - A = [0 I; V ω-H0] - B = [I 0; 0 V'] + A = [0 I; -h₊ ω-h₀] + B = [I 0; 0 h₋] -This is the matrix form of the problem `λ(ω-H0)φ - Vφ - λ²V'φ = 0`, where `φχ = [φ; λφ]`, +This is the matrix form of the problem `λ(ω-h₀)φ - h₊φ - λ²h₋φ = 0`, where `φχ = [φ; λφ]`, and `φ` are `ω`-energy eigenmodes, with (possibly complex) momentum `q`, and eigenvalues are `λ = exp(-iqa₀)`. The algorithm assumes the Hamiltonian has only `dn = (0,)` and `dn = (±1, -)` Bloch harmonics (`H0`, `V` and `V'`), so its unit cell will be enlarged before applying +)` Bloch harmonics (`h₀`, `h₊` and `h₋`), so its unit cell will be enlarged before applying the solver if needed. Bound states in the spectrum will yield delta functions in the density of states that can be resolved by adding a broadening in the form of a small positive imaginary part to `ω`. -To avoid singular solutions `λ=0,∞`, the nullspace of `V` is projected out of the problem. -This produces a new `A´` and `B´` with reduced dimensions. `B´` can often be inverted, -turning this into a standard eigenvalue problem, which is slightly faster to solve. This is -achieved with `direct = true`. However, `B´` sometimes is still non-invertible for some -values of `ω`. In this case use `direct = false` (the default). +For performace, the eigenvalue equation may be `deflated', i.e. singular solutions `λ=0,∞` +will be removed within the absolute tolerance specified by the keyword `deflation`. If no +deflation is desired, use `deflation = nothing`. # Examples ```jldoctest @@ -121,10 +135,10 @@ julia> using LinearAlgebra julia> h = LP.honeycomb() |> hamiltonian(hopping(1)) |> unitcell((1,-1), (10,10)) |> Quantica.wrap(2); -julia> g = greens(h, SingleShot1D(), boundaries = (0,)) -GreensFunction{SingleShot1DGreensSolver}: Green's function using the single-shot 1D method +julia> g = greens(h, Schur1D(), boundaries = (0,)) +GreensFunction{Schur1DGreensSolver}: Green's function using the Schur1D method Matrix size : 40 × 40 - Reduced size : 20 × 20 + Deflated size : 20 × 20 Element type : scalar (ComplexF64) Boundaries : (0,) @@ -135,414 +149,461 @@ julia> tr(g(0.3)) # See also `greens` """ -struct SingleShot1D - invert_B::Bool - cutoff::Int +struct Schur1D{R} + atol::R # could be missing for default_tol(T) end -SingleShot1D(; direct = false, cutoff = 1) = SingleShot1D(direct, cutoff) - -struct SingleShot1DTemporaries{T} - b2s::Matrix{T} # size = dim_b, 2dim_s - φ::Matrix{T} # size = dim_H0, 2dim_s - χ::Matrix{T} # size = dim_H0, 2dim_s - ss1::Matrix{T} # size = dim_s, dim_s - ss2::Matrix{T} # size = dim_s, dim_s - Hs1::Matrix{T} # size = dim_H0, dim_s - Hs2::Matrix{T} # size = dim_H0, dim_s - HH0::Matrix{T} # size = dim_H0, dim_H0 - HH1::Matrix{T} # size = dim_H0, dim_H0 - HH2::Matrix{T} # size = dim_H0, dim_H0 - HH3::Matrix{T} # size = dim_H0, dim_H0 - vH::Vector{T} # size = dim_H0 - v2s1::Vector{Int} # size = 2dim_s - v2s2::Vector{Int} # size = 2dim_s +Schur1D(; deflation = missing) = Schur1D(deflation) + +struct DeflatorWorkspace{T} + nn::Matrix{T} + nr1::Matrix{T} + nr2::Matrix{T} + rr0::Matrix{T} + rr1::Matrix{T} + rr2::Matrix{T} + rr3::Matrix{T} + rr4::Matrix{T} + mb::Matrix{T} end -function SingleShot1DTemporaries{T}(dim_H, dim_s, dim_b) where {T} - b2s = Matrix{T}(undef, dim_b, 2dim_s) - φ = Matrix{T}(undef, dim_H, 2dim_s) - χ = Matrix{T}(undef, dim_H, 2dim_s) - ss1 = Matrix{T}(undef, dim_s, dim_s) - ss2 = Matrix{T}(undef, dim_s, dim_s) - Hs1 = Matrix{T}(undef, dim_H, dim_s) - Hs2 = Matrix{T}(undef, dim_H, dim_s) - HH0 = Matrix{T}(undef, dim_H, dim_H) - HH1 = Matrix{T}(undef, dim_H, dim_H) - HH2 = Matrix{T}(undef, dim_H, dim_H) - HH3 = Matrix{T}(undef, dim_H, dim_H) - vH = Vector{T}(undef, dim_H) - v2s1 = Vector{Int}(undef, 2dim_s) - v2s2 = Vector{Int}(undef, 2dim_s) - return SingleShot1DTemporaries{T}(b2s, φ, χ, ss1, ss2, Hs1, Hs2, HH0, HH1, HH2, HH3, vH, v2s1, v2s2) +DeflatorWorkspace{T}(n, r) where {T} = + DeflatorWorkspace(Matrix{T}.(undef, + ((n, n), (n, r), (n, r), (r, r), (r, r), (r, r), (r, r), (r, r), (n+r, n-r)))...) + +struct Deflator{T,M<:AbstractMatrix{T},R<:Real,S} + hp::M + hm::M + hmQ0::M # h₋*Q0 where Q0 = [rowspace(A0) nullspace(A0)]. h₊ = [R' 0] Q0'. h₋ = Q0 [R; 0] + R::Matrix{T} # A0 = [-hRR 0; -hBR 0] * [R'; B' ]. R = orthogonal complement of nullspace(A0) === rowspace(A0) + Adef::Matrix{T} # deflated A + Bdef::Matrix{T} # deflated B + Ablock::Matrix{T} # Adef = Ablock * QR; Ablock = [0 I 0; -hRR gRR⁻¹ gRB⁻¹] + Bblock::Matrix{T} # Bdef = Bblock * QR; Bblock = [I 0 0; 0 hRR' hBR'] + Vblock´::Matrix{T} # # Vblock = [-hBR gBR⁻¹ gBB⁻¹] = [R 0] [QB'; QR']. Vblock´ = Matrix(Vblock') + ig0::Matrix{T} # Matrix(-h₀) + ωshifter::S # metadata to aid in ω-shifting the relevant A subblocks + atol::R # A0, A2 deflation tolerance + tmp::DeflatorWorkspace{T} end -struct SingleShot1DGreensSolver{T<:Complex,R<:Real,H<:Hessenberg{T}} <: AbstractGreensSolver - invert_B::Bool - λ2min::R - A::Matrix{T} - B::Matrix{T} - minusH0::SparseMatrixCSC{T,Int} - V::SparseMatrixCSC{T,Int} - Pb::Matrix{T} # size = dim_b, dim_H0 - Ps::Matrix{T} # size = dim_s, dim_H0 - Ps´::Matrix{T} # size = dim_s´, dim_H0 - H0ss::Matrix{T} - H0bs::Matrix{T} - Vss::Matrix{T} - Vbs::Matrix{T} - hessbb::H - velocities::Vector{R} # size = 2dim_s = 2num_modes - maxdn::Int - temps::SingleShot1DTemporaries{T} +struct Schur1DGreensSolver{D<:Union{Deflator,Missing},M} <: AbstractGreensSolver + h0::M + hp::M + hm::M + unitcells::Int + deflatorR::D + deflatorL::D end -function Base.show(io::IO, g::GreensFunction{<:SingleShot1DGreensSolver}) +function Base.show(io::IO, g::GreensFunction{<:Schur1DGreensSolver}) print(io, summary(g), "\n", -" Matrix size : $(size(g.solver.V, 1)) × $(size(g.solver.V, 2)) - Reduced size : $(size(g.solver.Ps, 1)) × $(size(g.solver.Ps, 1)) +" Matrix size : $(size(g.solver.h0, 1)) × $(size(g.solver.h0, 2)) + Deflated size : $(deflated_size_text(g)) Element type : $(displayelements(g.h)) Boundaries : $(g.boundaries)") end -Base.summary(g::GreensFunction{<:SingleShot1DGreensSolver}) = - "GreensFunction{SingleShot1DGreensSolver}: Green's function using the single-shot 1D method" +function deflated_size_text(g::GreensFunction) + text = hasdeflator(g.solver) <= 0 ? "No deflation" : + "$(deflated_size_text(g.solver.deflatorR))" + return text +end + +deflated_size_text(d::Deflator) = "$(size(d.Adef, 1) ÷ 2) × $(size(d.Adef, 2) ÷ 2)" -hasbulk(gs::SingleShot1DGreensSolver) = !iszero(size(gs.Pb, 1)) +Base.summary(g::GreensFunction{<:Schur1DGreensSolver}) = + "GreensFunction{Schur1DGreensSolver}: Green's function using the Schur1D method" -function greensolver(s::SingleShot1D, h) - latdim(h) == 1 || throw(ArgumentError("Cannot use a SingleShot1D Green function solver with an $(latdim(h))-dimensional Hamiltonian")) - return SingleShot1DGreensSolver(h, s.invert_B, s.cutoff) -end +Base.size(s::Schur1DGreensSolver, args...) = size(s.deflatorR, args...) +Base.size(d::Deflator, args...) = size(d.Adef, args...) -## Preparation +hasdeflator(::Schur1DGreensSolver{<:Deflator}) = true +hasdeflator(::Schur1DGreensSolver{Missing}) = false -function SingleShot1DGreensSolver(h::Hamiltonian, invert_B, cutoff) - latdim(h) == 1 || throw(ArgumentError("Cannot use a SingleShot1D Green function solver with an $(latdim(h))-dimensional Hamiltonian")) - maxdn = max(1, maximum(har -> abs(first(har.dn)), h.harmonics)) - H = flatten(maxdn == 1 ? h : unitcell(h, (maxdn,))) +function greensolver(s::Schur1D, h) + latdim(h) == 1 || throw(ArgumentError("Cannot use a Schur1D Green function solver with an $(latdim(h))-dimensional Hamiltonian")) + unitcells = max(1, maximum(har -> abs(first(har.dn)), h.harmonics)) + H = maybeflatten(unitcells == 1 ? h : unitcell(h, (unitcells,))) + h₊, h₀, h₋ = H[(1,)], H[(0,)], H[(-1,)] + n = size(H, 1) T = complex(blockeltype(H)) - λ2min = cutoff^2 * eps(real(T)) - H0, V, V´ = H[(0,)], H[(1,)], H[(-1,)] - Pb, Ps, Ps´ = bulk_surface_projectors(H0, V, V´, λ2min) - H0ss = Ps * H0 * Ps' - H0bs = Pb * H0 * Ps' - Vss = Ps * V * Ps' - Vbs = Pb * V * Ps' - hessbb = hessenberg!(Pb * (-H0) * Pb') - dim_s = size(Ps, 1) - A = zeros(T, 2dim_s, 2dim_s) - B = zeros(T, 2dim_s, 2dim_s) - dim_s, dim_b, dim_H = size(Ps, 1), size(Pb, 1), size(H0, 2) - velocities = fill(zero(real(T)), 2dim_s) - temps = SingleShot1DTemporaries{T}(dim_H, dim_s, dim_b) - return SingleShot1DGreensSolver(invert_B, λ2min, A, B, -H0, V, Pb, Ps, Ps´, H0ss, H0bs, Vss, Vbs, hessbb, velocities, maxdn, temps) + atol = s.atol === missing ? default_tol(T) : s.atol + deflatorR = Deflator(atol, h₊, h₀, h₋) + deflatorL = Deflator(atol, h₋, h₀, h₊) + return Schur1DGreensSolver(h₀, h₊, h₋, unitcells, deflatorR, deflatorL) end -function bulk_surface_projectors(H0::AbstractMatrix{T}, V, V´, cutoff2) where {T} - SVD = svd(Matrix(V), full = true) - W, S, U = SVD.U, SVD.S, SVD.V - dim_b = count(s -> abs2(s) < cutoff2, S) - dim_s = length(S) - dim_b - if iszero(dim_b) - Ps = Matrix{T}(I, dim_s, dim_s) - Ps´ = copy(Ps) - Pb = Ps[1:0, :] - else - Ps = U'[1:dim_s, :] - Pb = U'[dim_s+1:end, :] - Ps´ = W'[1:dim_s, :] - Pb´ = W'[dim_s+1:end, :] - end - return Pb, Ps, Ps´ - return Pb, Ps, Ps´ +Deflator(atol::Nothing, As...) = missing + +function Deflator(atol::Real, h₊::M, h₀::M, h₋::M) where {M} + rowspaceR, _, nullspaceR = fullrank_decomposition_qr(h₊, atol) + B = Matrix(nullspaceR) # nullspace(A0) + R = Matrix(rowspaceR) # orthogonal complement of nullspace(h₊) + h₋Q0 = h₋ * parent(rowspaceR) # h₋ * [R B] = h₋ * Q0, needed for Jordan chain + n = size(h₀, 2) + r = size(R, 2) + b = size(B, 2) + T = eltype(h₀) + h₊R = h₊*R + hRR = R'*h₊R + hBR = B'*h₊R + gRR⁻¹ = - R'*h₀*R + gBB⁻¹ = - B'*h₀*B + gRB⁻¹ = - R'*h₀*B + gBR⁻¹ = gRB⁻¹' + g0⁻¹ = Matrix(-h₀) + Adef = Matrix{T}(undef, 2r, 2r) # Needs to be dense for schur!(Adef, Bdef) + Bdef = Matrix{T}(undef, 2r, 2r) # Needs to be dense for schur!(Adef, Bdef) + Ablock = Matrix([0I I spzeros(r, b); -hRR gRR⁻¹ gRB⁻¹]) + Bblock = Matrix([I 0I spzeros(r, b); 0I hRR' hBR']) + Vblock´ = Matrix([-hBR gBR⁻¹ gBB⁻¹]') + ωshifter = diag(gRR⁻¹), (r+1:2r, r+1:2r), diag(gBB⁻¹), (2r+1:2r+b, 1:b), diag(g0⁻¹) + tmp = DeflatorWorkspace{T}(n, r) + return Deflator(h₊, h₋, h₋Q0, R, Adef, Bdef, Ablock, Bblock, Vblock´, g0⁻¹, ωshifter, atol, tmp) end -## Solver execution - -function (gs::SingleShot1DGreensSolver)(ω) - A, B = single_shot_surface_matrices(gs, ω) - λs, φχs = eigen_funcbarrier!(A, B, gs) - dim_s = size(φχs, 1) ÷ 2 - φs = view(φχs, 1:dim_s, :) - χs = view(φχs, dim_s+1:2dim_s, :) - φ = gs.temps.φ - χ = gs.temps.χ - if hasbulk(gs) - reconstruct_bulk!(φ, ω, λs, φs, gs) - reconstruct_bulk!(χ, ω, λs, χs, gs) - else - copy!(φ, φs) - copy!(χ, χs) +## Tools + +function shiftω!(d::Deflator, ω) + diagRR, rowcolA, diagBB, rowcolV, diagg0⁻¹ = d.ωshifter + for (v, row, col) in zip(diagRR, rowcolA...) + d.Ablock[row, col] = ω + v end - return classify_retarded_advanced(λs, φs, χs, φ, χ, gs) + for (v, row, col) in zip(diagBB, rowcolV...) + d.Vblock´[row, col] = conj(ω) + v + end + for (n, v) in enumerate(diagg0⁻¹) + d.ig0[n, n] = ω + v + end + return d end -function eigen_funcbarrier!(A::AbstractMatrix{T}, B, gs)::Tuple{Vector{T},Matrix{T}} where {T<:Complex} - if gs.invert_B - factB = lu!(B) - B⁻¹A = ldiv!(factB, A) - λs, φχs = eigen!(B⁻¹A; sortby = abs) - clean_λ!(λs, gs.λ2min) - else - alpha, beta, _, φχ´ = LAPACK.ggev!('N', 'V', A, B) - λ´ = clean_λ!(alpha, beta, gs.λ2min) - λs, φχs = GeneralizedEigen(LinearAlgebra.sorteig!(λ´, φχ´, abs)...) +function shiftω!(mat::AbstractMatrix, ω) + @inbounds for i in axes(mat, 1) + mat[i, i] += ω end - return λs, φχs + return mat end -function clean_λ!(λ::AbstractVector{T}, cutoff2) where {T} - λs .= (λ -> ifelse(isnan(λ), T(Inf), - ifelse(abs2(λ) < λ2min, zero(T), - ifelse(abs2(λ) > 1/λ2min, T(Inf), λ)))).(λs) - return λs +function orthobasis_decomposition_qr(mat, atol) + q = pqr(mat) + basis = getQ(q) + RP´ = getRP´(q) + n = size(basis, 2) + r = nonzero_rows(RP´, atol) + orthobasis = view(basis, :, 1:r) + complement = view(basis, :, r+1:n) + r = view(RP´, 1:r, :) + return orthobasis, r, complement end -function clean_λ!(αs::AbstractVector{T}, βs, cutoff2) where {T} - λs = αs - for i in eachindex(αs) - α2, β2 = abs2(αs[i]), abs2(βs[i]) - if α2 < cutoff2 - λs[i] = zero(T) - elseif β2 < cutoff2 - λs[i] = T(Inf) - else - λs[i] = αs[i] / βs[i] - end - end - return λs +function fullrank_decomposition_qr(mat, atol) + rowspace, r, nullspace = orthobasis_decomposition_qr(mat', atol) + return rowspace, r', nullspace end -function single_shot_surface_matrices(gs::SingleShot1DGreensSolver{T}, ω) where {T} - A, B = gs.A, gs.B - dim_s = size(gs.H0bs, 2) - fill!(A, 0) - fill!(B, 0) - for i in 1:dim_s - A[i, dim_s + i] = B[i, i] = one(T) - A[dim_s + i, dim_s + i] = ω - end - tmp = view(gs.temps.b2s, :, 1:dim_s) # gs.temps.b2s has 2dim_s cols - dim = size(A, 1) ÷ 2 - A21 = view(A, dim+1:2dim, 1:dim) - A22 = view(A, dim+1:2dim, dim+1:2dim) - B22 = view(B, dim+1:2dim, dim+1:2dim) - - # A22 = ωI - H₀ₛₛ - H₀ₛᵦ g₀ᵦᵦ H₀ᵦₛ - V'ₛᵦ g₀ᵦᵦ Vᵦₛ - A22 .-= gs.H0ss # ω already added to diagonal above - if hasbulk(gs) - copy!(tmp, gs.H0bs) - ldiv!(gs.hessbb + ω*I, tmp) - mul!(A22, gs.H0bs', tmp, -1, 1) - copy!(tmp, gs.Vbs) - ldiv!(gs.hessbb + ω*I, tmp) - mul!(A22, gs.Vbs', tmp, -1, 1) - end +nullspace_qr(mat, atol) = last(fullrank_decomposition_qr(mat, atol)) - # A21 = - Vₛₛ - H₀ₛᵦ g₀ᵦᵦ Vᵦₛ - A21 .= .- gs.Vss - if hasbulk(gs) - copy!(tmp, gs.Vbs) - ldiv!(gs.hessbb + ω*I, tmp) - mul!(A21, gs.H0bs', tmp, -1, 1) - end +rowspace_qr(mat, atol) = first(fullrank_decomposition_qr(mat, atol)) + +## Deflate - # B22 = -A21' - B22 .= .- A21' +function deflate(d::Deflator{T,<:SparseMatrixCSC}) where {T} + r = size(d.R, 2) + m = size(d.Vblock´, 1) # m = 2r+b + # Vblock´ = [RP´ 0] * Q' = b × 2r+b; Q = [rowspaceV nullspaceV] = 2r+b × 2r+b = m × m + # nullspaceV is 2r+b × 2r + Vblock´ = copy!(d.tmp.mb, d.Vblock´) # ::Matrix{T} + nullspaceV = getQ(qr!(Vblock´, Val(true)), m-2r+1:m) + Adef = mul!(d.Adef, d.Ablock, nullspaceV) + Bdef = mul!(d.Bdef, d.Bblock, nullspaceV) + Q1 = view(nullspaceV, 1:r, :) + Q2 = view(nullspaceV, r+1:m, :) + return Adef, Bdef, Q1, Q2 +end - chkfinite(A) - chkfinite(B) +## Solver execution: compute self-energy, with or without deflation +(s::Schur1DGreensSolver)(ω, which = Val{:R}) = schursolver(s, ensurecomplex(ω), which) - return A, B +ensurecomplex(ω::T) where {T<:Real} = ω + im * default_tol(T) +ensurecomplex(ω::Complex) = ω + +function schursolver(s::Schur1DGreensSolver{Missing}, ω, which) + A = Matrix([ω*I - s.h0 -s.hp; -I 0I]) + B = Matrix([s.hm 0I; 0I -I]) + sch = schur!(A, B) + Σ = nondeflated_selfenergy(which, s, sch) + # @show sum(abs.(Σ - s.hm * ((ω*I - s.h0 - Σ) \ Matrix(s.hp)))) + return Σ end -function chkfinite(A::AbstractMatrix) - for a in A - if !isfinite(a) - throw(ArgumentError("Matrix contains Infs or NaNs. This may happen when the energy ω exactly matches a bound state in the spectrum. Try adding a small positive imaginary part to ω.")) - end - end - return true +function nondeflated_selfenergy(::Type{Val{:R}}, s, sch) + n = size(s.h0, 1) + ordschur!(sch, abs.(sch.α ./ sch.β) .<= 1) + ϕΛR⁻¹ = view(sch.Z, 1:n, 1:n) + ϕR⁻¹ = view(sch.Z, n+1:2n, 1:n) + ΣR = s.hm * ϕΛR⁻¹ / ϕR⁻¹ + return ΣR end -# φ = [Pₛ' + Pᵦ' g₀ᵦᵦ (λ⁻¹Vᵦₛ + H₀ᵦₛ)] φₛ -function reconstruct_bulk!(φ, ω, λs, φs, gs) - tmp = gs.temps.b2s - mul!(tmp, gs.Vbs, φs) - tmp ./= transpose(λs) - mul!(tmp, gs.H0bs, φs, 1, 1) - ldiv!(gs.hessbb + ω*I, tmp) - mul!(φ, gs.Pb', tmp) - mul!(φ, gs.Ps', φs, 1, 1) - return φ +function nondeflated_selfenergy(::Type{Val{:L}}, s, sch) + n = size(s.h0, 1) + ordschur!(sch, abs.(sch.β ./ sch.α) .<= 1) + ϕΛ⁻¹R⁻¹ = view(sch.Z, 1:n, 1:n) + ϕR⁻¹ = view(sch.Z, n+1:2n, 1:n) + ΣL = s.hp * ϕR⁻¹ / ϕΛ⁻¹R⁻¹ + return ΣL end -# function classify_retarded_advanced(λs, φs, φ, χ, gs) -function classify_retarded_advanced(λs, φs, χs, φ, χ, gs) - vs = compute_velocities_and_normalize!(φ, χ, λs, gs) - # order for ret-evan, ret-prop, adv-prop, adv-evan - p = sortperm!(gs.temps.v2s1, vs; rev = true) - p´ = gs.temps.v2s2 - Base.permute!!(vs, copy!(p´, p)) - Base.permute!!(λs, copy!(p´, p)) - Base.permutecols!!(φ, copy!(p´, p)) - Base.permutecols!!(χ, copy!(p´, p)) - Base.permutecols!!(φs, copy!(p´, p)) - - ret, adv = nonsingular_ret_adv(λs, vs) - # @show vs - # @show abs.(λs) - # @show ret, adv, length(λs), length(vs) - - λR = view(λs, ret) - χR = view(χ, :, ret) # This resides in part of gs.temps.χ - φR = view(φ, :, ret) # This resides in part of gs.temps.φ - # overwrite output of eigen to preserve normalization of full φ - φRs = mul!(view(φs, :, ret), gs.Ps, φR) - iφRs = issquare(φRs) ? rdiv!(copyto!(gs.temps.ss1, I), lu!(φRs)) : pinv(φRs) - - iλA = view(λs, adv) - iλA .= inv.(iλA) - χA = view(χ, :, adv) # This resides in part of gs.temps.χ - φA = view(φ, :, adv) # This resides in part of gs.temps.φ - # overwrite output of eigen to preserve normalization of full χ - χAs´ = mul!(view(χs, :, adv), gs.Ps´, χA) - iχAs´ = issquare(χAs´) ? rdiv!(copyto!(gs.temps.ss2, I), lu!(χAs´)) : pinv(χAs´) - - return λR, χR, iφRs, iλA, φA, iχAs´ +nondeflated_selfenergy(::Type{Val{:RL}}, s, sch) = + nondeflated_selfenergy(Val{:R}, s, sch), nondeflated_selfenergy(Val{:L}, s, sch) + +schursolver(s::Schur1DGreensSolver{<:Deflator}, ω, ::Type{Val{:R}}) = + deflated_selfenergy(s.deflatorR, s, ω) + +schursolver(s::Schur1DGreensSolver{<:Deflator}, ω, ::Type{Val{:L}}) = + deflated_selfenergy(s.deflatorL, s, ω) + +schursolver(s::Schur1DGreensSolver{<:Deflator}, ω, ::Type{Val{:RL}}) = + deflated_selfenergy(s.deflatorR, s, ω), deflated_selfenergy(s.deflatorL, s, ω) + +function deflated_selfenergy(d::Deflator{T,M}, s::Schur1DGreensSolver, ω) where {T,M} + # if r = 0 (absolute deflation), Σ=0 + size(d) == (0, 0) && return fill!(d.tmp.nn, zero(T)) + shiftω!(d, ω) + A, B, Q1, Q2 = deflate(d) + # find right-moving eigenvectors with atol < |λ| < 1 + sch = schur!(A, B) + rmodes = retarded_modes(sch, d.atol) + nr = sum(rmodes) + ordschur!(sch, rmodes) + Zret = view(sch.Z, :, 1:nr) + R, h₋Q0 = d.R, d.hmQ0 + ## Qs = [Q1; Q2]; [φR; χR; χB] = Qs * Zret * R11 + ## R'φ = φR = R'Z11 * R11, where R'Z11 = Q1 * Zret + ## Q0'*χ = Q0'*φ*Λ = [χR; χB] + ## h₋χ = h₋ * Q0 * [χR; χB] = h₋ * Q0 * Q2 * Zret * R11 = Z21 * R11, where Z21 = h₋ * Q0 * Q2 * Zret + ## R´Z11 = Q1 * Zret + ## Z21 = h₋Q0 * Q2 * Zret + R´Z11 = Q1 * Zret + Z21 = h₋Q0 * Q2 * Zret + + ## add generalized eigenvectors until we span the full R space + R´source, target = add_jordan_chain(d, R´Z11, Z21) + + ΣR = mul!(d.tmp.nn, rdiv!(target, lu!(R´source)), R') + # @show sum(abs.(ΣR - s.hm * (((ω * I - s.h0) - ΣR) \ Matrix(s.hp)))) + + return ΣR end -# The velocity is given by vᵢ = im * φᵢ'(V'λᵢ-Vλᵢ')φᵢ / φᵢ'φᵢ = 2*imag(χᵢ'Vφᵢ)/φᵢ'φᵢ -function compute_velocities_and_normalize!(φ, χ, λs, gs) - vs = gs.velocities - tmp = gs.temps.vH - for (i, λ) in enumerate(λs) - abs2λ = abs2(λ) - if abs2λ ≈ 1 - φi = view(φ, :, i) - χi = view(χ, :, i) - norm2φi = dot(φi, φi) - mul!(tmp, gs.V, φi) - v = 2*imag(dot(χi, tmp))/norm2φi - φi .*= inv(sqrt(norm2φi * abs(v))) - χi .*= inv(sqrt(norm2φi * abs(v))) - vs[i] = v - else - vs[i] = abs2λ < 1 ? Inf : -Inf - end - end # sortperm(vs) would now give the order of adv-evan, adv-prop, ret-prop, ret-evan - return view(vs, 1:length(λs)) +# need this barrier for type-stability (sch.α and sch.β are finicky) +function retarded_modes(sch, atol) + rmodes = Vector{Bool}(undef, length(sch.α)) + rmodes .= atol .< abs.(sch.α ./ sch.β) .< 1 + return rmodes end -function nonsingular_ret_adv(λs::AbstractVector{T}, vs) where {T} - rmin, rmax = 1, 0 - amin, amax = 1, 0 - for (i, v) in enumerate(vs) - aλ2 = abs2(λs[i]) - if iszero(aλ2) - rmin = i + 1 - elseif v > 0 - rmax = i - amin = i + 1 - elseif v < 0 && isfinite(aλ2) - amax = i - end +function add_jordan_chain(d::Deflator, R´Z11, Z21) + local R´φg_candidates, source_rowspace + g, hg, gh, hgh = integrate_out_bulk!(d) + Σ11 = copy(hgh) + Σn1 = copy(Σ11) + invgΣ = similar(g) + hgΣ = similar(g) + R´source = similar(R´Z11, size(R´Z11, 1), 0) + maxiter = 10 + for n in 1:maxiter + # Σ11 <-- hgh + hg * Σ11 * inv(1 - g*Σ11) * gh + # Σn1 <-- Σn1 * inv(1 - g*Σ11) * gh + one!(invgΣ) + mul!(invgΣ, g, Σ11, -1, 1) + invgΣgh = ldiv!(lu!(invgΣ), copy!(d.tmp.rr0, gh)) + mul!(hgΣ, hg, Σ11) + Σn1´ = mul!(Σ11, Σn1, invgΣgh) # save allocations + Σ11´ = mul!(Σn1, hgΣ, invgΣgh) # save allocations + Σ11´ .+= hgh + Σn1, Σ11 = Σn1´, Σ11´ + + R´φg_candidates = nullspace_qr(Σn1, d.atol) + # iterate [R´Z11 R´φg_candidates] = [R´source 0] [source_rowspace'; ....] until R´source is full rank + # Then Σ = [Z21 φgJ_candidates] ([R´Z11 R´φg_candidates] \ R') = [R´Z11 φgJ_candidates] * source_rowspace * inv(R´source) * R' + # we have built a full-rank basis R´source of the space spanned by [R´Z11 φgJ_candidates], which we can invert + source_rowspace, R´source, _ = fullrank_decomposition_qr([R´Z11 R´φg_candidates], d.atol) + # when R´source is square, it will be full rank and invertible. + size(R´source, 1) == size(R´source, 2) && break end - return rmin:rmax, amin:amax + φgJ_candidates = d.R * Σ11 * R´φg_candidates + target = [Z21 φgJ_candidates] * source_rowspace + # R´source is an Adjoint, must covert to do lu! later + return copy(R´source), target end -## Greens execution - -(g::GreensFunction{<:SingleShot1DGreensSolver})(ω, cells) = g(ω, missing)(sanitize_cells(cells, g)) - -# Infinite: G∞_{N} = GVᴺ G∞_{0} -function (g::GreensFunction{<:SingleShot1DGreensSolver,1,Tuple{Missing}})(ω, ::Missing) - gs = g.solver - factors = gs(ω) - G∞⁻¹ = inverse_G∞!(gs.temps.HH0, ω, gs, factors) |> lu! - return cells -> G_infinite!(gs, factors, G∞⁻¹, cells) +# computes R'*(g0, h₋ g0, g0 h₊, h₋ g0 h₊) R +function integrate_out_bulk!(d::Deflator) + R = d.R + g0⁻¹ = copy!(d.tmp.nn, d.ig0) + # ig0⁻¹R, ig0⁻¹L = g0⁻¹ \ R, g0⁻¹ \ L + lug0⁻¹ = lu!(g0⁻¹) + g0R = ldiv!(lug0⁻¹, copy!(d.tmp.nr1, R)) + R´g0R = mul!(d.tmp.rr1, R', g0R) + h₋g0R = mul!(d.tmp.nr2, d.hm, g0R) + R´h₋g0R = mul!(d.tmp.rr2, R', h₋g0R) + h₊R = mul!(d.tmp.nr1, d.hp, R) + g0h₊R = ldiv!(lug0⁻¹, copy!(d.tmp.nr2, h₊R)) + R´g0h₊R = mul!(d.tmp.rr3, R', g0h₊R) + h₋g0h₊R = mul!(d.tmp.nr1, d.hm, g0h₊R) + R´h₋g0h₊R = mul!(d.tmp.rr4, R', h₋g0h₊R) + return R´g0R, R´h₋g0R, R´g0h₊R, R´h₋g0h₊R end -# Semiinifinite: G_{N,M} = (GVᴺ⁻ᴹ - GVᴺGV⁻ᴹ)G∞_{0} -function (g::GreensFunction{<:SingleShot1DGreensSolver,1,Tuple{Int}})(ω, ::Missing) - gs = g.solver - factors = gs(ω) - G∞⁻¹ = inverse_G∞!(gs.temps.HH0, ω, gs, factors) |> lu! - N0 = only(g.boundaries) - return cells -> G_semiinfinite!(gs, factors, G∞⁻¹, shift_cells(cells, N0)) +### Greens execution + +# Choose codepath +function (g::GreensFunction{<:Schur1DGreensSolver})(ω, cells; kw...) + cells´ = sanitize_cells(cells, g) + if is_infinite_local(cells´, g) + gω = local_fastpath(g, ω; kw...) + elseif is_across_boundary(cells´, g) + gω = Matrix(zero(g.solver.h0)) + elseif is_at_surface(cells´, g) + gω = surface_fastpath(g, ω, dist_to_boundary(cells´, g); kw...) + else # general form + gω = g(ω, missing; kw...)(cells´) + end + gω´ = slice_original_cell(gω, g.solver.unitcells) + gω´´ = maybeunflatten_blocks(gω´, orbitalstructure(g.h)) + return gω´´ end -function invG(g::GreensFunction{<:SingleShot1DGreensSolver}, ω) - gs = g.solver - factors = gs(ω) - onlyhalf = only(g.boundaries) === missing ? false : true - G∞⁻¹ = inverse_G∞!(gs.temps.HH0, ω, gs, factors, onlyhalf) - return G∞⁻¹ +slice_original_cell(mat, n) = n === 1 ? mat : mat[1:size(mat,1) ÷ n, 1:size(mat, 2)] + +is_infinite_local((src, dst), g::GreensFunction{<:Schur1DGreensSolver,1,Tuple{Missing}}) = + only(src) == only(dst) +is_infinite_local(cells, g) = false + +is_at_surface((src, dst), g::GreensFunction{<:Schur1DGreensSolver,1,Tuple{Int}}) = + abs(dist_to_boundary(src, g)) == 1 || abs(dist_to_boundary(dst, g)) == 1 +is_at_surface(cells, g) = false + +is_across_boundary((src, dst), g::GreensFunction{<:Schur1DGreensSolver,1,Tuple{Int}}) = + sign(dist_to_boundary(src, g)) != sign(dist_to_boundary(src, g)) || + dist_to_boundary(src, g) == 0 || dist_to_boundary(dst, g) == 0 +is_across_boundary(cells, g) = false + +dist_to_boundary(cell, g::GreensFunction{<:Schur1DGreensSolver,1,Tuple{Int}}) = + only(cell) - only(g.boundaries) +dist_to_boundary((src, dst)::Pair, g::GreensFunction{<:Schur1DGreensSolver,1,Tuple{Int}}) = + dist_to_boundary(src, g), dist_to_boundary(dst, g) + +## Fast-paths + +# Surface-bulk semi-infinite: +# G_{1,1} = (ω*I - h0 - ΣR)⁻¹, G_{-1,-1} = (ω*I - h0 - ΣL)⁻¹ +# G_{N,1} = (G_{1,1}h₁)ᴺ⁻¹G_{1,1}, where G_{1,1} = (ω*I - h0 - ΣR)⁻¹ +# G_{-N,-1} = (G_{-1,-1}h₋₁)ᴺ⁻¹G_{-1,-1}, where G_{-1,-1} = (ω*I - h0 - ΣL)⁻¹ +function surface_fastpath(g, ω, (dsrc, ddst); sources = I) + dist = ddst - dsrc + Σ = dsrc > 0 ? g.solver(ω, Val{:R}) : g.solver(ω, Val{:L}) + h0 = g.solver.h0 + sourcemat = flat_padded_source_matrix(sources, g) + # in-place optimization of luG⁻¹ = lu(ω*I - h0 - Σ) + G⁻¹ = Σ + @. G⁻¹ = - h0 - Σ + shiftω!(G⁻¹, ω) + luG⁻¹ = lu!(G⁻¹) + if dist == 0 + G = ldiv!(luG⁻¹, sourcemat) + else + h = dist > 0 ? g.solver.hp : g.solver.hm + Gh = ldiv!(luG⁻¹, Matrix(h)) + G = Gh^(abs(dist)-1) * ldiv!(luG⁻¹, sourcemat) + end + return G end -function G_infinite!(gs, factors, G∞⁻¹, (src, dst)) - src´ = div(only(src), gs.maxdn, RoundToZero) - dst´ = div(only(dst), gs.maxdn, RoundToZero) - N = dst´ - src´ - GVᴺ = GVᴺ!(gs.temps.HH1, N, gs, factors) - G∞ = rdiv!(GVᴺ, G∞⁻¹) +# Local infinite: G∞_{n,n} = (ω*I - h0 - ΣR - ΣL)⁻¹ +function local_fastpath(g, ω; sources = I) + ΣR, ΣL = g.solver(ω, Val{:RL}) + h0 = g.solver.h0 + sourcemat = flat_padded_source_matrix(sources, g) + # in-place optimization of luG∞⁻¹ = lu(ω*I - h0 - ΣL - ΣR) + G∞⁻¹ = ΣR + @. G∞⁻¹ = - h0 - ΣR - ΣL + shiftω!(G∞⁻¹, ω) + luG∞⁻¹ = lu!(G∞⁻¹) + G∞ = ldiv!(luG∞⁻¹, sourcemat) return G∞ end -function G_semiinfinite!(gs, factors, G∞⁻¹, (src, dst)) - M = div(only(src), gs.maxdn, RoundToZero) - N = div(only(dst), gs.maxdn, RoundToZero) - if sign(N) != sign(M) - G∞ = fill!(gs.temps.HH1, 0) - else - GVᴺ⁻ᴹ = GVᴺ!(gs.temps.HH1, N-M, gs, factors) - GVᴺ = GVᴺ!(gs.temps.HH2, N, gs, factors) - GV⁻ᴹ = GVᴺ!(gs.temps.HH3, -M, gs, factors) - mul!(GVᴺ⁻ᴹ, GVᴺ, GV⁻ᴹ, -1, 1) # (GVᴺ⁻ᴹ - GVᴺGV⁻ᴹ) - G∞ = rdiv!(GVᴺ⁻ᴹ , G∞⁻¹) +function flat_padded_source_matrix(::UniformScaling, g) + n = flatsize(g, 1) + m = n * g.solver.unitcells + return Matrix(one(blockeltype(g)) * I, m, n) +end + +function flat_padded_source_matrix(sources::Vector{Int}, g) + n = length(sources) + m = size(g, 1) + T = blocktype(g) + sourcemat = zeros(T, m, n) + for (i, s) in enumerate(sources) + sourcemat[s, i] = one(T) end - return G∞ + flatmat = flatten(sourcemat, orbitalstructure(g)) + T´ = eltype(flatmat) + nu = g.solver.unitcells + paddedmat = nu == 1 ? flatmat : vcat(flatmat, zeros(T´, (nu-1)*size(flatmat, 1), size(flatmat, 2))) + return paddedmat end -shift_cells((src, dst), N0) = (only(src) - N0, only(dst) - N0) +## General paths -# G∞⁻¹ = G₀⁻¹ - V´GrV - VGlV´ with V´GrV = V'*χR*φRs⁻¹Ps and VGlV´ = iχA*φAs´⁻¹Ps´ -function inverse_G∞!(matrix, ω, gs, (λR, χR, iφRs, iλA, φA, iχAs´), onlyhalf = false) - G0⁻¹ = inverse_G0!(matrix, ω, gs) - V´GrV = mul!(gs.temps.Hs2, gs.V', mul!(gs.temps.Hs1, χR, iφRs)) - mul!(G0⁻¹, V´GrV, gs.Ps, -1, 1) - if onlyhalf - G∞⁻¹ = G0⁻¹ - else - VGlV´ = mul!(gs.temps.Hs2, gs.V, mul!(gs.temps.Hs1, φA, iχAs´)) - G∞⁻¹ = mul!(G0⁻¹, VGlV´, gs.Ps´, -1, 1) - end - return G∞⁻¹ +# Infinite: G∞_{N} = GVᴺ G∞_{0} +function (g::GreensFunction{<:Schur1DGreensSolver,1,Tuple{Missing}})(ω, ::Missing) + G∞⁻¹, GRh₊, GLh₋ = Gfactors(g.solver, ω) + luG∞⁻¹ = lu(G∞⁻¹) + return cells -> G_infinite(luG∞⁻¹, GRh₊, GLh₋, cells) end -# G₀⁻¹ = ω - H₀ -function inverse_G0!(G0⁻¹, ω, gs) - copy!(G0⁻¹, gs.minusH0) - for i in axes(G0⁻¹, 1) - G0⁻¹[i, i] += ω - end - return G0⁻¹ +function G_infinite(luG∞⁻¹, GRh₊, GLh₋, (src, dst)) + N = only(dst) - only(src) + N == 0 && return inv(luG∞⁻¹) + Ghᴺ = Gh_power(GRh₊, GLh₋, N) + G∞ = rdiv!(Ghᴺ, luG∞⁻¹) + return G∞ end -function GVᴺ!(GVᴺ, N, gs, (λR, χR, iφRs, iλA, φA, iχAs´)) - if N == 0 - copyto!(GVᴺ, I) - elseif N > 0 - χᴿλᴿᴺ´ = N == 1 ? χR : (gs.temps.Hs1 .= χR .* transpose(λR) .^ (N-1)) - mul!(GVᴺ, mul!(gs.temps.Hs2, χᴿλᴿᴺ´, iφRs), gs.Ps) - else - φᴬλᴬᴺ´ = N == -1 ? φA : (gs.temps.Hs1 .= φA .* transpose(iλA) .^ (-N-1)) - mul!(GVᴺ, mul!(gs.temps.Hs2, φᴬλᴬᴺ´, iχAs´), gs.Ps´) - end - return GVᴺ +# Semiinifinite: G_{N,M} = (Ghᴺ⁻ᴹ - GhᴺGh⁻ᴹ)G∞_{0} +function (g::GreensFunction{<:Schur1DGreensSolver,1,Tuple{Int}})(ω, ::Missing) + G∞⁻¹, GRh₊, GLh₋ = Gfactors(g.solver, ω) + return cells -> G_semiinfinite(G∞⁻¹, GRh₊, GLh₋, dist_to_boundary.(cells, Ref(g))) end +function G_semiinfinite(G∞⁻¹, Gh₊, Gh₋, (N, M)) + Ghᴺ⁻ᴹ = Gh_power(Gh₊, Gh₋, N-M) + Ghᴺ = Gh_power(Gh₊, Gh₋, N) + Gh⁻ᴹ = Gh_power(Gh₊, Gh₋, -M) + mul!(Ghᴺ⁻ᴹ, Ghᴺ, Gh⁻ᴹ, -1, 1) # (Ghᴺ⁻ᴹ - GhᴺGh⁻ᴹ) + # G∞ = rdiv!(Ghᴺ⁻ᴹ , luG∞⁻¹) # This is not defined in Julia (v1.7) yet + G∞ = Ghᴺ⁻ᴹ / G∞⁻¹ + return G∞ +end + +function Gfactors(solver::Schur1DGreensSolver, ω) + ΣR, ΣL = solver(ω, Val{:RL}) + A1 = ω*I - solver.h0 + GR⁻¹ = A1 - ΣR + GRh₊ = GR⁻¹ \ Matrix(solver.hp) + GL⁻¹ = A1 - ΣL + GLh₋ = GL⁻¹ \ Matrix(solver.hm) + G∞⁻¹ = GL⁻¹ - ΣR + return G∞⁻¹, GRh₊, GLh₋ +end + +Gh_power(Gh₊, Gh₋, N) = N == 0 ? one(Gh₊) : N > 0 ? Gh₊^N : Gh₋^-N + ####################################################################### # BandGreensSolver ####################################################################### diff --git a/src/hamiltonian.jl b/src/hamiltonian.jl index aa1b3d95..8fd6d5cd 100644 --- a/src/hamiltonian.jl +++ b/src/hamiltonian.jl @@ -154,6 +154,11 @@ Rebuild `s` by flattening its basis to have a scalar eltype. Curried form equivalent to `flatten(h)` of `h |> flatten` (included for consistency with the rest of the API). + flatten(x, o::OrbitalStructure) + +Flatten object x, if applicable, using the orbital structure o, as obtained from a +Hamiltonian `h` with `orbitalstructure(h)` + # Examples ```jldoctest @@ -186,21 +191,49 @@ Hamiltonian{<:Lattice} : Hamiltonian on a 2D Lattice in 2D space flatten() = h -> flatten(h) function flatten(h::Hamiltonian) - all(isequal(1), norbitals(h)) && return copy(h) + isflat(h) && return copy(h) harmonics´ = [flatten(har, h.orbstruct) for har in h.harmonics] lattice´ = flatten(h.lattice, h.orbstruct) orbs´ = (_ -> (:flat, )).(orbitals(h)) return Hamiltonian(lattice´, harmonics´, orbs´) end +# special method: if already flat, don't make a copy +maybeflatten(h, args...) = isflat(h) ? h : flatten(h, args...) + +isflat(h::Hamiltonian) = all(isequal(1), norbitals(h)) + flatten(h::HamiltonianHarmonic, orbstruct) = - HamiltonianHarmonic(h.dn, _flatten(h.h, orbstruct)) + HamiltonianHarmonic(h.dn, flatten(h.h, orbstruct)) -_flatten(src::AbstractArray{<:Number}, orbstruct) = src -_flatten(src::AbstractArray{T}, orbstruct, ::Type{T}) where {T<:Number} = src -_flatten(src::AbstractArray{T}, orbstruct, ::Type{T´}) where {T<:Number,T´<:Number} = T´.(src) +function flatten(lat::Lattice, orbstruct) + length(orbitals(orbstruct)) == nsublats(lat) || throw(ArgumentError("Mismatch between sublattices and orbitals")) + unitcell´ = flatten(lat.unitcell, orbstruct) + bravais´ = lat.bravais + lat´ = Lattice(bravais´, unitcell´) +end -function _flatten(src::SparseMatrixCSC{<:SMatrix{N,N,T}}, orbstruct, ::Type{T´} = T) where {N,T,T´} +function flatten(unitcell::Unitcell, orbstruct) # orbs::NTuple{S,Any}) where {S} + norbs = length.(orbitals(orbstruct)) + nsublats = length(sublats(orbstruct)) + offsets´ = orbstruct.flatoffsets + ns´ = last(offsets´) + sites´ = similar(unitcell.sites, ns´) + i = 1 + for sl in 1:nsublats, site in sitepositions(unitcell, sl), rep in 1:norbs[sl] + sites´[i] = site + i += 1 + end + names´ = unitcell.names + unitcell´ = Unitcell(sites´, names´, offsets´) + return unitcell´ +end + +flatten(src::AbstractArray{<:Number}, orbstruct) = src +flatten(src::AbstractArray{T}, orbstruct, ::Type{T}) where {T<:Number} = src +flatten(src::AbstractArray{T}, orbstruct, ::Type{T´}) where {T<:Number,T´<:Number} = T´.(src) + +function flatten(src::SparseMatrixCSC{<:SMatrix{N,N,T}}, orbstruct, ::Type{T´} = T) where {N,T,T´} norbs = length.(orbitals(orbstruct)) offsets´ = orbstruct.flatoffsets dim´ = last(offsets´) @@ -226,7 +259,7 @@ function _flatten(src::SparseMatrixCSC{<:SMatrix{N,N,T}}, orbstruct, ::Type{T´} return matrix end -function _flatten(src::StridedMatrix{<:SMatrix{N,N,T}}, orbstruct, ::Type{T´} = T) where {N,T,T´} +function flatten(src::StridedMatrix{<:SMatrix{N,N,T}}, orbstruct, ::Type{T´} = T) where {N,T,T´} norbs = length.(orbitals(orbstruct)) offsets´ = orbstruct.flatoffsets dim´ = last(offsets´) @@ -246,7 +279,7 @@ function _flatten(src::StridedMatrix{<:SMatrix{N,N,T}}, orbstruct, ::Type{T´} = end # for Subspace bases -function _flatten(src::StridedMatrix{<:SVector{N,T}}, orbstruct, ::Type{T´} = T) where {N,T,T´} +function flatten(src::StridedMatrix{<:SVector{N,T}}, orbstruct, ::Type{T´} = T) where {N,T,T´} norbs = length.(orbitals(orbstruct)) offsets´ = orbstruct.flatoffsets dim´ = last(offsets´) @@ -264,29 +297,6 @@ function _flatten(src::StridedMatrix{<:SVector{N,T}}, orbstruct, ::Type{T´} = T return matrix end -function flatten(lat::Lattice, orbstruct) - length(orbitals(orbstruct)) == nsublats(lat) || throw(ArgumentError("Mismatch between sublattices and orbitals")) - unitcell´ = flatten(lat.unitcell, orbstruct) - bravais´ = lat.bravais - lat´ = Lattice(bravais´, unitcell´) -end - -function flatten(unitcell::Unitcell, orbstruct) # orbs::NTuple{S,Any}) where {S} - norbs = length.(orbitals(orbstruct)) - nsublats = length(sublats(orbstruct)) - offsets´ = orbstruct.flatoffsets - ns´ = last(offsets´) - sites´ = similar(unitcell.sites, ns´) - i = 1 - for sl in 1:nsublats, site in sitepositions(unitcell, sl), rep in 1:norbs[sl] - sites´[i] = site - i += 1 - end - names´ = unitcell.names - unitcell´ = Unitcell(sites´, names´, offsets´) - return unitcell´ -end - function flatten_sparse_copy!(dst, src, o::OrbitalStructure) fill!(dst, zero(eltype(dst))) norbs = length.(orbitals(o)) @@ -308,6 +318,7 @@ function flatten_sparse_copy!(dst, src, o::OrbitalStructure) return dst end +# Specialized flattening copy! and muladd! function flatten_sparse_muladd!(dst, src, o::OrbitalStructure, α = I) norbs = length.(orbitals(o)) coloffset = 0 @@ -358,10 +369,9 @@ end ## unflatten ## """ - unflatten(v::AbstractArray, o::OrbitalStructure{T}) + unflatten(x, o::OrbitalStructure) -Rebuild `v` to have element type `T` and orbital structure `o` by performing the inverse of -`flatten(v)`. +Rebuild object `x` performing the inverse of `flatten(x)` or `flatten(x, o)`. # Examples ```jldoctest @@ -392,10 +402,22 @@ Subspace{0}: eigenenergy subspace on a 0D manifold # See also `flatten`, `orbitalstructure` """ -unflatten(vflat::AbstractVector, o::OrbitalStructure{T}) where {T} = - unflatten!(similar(vflat, T, dimh(o)), vflat, o) -unflatten(vflat::AbstractMatrix, o::OrbitalStructure{T}) where {T} = - unflatten!(similar(vflat, T, dimh(o), size(vflat, 2)), vflat, o) +unflatten + +# Pending implementation for Hamiltonian, Lattice, Unitcell + +unflatten_orbitals(v::AbstractVector, o::OrbitalStructure) = + isunflat_orbitals(v, o) ? copy(v) : unflatten!(similar(v, orbitaltype(o), dimh(o)), v, o) +unflatten_orbitals(m::AbstractMatrix, o::OrbitalStructure) = + isunflat_orbitals(m, o) ? copy(m) : unflatten!(similar(m, orbitaltype(o), dimh(o), size(m, 2)), m, o) +unflatten_blocks(m::AbstractMatrix, o::OrbitalStructure) = + isunflat_blocks(m, o) ? copy(m) : unflatten!(similar(m, blocktype(o), dimh(o), dimh(o)), m, o) + +maybeunflatten_orbitals(x, o) = isunflat_orbitals(x, o) ? x : unflatten_orbitals(x, o) +maybeunflatten_blocks(x, o) = isunflat_blocks(x, o) ? x : unflatten_blocks(x, o) + +isunflat_orbitals(m, o) = orbitaltype(o) == eltype(m) && size(m, 1) == dimh(o) +isunflat_blocks(m, o) = blocktype(o) == eltype(m) && size(m, 1) == dimh(o) function unflatten!(v::AbstractArray{T}, vflat::AbstractArray, o::OrbitalStructure) where {T} norbs = length.(orbitals(o)) @@ -404,51 +426,68 @@ function unflatten!(v::AbstractArray{T}, vflat::AbstractArray, o::OrbitalStructu check_unflatten_dst_dims(v, dimh(o)) check_unflatten_src_dims(vflat, dimflat) check_unflatten_eltypes(v, o) - for col in 1:size(v, 2) - row = 0 - for s in sublats(o) - N = norbs[s] - for i in flatoffsets[s]+1:N:flatoffsets[s+1] - row += 1 - v[row, col] = padright(view(vflat, i:i+N-1, col:col), T) + col = 0 + for scol in sublats(o) + M = colstride(norbs, scol, T) + for j in flatoffsets[scol]+1:M:flatoffsets[scol+1] + col +=1 + row = 0 + for srow in sublats(o) + N = rowstride(norbs, srow) + for i in flatoffsets[srow]+1:N:flatoffsets[srow+1] + row += 1 + val = view(vflat, i:i+N-1, j:j+M-1) + v[row, col] = padtotype(val, T) + end end end end return v end -## unflatten_or_reinterpret: call unflatten but only if we cannot unflatten via reinterpret -unflatten_or_reinterpret(vflat, ::Missing) = vflat +rowstride(norbs, s) = norbs[s] +colstride(norbs, s, ::Type{<:SVector}) = 1 +colstride(norbs, s, ::Type{<:SMatrix}) = rowstride(norbs, s) + +# dest v should be a vector or matrix such that H*v is possible +check_unflatten_dst_dims(v, dimh) = + size(v, 1) == dimh || + throw(ArgumentError("Dimension of destination array is inconsistent with orbital structure")) + +# dest v should be a square matrix like the Hamiltonian +check_unflatten_dst_dims(v::AbstractArray{<:SMatrix}, dimh) = + size(v, 1) == dimh && size(v, 2) == dimh || + throw(ArgumentError("Dimension of destination array is inconsistent with orbital structure")) + +check_unflatten_src_dims(vflat, dimflat) = + size(vflat, 1) == dimflat || + throw(ArgumentError("Dimension of source array is inconsistent with orbital structure")) + +check_unflatten_eltypes(::AbstractArray{T}, o::OrbitalStructure) where {T} = + T === orbitaltype(o) || T === blocktype(o) || + throw(ArgumentError("Eltype of desination array is inconsistent with orbital structure")) + +## unflatten_orbitals_or_reinterpret: call unflatten_orbitals but only if we cannot unflatten via reinterpret +unflatten_orbitals_or_reinterpret(vflat, ::Missing) = vflat # source is already of the correct orbitaltype(h) -function unflatten_or_reinterpret(vflat::AbstractArray{T}, o::OrbitalStructure{T}) where {T} +function unflatten_orbitals_or_reinterpret(vflat::AbstractArray{T}, o::OrbitalStructure{T}) where {T} check_unflatten_dst_dims(vflat, dimh(o)) return vflat end -function unflatten_or_reinterpret(vflat::AbstractArray{<:Number}, o::OrbitalStructure{<:Number}) +function unflatten_orbitals_or_reinterpret(vflat::AbstractArray{<:Number}, o::OrbitalStructure{<:Number}) check_unflatten_dst_dims(vflat, dimh(o)) return vflat end # source can be reinterpreted, because the number of orbitals is the same M for all N sublattices -unflatten_or_reinterpret(vflat::AbstractArray{T}, o::OrbitalStructure{S,N,<:NTuple{N,NTuple{M}}}) where {T,S,N,M} = +unflatten_orbitals_or_reinterpret(vflat::AbstractArray{T}, o::OrbitalStructure{S,N,<:NTuple{N,NTuple{M}}}) where {T,S,N,M} = reinterpret(SVector{M,T}, vflat) -# otherwise call unflatten -unflatten_or_reinterpret(vflat, o) = unflatten(vflat, o) +# otherwise call unflatten_orbitals +unflatten_orbitals_or_reinterpret(vflat, o) = unflatten_orbitals(vflat, o) -check_unflatten_dst_dims(v, dimh) = - size(v, 1) == dimh || - throw(ArgumentError("Dimension of destination array is inconsistent with Hamiltonian")) - -check_unflatten_src_dims(vflat, dimflat) = - size(vflat, 1) == dimflat || - throw(ArgumentError("Dimension of source array is inconsistent with Hamiltonian")) - -check_unflatten_eltypes(::AbstractArray{T}, ::OrbitalStructure{S}) where {T,S} = - T === S || throw(ArgumentError("Eltype of desination array is inconsistent with Hamiltonian")) - -valdim(::Type{<:Number}) = Val(1) -valdim(::Type{S}) where {N,S<:SVector{N}} = Val(N) +# valdim(::Type{<:Number}) = Val(1) +# valdim(::Type{S}) where {N,S<:SVector{N}} = Val(N) ####################################################################### # similarmatrix @@ -506,7 +545,7 @@ _similarmatrix(::Type{A}, ::Type{A´}, h) where {T<:Number,A<:AbstractSparseMatr _similarmatrix(::Type{A}, ::Type{A´}, h) where {N,T<:SMatrix{N,N},A<:AbstractSparseMatrix{T},T´<:SMatrix{N,N},A´<:AbstractSparseMatrix{T´}} = similar_merged(h.harmonics, T) _similarmatrix(::Type{A}, ::Type{A´}, h) where {N,T<:Number,A<:AbstractSparseMatrix{T},T´<:SMatrix{N,N},A´<:AbstractSparseMatrix{T´}} = - _flatten(similar_merged(h.harmonics), h.orbstruct, T) + flatten(similar_merged(h.harmonics), h.orbstruct, T) _similarmatrix(::Type{A}, ::Type{A´}, h) where {T<:Number,A<:Matrix{T},T´<:Number,A´<:AbstractMatrix{T´}} = similar(A, size(h)) _similarmatrix(::Type{A}, ::Type{A´}, h) where {N,T<:SMatrix{N,N},A<:Matrix{T},T´<:SMatrix{N,N},A´<:AbstractMatrix{T´}} = @@ -744,6 +783,7 @@ _blocktype(::Type{S}) where {N,Tv,S<:SVector{N,Tv}} = SMatrix{N,N,Tv,N*N} _blocktype(::Type{S}) where {S<:Number} = S blocktype(h::Hamiltonian{LA,L,M}) where {LA,L,M} = M +blocktype(o::OrbitalStructure) = _blocktype(orbitaltype(o)) promote_blocktype(hs::Hamiltonian...) = promote_blocktype(blocktype.(hs)...) promote_blocktype(s1::Type, s2::Type, ss::Type...) = @@ -758,9 +798,7 @@ blockdim(::Type{S}) where {N,S<:SMatrix{N,N}} = N blockdim(::Type{T}) where {T<:Number} = 1 orbitaltype(h::Hamiltonian) = orbitaltype(h.orbstruct) - -orbitaltype(o::OrbitalStructure{T}) where {T} = T - +orbitaltype(o::OrbitalStructure) = o.orbtype orbitaltype(::Type{M}) where {M<:Number} = M orbitaltype(::Type{S}) where {N,T,S<:SMatrix{N,N,T}} = SVector{N,T} diff --git a/src/lattice.jl b/src/lattice.jl index b9c2b4ca..efc84c19 100644 --- a/src/lattice.jl +++ b/src/lattice.jl @@ -676,8 +676,8 @@ supercell_center(lat::AbstractLattice{E,L,T}) where {E,L,T} = # supplements supercell with most orthogonal bravais axes function extended_supercell(bravais, supercell::SMatrix{L,L´}) where {L,L´} L == L´ && return supercell - bravais_new_norm = normalize_cols(bravais * supercell) - bravais_norm = normalize_cols(bravais) + bravais_new_norm = normalize_columns(bravais * supercell) + bravais_norm = normalize_columns(bravais) # νnorm are the L projections of old bravais on new bravais axis subspace ν = bravais_norm' * bravais_new_norm # L×L´ νnorm = SVector(ntuple(row -> norm(ν[row,:]), Val(L))) @@ -687,10 +687,6 @@ function extended_supercell(bravais, supercell::SMatrix{L,L´}) where {L,L´} return ext_supercell end -normalize_cols(s::SMatrix{L,0}) where {L} = s -normalize_cols(s::SMatrix{0,L}) where {L} = s -normalize_cols(s) = mapslices(v -> v/norm(v), s, dims = 1) - function _superlat(lat, scmatrix, pararegion, selector_perp, seed) br = bravais(lat) rsel = resolve(selector_perp, lat) diff --git a/src/plot_vegalite.jl b/src/plot_vegalite.jl index fc651b5a..069596fa 100644 --- a/src/plot_vegalite.jl +++ b/src/plot_vegalite.jl @@ -300,9 +300,9 @@ sanitize_colorrange(r::Tuple, table) = [first(r), last(r)] needslegend(x::Number) = nothing needslegend(x) = true -unflatten_or_reinterpret_or_missing(psi::Missing, h) = missing -unflatten_or_reinterpret_or_missing(psi::AbstractArray, h) = unflatten_or_reinterpret(psi, h.orbstruct) -unflatten_or_reinterpret_or_missing(s::Subspace, h) = unflatten_or_reinterpret(s.basis, h.orbstruct) +unflatten_orbitals_or_reinterpret_or_missing(psi::Missing, h) = missing +unflatten_orbitals_or_reinterpret_or_missing(psi::AbstractArray, h) = unflatten_orbitals_or_reinterpret(psi, h.orbstruct) +unflatten_orbitals_or_reinterpret_or_missing(s::Subspace, h) = unflatten_orbitals_or_reinterpret(s.basis, h.orbstruct) checkdims_psi(h, psi) = size(h, 2) == size(psi, 1) || throw(ArgumentError("The eigenstate length $(size(psi,1)) must match the Hamiltonian dimension $(size(h, 2))")) @@ -326,7 +326,7 @@ end function linkstable!(table, h, psi, cell = missing; axes, digits, onlyonecell, plotsites, plotlinks, sitesize, siteopacity, sitecolor, linksize, linkopacity, linkcolor) where {LA,L} - psi´ = unflatten_or_reinterpret_or_missing(psi, h) + psi´ = unflatten_orbitals_or_reinterpret_or_missing(psi, h) checkdims_psi(h, psi´) d = (; axes = axes, digits = digits, sitesize_shader = site_shader(sitesize, h, psi´), diff --git a/src/tools.jl b/src/tools.jl index 3d5c49c1..89edae5f 100644 --- a/src/tools.jl +++ b/src/tools.jl @@ -106,22 +106,31 @@ padright(t::NTuple{N´,Any}, x, ::Val{N}) where {N´,N} = ntuple(i -> i > N´ ? padright(t::NTuple{N´,Any}, ::Val{N}) where {N´,N} = ntuple(i -> i > N´ ? 0 : t[i], Val(N)) padright(v, ::Type{<:Number}) = first(v) -padright(v, ::Type{S}) where {E,T,S<:SVector{E,T}} = padright(v, zero(T), S) +padright(v, ::Type{S}) where {T,S<:SVector{<:Any,T}} = padright(v, zero(T), S) padright(v, x::T, ::Type{S}) where {E,T,S<:SVector{E,T}} = SVector{E,T}(ntuple(i -> i > length(v) ? x : convert(T, v[i]), Val(E))) # Pad element type to a "larger" type -@inline padtotype(s::SMatrix{E,L}, ::Type{S}) where {E,L,E2,L2,S<:SMatrix{E2,L2}} = +padtotype(s::SMatrix{E,L}, ::Type{S}) where {E,L,E2,L2,T,S<:SMatrix{E2,L2,T}} = S(SMatrix{E2,E}(I) * s * SMatrix{L,L2}(I)) -@inline padtotype(s::StaticVector, ::Type{S}) where {N,T,S<:SVector{N,T}} = +padtotype(s::StaticVector, ::Type{S}) where {N,T,S<:SVector{N,T}} = padright(T.(s), Val(N)) -@inline padtotype(x::Number, ::Type{S}) where {E,L,S<:SMatrix{E,L}} = +padtotype(x::Number, ::Type{S}) where {E,L,S<:SMatrix{E,L}} = S(x * (SMatrix{E,1}(I) * SMatrix{1,L}(I))) -@inline padtotype(s::Number, ::Type{S}) where {N,T,S<:SVector{N,T}} = +padtotype(s::Number, ::Type{S}) where {N,T,S<:SVector{N,T}} = padright(SA[T(s)], Val(N)) -@inline padtotype(x::Number, ::Type{T}) where {T<:Number} = T(x) -@inline padtotype(u::UniformScaling, ::Type{T}) where {T<:Number} = T(u.λ) -@inline padtotype(u::UniformScaling, ::Type{S}) where {S<:SMatrix} = S(u) +padtotype(x::Number, ::Type{T}) where {T<:Number} = T(x) +padtotype(u::UniformScaling, ::Type{T}) where {T<:Number} = T(u.λ) +padtotype(u::UniformScaling, ::Type{S}) where {S<:SMatrix} = S(u) +padtotype(a::AbstractVector, t::SVector) = padright(a, t) + +function padtotype(a::AbstractMatrix, ::Type{S}) where {N,M,T,S<:SMatrix{N,M,T}} + t = ntuple(Val(N*M)) do i + n, m = mod1(i, N), fld1(i, N) + @inbounds n > size(a, 1) || m > size(a, 2) ? zero(T) : T(a[n,m]) + end + return S(t) +end ## Work around BUG: -SVector{0,Int}() isa SVector{0,Union{}} negative(s::SVector{L,<:Number}) where {L} = -s @@ -187,9 +196,24 @@ function ispositive(ndist) return result end -chop(x::T, x0 = one(T)) where {T<:Real} = ifelse(abs(x) < √eps(T(x0)), zero(T), x) +chop(x::T, x0 = one(T)) where {T<:Real} = ifelse(abs2(x) < eps(T(x0)), zero(T), x) chop(x::C, x0 = one(R)) where {R<:Real,C<:Complex{R}} = chop(real(x), x0) + im*chop(imag(x), x0) +# function chop!(A::AbstractArray{T}, atol = default_tol(T)) where {T} +# for (i, a) in enumerate(A) +# # if abs(a) < atol +# # A[i] = zero(T) +# if !iszero(a) #&& (abs(real(a)) < atol || abs(imag(a)) < atol) +# A[i] = chop(a) +# elseif abs(a) > 1/atol || isnan(a) +# A[i] = T(Inf) +# end +# end +# return A +# end + +default_tol(::Type{T}) where {T} = sqrt(eps(real(float(T)))) + function unique_sorted_approx!(v::AbstractVector{T}) where {T} i = 1 xprev = first(v) @@ -214,6 +238,10 @@ function normalize_columns!(kmat::AbstractMatrix, cols) return kmat end +normalize_columns(s::SMatrix{L,0}) where {L} = s +normalize_columns(s::SMatrix{0,L}) where {L} = s +normalize_columns(s) = mapslices(v -> v/norm(v), s, dims = 1) + eltypevec(::AbstractMatrix{T}) where {T<:Number} = T eltypevec(::AbstractMatrix{S}) where {N,T<:Number,S<:SMatrix{N,N,T}} = SVector{N,T} @@ -282,6 +310,14 @@ permutations!(p, ::Tuple{}, s2) = push!(p, s2) delete(t::NTuple{N,Any}, i) where {N} = ntuple(j -> j < i ? t[j] : t[j+1], Val(N-1)) +function one!(m::StridedMatrix{T}) where {T} + @inbounds for c in CartesianIndices(m) + m[c] = ifelse(c[1] == c[2], one(T), zero(T)) + end + return m +end + + ###################################################################### # convert a matrix/number block to a matrix/inlinematrix string ###################################################################### @@ -349,4 +385,82 @@ end rclamp(r1::UnitRange, r2::UnitRange) = isempty(r1) ? r1 : clamp(minimum(r1), extrema(r2)...):clamp(maximum(r1), extrema(r2)...) iclamp(minmax, r::Missing) = minmax -iclamp((x1, x2), (xmin, xmax)) = (max(x1, xmin), min(x2, xmax)) \ No newline at end of file +iclamp((x1, x2), (xmin, xmax)) = (max(x1, xmin), min(x2, xmax)) + +############################################################################################ +# QR utils +############################################################################################ +# Sparse QR from SparseSuite is also pivoted +pqr(a::SparseMatrixCSC) = qr(a) +pqr(a::Adjoint{T,<:SparseMatrixCSC}) where {T} = qr(copy(a)) +pqr(a::Transpose{T,<:SparseMatrixCSC}) where {T} = qr(copy(a)) +pqr(a) = qr!(copy(a), Val(true)) + +getQ(qr::Factorization, cols = :) = qr.Q * Idense(size(qr, 1), cols) +getQ(qr::SuiteSparse.SPQR.QRSparse, cols = :) = Isparse(size(qr, 1), :, qr.prow) * sparse(qr.Q * Idense(size(qr, 1), cols)) +getQ_dense(qr::SuiteSparse.SPQR.QRSparse, cols = :) = Isparse(size(qr, 1), :, qr.prow) * (qr.Q * Idense(size(qr, 1), cols)) + +getQ´(qr::Factorization, cols = :) = qr.Q' * Idense(size(qr, 1), cols) +getQ´(qr::SuiteSparse.SPQR.QRSparse, cols = :) = sparse((qr.Q * Idense(size(qr, 1), cols))') * Isparse(size(qr,1), qr.prow, :) + +getRP´(qr::Factorization) = qr.R * qr.P' +getRP´(qr::SuiteSparse.SPQR.QRSparse) = qr.R * Isparse(size(qr, 2), qr.pcol, :) + +getPR´(qr::Factorization) = qr.P * qr.R' +getPR´(qr::SuiteSparse.SPQR.QRSparse) = Isparse(size(qr, 2), :, qr.pcol) * qr.R' + +Idense(n, ::Colon) = Matrix(I, n, n) + +function Idense(n, cols) + m = zeros(Bool, n, length(cols)) + for (j, col) in enumerate(cols) + m[col, j] = true + end + return m +end + +# Equivalent to I(n)[rows, cols], but faster +Isparse(n, rows, cols) = Isparse(inds(rows, n), inds(cols, n)) + +function Isparse(rows, cols) + rowval = Int[] + nzval = Bool[] + colptr = Vector{Int}(undef, length(cols) + 1) + colptr[1] = 1 + for (j, col) in enumerate(cols) + push_rows!(rowval, nzval, rows, col) + colptr[j+1] = length(rowval) + 1 + end + return SparseMatrixCSC(length(rows), length(cols), colptr, rowval, nzval) +end + +inds(::Colon, n) = 1:n +inds(is, n) = is + +function push_rows!(rowval, nzval, rows::AbstractUnitRange, col) + if col in rows + push!(rowval, 1 + rows[1 + col - first(rows)] - first(rows)) + push!(nzval, true) + end + return nothing +end + +function push_rows!(rowval, nzval, rows, col) + for (j, row) in enumerate(rows) + if row == col + push!(rowval, j) + push!(nzval, true) + end + end + return nothing +end + +function nonzero_rows(m::AbstractMatrix{T}, atol = default_tol(T)) where {T} + nrows = size(m, 1) + zerorows = 0 + for j in 0:nrows-1 + maximum(abs, view(m, nrows - j, :)) > atol && break + zerorows += 1 + end + return nrows - zerorows +end \ No newline at end of file diff --git a/test/test_bandstructure.jl b/test/test_bandstructure.jl index d7149a7c..b6122922 100644 --- a/test/test_bandstructure.jl +++ b/test/test_bandstructure.jl @@ -112,7 +112,7 @@ end @testset "unflatten" begin h = LatticePresets.honeycomb() |> hamiltonian(onsite(2I) + hopping(I, range = 1), orbitals = (Val(2), Val(1))) |> unitcell(2) |> unitcell sp = spectrum(h).states[:,1] - sp´ = Quantica.unflatten_or_reinterpret(sp, h.orbstruct) + sp´ = Quantica.unflatten_orbitals_or_reinterpret(sp, h.orbstruct) l = size(h, 1) @test length(sp) == 1.5 * l @test length(sp´) == l @@ -131,7 +131,7 @@ end h = LatticePresets.honeycomb() |> hamiltonian(onsite(2I) + hopping(I, range = 1), orbitals = Val(2)) |> unitcell(2) |> unitcell sp = spectrum(h).states[:,1] - sp´ = Quantica.unflatten_or_reinterpret(sp, h.orbstruct) + sp´ = Quantica.unflatten_orbitals_or_reinterpret(sp, h.orbstruct) l = size(h, 1) @test length(sp) == 2 * l @test length(sp´) == l @@ -139,7 +139,7 @@ end h = LatticePresets.honeycomb() |> hamiltonian(onsite(2I) + hopping(I, range = 1), orbitals = Val(2)) |> unitcell(2) |> unitcell sp = spectrum(h).states[:,1] - sp´ = Quantica.unflatten_or_reinterpret(sp, h.orbstruct) + sp´ = Quantica.unflatten_orbitals_or_reinterpret(sp, h.orbstruct) @test sp === sp end diff --git a/test/test_greens.jl b/test/test_greens.jl index e2b6ea7d..aa81e8a2 100644 --- a/test/test_greens.jl +++ b/test/test_greens.jl @@ -1,51 +1,62 @@ using LinearAlgebra: tr -# @testset "greens bandstructure" begin -# g = LatticePresets.square() |> hamiltonian(hopping(-1)) |> unitcell((2,0), (0,1)) |> -# greens(bandstructure(subticks = 17)) -# @test g(0.2) ≈ transpose(g(0.2)) -# @test imag(tr(g(1))) < 0 -# @test Quantica.chop(imag(tr(g(5)))) == 0 -# end +@testset "greens singleshot selfenergy" begin + h1 = LP.honeycomb() |> hamiltonian(hopping(1)) |> unitcell((1,-1), region = r->abs(r[2])<3) + h2 = LP.honeycomb() |> hamiltonian(hopping(1) + onsite(0.02, sublats = :A) - onsite(0.02, sublats = :B)) |> unitcell((4,4), region = r -> -5 < r[1] < 5) + res(g, ω) = (s = g.solver; Σ = s(ω); sum(abs.(Σ - s.hm * (((ω * I - s.h0) - Σ) \ Matrix(s.hp))))) + for h in (h1, h2) + g = greens(h, Schur1D(), boundaries = (0,)) + residual = maximum(res(g, ω) for ω in range(-3.2, 3.2, length = 101)) + g = greens(h, Schur1D(deflation = nothing), boundaries = (0,)) + residual0 = maximum(res(g, ω) for ω in range(-3.2, 3.2, length = 101)) + @test residual < 2*residual0 + end +end + +@testset "greens multiorbital" begin + h = LP.honeycomb() |> hamiltonian(hopping(SA[0 1; 1 0], range = 2), orbitals = Val(2)) |> unitcell((2,-1), region = r->abs(r[2])<3) + g = greens(h, Schur1D(), boundaries = (0,)) + g0 = greens(flatten(h), Schur1D(), boundaries = (0,)) + @test tr(tr(g(0.2))) ≈ tr(g0(0.2)) +end @testset "greens singleshot spectra" begin # Broken until https://github.com/Reference-LAPACK/lapack/issues/477 comes to julia - # h = LP.honeycomb() |> hamiltonian(hopping(1)) |> unitcell((1,-1), region = r->abs(r[2])<3) - # g = greens(h, SingleShot1D()) - # @test_throws ArgumentError g(0) - # dos = [-imag(tr(g(w + 1.0e-6im))) for w in range(-3.2, 3.2, length = 1001)] - # @test all(x -> Quantica.chop(x) >= 0, dos) + h = LP.honeycomb() |> hamiltonian(hopping(1)) |> unitcell((1,-1), region = r->abs(r[2])<3) + g = greens(h, Schur1D()) + dos = [-imag(tr(g(w))) for w in range(-3.2, 3.2, length = 101)] + @test all(x -> Quantica.chop(x) >= 0, dos) # Broken until https://github.com/Reference-LAPACK/lapack/issues/477 comes to julia - # h = LP.honeycomb() |> hamiltonian(hopping(1, range = 2)) |> unitcell((1,-1), region = r->abs(r[2])<3) - # g = greens(h, SingleShot1D()) - # dos = [-imag(tr(g(w))) for w in range(-6, 6, length = 1001)] - # @test all(x -> Quantica.chop(x) >= 0, dos) + h = LP.honeycomb() |> hamiltonian(hopping(1, range = 2)) |> unitcell((1,-1), region = r->abs(r[2])<3) + g = greens(h, Schur1D()) + dos = [-imag(tr(g(w))) for w in range(-6, 6, length = 101)] + @test all(x -> Quantica.chop(x) >= 0, dos) # Broken until https://github.com/Reference-LAPACK/lapack/issues/477 comes to julia - # h = LatticePresets.honeycomb() |> hamiltonian(hopping(1)) |> unitcell((1,-1),(2,2)) |> wrap(2) - # g = greens(h, SingleShot1D()) - # dos = [-imag(tr(g(w))) for w in range(-3.2, 3.2, length = 1001)] - # @test all(x -> Quantica.chop(x) >= 0, dos) + h = LatticePresets.honeycomb() |> hamiltonian(hopping(1)) |> unitcell((1,-1),(2,2)) |> wrap(2) + g = greens(h, Schur1D()) + dos = [-imag(tr(g(w))) for w in range(-3.2, 3.2, length = 101)] + @test all(x -> Quantica.chop(x) >= 0, dos) # Broken until https://github.com/Reference-LAPACK/lapack/issues/477 comes to julia - # h = LatticePresets.honeycomb() |> hamiltonian(hopping((r, dr) -> 1/(dr'*dr), range = 1)) |> unitcell((1,-1),(2,2)) |> wrap(2) - # g = greens(h, SingleShot1D()) - # dos = [-imag(tr(g(w))) for w in range(-3.2, 3.2, length = 1001)] - # @test all(x -> Quantica.chop(x) >= 0, dos) + h = LatticePresets.honeycomb() |> hamiltonian(hopping((r, dr) -> 1/(dr'*dr), range = 1)) |> unitcell((1,-1),(2,2)) |> wrap(2) + g = greens(h, Schur1D()) + dos = [-imag(tr(g(w))) for w in range(-3.2, 3.2, length = 101)] + @test all(x -> Quantica.chop(x) >= 0, dos) h = LP.honeycomb() |> hamiltonian(hopping(1)) |> unitcell((2,-2),(3,3)) |> unitcell(1, 1, indices = not(1)) |> wrap(2) - g = greens(h, SingleShot1D(direct = true)) + g = greens(h, Schur1D()) # @test_broken Quantica.chop(imag(tr(g(0.2)))) <= 0 - g = greens(h, SingleShot1D()) + g = greens(h, Schur1D()) @test Quantica.chop(imag(tr(g(0.2)))) <= 0 - dos = [-imag(tr(g(w + 1.0e-6im))) for w in range(-3.2, 3.2, length = 1001)] + dos = [-imag(tr(g(w))) for w in range(-3.2, 3.2, length = 101)] @test all(x -> Quantica.chop(x) >= 0, dos) end @testset "greens singleshot spatial" begin h = LatticePresets.honeycomb() |> hamiltonian(hopping(1)) |> unitcell((1,-1),(3,3)) |> unitcell((1,0)) - g = greens(h, SingleShot1D(), boundaries = (0,)) + g = greens(h, Schur1D(), boundaries = (0,)) g0 = g(0.01, missing) ldos = [-imag(tr(g0(n=>n))) for n in 1:200] @test all(x -> Quantica.chop(x) >= 0, ldos) From 647f56123f90ee0c26ffebe314d22d2c3e504518 Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Fri, 14 May 2021 19:18:38 +0200 Subject: [PATCH 2/4] fix unflatten tests --- src/hamiltonian.jl | 33 +++++++++++++++++++++++++++++---- src/tools.jl | 2 +- 2 files changed, 30 insertions(+), 5 deletions(-) diff --git a/src/hamiltonian.jl b/src/hamiltonian.jl index 8fd6d5cd..ff8c47bb 100644 --- a/src/hamiltonian.jl +++ b/src/hamiltonian.jl @@ -419,13 +419,19 @@ maybeunflatten_blocks(x, o) = isunflat_blocks(x, o) ? x : unflatten_blocks(x, o) isunflat_orbitals(m, o) = orbitaltype(o) == eltype(m) && size(m, 1) == dimh(o) isunflat_blocks(m, o) = blocktype(o) == eltype(m) && size(m, 1) == dimh(o) -function unflatten!(v::AbstractArray{T}, vflat::AbstractArray, o::OrbitalStructure) where {T} - norbs = length.(orbitals(o)) - flatoffsets = o.flatoffsets - dimflat = last(flatoffsets) +function unflatten!(v::AbstractArray, vflat::AbstractArray, o::OrbitalStructure) + dimflat = last(o.flatoffsets) check_unflatten_dst_dims(v, dimh(o)) check_unflatten_src_dims(vflat, dimflat) check_unflatten_eltypes(v, o) + v = _unflatten!(v, vflat, o) + return v +end + +# unflatten into SMatrix blocks +function _unflatten!(v::AbstractArray{T}, vflat::AbstractArray, o::OrbitalStructure) where {T<:SMatrix} + norbs = length.(orbitals(o)) + flatoffsets = o.flatoffsets col = 0 for scol in sublats(o) M = colstride(norbs, scol, T) @@ -445,6 +451,25 @@ function unflatten!(v::AbstractArray{T}, vflat::AbstractArray, o::OrbitalStructu return v end +# unflatten into SVector orbitals +function _unflatten!(v::AbstractArray{T}, vflat::AbstractArray, o::OrbitalStructure) where {T<:SVector} + norbs = length.(orbitals(o)) + flatoffsets = o.flatoffsets + col = 0 + for col in 1:size(v, 2) + row = 0 + for srow in sublats(o) + N = rowstride(norbs, srow) + for i in flatoffsets[srow]+1:N:flatoffsets[srow+1] + row += 1 + val = view(vflat, i:i+N-1, col:col) + v[row, col] = padtotype(val, T) + end + end + end + return v +end + rowstride(norbs, s) = norbs[s] colstride(norbs, s, ::Type{<:SVector}) = 1 colstride(norbs, s, ::Type{<:SMatrix}) = rowstride(norbs, s) diff --git a/src/tools.jl b/src/tools.jl index 89edae5f..f599f6b1 100644 --- a/src/tools.jl +++ b/src/tools.jl @@ -122,7 +122,7 @@ padtotype(s::Number, ::Type{S}) where {N,T,S<:SVector{N,T}} = padtotype(x::Number, ::Type{T}) where {T<:Number} = T(x) padtotype(u::UniformScaling, ::Type{T}) where {T<:Number} = T(u.λ) padtotype(u::UniformScaling, ::Type{S}) where {S<:SMatrix} = S(u) -padtotype(a::AbstractVector, t::SVector) = padright(a, t) +padtotype(a::AbstractArray, t::Type{<:SVector}) = padright(a, t) function padtotype(a::AbstractMatrix, ::Type{S}) where {N,M,T,S<:SMatrix{N,M,T}} t = ntuple(Val(N*M)) do i From db363c9d800eec81ac0ddbf77dcdf703f270a7d2 Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Fri, 14 May 2021 19:25:26 +0200 Subject: [PATCH 3/4] bypass LAPACK bug in tests --- test/test_greens.jl | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/test/test_greens.jl b/test/test_greens.jl index aa81e8a2..28ec48e0 100644 --- a/test/test_greens.jl +++ b/test/test_greens.jl @@ -21,33 +21,31 @@ end end @testset "greens singleshot spectra" begin - # Broken until https://github.com/Reference-LAPACK/lapack/issues/477 comes to julia h = LP.honeycomb() |> hamiltonian(hopping(1)) |> unitcell((1,-1), region = r->abs(r[2])<3) g = greens(h, Schur1D()) dos = [-imag(tr(g(w))) for w in range(-3.2, 3.2, length = 101)] @test all(x -> Quantica.chop(x) >= 0, dos) - # Broken until https://github.com/Reference-LAPACK/lapack/issues/477 comes to julia h = LP.honeycomb() |> hamiltonian(hopping(1, range = 2)) |> unitcell((1,-1), region = r->abs(r[2])<3) g = greens(h, Schur1D()) dos = [-imag(tr(g(w))) for w in range(-6, 6, length = 101)] @test all(x -> Quantica.chop(x) >= 0, dos) - # Broken until https://github.com/Reference-LAPACK/lapack/issues/477 comes to julia h = LatticePresets.honeycomb() |> hamiltonian(hopping(1)) |> unitcell((1,-1),(2,2)) |> wrap(2) g = greens(h, Schur1D()) dos = [-imag(tr(g(w))) for w in range(-3.2, 3.2, length = 101)] @test all(x -> Quantica.chop(x) >= 0, dos) # Broken until https://github.com/Reference-LAPACK/lapack/issues/477 comes to julia - h = LatticePresets.honeycomb() |> hamiltonian(hopping((r, dr) -> 1/(dr'*dr), range = 1)) |> unitcell((1,-1),(2,2)) |> wrap(2) - g = greens(h, Schur1D()) - dos = [-imag(tr(g(w))) for w in range(-3.2, 3.2, length = 101)] - @test all(x -> Quantica.chop(x) >= 0, dos) + if VERSION >= v"1.7.0-DEV" + h = LatticePresets.honeycomb() |> hamiltonian(hopping((r, dr) -> 1/(dr'*dr), range = 1)) |> unitcell((1,-1),(2,2)) |> wrap(2) + g = greens(h, Schur1D()) + dos = [-imag(tr(g(w))) for w in range(-3.2, 3.2, length = 101)] + @test all(x -> Quantica.chop(x) >= 0, dos) + end h = LP.honeycomb() |> hamiltonian(hopping(1)) |> unitcell((2,-2),(3,3)) |> unitcell(1, 1, indices = not(1)) |> wrap(2) g = greens(h, Schur1D()) - # @test_broken Quantica.chop(imag(tr(g(0.2)))) <= 0 g = greens(h, Schur1D()) @test Quantica.chop(imag(tr(g(0.2)))) <= 0 dos = [-imag(tr(g(w))) for w in range(-3.2, 3.2, length = 101)] From 4825d77d017af183bd198d1fb2384c9ff0ed1b91 Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Fri, 14 May 2021 19:43:11 +0200 Subject: [PATCH 4/4] more LAPACK bug workarounds --- test/test_greens.jl | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/test/test_greens.jl b/test/test_greens.jl index 28ec48e0..75a2cbf3 100644 --- a/test/test_greens.jl +++ b/test/test_greens.jl @@ -21,23 +21,23 @@ end end @testset "greens singleshot spectra" begin - h = LP.honeycomb() |> hamiltonian(hopping(1)) |> unitcell((1,-1), region = r->abs(r[2])<3) - g = greens(h, Schur1D()) - dos = [-imag(tr(g(w))) for w in range(-3.2, 3.2, length = 101)] - @test all(x -> Quantica.chop(x) >= 0, dos) + # Broken until https://github.com/Reference-LAPACK/lapack/issues/477 comes to julia + if VERSION >= v"1.7.0-DEV" + h = LP.honeycomb() |> hamiltonian(hopping(1)) |> unitcell((1,-1), region = r->abs(r[2])<3) + g = greens(h, Schur1D()) + dos = [-imag(tr(g(w))) for w in range(-3.2, 3.2, length = 101)] + @test all(x -> Quantica.chop(x) >= 0, dos) - h = LP.honeycomb() |> hamiltonian(hopping(1, range = 2)) |> unitcell((1,-1), region = r->abs(r[2])<3) - g = greens(h, Schur1D()) - dos = [-imag(tr(g(w))) for w in range(-6, 6, length = 101)] - @test all(x -> Quantica.chop(x) >= 0, dos) + h = LP.honeycomb() |> hamiltonian(hopping(1, range = 2)) |> unitcell((1,-1), region = r->abs(r[2])<3) + g = greens(h, Schur1D()) + dos = [-imag(tr(g(w))) for w in range(-6, 6, length = 101)] + @test all(x -> Quantica.chop(x) >= 0, dos) - h = LatticePresets.honeycomb() |> hamiltonian(hopping(1)) |> unitcell((1,-1),(2,2)) |> wrap(2) - g = greens(h, Schur1D()) - dos = [-imag(tr(g(w))) for w in range(-3.2, 3.2, length = 101)] - @test all(x -> Quantica.chop(x) >= 0, dos) + h = LatticePresets.honeycomb() |> hamiltonian(hopping(1)) |> unitcell((1,-1),(2,2)) |> wrap(2) + g = greens(h, Schur1D()) + dos = [-imag(tr(g(w))) for w in range(-3.2, 3.2, length = 101)] + @test all(x -> Quantica.chop(x) >= 0, dos) - # Broken until https://github.com/Reference-LAPACK/lapack/issues/477 comes to julia - if VERSION >= v"1.7.0-DEV" h = LatticePresets.honeycomb() |> hamiltonian(hopping((r, dr) -> 1/(dr'*dr), range = 1)) |> unitcell((1,-1),(2,2)) |> wrap(2) g = greens(h, Schur1D()) dos = [-imag(tr(g(w))) for w in range(-3.2, 3.2, length = 101)]