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.gav...@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@scipy.org > http://mail.scipy.org/mailman/listinfo/numpy-discussion > >
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
