[Numpy-discussion] Numpy/Fortran puzzle (?)
Andrea Gavana
andrea.gavana at gmail.com
Tue Aug 19 10:23:10 EDT 2014
Hi All,
I have the following (very ugly) line of code:
all_results = np.asarray([transm_hist[date_idx, :, idx_main_set[date_idx]
]*main_flow[date_idx, 0:n_fluids] for date_idx in xrange(n_dates)])
where transm_hist.shape = (n_dates, n_fluids, n_nodes), main_flow.shape =
(n_dates, n_fluids) and idx_main_set is an array containing integer indices
with idx_main_set.shape = (n_dates, ) . The resulting variable
all_results.shape = (n_dates, n_fluids)
Since that line of code is relatively slow if done repeatedly, I thought
I'd be smart to rewrite it in Fortran and then use f2py to wrap the
subroutine. So I wrote this:
subroutine matmul(transm_hist, idx_main_set, main_flow, all_results, &
n_dates, n_fluids, n_nodes)
implicit none
integer ( kind = 4 ), intent(in) :: n_dates, n_fluids, n_nodes
real ( kind = 4 ), intent(in) :: transm_hist(n_dates, n_fluids,
n_nodes)
real ( kind = 4 ), intent(in) :: main_flow(n_dates, n_fluids)
integer ( kind = 4 ), intent(in) :: idx_main_set(n_dates)
real ( kind = 4 ), intent(out):: all_results(n_dates, n_fluids)
integer (kind = 4) i, node
do i = 1, n_dates
node = int(idx_main_set(i))
all_results(i, :) = transm_hist(i, 1:n_fluids, node)*main_flow(i,
1:n_fluids)
enddo
end
Unfortunately, it appears that I am not getting out quite the same
results... I know it's a bit of a stretch with so little information, but
does anyone have a suggestion on where the culprit might be? Maybe the
elementwise multiplication is done differently in Numpy and Fortran, or I
am misunderstanding what the np.asarray is doing with the list
comprehension above?
I appreciate any suggestion, which can also be related to improvement in
the code. Thank you in advance.
Andrea.
"Imagination Is The Only Weapon In The War Against Reality."
http://www.infinity77.net
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20140819/58e99dff/attachment.html>
More information about the NumPy-Discussion
mailing list