%% test if eig exploits block-diagonal structure [NO, block one is better]
nr_blks = 20;
blk_sizes = randi(40,nr_blks,1);
M = cell(nr_blks,1);
for i = 1:nr_blks
    M{i} = rand(blk_sizes(i));
    M{i} = (M{i}+M{i}')/2;
end


tic
for exp = 1:1e3
    for i = 1:nr_blks
        [V,D] = eig(M{i});
        IND =  find(diag(D) > 1e-5);
        M_plus = V(:,IND)*D(IND,IND)*V(:,IND)';
        M_plus = (M_plus+M_plus')/2;
    end
end
toc


M_full = blkdiag(M{:});
tic
for exp = 1:1e3
    [V,D] = eig(M_full);
    IND =  find(diag(D) > 1e-5);
    M_plus = V(:,IND)*D(IND,IND)*V(:,IND)';
    M_plus = (M_plus+M_plus')/2;
end
toc

%% difference between ~logical(m) and false(m)
m = 1e4;

tic
for i = 1:m
    blk_nz = ~logical(m);
end
toc

m = 1e4;

tic
for i = 1:m
    blk_nz = false(m);
end
toc

%% dense * sparse matrix, [it helps if one of them is sparse]
n = 1000;
A = sprand(1000,1000,0.005);
B = sprand(1000,1000,0.005);
Af = full(A);
Bf = full(B);

timeit(@() Af*Bf)
timeit(@() A*B)
timeit(@() Af*B)
timeit(@() A*Bf)

% %%
% C = rand(n,n);
% timeit(@() C*B*C')
% timeit(@() C*B*C')