% Block Matrix Inversion
% Minas Spetsakis, Sep. 2023

N = 5;
A = randn(N,N);
A = A*A';
B = randn(N,1);
D = 1;
S = [A, B; B', 1];

Ai = inv(A);
tmp = 1/(D-B'*Ai*B);

Si = [ Ai+Ai*B*tmp*B'*Ai, -Ai*B*tmp; -tmp*B'*Ai, tmp];
norm(Si-inv(S))
