@@ -1,3 +1,4 @@ | |||||
# Fast functions for LOQC | # Fast functions for LOQC | ||||
- Permanent (hardcoded up to 4x4 and Ryser) | - Permanent (hardcoded up to 4x4 and Ryser) | ||||
- Parallel permanent | - Parallel permanent | ||||
- No Cython dependency |
@@ -1 +1 @@ | |||||
from permanent import permanent |
@@ -4,32 +4,8 @@ import multiprocessing as mp | |||||
import numpy as np | import numpy as np | ||||
import lib | 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) | |||||
dimension=8 | |||||
real=np.random.uniform(0,1,(dimension, dimension)) | |||||
imag=np.random.uniform(0,1,(dimension, dimension)) | |||||
submatrix=real+1j*imag | |||||
print "OUTPUT: ", lib.permanent(submatrix) |
@@ -1,5 +1,5 @@ | |||||
#!/bin/bash | #!/bin/bash | ||||
rm lib/*.so | |||||
#rm lib/*.so | |||||
python ./setup.py build_ext --inplace && | python ./setup.py build_ext --inplace && | ||||
python ./run-tests.py | python ./run-tests.py |
@@ -11,7 +11,7 @@ setup(name = "loqcmath", | |||||
url = "https://www.peteshadbolt.co.uk", | url = "https://www.peteshadbolt.co.uk", | ||||
ext_modules = [ | ext_modules = [ | ||||
Extension( | Extension( | ||||
'perm_c', ['src/perm_c'], | |||||
'permanent', ['src/permanent.c'], | |||||
extra_compile_args=["-Ofast", "-march=native"]), | extra_compile_args=["-Ofast", "-march=native"]), | ||||
], | ], | ||||
@@ -41,57 +41,14 @@ PyMODINIT_FUNC initpermanent(void) { | |||||
(x1) * PyArray_STRIDES(py_F)[1]))) | (x1) * PyArray_STRIDES(py_F)[1]))) | ||||
#define F_shape(i) (py_F->dimensions[(i)]) | #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) { | static PyObject *permanent(PyObject *self, PyObject *args) { | ||||
// Declare variables. | // Declare variables. | ||||
npy_int64 d, i, j; | npy_int64 d, i, j; | ||||
npy_float64 output; | |||||
Py_complex sum; | |||||
PyArrayObject *submatrix; | PyArrayObject *submatrix; | ||||
// Parse variables. | // Parse variables. | ||||
if (!PyArg_ParseTuple(args, "ldllO!O!O!O!", | |||||
&threads, | |||||
&dt, | |||||
&steps, | |||||
&N, | |||||
&PyArray_Type, &submatrix)) { | |||||
if (!PyArg_ParseTuple(args, "O!", &PyArray_Type, &submatrix)) { | |||||
return NULL; | return NULL; | ||||
} | } | ||||
@@ -102,5 +59,6 @@ static PyObject *permanent(PyObject *self, PyObject *args) { | |||||
} | } | ||||
} | } | ||||
Py_RETURN_NONE; | |||||
PyObject *output=PyComplex_FromCComplex(sum); | |||||
return output; | |||||
} | } |