[SciPy-dev] Performance Python

[Ramon Crehuet]

Hi all,
A few days ago I posted a contribution to the Performace Python page at
the Scipy-users list. This list is probably a better place to send it.
I would like to add that to the wiki. Can I be authorized for that?
Best regards,
Ramon



Subject: Contribution to Performance Python
From: Ramon Crehuet <rcsqtc at iqac.csic.es>
Date: Fri, 06 Nov 2009 11:53:25 +0100
To: scipy-user at scipy.net

Hi all,
After reading the "Performance Python" page at:
http://www.scipy.org/PerformancePython?action=show
I thought some code with Fortran 90/95 was missing, in partuclar
considering its useful array features. So I have written a couple of
examples with Fortran 90 arrays and with the Fortran 95 forall
construct. Both are nicer (to me  :-)  ) than the FORTRAN77 loops and also
faster! In my PC a 1000x1000 array gives:

Doing 100 iterations on a 1000x1000 grid
numeric took 5.86 seconds
fortran77 took 3.53 seconds
fortran90-arrays took 1.58 seconds
fortran95-forall took 1.58 seconds
slow (1 iteration) took 9.13 seconds
100 iterations should take about 913.000000 seconds

If this is interesting to the community, who should I contact to have
this included in the scipy web page?
Cheers,
Ramon


This are the two new subroutines. I can send the modified laplace.py to
whoever wants it.
******************************************************
! File flaplace90_arrays.f90

subroutine timestep(u,n,m,dx,dy,error)

implicit none

real (kind=8), dimension(0:n-1,0:m-1), intent(inout):: u


real (kind=8), intent(in) :: dx,dy
real (kind=8), intent(out) :: error
integer, intent(in) :: n,m
real (kind=8), dimension(0:n-1,0:m-1) :: diff
real (kind=8) :: dx2,dy2,dnr_inv

!f2py intent(in) :: dx,dy
!f2py intent(in,out) :: u
!f2py intent(out) :: error
!f2py intent(hide) :: n,m

dx2 = dx*dx
dy2 = dy*dy
dnr_inv = 0.5d0 / (dx2+dy2)

diff=u

u(1:n-2, 1:m-2) = ((u(0:n-3, 1:m-2) + u(2:n-1, 1:m-2))*dy2 + &
                         (u(1:n-2,0:m-3) + u(1:n-2, 2:m-1))*dx2)*dnr_inv

error=sqrt(sum((u-diff)**2))

end subroutine
******************************************************
! File flaplace95_forall.f90

subroutine timestep(u,n,m,dx,dy,error)

implicit none

real (kind=8), dimension(0:n-1,0:m-1), intent(inout):: u


real (kind=8), intent(in) :: dx,dy
real (kind=8), intent(out) :: error
integer, intent(in) :: n,m
real (kind=8), dimension(0:n-1,0:m-1) :: diff
real (kind=8) :: dx2,dy2,dnr_inv
integer :: i,j

!f2py intent(in) :: dx,dy
!f2py intent(in,out) :: u
!f2py intent(out) :: error
!f2py intent(hide) :: n,m

dx2 = dx*dx
dy2 = dy*dy
dnr_inv = 0.5d0 / (dx2+dy2)

diff=u

forall (i=1:n-2,j=1:m-2)
    u(i,j) = ((u(i-1,j) + u(i+1,j))*dy2+(u(i,j-1) + u(i,j+1))*dx2)*dnr_inv
end forall

error=sqrt(sum((u-diff)**2))

end subroutine

******************************************************