Author: Hakan Ardo <ha...@debian.org> Branch: extradoc Changeset: r4533:a6dc6af7135f Date: 2012-08-13 09:35 +0200 http://bitbucket.org/pypy/extradoc/changeset/a6dc6af7135f/
Log: python version of LU diff --git a/talk/iwtc11/benchmarks/benchmark.sh b/talk/iwtc11/benchmarks/benchmark.sh --- a/talk/iwtc11/benchmarks/benchmark.sh +++ b/talk/iwtc11/benchmarks/benchmark.sh @@ -57,4 +57,5 @@ $* ./runner.py $EXTRA_OPTS scimark.py SparseMatMult 1000 5000 262144 $* ./runner.py $EXTRA_OPTS scimark.py SparseMatMult 100000 1000000 1024 $* ./runner.py $EXTRA_OPTS scimark.py MonteCarlo 268435456 + $* ./runner.py $EXTRA_OPTS scimark.py LU 100 4096 fi diff --git a/talk/iwtc11/benchmarks/convolution/convolution.py b/talk/iwtc11/benchmarks/convolution/convolution.py --- a/talk/iwtc11/benchmarks/convolution/convolution.py +++ b/talk/iwtc11/benchmarks/convolution/convolution.py @@ -64,6 +64,9 @@ for x in xrange(self.width): yield x, y + def copy_data_from(self, other): + self.data[:] = other.data[:] + class NumpyArray(Array2D): def __init__(self, w, h): self.width = w diff --git a/talk/iwtc11/benchmarks/scimark.py b/talk/iwtc11/benchmarks/scimark.py --- a/talk/iwtc11/benchmarks/scimark.py +++ b/talk/iwtc11/benchmarks/scimark.py @@ -59,6 +59,10 @@ else: return self.dm1 * float(k); + def RandomMatrix(self, a): + for x, y in a.indexes(): + a[x, y] = self.nextDouble() + return a class ArrayList(Array2D): def __init__(self, w, h, data=None): @@ -80,6 +84,10 @@ else: self.data[idx] = val + def copy_data_from(self, other): + for l1, l2 in zip(self.data, other.data): + l1[:] = l2 + def SOR_execute(omega, G, num_iterations): for p in xrange(num_iterations): for y in xrange(1, G.height - 1): @@ -138,3 +146,42 @@ MonteCarlo_integrate(n) return 'MonteCarlo(%d)' % n +def LU_factor(A, pivot): + M, N = A.height, A.width + minMN = min(M, N) + for j in xrange(minMN): + jp = j + t = abs(A[j][j]) + for i in xrange(j + 1, M): + ab = abs(A[i][j]) + if ab > t: + jp = i + t = ab + pivot[j] = jp + + if A[jp][j] == 0: + raise Exception("factorization failed because of zero pivot") + + if jp != j: + A[j], A[jp] = A[jp], A[j] + + if j < M-1: + recp = 1.0 / A[j][j] + for k in xrange(j + 1, M): + A[k][j] *= recp + + if j < minMN-1: + for ii in xrange(j + 1, M): + for jj in xrange(j + 1, N): + A[ii][jj] -= A[ii][j] * A[j][jj] + +def LU(args): + N, cycles = map(int, args) + rnd = Random(7) + A = rnd.RandomMatrix(ArrayList(N, N)) + lu = ArrayList(N, N) + pivot = array('i', [0]) * N + for i in xrange(cycles): + lu.copy_data_from(A) + LU_factor(lu, pivot) + return 'LU(%d, %d)' % (N, cycles) diff --git a/talk/iwtc11/benchmarks/scimark/kernel.c b/talk/iwtc11/benchmarks/scimark/kernel.c --- a/talk/iwtc11/benchmarks/scimark/kernel.c +++ b/talk/iwtc11/benchmarks/scimark/kernel.c @@ -220,6 +220,7 @@ cycles *= 2; } /* approx Mflops */ + printf("LU: N=%d, cycles=%d\n", N, cycles); result = LU_num_flops(N) * cycles / Stopwatch_read(Q) * 1.0e-6; Stopwatch_delete(Q); diff --git a/talk/iwtc11/benchmarks/test_scimark.py b/talk/iwtc11/benchmarks/test_scimark.py --- a/talk/iwtc11/benchmarks/test_scimark.py +++ b/talk/iwtc11/benchmarks/test_scimark.py @@ -1,4 +1,5 @@ -from scimark import SOR_execute, Array2D, ArrayList, Random, MonteCarlo_integrate +from scimark import SOR_execute, Array2D, ArrayList, Random, MonteCarlo_integrate, LU_factor +from array import array from cffi import FFI import os @@ -11,16 +12,18 @@ void SOR_execute(int M, int N,double omega, double **G, int num_iterations); double MonteCarlo_integrate(int Num_samples); + int LU_factor(int M, int N, double **A, int *pivot); """) C = ffi.verify(""" #include <SOR.h> #include <Random.h> #include <MonteCarlo.h> + #include <LU.h> """, extra_compile_args=['-I' + os.path.join(os.getcwd(), 'scimark')], extra_link_args=['-fPIC'], extra_objects=[os.path.join(os.getcwd(), 'scimark', f) - for f in ['SOR.c', 'Random.c', 'MonteCarlo.c']]) + for f in ['SOR.c', 'Random.c', 'MonteCarlo.c', 'LU.c']]) class TestWithArray2D(object): Array = Array2D @@ -35,9 +38,40 @@ for x, y in b.indexes(): assert a[y][x] == b[x, y] + def test_copy_random_matrix(self): + rnd_C = C.new_Random_seed(7) + rnd_py = Random(7) + c_mat = C.RandomMatrix(20, 10, rnd_C) + py_mat = rnd_py.RandomMatrix(self.Array(10, 20)) + py_mat_cpy = self.Array(10, 20) + py_mat_cpy.copy_data_from(py_mat) + for x, y in py_mat.indexes(): + assert c_mat[y][x] == py_mat[x, y] == py_mat_cpy[x, y] + + class TestWithArrayList(TestWithArray2D): Array = ArrayList + def test_LU(self): + rnd = C.new_Random_seed(7) + for height in [10, 20, 30]: + for width in [10, 20, 30]: + c_mat = C.RandomMatrix(height, width, rnd) + c_pivot = ffi.new('int []', min(width, height)) + py_mat = self.Array(width, height, data=c_mat) + py_pivot = array('i', [0]) * min(width, height) + C.LU_factor(height, width, c_mat, c_pivot) + LU_factor(py_mat, py_pivot) + + for a, b in zip(c_pivot, py_pivot): + assert a == b + for x, y in py_mat.indexes(): + assert c_mat[y][x] == py_mat[x, y] + + + + + def test_random(): rnd_C = C.new_Random_seed(7) rnd_py = Random(7) _______________________________________________ pypy-commit mailing list pypy-commit@python.org http://mail.python.org/mailman/listinfo/pypy-commit