[Numpy-discussion] Numpy/Fortran puzzle (?)

Xavier Barthelemy xabart at gmail.com
Thu Aug 21 15:26:20 EDT 2014


Hi Andrea
You should add a dimension argument in your Fortran code, also you should
write a f2py header in the same Fortran code.
Remember, numpy memory is C order wise.
You can specify in numpy the ordering of the matrices you pass when you
create them.
F2py automatically deals with matrices , but tends to mix dimensions when
there are too many matrices.
Manual declaration of dimensions should do the trick
Xavier
On 21/08/2014 2:07 am, "Andrea Gavana" <andrea.gavana at gmail.com> wrote:

> 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
>
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20140821/4cafce65/attachment.html>


More information about the NumPy-Discussion mailing list