Simplify similarmatrix, remove optimize!

README.md

@@ -18,6 +18,7 @@ The Quantica.jl package provides an expressive API to build arbitrary quantum sy
 - `lattice`, `sublat`: build lattices
 - `hopping`, `onsite`, `siteselector`, `hopselector`, `nrange`, `not`: build tightbinding models
 - `hamiltonian`: build a Hamiltonian from tightbinding model and a lattice
+- `bloch`, `bloch!`, `similarmatrix`: build the Bloch matrix of a Hamiltonian
 - `parametric`, `@onsite!`, `@hopping!`, `parameters`: build a parametric Hamiltonian
 - `dims`, `sitepositions`, `siteindices`, `bravais`: inspect lattices and Hamiltonians
 - `supercell`, `unitcell`, `flatten`, `wrap`, `transform!`, `combine`: build derived lattices or Hamiltonians

src/KPM.jl

@@ -31,10 +31,10 @@ function KPMBuilder(h, A, kets, order, bandrange)
     return builder
 end
 
-function matrixKPM(h::Hamiltonian{<:Lattice,L}, method = missing) where {L}
+function matrixKPM(h::Hamiltonian{<:Lattice,L}, matrixtype = missing) where {L}
     iszero(L) ||
         throw(ArgumentError("Hamiltonian is defined on an infinite lattice. Reduce it to zero-dimensions with `wrap` or `unitcell`."))
-    m = similarmatrix(h, method)
+    m = similarmatrix(h, matrixtype)
     return bloch!(m, h)
 end
 
@@ -197,7 +197,7 @@ function bandbracketKPM(h, ::Missing)
 end
 bandbracketKPM(h, (ϵmin, ϵmax)::Tuple{T,T}, pad = float(T)(0.01)) where {T} = ((ϵmax + ϵmin) / 2, (ϵmax - ϵmin) / (2 - pad))
 
-bandrangeKPM(h::Hamiltonian) = bandrangeKPM(matrixKPM(h, ArnoldiMethodPackage()))
+bandrangeKPM(h::Hamiltonian) = bandrangeKPM(matrixKPM(h, flatten))
 
 function bandrangeKPM(h::AbstractMatrix{T}) where {T}
     @warn "Computing spectrum bounds... Consider using the `bandrange` option for faster performance."

src/Quantica.jl

@@ -21,10 +21,10 @@ export sublat, bravais, lattice, dims, supercell, unitcell,
        hopping, onsite, @onsite!, @hopping!, parameters, siteselector, hopselector, nrange,
        sitepositions, siteindices, not,
        ket, randomkets,
-       hamiltonian, parametric, bloch, bloch!, optimize!, similarmatrix,
+       hamiltonian, parametric, bloch, bloch!, similarmatrix,
        flatten, wrap, transform!, combine,
        spectrum, bandstructure, marchingmesh, linearmesh, buildmesh, buildlift,
-       defaultmethod, bands, vertices,
+       bands, vertices,
        energies, states,
        momentaKPM, dosKPM, averageKPM, densityKPM, bandrangeKPM,
        greens, greensolver

src/bandstructure.jl

@@ -7,7 +7,7 @@ struct Spectrum{E,T,A<:AbstractMatrix{T}}
 end
 
 """
-    spectrum(h; method = defaultmethod(h), transform = missing)
+    spectrum(h; method = LinearAlgebraPackage(), transform = missing)
 
 Compute the spectrum of a 0D Hamiltonian `h` (or alternatively of the bounded unit cell of a
 finite dimensional `h`) using one of the following `method`s
@@ -27,8 +27,8 @@ The energies and eigenstates in the resulting `s::Spectrum` object can be access
     `energies`, `states`, `bandstructure`
 
 """
-function spectrum(h; method = defaultmethod(h), transform = missing)
-    matrix = similarmatrix(h, method)
+function spectrum(h; method = LinearAlgebraPackage(), transform = missing)
+    matrix = similarmatrix(h, method_matrixtype(method, h))
     bloch!(matrix, h)
     (ϵk, ψk) = diagonalize(matrix, method)
     s = Spectrum(ϵk, ψk)
@@ -193,7 +193,7 @@ Curried form of the above equivalent to `bandstructure(h, [mesh]; kw...)`.
 
 The default options are
 
-    (lift = missing, minoverlap = 0, method = defaultmethod(h), transform = missing)
+    (lift = missing, minoverlap = 0, method = LinearAlgebraPackage(), transform = missing)
 
 `lift`: when not `missing`, `lift` is a function `lift = (vs...) -> ϕ`, where `vs` are the
 coordinates of a mesh vertex and `ϕ` are Bloch phases if sampling a `h::Hamiltonian`, or
@@ -264,9 +264,9 @@ function bandstructure(h::Union{Hamiltonian,ParametricHamiltonian}, spec::MeshSp
 end
 
 function bandstructure(h::Union{Hamiltonian,ParametricHamiltonian}, mesh::Mesh;
-                       method = defaultmethod(h), lift = missing, minoverlap = 0, transform = missing)
+                       method = LinearAlgebraPackage(), lift = missing, minoverlap = 0, transform = missing)
     # ishermitian(h) || throw(ArgumentError("Hamiltonian must be hermitian"))
