From 00f4f0a00a981f13ec4857946ab30f79050bfc8a Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Fri, 13 Jun 2014 17:11:40 +0200 Subject: [PATCH] More work on performance but broke the algorithm --- python/_project.pyx | 40 ++++++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/python/_project.pyx b/python/_project.pyx index f2c31d9..2d85d32 100644 --- a/python/_project.pyx +++ b/python/_project.pyx @@ -535,18 +535,6 @@ cdef DTYPE_t cube_integral(DTYPE_t u[3], DTYPE_t u0[3], int r[1]) nogil: return alpha_max -@cython.boundscheck(False) -@cython.cdivision(True) -cdef DTYPE_t mysum(DTYPE_t *v, int q) nogil: - cdef int i - cdef DTYPE_t s - - s = 0 - - for i in xrange(q): - s += v[i] - return s - @cython.boundscheck(False) @cython.cdivision(True) cdef DTYPE_t cube_integral_trilin(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t vertex_value[8], int err[1]) nogil: @@ -572,8 +560,10 @@ cdef DTYPE_t cube_integral_trilin(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t alpha_max = tmp_a j = i - if (u0[i] >= 0 || u0[i] <= 1) + if (u0[i] < 0 or u0[i] > 1): err[0] = 1 + with gil: + print("Bouh! u0[%d] = %lg" % (i, u0[i])) return 0 I = compute_projection(vertex_value, u, u0, alpha_max) @@ -589,7 +579,7 @@ cdef DTYPE_t cube_integral_trilin(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t @cython.boundscheck(False) cdef DTYPE_t integrator0(DTYPE_t[:,:,:] density, - DTYPE_t u[3], DTYPE_t u0[3], int u_delta[3], int iu0[3], int jumper[1], int err[1]) nogil: + DTYPE_t u[3], DTYPE_t u0[3], int u_delta[3], int iu0[3], int jumper[1], int err[1]) nogil except? 0: cdef DTYPE_t d d = density[iu0[0], iu0[1], iu0[2]] @@ -598,7 +588,7 @@ cdef DTYPE_t integrator0(DTYPE_t[:,:,:] density, @cython.boundscheck(False) cdef DTYPE_t integrator1(DTYPE_t[:,:,:] density, - DTYPE_t u[3], DTYPE_t u0[3], int u_delta[3], int iu0[3], int jumper[1], int err[1]) nogil: + DTYPE_t u[3], DTYPE_t u0[3], int u_delta[3], int iu0[3], int jumper[1], int err[1]) nogil except? 0: cdef DTYPE_t vertex_value[8] cdef DTYPE_t d cdef int a[3][2] @@ -635,7 +625,7 @@ cdef DTYPE_t integrator1(DTYPE_t[:,:,:] density, cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density, DTYPE_t a_u[3], DTYPE_t min_distance, - DTYPE_t max_distance, DTYPE_t[:] shifter, int integrator_id, int out_err[1]) nogil: + DTYPE_t max_distance, DTYPE_t[:] shifter, int integrator_id, int out_err[1]) nogil except? 0: cdef DTYPE_t u[3] cdef DTYPE_t ifu0[3] @@ -652,7 +642,7 @@ cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density, cdef int jumper[1] cdef DTYPE_t (*integrator)(DTYPE_t[:,:,:], - DTYPE_t u[3], DTYPE_t u0[3], int u_delta[3], int iu0[3], int jumper[1], int err[1]) nogil + DTYPE_t u[3], DTYPE_t u0[3], int u_delta[3], int iu0[3], int jumper[1], int err[1]) nogil except? 0 out_err[0] = 0 @@ -675,7 +665,7 @@ cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density, if (not ((iu0[i]>= 0) and (iu0[i] < N))): out_err[0] = 1 return 0 - + if (not (u0[i]>=0 and u0[i]<=1)): out_err[0] = 1 return 0 @@ -693,6 +683,8 @@ cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density, I0 += integrator(density, u, u0, u_delta, iu0, jumper, err) if err[0] == 1: + with gil: + print("Bah! error!") out_err[0] = 1 return 0 @@ -734,6 +726,7 @@ def line_of_sight_projection(DTYPE_t[:,:,:] density not None, u[1] = a_u[1] u[2] = a_u[2] + out_err[0] = 0 v = C_line_of_sight_projection(density, u, min_distance, @@ -754,9 +747,15 @@ cdef double _spherical_projloop(double theta, double phi, DTYPE_t[:,:,:] density return C_line_of_sight_projection(density, u0, min_distance, max_distance, shifter, integrator_id, out_err) +ctypedef fused sum_type: + cython.int + cython.float + npx.uint64_t + npx.uint32_t + @cython.boundscheck(False) -cdef npx.uint64_t _mysum(int[:] jobs) nogil: - cdef npx.uint64_t s +cdef sum_type _mysum(sum_type[:] jobs) nogil: + cdef sum_type s cdef npx.uint64_t N cdef int i @@ -809,6 +808,7 @@ def spherical_projection(int Nside, if omp_get_thread_num() == 0 and progress != 0 and (i%booster) == 0: with gil: p.update(_mysum(job_done)) + out_err[0] = 0 outm[i] = _spherical_projloop(theta[i], phi[i], density_view, min_distance, max_distance, shifter, integrator_id, out_err) job_done[omp_get_thread_num()] = i