[julia-users] Modified Gram-Schmidt: getting the best out of Julia

[Ján Dolinský]

Hello,

I have been using Julia for couple of weeks now. I did some benchmarking 
like computing a dot product etc. in Julia and Octave but I soon realized 
that I am benchmarking BLAS/LAPACK implementations of elementary algebra 
functions rather then Julia or Octave itself.

To learn Julia further, I tried to write Modified Gram-Schmidt procedure 
(MGS) in Julia and Octave. First, I would like to write MGS in a serial 
fashion properly and then perhaps move on to a parallelized version (and 
release it as open source learning material).

My attempt in Octave is as follows:
function [X, R] = mgs(X)

[n p] = size(X);

R = eye(p);

for l=1:p
    R(l,l+1:p)=(X(:,l)'*X(:,l+1:p))/(X(:,l)'*X(:,l));
    X(:,l+1:p)=X(:,l+1:p)-X(:,l)*R(l,l+1:p);
end

X = rand(1000);
octave:18> tic; [Q, R] = mgs(X); toc;
Elapsed time is 3.99396 seconds.

I rewrite the above in Julia:
using ArrayViews

function mgs(X)

  n, p = size(X)
  R = eye(nvars)

  for l=1:nvars-1
    R[l,l+1:p]=(view(X,:,l)'*view(X,:,l+1:nvars))/sumabs2(view(X,:,l))
    X[:,l+1:p]=view(X,:,l+1:nvars)-view(X,:,l)*view(R,l,l+1:p)
  end

  X, R

end

X = rand(1000,1000);
@time Q,R = mgs(X);
elapsed time: 7.986772918 seconds (8085500312 bytes allocated, 34.69% gc 
time)

At the moment, my Julia implementation is a lot slower but I assume my 
understanding of Julia and julian writing style can be improved a lot :). 

Thanks in advance for any hints.

Best Regards,
Jan