-    matrix = similarmatrix(h, method)
+    matrix = similarmatrix(h, method_matrixtype(method, h))
     codiag = codiagonalizer(h, matrix, mesh, lift)
     diag = diagonalizer(method, codiag, minoverlap)
     matrixf(ϕs) = bloch!(matrix, h, applylift(lift, ϕs))
@@ -276,12 +276,11 @@ function bandstructure(h::Union{Hamiltonian,ParametricHamiltonian}, mesh::Mesh;
 end
 
 function bandstructure(matrixf::Function, mesh::Mesh;
-                       method = missing, lift = missing, minoverlap = 0, transform = missing)
+                       method = LinearAlgebraPackage(), lift = missing, minoverlap = 0, transform = missing)
     matrixf´ = _wraplift(matrixf, lift)
     matrix = _samplematrix(matrixf´, mesh)
-    method´ = method === missing ? defaultmethod(matrix) : method
     codiag = codiagonalizer(matrixf´, matrix, mesh, missing)
-    diag = diagonalizer(method´, codiag, minoverlap)
+    diag = diagonalizer(method, codiag, minoverlap)
     b = _bandstructure(matrixf´, matrix, mesh, diag)
     transform === missing || transform!(transform, b)
     return b

src/diagonalizer.jl

@@ -14,8 +14,6 @@ diagonalizer(method, codiag, minoverlap) = Diagonalizer(method, codiag, Float64(
 
 ## Diagonalize methods ##
 
-defaultmethod(h::Union{Hamiltonian,ParametricHamiltonian,AbstractMatrix}) = LinearAlgebraPackage()
-
 checkloaded(package::Symbol) = isdefined(Main, package) ||
     throw(ArgumentError("Package $package not loaded, need to be `using $package`."))
 
@@ -31,8 +29,6 @@ function diagonalize(matrix, method::LinearAlgebraPackage)
     return ϵ, ψ
 end
 
-similarmatrix(h, ::LinearAlgebraPackage) = similarmatrix(h, Matrix{blockeltype(h)})
-
 ## Arpack ##
 struct ArpackPackage{K<:NamedTuple} <: AbstractDiagonalizeMethod
     kw::K
@@ -45,8 +41,6 @@ function diagonalize(matrix, method::ArpackPackage)
     return ϵ, ψ
 end
 
-similarmatrix(h, ::ArpackPackage) = similarmatrix(h, SparseMatrixCSC{blockeltype(h)})
-
 ## ArnoldiMethod ##
 struct ArnoldiMethodPackage{K<:NamedTuple} <: AbstractDiagonalizeMethod
     kw::K
@@ -59,8 +53,6 @@ function diagonalize(matrix, method::ArnoldiMethodPackage)
     return ϵ, ψ
 end
 
-similarmatrix(h, ::ArnoldiMethodPackage) = similarmatrix(h, SparseMatrixCSC{blockeltype(h)})
-
 ## IterativeSolvers ##
 
 struct KrylovKitPackage{K<:NamedTuple} <: AbstractDiagonalizeMethod
@@ -92,7 +84,12 @@ function diagonalize(matrix::AbstractMatrix{M}, method::KrylovKitPackage) where
     return ϵ´, ψ´
 end
 
-similarmatrix(h, ::KrylovKitPackage) = similarmatrix(h, SparseMatrixCSC{blockeltype(h)})
+### matrix types
+
+similarmatrix(h, method::AbstractDiagonalizeMethod) = similarmatrix(h, method_matrixtype(method, h))
+
+method_matrixtype(::LinearAlgebraPackage, h) = Matrix{blockeltype(h)}
+method_matrixtype(::AbstractDiagonalizeMethod, h) = flatten
 
 #######################################################################
 # shift and invert methods
@@ -198,7 +195,7 @@ function codiagonalizer(h, matrix::AbstractMatrix{T}, mesh, lift) where {T}
 end
 
 function codiag_function(h::Union{Hamiltonian,ParametricHamiltonian}, matrix, lift, dirs, delta)
-    hdual = Dual(h)
+    hdual = dual_if_parametric(h)
     matrixdual = dualarray(matrix)
     anyold = anyoldmatrix(matrix)
     ndirs = length(dirs)
@@ -216,6 +213,9 @@ function codiag_function(h::Union{Hamiltonian,ParametricHamiltonian}, matrix, li
     return comatrix, matrixindices
 end
 
+dual_if_parametric(ph::ParametricHamiltonian) = Dual(ph)
+dual_if_parametric(h::Hamiltonian) = h
+
 # In the Function case we cannot know what directions to scan (arguments of matrixf). Also,
 # we cannot be sure that dual numbers propagate. We thus restrict to finite differences in the mesh
 function codiag_function(matrixf::Function, matrix, lift, meshdirs, delta)

src/greens.jl

@@ -56,7 +56,6 @@ greens(solver::Function) = h -> greens(h, solver(h))
 # Needed to make similarmatrix work with GreensFunction
 matrixtype(g::GreensFunction) = Matrix{eltype(g.h)}
 Base.parent(g::GreensFunction) = g.h
-optimize!(g::GreensFunction) = g
 
 #######################################################################
 # BandGreenSolver