Fixes for local updates (#162)

* fix SFD * Skip init! on creation * add fallback exponentiation to CMat and BlockDiagonal * use bigfloats to be extra safe * reduce precision to 128Bit * fix tests * add missing detratio methods * fix typo * fix out of bounds error when model flavors > field flavor * add more missing methods * more nflavor fixes * fix missing vmuladd! * allow global check to fail * log filename of failed load * allow adjustments of dq
carstenbauer · Feb 24, 2022 · a4dbd32 · a4dbd32
1 parent 177cc1e
commit a4dbd32
Show file tree

Hide file tree

Showing 10 changed files with 142 additions and 44 deletions.
diff --git a/src/FileIO.jl b/src/FileIO.jl
@@ -172,17 +172,23 @@ will be attempted. (This is useful to reduce the memory requirements.)
 """
 function load(filename::String, groups::String...; kwargs...)
     if isfile(filename)
-        data = if endswith(filename, "jld2")
-            FileWrapper(JLD2.jldopen(filename, "r"), filename)
-        else 
-            FileWrapper(JLD.load(filename), filename)
-        end
-        output = try 
-            if haskey(data, "MC") && !("MC" in groups)
-                _load(data, "MC", groups...) else _load(data, groups...)
+        output = try
+            data = if endswith(filename, "jld2")
+                FileWrapper(JLD2.jldopen(filename, "r"), filename)
+            else 
+                FileWrapper(JLD.load(filename), filename)
+            end
+            output = try 
+                if haskey(data, "MC") && !("MC" in groups)
+                    _load(data, "MC", groups...) else _load(data, groups...)
+                end
+            finally
+                endswith(filename, "jld2") && close(data.file)
             end
-        finally
-            endswith(filename, "jld2") && close(data.file)
+            output
+        catch e
+            println("Error loading file $filename:")
+            rethrow()
         end
 
         return output

diff --git a/src/flavors/DQMC/DQMC.jl b/src/flavors/DQMC/DQMC.jl
@@ -184,12 +184,16 @@ See also: [`resume!`](@ref)
     propagate(mc)
 
     # Check assumptions for global updates
-    copyto!(mc.stack.tmp2, mc.stack.greens)
-    udt_AVX_pivot!(mc.stack.tmp1, mc.stack.tempvf, mc.stack.tmp2, mc.stack.pivot, mc.stack.tempv)
-    ud = det(Matrix(mc.stack.tmp1))
-    td = det(Matrix(mc.stack.tmp2))
-    if !(0.9999999 <= abs(td) <= 1.0000001) || !(0.9999999 <= abs(ud) <= 1.0000001)
-        @error("Assumptions for global updates broken! ($td, $ud should be 1)")
+    try
+        copyto!(mc.stack.tmp2, mc.stack.greens)
+        udt_AVX_pivot!(mc.stack.tmp1, mc.stack.tempvf, mc.stack.tmp2, mc.stack.pivot, mc.stack.tempv)
+        ud = det(Matrix(mc.stack.tmp1))
+        td = det(Matrix(mc.stack.tmp2))
+        if !(0.9999999 <= abs(td) <= 1.0000001) || !(0.9999999 <= abs(ud) <= 1.0000001)
+            @error("Assumptions for global updates broken! ($td, $ud should be 1)")
+        end
+    catch e
+        @warn "Could not verify global update" exception = e
     end
 
     min_sweeps = round(Int, 1 / min_update_rate)

diff --git a/src/flavors/DQMC/fields.jl b/src/flavors/DQMC/fields.jl
@@ -65,6 +65,14 @@ function calculate_detratio!(cache::StandardFieldCache, Δ::Number, G::Union{FMa
     cache.detratio = cache.R * cache.R
 end
 
+function calculate_detratio!(cache::StandardFieldCache, Δ::Number, G::BlockDiagonal{Float64}, i)
+    # Unrolled R = I + Δ * (I - G)
+    @inbounds for b in eachindex(G.blocks)
+        cache.R[b] = 1.0 + Δ * (1.0 - G.blocks[b][i, i])
+    end
+    # determinant of Diagonal
+    cache.detratio = prod(cache.R)
+end
 
 function calculate_detratio!(cache::StandardFieldCache, Δ::FVec64, G::BlockDiagonal{Float64}, i)
     # Unrolled R = I + Δ * (I - G)
@@ -74,17 +82,49 @@ function calculate_detratio!(cache::StandardFieldCache, Δ::FVec64, G::BlockDiag
     # determinant of Diagonal
     cache.detratio = prod(cache.R)
 end
-# function calculate_detratio!(cache::StandardFieldCache, Δ::FVec64, G::BlockDiagonal{ComplexF64}, i)
-#     # Unrolled R = I + Δ * (I - G)
-#     @inbounds for b in eachindex(G.blocks)
-#         cache.R.re[b] = 1.0 + Δ[b] * (1.0 - G.blocks[b].re[i, i])
-#     end
-#     @inbounds for b in eachindex(G.blocks)
-#         cache.R.im[b] = - Δ[b] * G.blocks.im[b][i, i]
-#     end
-#     # determinant of Diagonal (don't think this is worth optimizing)
-#     cache.detratio = prod(cache.R)
-# end
+
+function calculate_detratio!(cache::StandardFieldCache, Δ::Float64, G::BlockDiagonal{ComplexF64}, i)
+    # Unrolled R = I + Δ * (I - G)
+    @inbounds for b in eachindex(G.blocks)
+        cache.R.re[b] = 1.0 + Δ * (1.0 - G.blocks[b].re[i, i])
+    end
+    @inbounds for b in eachindex(G.blocks)
+        cache.R.im[b] = - Δ * G.blocks[b].im[i, i]
+    end
+    # determinant of Diagonal (don't think this is worth optimizing)
+    cache.detratio = prod(cache.R)
+end
+
+function calculate_detratio!(cache::StandardFieldCache, Δ::FVec64, G::BlockDiagonal{ComplexF64}, i)
+    # Unrolled R = I + Δ * (I - G)
+    @inbounds for b in eachindex(G.blocks)
+        cache.R.re[b] = 1.0 + Δ[b] * (1.0 - G.blocks[b].re[i, i])
+    end
+    @inbounds for b in eachindex(G.blocks)
+        cache.R.im[b] = - Δ[b] * G.blocks[b].im[i, i]
+    end
+    # determinant of Diagonal (don't think this is worth optimizing)
+    cache.detratio = prod(cache.R)
+end
+
+function calculate_detratio!(cache::StandardFieldCache, Δ::ComplexF64, G::BlockDiagonal{ComplexF64}, i)
+    # Unrolled R = I + Δ * (I - G)
+    @inbounds for b in eachindex(G.blocks)
+        cache.R.re[b] = 1.0 + real(Δ) * (1.0 - G.blocks[b].re[i, i])
+    end
+    @inbounds for b in eachindex(G.blocks)
+        cache.R.re[b] += imag(Δ) * G.blocks[b].im[i, i]
+    end
+    @inbounds for b in eachindex(G.blocks)
+        cache.R.im[b] = imag(Δ) * (1.0 - G.blocks[b].re[i, i])
+    end
+    @inbounds for b in eachindex(G.blocks)
+        cache.R.im[b] -= real(Δ) * G.blocks[b].im[i, i]
+    end
+    # determinant of Diagonal
+    cache.detratio = prod(cache.R)
+end
+
 function calculate_detratio!(cache::StandardFieldCache, Δ::CVec64, G::BlockDiagonal{ComplexF64}, i)
     # Unrolled R = I + Δ * (I - G)
     @inbounds for b in eachindex(G.blocks)
@@ -149,13 +189,31 @@ function vldiv22!(cache::StandardFieldCache, R::FMat64, Δ::FVec64)
     nothing
 end
 
+function vldiv22!(cache::StandardFieldCache, R::FMat64, Δ::Float64)
+    @inbounds begin
+        inv_div = 1.0 / cache.detratio
+        cache.invRΔ[1, 1] =   R[2, 2] * Δ * inv_div
+        cache.invRΔ[1, 2] = - R[1, 2] * Δ * inv_div
+        cache.invRΔ[2, 1] = - R[2, 1] * Δ * inv_div
+        cache.invRΔ[2, 2] =   R[1, 1] * Δ * inv_div
+    end
+    nothing
+end
+
 function vldiv22!(cache::StandardFieldCache, R::FVec64, Δ::FVec64)
     @inbounds begin
         cache.invRΔ[1] = Δ[1] / R[1]
         cache.invRΔ[2] = Δ[2] / R[2]
     end
     nothing
 end
+function vldiv22!(cache::StandardFieldCache, R::FVec64, Δ::Float64)
+    @inbounds begin
+        cache.invRΔ[1] = Δ / R[1]
+        cache.invRΔ[2] = Δ / R[2]
+    end
+    nothing
+end
 
 function vldiv22!(cache::StandardFieldCache, R::CVec64, Δ::CVec64)
     @inbounds begin
@@ -170,6 +228,31 @@ function vldiv22!(cache::StandardFieldCache, R::CVec64, Δ::CVec64)
     nothing
 end
 
+function vldiv22!(cache::StandardFieldCache, R::CVec64, Δ::ComplexF64)
+    @inbounds begin
+        # Reminder: 1/c = c* / (cc*) (complex conjugate)
+        f1 = 1.0 / (R.re[1] * R.re[1] + R.im[1] * R.im[1])
+        f2 = 1.0 / (R.re[2] * R.re[2] + R.im[2] * R.im[2])
+        cache.invRΔ.re[1] = f1 * (real(Δ) * R.re[1] + imag(Δ) * R.im[1])
+        cache.invRΔ.re[2] = f2 * (real(Δ) * R.re[2] + imag(Δ) * R.im[2])
+        cache.invRΔ.im[1] = f1 * (imag(Δ) * R.re[1] - real(Δ) * R.im[1])
+        cache.invRΔ.im[2] = f2 * (imag(Δ) * R.re[2] - real(Δ) * R.im[2])
+    end
+    nothing
+end
+function vldiv22!(cache::StandardFieldCache, R::CVec64, Δ::Float64)
+    @inbounds begin
+        # Reminder: 1/c = c* / (cc*) (complex conjugate)
+        f1 = 1.0 / (R.re[1] * R.re[1] + R.im[1] * R.im[1])
+        f2 = 1.0 / (R.re[2] * R.re[2] + R.im[2] * R.im[2])
+        cache.invRΔ.re[1] = f1 * Δ * R.re[1]
+        cache.invRΔ.re[2] = f2 * Δ * R.re[2]
+        cache.invRΔ.im[1] = - f1 * Δ * R.im[1]
+        cache.invRΔ.im[2] = - f2 * Δ * R.im[2]
+    end
+    nothing
+end
+
 
 function update_greens!(cache::StandardFieldCache, G, i, N)
     # calculate Δ R⁻¹
@@ -280,9 +363,9 @@ end
 nflavors(::DensityHirschField) = 1
 
 @inline function interaction_matrix_exp!(f::DensityHirschField, result::Diagonal, slice, power)
-    # TODO this will index out of bound with 2 flavors
+    N = size(f.conf, 1)
     @inbounds for i in eachindex(result.diag)
-        result.diag[i] = exp(power * f.α * f.conf[i, slice])
+        result.diag[i] = exp(power * f.α * f.conf[mod1(i, N), slice])
     end
     nothing
 end
@@ -518,8 +601,9 @@ energy_boson(f::DensityGHQField, conf = f.conf) = f.α * sum(conf)
 
 # TODO: Maybe worth adding a complex method?
 @inline @bm function interaction_matrix_exp!(f::DensityGHQField, result::Diagonal, slice, power)
+    N = size(f.conf, 1)
     @inbounds for i in eachindex(result.diag)
-        result.diag[i] = exp(+power * f.α * f.η[f.conf[i, slice]])
+        result.diag[i] = exp(+power * f.α * f.η[f.conf[mod1(i, N), slice]])
     end
     return nothing
 end

diff --git a/src/flavors/DQMC/linalg/complex.jl b/src/flavors/DQMC/linalg/complex.jl
@@ -22,6 +22,11 @@ function vmuladd!(C::Matrix{T}, A::Matrix{T}, B::Diagonal{T}, factor::T = T(1))
         C[m,n] += factor * A[m,n] * B.diag[n]
     end
 end
+function vmuladd!(C::Matrix{T}, A::Diagonal{T}, B::Matrix{T}, factor::T = T(1)) where {T <: Real}
+    @turbo for m in 1:size(A, 1), n in 1:size(A, 2)
+        C[m,n] += factor * A.diag[m] * B[m,n]
+    end
+end
 function vmuladd!(C::Matrix{T}, A::Matrix{T}, X::Adjoint{T}, factor::T = T(1)) where {T <: Real}
     B = X.parent
     @turbo for m in 1:size(A, 1), n in 1:size(B, 2)

diff --git a/src/flavors/DQMC/measurements/extensions.jl b/src/flavors/DQMC/measurements/extensions.jl
@@ -188,9 +188,8 @@ function dia_K_x(mc, G, idxs)
 end
 
 # uses directions with filter > 0
-function para_ccc(mc, ccs, ind_dir)
+function para_ccc(mc, ccs, ind_dir, dq = [0.0, 0.0])
     #dq = [2pi / size(lattice(mc))[2], 0.0]
-    dq = [0.0, 0.0]
     dirs = directions(lattice(mc))
     Λxx = ComplexF64(0)
 

diff --git a/src/flavors/DQMC/measurements/generic.jl b/src/flavors/DQMC/measurements/generic.jl
@@ -147,7 +147,7 @@ _simple_buffer(mc, t::Type, T) = _simple_buffer(mc, t(), T)
 _simple_buffer(mc, ::Type{Nothing}, T) = nothing
 _simple_buffer(mc, ::Nothing, T) = nothing
 _simple_buffer(mc, ::EachSite, T) = zeros(T, length(lattice(mc)))
-_simple_buffer(mc, ::EachSiteAndFlavor, T) = zeros(T, nflavors(field(mc)) * length(lattice(mc)))
+_simple_buffer(mc, ::EachSiteAndFlavor, T) = zeros(T, nflavors(mc) * length(lattice(mc)))
 _simple_buffer(mc, ::EachSitePair, T) = zeros(T, length(lattice(mc)), length(lattice(mc)))
 
 
@@ -364,7 +364,7 @@ end
 
 # Lattice irrelevant
 @bm function measure!(::Nothing, measurement, mc::DQMC, model, sweep, packed_greens)
-    flv = Val(max(nflavors(field(mc)), nflavors(model)))
+    flv = Val(nflavors(mc))
     push!(measurement.observable, measurement.kernel(mc, model, packed_greens, flv))
     nothing
 end
@@ -378,7 +378,7 @@ end
 
 # Call kernel for each site (linear index)
 @bm function apply!(temp::Array, iter::DirectLatticeIterator, measurement, mc::DQMC, model, packed_greens)
-    flv = Val(max(nflavors(field(mc)), nflavors(model)))
+    flv = Val(nflavors(mc))
     for i in iter
         temp[i] += measurement.kernel(mc, model, i, packed_greens, flv)
     end
@@ -387,7 +387,7 @@ end
 
 # Call kernel for each pair (src, trg) (Nsties² total)
 @bm function apply!(temp::Array, iter::EachSitePair, measurement, mc::DQMC, model, packed_greens)
-    flv = Val(max(nflavors(field(mc)), nflavors(model)))
+    flv = Val(nflavors(mc))
     for (i, j) in iter
         temp[i, j] += measurement.kernel(mc, model, (i, j), packed_greens, flv)
     end
@@ -396,15 +396,15 @@ end
 
 # Call kernel for each pair (site, site) (i.e. on-site) 
 @bm function apply!(temp::Array, iter::_OnSite, measurement, mc::DQMC, model, packed_greens)
-    flv = Val(max(nflavors(field(mc)), nflavors(model)))
+    flv = Val(nflavors(mc))
     for (i, j) in iter
         temp[i] += measurement.kernel(mc, model, (i, j), packed_greens, flv)
     end
     nothing
 end
 
 @bm function apply!(temp::Array, iter::DeferredLatticeIterator, measurement, mc::DQMC, model, packed_greens)
-    flv = Val(max(nflavors(field(mc)), nflavors(model)))
+    flv = Val(nflavors(mc))
     @inbounds for idxs in iter
         temp[first(idxs)] += measurement.kernel(mc, model, idxs[2:end], packed_greens, flv)
     end
@@ -414,14 +414,14 @@ end
 
 # Sums
 @bm function apply!(temp::Array, iter::_Sum{<: DirectLatticeIterator}, measurement, mc::DQMC, model, packed_greens)
-    flv = Val(max(nflavors(field(mc)), nflavors(model)))
+    flv = Val(nflavors(mc))
     @inbounds for idxs in iter
         temp[1] += measurement.kernel(mc, model, idxs, packed_greens, flv)
     end
     nothing
 end
 @bm function apply!(temp::Array, iter::_Sum{<: DeferredLatticeIterator}, measurement, mc::DQMC, model, packed_greens)
-    flv = Val(max(nflavors(field(mc)), nflavors(model)))
+    flv = Val(nflavors(mc))
     @inbounds for idxs in iter
         temp[1] += measurement.kernel(mc, model, idxs[2:end], packed_greens, flv)
     end

diff --git a/src/flavors/DQMC/measurements/measurements.jl b/src/flavors/DQMC/measurements/measurements.jl
@@ -52,7 +52,7 @@ function greens_measurement(
         mc::DQMC, model::Model, greens_iterator = Greens(); 
         capacity = _default_capacity(mc), eltype = geltype(mc),
         obs = let
-            N = length(lattice(model)) * nflavors(field(mc))
+            N = length(lattice(model)) * nflavors(mc)
             LogBinner(zeros(eltype, (N, N)), capacity=capacity)
         end, kwargs...
     )

diff --git a/src/flavors/DQMC/stack.jl b/src/flavors/DQMC/stack.jl
@@ -150,7 +150,7 @@ function initialize_stack(mc::DQMC, ::DQMCStack)
     GreensMatType = gmattype(mc)
     HoppingElType = heltype(mc)
     N = length(lattice(mc))
-    flv = nflavors(field(mc))
+    flv = nflavors(mc)
 
     # Generate safe multiplication chunks
     # - every chunk must have ≤ safe_mult elements

diff --git a/src/flavors/DQMC/unequal_time_stack.jl b/src/flavors/DQMC/unequal_time_stack.jl
@@ -61,7 +61,7 @@ function initialize_stack(mc::DQMC, s::UnequalTimeStack)
     GreensElType  = geltype(mc)
     GreensMatType = gmattype(mc)
     N = length(lattice(mc))
-    flv = nflavors(field(mc))
+    flv = nflavors(mc)
     # mirror stack
     M = mc.stack.n_elements
 

diff --git a/src/lattices/lattice_iterators.jl b/src/lattices/lattice_iterators.jl
@@ -179,7 +179,7 @@ Creates an iterator template which iterates through the diagonal of the Greensfu
 """
 struct EachSiteAndFlavor <: DirectLatticeIteratorTemplate end
 function (::EachSiteAndFlavor)(mc::MonteCarloFlavor, model::Model)
-    RangeIterator(1 : length(lattice(mc)) * nflavors(field(mc)))
+    RangeIterator(1 : length(lattice(mc)) * nflavors(mc))
 end
 
 """