%----------------------------------------------
% block operations are faster
m = 1000;
n = 2*m;
p = m;

A = rand(m,n);
B = rand(n,p);

tic; C1 = A*B; toc

C2 = zeros(m,p);
tic
for di = 1:p
    C2(:,di) = A*B(:,di);
end
toc

norm(C1-C2,'fro')

C3 = zeros(m,p);
tic
for di = 1:m
    for dj = 1:p
        C3(di,dj) = A(di,:)*B(:,dj);
    end
end
toc
norm(C1-C3,'fro')


% C3 = zeros(m,p);
% tic
% for di = 1:n
%     C3 = C3 + A(:,di)*B(di,:);
% end
% toc
% norm(C1-C3,'fro')



%----------------------------------------------

%----------------------------------------------
% chol
% n = 100;
%A = gallery('poisson', n); 
nx = 50; ny = nx;
A = gallery('wathen',nx,ny);

figure(1); 
subplot(1,3,1); spy(A);

% A = L*L'
L1 = chol(A, 'lower');
subplot(1,3,2); spy(L1);

% L*L'=S'*A*S
[L2, pp, ss] = chol(A, 'lower','vector');
subplot(1,3,3); spy(L2);

[mm, nn] = size(A);
uu = rand(mm,1);
bb = A*uu;

tic; 
x1 = L1'\(L1\bb);
toc
norm(x1-uu)

tic; 
x2 = L2'\(L2\bb(ss));
toc
norm(x2-uu(ss))


%----------------------------------------------

A = full(delsq(numgrid('L', 10)));
B = gallery('uniformdata',10,0);
M = [eye(10) B; B' zeros(10)]; 

[LA,DA] = ldl(A); 
fprintf('The factorization error ||A - LA*DA*LA''|| is %g\n', ...
norm(A - LA*DA*LA'));
neginds = find(diag(DA) < 0)

figure; spy(DA); title('Structure of D from ldl(A)');
[Las, Das] = ldl(A - 4*eye(size(A)));
figure; spy(Das); 
title('Structure of D from ldl(A - 4*eye(size(A)))');

[Lm, Dm, Pm] = ldl(M);
fprintf(1, ...
'The error norm ||Pm''*M*Pm - Lm*Dm*Lm''|| is %g\n', ...
norm(Pm'*M*Pm - Lm*Dm*Lm'));
fprintf(1, ...
'The difference between Lm and tril(Lm) is %g\n', ...
norm(Lm - tril(Lm)));

LA = chol(M)
[LA, pp] = chol(M)

%----------------------------------------------
% SMW formula

n = 1000;
k = max(1, round(0.1*n));
B = rand(n,k);

A = eye(n) + B*B';
uu = rand(n,1);
bb = A*uu;

tic;
A = eye(n) + B*B';
x1 = A\bb;
toc;
norm(x1-uu)

tic;
% A = eye(n) + B*B';
x1 = A\bb;
toc;
norm(x1-uu)

tic;
x2 = bb - B* ((eye(k) + B'*B)\(B'*bb));
toc
norm(x2-uu)