[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,
  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,


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.


"Imagination Is The Only Weapon In The War Against Reality."
-------------- 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