Sparse Toeplitz matrices
Marco Caliari
marco.caliari at univr.it
Fri Apr 10 10:33:13 CDT 2009
Dear maintainers,
I usually use finite differences (sparse) matrices and I find quite
convenient to use the function toeplitz.m to generate them. However, I
discovered that
toeplitz(sparse([-2,1,zeros(1,N-2)]))
can be quite slow if N is big.
I wrote the attached function which should be used only to generate sparse
Toeplitz matrices. Against Octave 3.1.54 I can see the following speedups:
octave-3.1.54:1> N = 5000;
octave-3.1.54:2> tic,A = toeplitz(sparse([-2,1,zeros(1,N-2)]));,toc
Elapsed time is 4.78363 seconds.
octave-3.1.54:3> tic,B = mytoeplitz(sparse([-2,1,zeros(1,N-2)]));,toc
Elapsed time is 0.095167 seconds.
octave-3.1.54:4> A-B
ans =
Compressed Column Sparse (rows = 5000, cols = 5000, nnz = 0 [0%])
octave-3.1.54:5>
tic,A=toeplitz(sparse([0,-1,zeros(1,N-2)]),sparse([0,1,zeros(1,N-2)]));,toc
Elapsed time is 3.8 seconds.
octave-3.1.54:6>
tic,B=mytoeplitz(sparse([0,-1,zeros(1,N-2)]),sparse([0,1,zeros(1,N-2)]));,toc
Elapsed time is 0.021 seconds.
octave-3.1.54:7> A-B
ans =
Compressed Column Sparse (rows = 5000, cols = 5000, nnz = 0 [0%])
If interested, I can try to make a patch against toeplitz.m.
Best regards,
Marco
-------------- next part --------------
function A = mytoeplitz(c,r)
if (nargin == 1)
r = c;
c = conj (c);
c(1) = conj (c(1));
endif
if (r (1) != c (1))
warning ("toeplitz: column wins diagonal conflict");
r(1) = 0;
endif
n = length(c);
c = c(:).';
r = r(:).';
cidx = find(c);
ridx = find(r);
A = spdiags(repmat(full(c(cidx)),n,1),1-cidx,n,n)+...
spdiags(repmat(full(r(ridx)),n,1),ridx-1,n,n);
More information about the Octave-maintainers
mailing list