[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: sparse matrix operations
From: |
Miroslaw Kwasniak |
Subject: |
Re: sparse matrix operations |
Date: |
Tue, 16 Nov 2004 21:48:11 +0100 |
User-agent: |
Mutt/1.5.6+20040907i |
On Tue, Nov 16, 2004 at 11:32:28AM -0700, E. Joshua Rigler wrote:
> For example, suppose I only want the diagonal of the product of two
> square matrices. As it stands right now in Octave (with octave-forge),
> C=A*B will result in a full C matrix, and may potentially use up more
> memory than I wish to allow.
Yes
octave:20> a=[1 3 1;0 0 0;1 1 2],b=a',c=b*a,
a =
1 3 1
0 0 0
1 1 2
b =
1 0 1
3 0 1
1 0 2
c =
2 4 3
4 10 5
3 5 5
> Since I only want the diagonal of C, I may
> as well only perform those operations that place a vector product at
> C(i,i), thereby saving CPU cycles and memory.
>
> I have a simple .m file that does this, but of course it is very slow in
> the Octave interpreter since it must loop over each of the desired
> non-zero elements of C. I was hoping somebody out there might have
> something similar, but in a .oct format, or could show me a trick for
> doing this with standard Octave operators, before I go and try to
> reinvent this wheel.
I don't know how smart is sparse implementation in octave but you can try
octave:21> sum(a'.*b)
ans =
2 10 5
In which c=a'.*b should't use more memory then smallest of one of a,b. It
doesn't save so mach memory as in place operation but save more cpu the
m-file loop.
My first test gives resonable results:
octave:49> sp_diag
N = 1000
mult sp=0 t=0.075325 e=0
sum sp=0 t=0.046321 e=0
loop sp=0 t=1.667849 e=0
mult sp=1 t=0.116958 e=0
sum sp=1 t=0.157900 e=0
loop sp=1 t=5.076967 e=0
sum metod is better then mult on full matrix, but on sparse execution time
isn't to bad :)
Mirek
file: sp_diag.m
-----------------------------------
N=1000
clear a0 a b c c0;
a0=[1 0 7;0 0 3;1 0 2];
a=a0;
a=[a a a a a a a a a a]';
a=[a a a a a a a a a a];
L=size(a,1);
for sp=0:1
if sp; a=sparse(a); end
b=a';
tic
for n=1:N;
c=diag(a*b);
end
printf('mult sp=%d t=%f e=%d\n',sp,toc,0);
c0=c;
tic
for n=1:N;
c=sum(a'.*b);
end
printf('sum sp=%d t=%f e=%d\n',sp,toc,any(c'~=c0));
tic
for n=1:N;
c=zeros(L,1);
for l=1:L
c(l)=sum(a(l,:).*b(:,l)');
end
end
printf('loop sp=%d t=%f e=%d\n',sp,toc,any(c~=c0));
end
-----------------------------------
-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.
Octave's home on the web: http://www.octave.org
How to fund new projects: http://www.octave.org/funding.html
Subscription information: http://www.octave.org/archive.html
-------------------------------------------------------------