Skip to content

Commit

Permalink
#RFC: Updated spectralnorm.jl for ~30x speedup
Browse files Browse the repository at this point in the history
The previous version based on Scala ran in just under 5 minutes (~293s), while this version,
based on Javascript with tweaks, runs consistently in 10-12s. This puts Julia ahead of all
other dynamic languages and only behind C/Ada/Fortran using extremely low-level and in some
cases parallelization (which maybe we should explore). Julia FTW!
Let me know if we can make any other tweaks/optimizations. 
  • Loading branch information
quinnj committed Apr 9, 2013
1 parent 1230762 commit 06b8f79
Showing 1 changed file with 41 additions and 34 deletions.
75 changes: 41 additions & 34 deletions examples/shootout/spectralnorm.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,43 +3,45 @@
# spectral-norm benchmark
# https://shootout.alioth.debian.org/u32/performance.php?test=spectralnorm
#
# Based on the Scala program
#
# Based on the Javascript program
#

include("timing.jl")

A(i,j) = 1.0 / ((i+j)*(i+j+1.0)/2.0+i+1.0)

# infinite matrix function and its transpose
A(i, j) = 1.0 / ((i+j)*(i+j+1)/2 +i+1)
At(j, i) = 1.0 / ((i+j)*(i+j+1)/2 +i+1)

# w <- M * v
function mult(N, v, w, M)
for i = 1:N
s = 0.0
for j = 1:N
s += M(i-1, j-1) * v[j]
end
w[i] = s
end
function Au(u,w)
n = length(u)
for i = 1:n, j = 1:n
j == 1 && (w[i] = 0)
w[i] += A(i-1,j-1) * u[j]
end
end

function approximate(N)
u = ones(N)
v = ones(N)
w = ones(N)
for i = 1:10
mult(N, u, w, A)
mult(N, w, v, At)
mult(N, v, w, A)
mult(N, w, u, At)
end

vbv = 0.0
vv = 0.0
for i = 1:N
vbv += u[i] * v[i]
vv += v[i] * v[i]
end

sqrt(vbv / vv)
function Atu(w,v)
n = length(w)
for i = 1:n, j = 1:n
j == 1 && (v[i] = 0)
v[i] += A(j-1,i-1) * w[j]
end
end

function approximate(n)
u = ones(Float64,n)
v = zeros(Float64,n)
w = zeros(Float64,n)
vv = vBv = 0
for i = 1:10
Au(u,w)
Atu(w,v)
Au(v,w)
Atu(w,u)
end
for i = 1:n
vBv += u[i]*v[i]
vv += v[i]*v[i]
end
return sqrt(vBv/vv)
end

function spectralnorm(N)
Expand All @@ -52,3 +54,8 @@ else
N = 100
end
spectralnorm(N)

# @assert spectralnorm(100) == 1.274219991
# @timeit spectralnorm(500) "spectralnorm(n=500)"
# @timeit spectralnorm(3000) "spectralnorm(n=3000)"
# @timeit spectralnorm(5500) "spectralnorm(n=5500)"

0 comments on commit 06b8f79

Please sign in to comment.