help-octave
[Top][All Lists]
Advanced

[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



reply via email to

[Prev in Thread] Current Thread [Next in Thread]