Thread-safety fixes in bands

docs/src/tutorial/bandstructures.md

@@ -17,6 +17,9 @@ States:
 ```
 The above destructuring syntax assigns eigenvalues and eigenvectors to `eᵢ` and `ψᵢ`, respectively. The available eigensolvers and their options can be checked in the `EigenSolvers` docstrings.
 
+!!! warning "Arpack solver is not thread-safe"
+    `EigenSolver.Arpack` relies on a Fortran library that is not currently thread-safe. If you launch Julia with multiple threads, they will not be used with this specific solver. Otherwise Arpack would segfault.
+
 We define a "bandstructure" of an `h::AbstractHamiltonian` as a linear interpolation of its eigenpairs over a portion of the Brillouin zone, which is discretized with a base mesh of `ϕᵢ` values. At each `ϕᵢ` of the base mesh, the Bloch matrix `h(ϕᵢ)` is diagonalized with `spectrum`. The adjacent eigenpairs `eⱼ(ϕᵢ), ψⱼ(ϕᵢ)` are then connected ("stitched") together into a number of band meshes with vertices `(ϕᵢ..., eⱼ(ϕᵢ))` by maximizing the overlap of adjacent `ψⱼ(ϕᵢ)` (since the bands should be continuuous). Degenerate eigenpairs are collected into a single node of the band mesh.
 
 The bandstructure of an `h::AbstractHamiltonian` is computed using `bands`:

src/apply.jl

@@ -219,12 +219,13 @@ function apply(solver::AbstractEigenSolver, h::AbstractHamiltonian, ::Type{S}, m
     h´ = minimal_callsafe_copy(h)
     # Some solvers (e.g. ES.LinearAlgebra) only accept certain matrix types
     # so this mat´ could be an alias of the call! output, or an unaliased conversion
-    mat´ = ES.input_matrix(solver, h)
+    mat´ = ES.input_matrix(solver, h´)
     function sfunc(φs)
-        φs´ = apply_map(mapping, φs)    # this can be a FrankenTuple
+        φs´ = apply_map(mapping, φs)      # this can be a FrankenTuple
         mat = call!(h´, φs´)
         mat´ === mat || copy!(mat´, mat)
-        # the solver always receives the matrix type declared by ES.input_matrix
+        # mat´ could be dense, while mat is sparse, so if not egal, we copy
+        # the solver always receives the type of matrix mat´ declared by ES.input_matrix
         eigen = solver(mat´)
         apply_transform!(eigen, transform)
         return eigen

src/bands.jl

@@ -88,8 +88,12 @@ function bands(h::Function, mesh::Mesh{S};
     return ss
 end
 
-eigensolvers_thread_pool(solver, h, S, mapping, transform) =
-    [apply(solver, h, S, mapping, transform) for _ in 1:Threads.nthreads()]
+function eigensolvers_thread_pool(solver, h, S, mapping, transform)
+    # if h::Function we cannot be sure it is thread-safe
+    nsolvers = ES.is_thread_safe(solver) && h isa AbstractHamiltonian ? Threads.nthreads() : 1
+    solvers = [apply(solver, h, S, mapping, transform) for _ in 1:nsolvers]
+    return solvers
+end
 
 function subbands(hf, solvers, basemesh::Mesh{SVector{L,T}};
          showprogress = true, defects = (), patches = 0, degtol = missing, split = true, warn = true, projectors = false) where {T,L}
@@ -250,11 +254,20 @@ function subbands_diagonalize!(data)
     baseverts = vertices(data.basemesh)
     meter = Progress(length(baseverts), "Step 1 - Diagonalizing: ")
     push!(data.coloffsets, 0) # first element
-    Threads.@threads :static for i in eachindex(baseverts)
-        vert = baseverts[i]
-        solver = data.solvers[Threads.threadid()]
-        data.eigens[i] = solver(vert)
-        data.showprogress && ProgressMeter.next!(meter)
+    if length(data.solvers) > 1
+        Threads.@threads :static for i in eachindex(baseverts)
+            vert = baseverts[i]
+            solver = data.solvers[Threads.threadid()]
+            data.eigens[i] = solver(vert)
+            data.showprogress && ProgressMeter.next!(meter)
+        end
+    else
+        solver = first(data.solvers)
+        for i in eachindex(baseverts)
+            vert = baseverts[i]
+            data.eigens[i] = solver(vert)
+            data.showprogress && ProgressMeter.next!(meter)
+        end
     end
     # Collect band vertices and store column offsets
     for (basevert, eigen) in zip(baseverts, data.eigens)

src/docstrings.jl

@@ -1088,7 +1088,7 @@ unflat
 alias `ES` can be used in place of `EigenSolvers`. Currently supported solvers are
 
     ES.LinearAlgebra(; kw...)       # Uses `eigen(mat; kw...)` from the `LinearAlgebra` package
-    ES.Arpack(; kw...)              # Uses `eigs(mat; kw...)` from the `Arpack` package
+    ES.Arpack(; kw...)              # Uses `eigs(mat; kw...)` from the `Arpack` package (WARNING: Arpack is not thread-safe)
     ES.KrylovKit(params...; kw...)  # Uses `eigsolve(mat, params...; kw...)` from the `KrylovKit` package
     ES.ArnoldiMethod(; kw...)       # Uses `partialschur(mat; kw...)` from the `ArnoldiMethod` package
 

src/solvers/eigen.jl

@@ -26,6 +26,8 @@ using Quantica: Quantica, AbstractEigenSolver, ensureloaded, SVector, SMatrix,
 # an alias of h's call! output makes apply call! conversion a no-op, see apply.jl
 input_matrix(::AbstractEigenSolver, h) = call!_output(h)
 
+is_thread_safe(::AbstractEigenSolver) = true
+
 #### LinearAlgebra #####
 
 struct LinearAlgebra{K} <: AbstractEigenSolver
@@ -60,16 +62,19 @@ function (solver::Arpack)(mat::AbstractMatrix{<:Number})
     return sanitize_eigen(ε, Ψ)
 end
 
+# See https://github.com/JuliaLinearAlgebra/Arpack.jl/issues/86
+is_thread_safe(::Arpack) = false
+
 #### KrylovKit #####
 
-struct KrylovKit{P,K} <: AbstractEigenSolver
+struct KrylovKit{P<:Tuple,K<:NamedTuple} <: AbstractEigenSolver
     params::P
     kwargs::K
 end
 
 function KrylovKit(params...; kw...)
     ensureloaded(:KrylovKit)
-    return KrylovKit(params, kw)
+    return KrylovKit(params, NamedTuple(kw))
 end
 
 function (solver::KrylovKit)(mat)

test/test_bandstructure.jl

@@ -1,4 +1,5 @@
 using Quantica: nsubbands, nvertices, nedges, nsimplices
+using Random
 
 @testset "basic bandstructures" begin
     h = LatticePresets.honeycomb() |> hamiltonian(hopping(-1))
@@ -57,6 +58,12 @@ end
     b = bands(hf, m, showprogress = false, mapping = x -> 2π * x)
     @test nsubbands(b) == 1
     @test nsimplices(b)  == 36
+    # teting thread safety - we should fall back to a single thread for hf::Function
+    hf((x,)) = Quantica.call!(hc, (x, -x))
+    m = subdiv(0,2π,40)
+    Random.seed!(1) # to have ArnoldiMethod be deterministic
+    b = bands(hf, m, showprogress = false, solver = ES.ArnoldiMethod(nev = 18))
+    @test nsubbands(b) <= 2    # there is a random, platform-dependent component to this
 
     hp2 = LatticePresets.honeycomb() |> hamiltonian(hopping(-1), @hopping!((t; s) -> s*t))
     hf2((s, x)) = Matrix(Quantica.call!(hp2, (x, x); s))
@@ -82,8 +89,9 @@ end
         supercell |> hamiltonian(@onsite!((o; k) -> o + k*I), @hopping!((t; k = 2, p = [1,2])-> t - k*I .+ p'p))
     b = bands(ph, mesh2D..., mapping = (k, φ) -> ftuple(; k = k, p = SA[1, φ]), showprogress = false)
     @test nsubbands(b)  == 1
-    # multithreading loop throws a CompositeException
-    @test_throws CompositeException bands(ph, mesh2D..., mapping = (k, φ) -> ftuple(; p = SA[1, φ]), showprogress = false)
+    # multithreading loop does not throw error
+    b = bands(ph, mesh2D..., mapping = (k, φ) -> ftuple(; k, p = SA[1, φ]), showprogress = false)
+    @test nsubbands(b) == 1
 end
 
 @testset "spectrum" begin