guilhem_lavaux/cosmotool
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

More work on performance but broke the algorithm

Browse Source
This commit is contained in:
Guilhem Lavaux 2014-06-13 17:11:40 +02:00
parent 7a61ad7519
commit 00f4f0a00a
1 changed files with 20 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

40
python/_project.pyx
View File

@ -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