Andreas Stahel wrote:
Dear Octave users
when testing code for a class a timing result on Octave 3.0.3 puzzled
me.
generate a very sparse, symmetric, positive definite matrix Anxny
(size 62500x62500) and time a few commands
x=Anxny\b -> 0.8 sec
R=chol(Anxny) -> 7.3 sec
x=R\(R'\b) -> 2.3 sec
[L,U,P]=splu(Anxny) -> 12 sec
I would expect the Cholesky back-substitution to be fastest and
cho(Anxny) to be comparable to Anxny\b !!
Would you happen to have hints on why this occurs
Try instead
[R, P, q] = chol (Anxny,'vector');
x(q) = R \ (R' \ b(q));
without the Q return value, chol can't use the sparsity preserving
column transformations.
With the script below I see
SolveTime = 0.19601
CholTime = 0.66004
CholSolveTime = 0.30802
Chol2Time = 0.20801
Chol2SolveTime = 0.028001
LUTime = 1.3121
and norm (x4.'-x3) equal to 1.5781e-11
D.
nx=100; ny=100;
hx=1/(nx+1); hy=1/(ny+1);
Dxx=spdiags([-ones(nx,1) 2*ones(nx,1) -ones(nx,1)],[-1 0 1],nx,nx)/(hx);
Dyy=spdiags([-ones(ny,1) 2*ones(ny,1) -ones(ny,1)],[-1 0 1],ny,ny)/(hy);
Anxny=kron(speye(ny),Dxx) + kron(Dyy,speye(nx));
b=ones(nx*ny,1);
t0=cputime;
x2=Anxny\b;
SolveTime=cputime()-t0
t0=cputime;
R=chol(Anxny);
CholTime=cputime()-t0
t0=cputime;
x3=R\(R'\b);
CholSolveTime=cputime()-t0
t0=cputime;
[R,P,q] = chol(Anxny,'vector');
Chol2Time=cputime()-t0
t0=cputime;
x4(q) = R\(R'\b(q));
Chol2SolveTime=cputime()-t0
t0=cputime;
[L,U,P]=splu(Anxny);
LUTime=cputime()-t0