[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Condensing indices in a matrix
From: |
Kai Torben Ohlhus |
Subject: |
Re: Condensing indices in a matrix |
Date: |
Tue, 22 Dec 2020 14:49:52 +0900 |
User-agent: |
Mozilla/5.0 (X11; Linux x86_64; rv:78.0) Gecko/20100101 Thunderbird/78.5.0 |
On 12/21/20 4:22 AM, Brett Green wrote:
On Tue, Dec 8, 2020 at 1:19 AM Kai Torben Ohlhus <k.ohlhus@gmail.com
<mailto:k.ohlhus@gmail.com>> wrote:
Was it possible to send some minimal working example how array M is
created? For u(x,y,z), maybe `sub2ind` [1] can be helpful.
HTH,
Kai
[1] https://octave.org/doc/v6.1.0/Advanced-Indexing.html
<https://octave.org/doc/v6.1.0/Advanced-Indexing.html>
Thank you and please excuse the slow response. Here is an MWE:
ind = @(T,n,m) 2*N*(m-1) + 2*(n-1) + T + 1;
M = zeros(2*N^2, 2*N^2);
for n=1:N
for m=1:N
M( ind(0,n,m) , ind(1,n,m) ) = 1;
M( ind(1,n,m) , ind(0,n,m) ) = 1;
M( ind(0,n,m) , ind(1,n+1,m) ) = 1;
M( ind(1,n+1,m) , ind(0,n,m) ) = 1;
M( ind(0,n,m) , ind(1,n,m+1) ) = 1;
M( ind(1,n,m+1) , ind(0,n,m) ) = 1;
end
end
While it's nice that sub2ind exists, I prefer to do the indexing myself
so that it is explicitly readable from the code, and so that I can
control it directly. The key problem for me is the assignment, which I
would like to do without for loops if possible.
It has been a longer project, but I hope it was worth the outcome.
Looking at your matrices M with the "spy" command, you can clearly see a
sparse band structure and create them more efficiently with my function
"createMfast" below.
Your function needs for dimension N = 50 (5100 x 5100 dense matrix)
about 4 seconds and about 200 MB on my machine. From N=100 the
computational cost become too huge.
The function "createMfast" uses sparse banded and boolean matrices (as
your values only deal with {0,1} and for N=500 (501000 x 501000 sparse
banded matrix) only 0.1 seconds creation time are required and it uses
only 17 MB of main memory. Thus you can grow larger if you want.
The following experiment shows a comparison up to dimension 50 and a
small timing benchmark.
````````````````````````````
clear all; clc;
function M = createM (N)
ind = @(T,n,m) 2*N*(m-1) + 2*(n-1) + T + 1;
M = zeros(2*N^2, 2*N^2);
for n=1:N
for m=1:N
M( ind(0,n,m) , ind(1,n,m) ) = 1;
M( ind(1,n,m) , ind(0,n,m) ) = 1;
M( ind(0,n,m) , ind(1,n+1,m) ) = 1;
M( ind(1,n+1,m) , ind(0,n,m) ) = 1;
M( ind(0,n,m) , ind(1,n,m+1) ) = 1;
M( ind(1,n,m+1) , ind(0,n,m) ) = 1;
end
end
endfunction
function M = createMfast (N)
d = false (2*N^2, 1);
d(1:2:end) = true;
dd = false (2*(N^2 + N), 6);
dd(1:length(d),1:3) = [d d d];
dd(2:length(d)+1,4) = d;
dd(4:length(d)+3,5) = d;
dd(end-length(d)+2:end,6) = d(1:end-1);
M = spdiags (dd, [-2*N-1 -3 -1 1 3 2*N+1], 2*(N^2 + N), 2*(N^2 + N));
endfunction
NN = [1:10, 15, 20, 30, 40, 50, 100, 200, 500];
T1 = [];
T2 = [];
for N = NN
if (N <= 50)
tic ();
M = createM (N);
T1(end+1) = toc ();
else
T1(end+1) = inf;
end
## Display the matrix structure.
#if (N == 10)
# figure ();
# spy (M);
#end
## Get diagonals for inspection.
#[a, b] = spdiags (M)
tic ();
MM = createMfast (N);
T2(end+1) = toc ();
if (N <= 50)
Mdiff = MM - sparse (M);
if (any (any (Mdiff)))
figure ();
spy (Mdiff);
end
end
end
figure ();
plot (NN, T1, 'b', NN, T2, 'r')
````````````````````````````
HTH,
Kai