Understanding static condensation

StaticCondensation.jl

using LinearAlgebra, Test

# Initialize matrices with proper dimensions
A11 = rand(2, 2)
A12 = rand(2, 2)
A21 = rand(2, 2)  # Doesn't necessarily need symmetry
A22 = rand(2, 2)
A23 = rand(2, 2)
A32 = rand(2, 2)  # Doesn't necessarily need symmetry
A33 = rand(2, 2)
A34 = rand(2, 2)
A43 = rand(2, 2)  # Doesn't necessarily need symmetry
A44 = rand(2, 2)
Z00 = zeros(2, 2)

# Assemble global matrix
A = [A11 A12 Z00 Z00;
     A21 A22 A23 Z00;
     Z00 A32 A33 A34;
     Z00 Z00 A43 A44]

# Initialize right-hand side
a = rand(2)
b = rand(2)
c = rand(2)
d = rand(2)
rhs = [a; b; c; d]
x = A \ rhs

# Static condensation steps
luA44 = lu(A44)
S34 = A33 - A34 * (luA44 \ A43)
luS34 = lu(S34)
S23 = A22 - A23 * (luS34 \ A32)
luS23 = lu(S23)
S12 = A11 - A12 * (luS23 \ A21)
luS12 = lu(S12)

# Solve system (corrected ordering) to get α with intermediate β, γ, and δ 
δ = luA44 \ d
γ = luS34 \ (c - A34 * δ)
β = luS23 \ (b - A23 * γ)
α = luS12 \ (a - A12 * β)

# Now apply corrections
β -= luS23 \ A21 * α
γ -= luS34 \ A32 * β
δ -= luA44 \ A43 * γ

# Combine solution
y = [α; β; γ; δ]

# Verify solution
@testset "Static Condenstation" begin
  @test x[1:2] ≈ α
  @test x[3:4] ≈ β
  @test x[5:6] ≈ γ
  @test x[7:8] ≈ δ
  @test A * y ≈ rhs # Check residual
end