-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathtest_speedcomparison.m
53 lines (43 loc) · 1.56 KB
/
test_speedcomparison.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
% test_speedcomparison.m
%
% Compare non-symmetric and symmetric versions of code
N = 50; % size of blocks
T = 1e4; % number of blocks
% ===== Build sparse block-tridiagonal matrix =============
fprintf('Building sparse block-tridiagonal matrix...\n');
tic;
% Build a random bi-diagonal matrix
[jj,ii] = meshgrid(1:N);
ii = ii(:); % convert to column vector
jj = jj(:); % convert to column vector
% First, compute indices for block bi-diagonal matrix
nelts = N^2*(T*2-1); % number of elements
Iinds = zeros(nelts,1); % row indices
Jinds = zeros(nelts,1); % column indices
irow = 0;
for k = 1:T-1
Iinds(irow+1:irow+(N^2*2)) = [ii+(k-1)*N; ii+k*N];
Jinds(irow+1:irow+(N^2*2)) = [jj+(k-1)*N; jj+(k-1)*N];
irow = irow+N^2*2;
end
Iinds(irow+1:irow+N^2) = ii+(T-1)*N; % last block
Jinds(irow+1:irow+N^2) = jj+(T-1)*N; % last block
D = sparse(Iinds,Jinds,randn(length(Iinds),1),T*N,T*N); % fill block bi-diag w/ random entries
% Now make symmetric sparse tri-diagonal matrix
B = D'*D+speye(N*T,N*T)*.1;
% Compute total build time
buildTime = toc
% Make plot showing block structure of L
ni = min(1000,nelts); % # of rows+cols to plot
clf; spy(B(1:ni,1:ni));
title('symmetric block-tridiagonal matrix');
drawnow;
%% ==== Compare timing of full and symmetric versions ============
tic;
[vdiag1,Dblocks1,OffDblocks1] = invblktridiag(B,N);
computeTime1 = toc;
tic;
[vdiag2,Dblocks2,OffDblocks2] = invblktridiag_sym(B,N);
computeTime2 = toc;
fprintf('Compute time using invblktridiag: %.2f\n',computeTime1);
fprintf('Compute time using invblktridiag_sym: %.2f\n',computeTime2);