# Usage

This section demonstrates how to use the Checkerboard.jl package.

## The Square Lattice

As an example, we will demonstrate how to construct and apply the checkerboard decomposition method in the case of a square lattice with isotropic nearest-neighbor hopping. We begin by importing all relevant packages.

using LinearAlgebra
using BenchmarkTools
using LatticeUtilities
using Checkerboard

# set number of BLAS threads to 1 for fair benchmarking purposes.
BLAS.set_num_threads(1);

The next step is to construct the neighbor table, and record the corresponding hopping amplitudes, for a square lattice with isotropic nearest-neighbor hopping. We also define a discretization in imaginary time $\Delta\tau$, which is used as the small parameter in the checkerboard decomposition approximation. We will use the package LatticeUtilities.jl to assist with constructing the neighbor table.

# construct square lattice neighbor table
L       = 12 # lattice size
square  = UnitCell(lattice_vecs = [[1.,0.],[0.,1.]], basis_vecs = [[0.,0.]])
lattice = Lattice(L = [L,L], periodic = [true,true])
bond_x  = Bond(orbitals = (1,1), displacement = [1,0])
bond_y  = Bond(orbitals = (1,1), displacement = [0,1])
nt      = build_neighbor_table([bond_x,bond_y], square, lattice)
N       = nsites(square, lattice)

# corresponding hopping for each pair of neighbors in the neighbor table
t = fill(1.0, size(nt,2))

# discretization in imaginary time i.e. the small parameter
# used in the checkerboard approximation
Δτ = 0.1

nt
2×288 Matrix{Int64}:
1  2  3  4  5  6  7  8   9  10  11  12  …  138  139  140  141  142  143  144
2  3  4  5  6  7  8  9  10  11  12   1       6    7    8    9   10   11   12

Next, for comparison purposes, we explicitly construct the hopping matrix $K$ and exponentiate it, giving us the exact matrix $e^{-\Delta\tau K}$ that the checkerboard decomposition matrix $\Gamma$ is intended to approximate.

K = zeros(Float64, N, N)
for c in 1:size(nt,2)
i      = nt[1,c]
j      = nt[2,c]
K[i,j] = -t[c]
K[j,i] = -conj(t[c])
end

expnΔτK = Hermitian(exp(-Δτ*K))

An instance of the CheckerboardMatrix type representing the checkerboard decomposition matrix $\Gamma$ can be instantiated as follows.

Γ = CheckerboardMatrix(nt, t, Δτ)

It is also straight forward to efficiently construct representations of both the transpose and inverse of the checkerboard matrix $\Gamma$ using the transpose and inv methods.

Γᵀ  = transpose(Γ)
Γ⁻¹ = inv(Γ)

Matrix-matrix multiplications involving $\Gamma$ can be performed using the mul! method just as you would with a standard matrix. Here we benchmark the matrix-matrix multiplication using the checkerboard decomposition as compared to multiplication by the dense matrix $e^{-\Delta\tau K}$.

I_dense = Matrix{Float64}(I,N,N)
B       = similar(I_dense)
@btime mul!(B, expnΔτK, I_dense); # 131.958 μs (0 allocations: 0 bytes)
@btime mul!(B, I_dense, expnΔτK); # 129.375 μs (0 allocations: 0 bytes)
@btime mul!(B, Γ,       I_dense); # 38.334 μs (0 allocations: 0 bytes)
@btime mul!(B, I_dense, Γ); # 17.500 μs (0 allocations: 0 bytes)

Using the @btime macro exported by the BenchmarkTools.jl, we see that performing both left and right matrix multiplies by the checkerboard matrix $\Gamma$ is faster than multiplying by the $e^{-\Delta\tau K}$. However, we see that right matrix multiplies are significantly faster than left multiplies by $\Gamma$, a result of the access order into I_dense being unavoidably more efficient for the right multiply. Therefore, it is important to keep this in mind when using this package to develop QMC codes: wherever possible do right multiplies by the checkerboard matrix instead of left multiplies.

Let us quickly verify that $\Gamma$ is a good approximation of $e^{-\Delta\tau K}$.

I_dense = Matrix{Float64}(I,N,N)
Γ_dense = similar(I_dense)
mul!(Γ_dense, Γ, I_dense)
norm(Γ_dense - expnΔτK) / norm(expnΔτK)
0.009848239506872646

It is also possible to do in-place left and right multiplications by $\Gamma$ using the lmul! and rmul! methods respectively.

A = Matrix{Float64}(I,N,N)
B = Matrix{Float64}(I,N,N)
rmul!(A,Γ) # A = A⋅Γ
lmul!(Γ,B) # B = Γ⋅B
A ≈ B
true

The checkerboard matrix is in reality a product of several individually sparse matrices that we will refer to as checkerboard color matrices, such that $\Gamma = \prod_{c} \Gamma_c = \Gamma_C \dots \Gamma_1$. Let us check to see how many checkerboard color matrices there are in our checkerboard decomposition.

Γ.Ncolors
4

It is also possible to multiply by a single one of the $\Gamma_c$ matrices with the mul!, lmul! and rmul! methods.

Γ₁ = Matrix{Float64}(I, N, N)
Γ₂ = Matrix{Float64}(I, N, N)
Γ₃ = Matrix{Float64}(I, N, N)
Γ₄ = Matrix{Float64}(I, N, N)

mul!(Γ₁, Γ, I_dense, 1)
mul!(Γ₂, I_dense, Γ, 2)
lmul!(Γ, Γ₃, 3)
rmul!(Γ₄, Γ, 4)

Γ_dense ≈ (Γ₄*Γ₃*Γ₂*Γ₁)
true

Lastly, it is easy to update the checkerboard decomposition with new hopping parameter values using the update! method, including disordered values.

@. t = 1.0 + 0.1*randn()
update!(Γ, t, Δτ)

This is important for QMC simulations of models like the Su-Schrieffer-Heeger model, where an electron-phonon coupling modulates the hopping amplitude, resulting in the hopping amplitudes changing over the course of a simulation.