[SciPy-User] Contribution to Performance Python
Ramon Crehuet
rcsqtc@iqac.csic...
Fri Nov 6 04:53:25 CST 2009
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
******************************************************
More information about the SciPy-User
mailing list