Add benchmark to track Julia's array allocation performance relative …

…to Python

test/perf/arrayalloc/Makefile

@@ -0,0 +1,14 @@
+ifndef JULIABIN
+ JULIABIN = julia
+endif
+ifndef PYTHONBIN
+ PYTHONBIN = python3
+endif
+
+python:
+ $(PYTHONBIN) lucompletepiv.py
+
+julia:
+ $(JULIABIN) lucompletepiv.jl
+
+all: python julia

test/perf/arrayalloc/lucompletepiv.jl

@@ -0,0 +1,53 @@
+function lucompletepivCopy!(A)
+ n = size(A, 1)
+ rowpiv=zeros(Int, n-1)
+ colpiv=zeros(Int, n-1)
+ for k = 1:n-1
+ As = abs(A[k:n, k:n])
+ μ, λ = ind2sub(size(As), indmax(As))
+ μ += k-1; λ += k-1
+ rowpiv[k] = μ
+ A[[k,μ], 1:n] = A[[μ,k], 1:n]
+ colpiv[k] = λ
+ A[1:n, [k,λ]] = A[1:n, [λ,k]]
+ if A[k,k] ≠ 0
+ ρ = k+1:n
+ A[ρ, k] = A[ρ, k]/A[k, k]
+ A[ρ, ρ] = A[ρ, ρ] - A[ρ, k] * A[k, ρ]
+ end
+ end
+ return (A, rowpiv, colpiv)
+end
+
+function lucompletepivSub!(A)
+ n = size(A, 1)
+ rowpiv=zeros(Int, n-1)
+ colpiv=zeros(Int, n-1)
+ for k = 1:n-1
+ As = abs(sub(A, k:n, k:n))
+ μ, λ = ind2sub(size(As), indmax(As))
+ μ += k-1; λ += k-1
+ rowpiv[k] = μ
+ A[[k,μ], 1:n] = sub(A, [μ,k], 1:n)
+ colpiv[k] = λ
+ A[1:n, [k,λ]] = sub(A, 1:n, [λ,k])
+ if A[k,k] ≠ 0
+ ρ = k+1:n
+ A[ρ, k] = sub(A, ρ, k)/A[k, k]
+ A[ρ, ρ] = sub(A, ρ, ρ) - sub(A, ρ, k) * sub(A, k, ρ)
+ end
+ end
+ return (A, rowpiv, colpiv)
+end
+
+println("Julia version with copy slices")
+for n = [100, 250, 500, 1000]
+ A = randn(n,n)
+ @printf("size %4d matrix factorized in %6.3f seconds\n", n, mean([@elapsed lucompletepivCopy!(copy(A)) for i = 1:3]))
+end
+
+println("\nJulia version with view slices")
+for n = [100, 250, 500, 1000]
+ A = randn(n,n)
+ @printf("size %4d matrix factorized in %6.3f seconds\n", n, mean([@elapsed lucompletepivSub!(copy(A)) for i = 1:3]))
+end

test/perf/arrayalloc/lucompletepiv.py

@@ -0,0 +1,35 @@
+import numpy as np
+import time
+
+def lucompletepiv(A):
+ assert np.size(A, 0) == np.size(A, 1)
+ n = np.size(A, 1)
+ rowpiv = np.zeros(n-1, dtype=int)
+ colpiv = np.zeros(n-1, dtype=int)
+ for k in range(n-1):
+ Asub = abs(A[k:n, k:n])
+ mu, lam = np.unravel_index(np.argmax(Asub), np.shape(Asub))
+ mu, lam = mu + k, lam + k
+ rowpiv[k] = mu
+ A[[k, mu], :n] = A[[mu, k], :n]
+ colpiv[k] = lam
+ A[:n, [k, lam]] = A[:n, [lam, k]]
+ if A[k, k] != 0:
+ rho = slice(k+1, n)
+ A[rho, k] = A[rho, k]/A[k, k]
+ A[rho, rho] -= np.dot(np.reshape(A[rho, k], (n - (k + 1), 1)),
+ np.reshape(A[k, rho], (1, n - (k + 1))))
+ return (A, rowpiv, colpiv)
+
+
+for n in [100, 250, 500, 1000]:
+ tt = []
+ for rep in range(3):
+ A = np.random.randn(n, n)
+ t0 = time.time()
+ lucompletepiv(A.copy())
+ tt.append(time.time() - t0)
+
+ print("size %4d matrix factorized in %6.3f seconds" % (n, sum(tt)/3))
+
+print("")