From 6f42ae28978af2ba75258206104a4b776b29a82a Mon Sep 17 00:00:00 2001 From: Pete Shadbolt Date: Wed, 22 Oct 2014 06:14:47 +1100 Subject: [PATCH] First commit --- .gitignore | 57 ++++++++++++++++++++++++++ LICENSE | 22 ++++++++++ README.md | 3 ++ run-tests.py | 35 ++++++++++++++++ run-tests.sh | 5 +++ setup.py | 18 ++++++++ src/permanent.c | 106 ++++++++++++++++++++++++++++++++++++++++++++++++ 7 files changed, 246 insertions(+) create mode 100644 .gitignore create mode 100644 LICENSE create mode 100644 README.md create mode 100644 run-tests.py create mode 100755 run-tests.sh create mode 100644 setup.py create mode 100644 src/permanent.c diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..585129b --- /dev/null +++ b/.gitignore @@ -0,0 +1,57 @@ +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] + +# C extensions +*.so + +# Distribution / packaging +.Python +env/ +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +*.egg-info/ +.installed.cfg +*.egg + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.coverage +.cache +nosetests.xml +coverage.xml + +# Translations +*.mo +*.pot + +# Django stuff: +*.log + +# Sphinx documentation +docs/_build/ + +# PyBuilder +target/ + +# PYC +*.pyc diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..c63dcbc --- /dev/null +++ b/LICENSE @@ -0,0 +1,22 @@ +The MIT License (MIT) + +Copyright (c) 2014 J. Pete Shadbolt + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. + diff --git a/README.md b/README.md new file mode 100644 index 0000000..55a8bfe --- /dev/null +++ b/README.md @@ -0,0 +1,3 @@ +# Fast functions for LOQC +- Permanent (hardcoded up to 4x4 and Ryser) +- Parallel permanent diff --git a/run-tests.py b/run-tests.py new file mode 100644 index 0000000..cce6d60 --- /dev/null +++ b/run-tests.py @@ -0,0 +1,35 @@ +import os +import time +import multiprocessing as mp +import numpy as np +import lib + + +def test_fn(perm_fn, name, dimension=6, steps=1000, dt=1e-3, bodies=101, threads=1): + # Test the speed of the evolution function. + sreal=np.random.uniform(0,1,(dimension, dimension)) + scomp=np.random.uniform(0,1,(dimension, dimension)) + submatrix=sreal+1j*scomp + + t0 = time.time() + for i in range(steps): + p=perm_fn(submatrix) + t1 = time.time() + + print "{0} ({1}): {2} steps/sec".format(name, threads, int(steps / (t1 - t0))) + + +if __name__ == "__main__": + d=4 + steps=1000 + + # Test on one core + test_fn(lib.perm_python, "1 core, C", steps=steps, threads=1) + + # Test on one core + test_fn(lib.perm_c, "1 core, C", steps=steps, threads=1) + + # Test on four cores + test_fn(lib.perm_c, "All cores, C", steps=steps/threads, threads=threads) + + diff --git a/run-tests.sh b/run-tests.sh new file mode 100755 index 0000000..add786a --- /dev/null +++ b/run-tests.sh @@ -0,0 +1,5 @@ +#!/bin/bash + +rm lib/*.so +python ./setup.py build_ext --inplace && +python ./run-tests.py diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..53caa55 --- /dev/null +++ b/setup.py @@ -0,0 +1,18 @@ +#!/usr/bin/env python + +from distutils.core import setup, Extension + +setup(name = "loqcmath", + version = "1.0", + description = "Fast maths for LOQC", + author = "Pete Shadbolt", + author_email = "pete.shadbolt@gmail.com", + maintainer = "pete.shadbolt@gmail.com", + url = "https://www.peteshadbolt.co.uk", + ext_modules = [ + Extension( + 'perm_c', ['src/perm_c'], + extra_compile_args=["-Ofast", "-march=native"]), + ], + +) diff --git a/src/permanent.c b/src/permanent.c new file mode 100644 index 0000000..334411e --- /dev/null +++ b/src/permanent.c @@ -0,0 +1,106 @@ +/* Computes the permanent, given a numpy array */ + +#include +#include +#include + +// Forward function declaration. +static PyObject *permanent(PyObject *self, PyObject *args); + +// Boilerplate: method list. +static PyMethodDef methods[] = { + { "permanent", permanent, METH_VARARGS, "Compute the permanent"}, + { NULL, NULL, 0, NULL } /* Sentinel */ +}; + +// Boilerplate: Module initialization. +PyMODINIT_FUNC initpermanent(void) { + (void) Py_InitModule("permanent", methods); + import_array(); +} + +/***************************************************************************** + * Array access macros. * + *****************************************************************************/ +#define m(x0) (*(npy_float64*)((PyArray_DATA(py_m) + \ + (x0) * PyArray_STRIDES(py_m)[0]))) +#define m_shape(i) (py_m->dimensions[(i)]) + +#define r(x0, x1) (*(npy_float64*)((PyArray_DATA(py_r) + \ + (x0) * PyArray_STRIDES(py_r)[0] + \ + (x1) * PyArray_STRIDES(py_r)[1]))) +#define r_shape(i) (py_r->dimensions[(i)]) + +#define v(x0, x1) (*(npy_float64*)((PyArray_DATA(py_v) + \ + (x0) * PyArray_STRIDES(py_v)[0] + \ + (x1) * PyArray_STRIDES(py_v)[1]))) +#define v_shape(i) (py_v->dimensions[(i)]) + +#define F(x0, x1) (*(npy_float64*)((PyArray_DATA(py_F) + \ + (x0) * PyArray_STRIDES(py_F)[0] + \ + (x1) * PyArray_STRIDES(py_F)[1]))) +#define F_shape(i) (py_F->dimensions[(i)]) + + +/***************************************************************************** + * compute_F * + *****************************************************************************/ +static inline void compute_F(npy_int64 N, + PyArrayObject *py_m, + PyArrayObject *py_r, + PyArrayObject *py_F) { + npy_int64 i, j; + npy_float64 sx, sy, Fx, Fy, s3, tmp; + + // Set all forces to zero. + for(i = 0; i < N; ++i) { + F(i, 0) = F(i, 1) = 0; + } + + // Compute forces between pairs of bodies. + for(i = 0; i < N; ++i) { + for(j = i + 1; j < N; ++j) { + sx = r(j, 0) - r(i, 0); + sy = r(j, 1) - r(i, 1); + + s3 = sqrt(sx*sx + sy*sy); + s3 *= s3 * s3; + + tmp = m(i) * m(j) / s3; + Fx = tmp * sx; + Fy = tmp * sy; + + F(i, 0) += Fx; + F(i, 1) += Fy; + + F(j, 0) -= Fx; + F(j, 1) -= Fy; + } + } +} + +static PyObject *permanent(PyObject *self, PyObject *args) { + // Declare variables. + npy_int64 d, i, j; + npy_float64 output; + PyArrayObject *submatrix; + + // Parse variables. + if (!PyArg_ParseTuple(args, "ldllO!O!O!O!", + &threads, + &dt, + &steps, + &N, + &PyArray_Type, &submatrix)) { + return NULL; + } + + // Compute the permanent + for(i = 0; i < d; ++i) { + for(j = 0; j