Fortran comparison Re: [SciPy-dev] scipy.weave versus simple C++.

[Pearu Peterson]

Hi Prabhu and Eric,

Just from a curiousity I tried your test also using Fortran and here are
the results (now comparing only cxx and Fortran):

500x500, 100 iterations, PII 400 laptop, Linux, gcc-2.95.4

Iterations took  6.12530398369  seconds.
Fortran iterations took  3.15447306633  seconds.

So, using Fortran you get speed up approx. _two_ times compared to cxx!!!
If one would use optimized blas and laplack functions, the speed up could
be even greater.

Below is the testing code that I used. Here is how I used it, step by
step, to illustrate that using Fortran from python is very simple:

1) Build fortran wrapper:
  f2py -c flaplace.f -m flaplace          # here you need f2py from CVS

As a result you'll get flaplace.so in the current directory.

2) Run tests:
$ python laplace.py 
Iterations took  6.51088702679  seconds.
Error:  38047.7567709
Fortran iterations took  3.22211503983  seconds.
Error:  21778.5117188

Regards,
	Pearu

--------------------
c File flaplace.f
      subroutine timestep(u,n,m,dx,dy,error)
      double precision u(n,m)
      double precision dx,dy,dx2,dy2,dnr_inv,tmp,err,diff
      integer n,m,i,j
cf2py intent(in) :: dx,dy
cf2py intent(in,out) :: u
cf2py intent(out) :: error
cf2py intent(hide) :: n,m
      dx2 = dx*dx
      dy2 = dy*dy
      dnr_inv = 0.5d0 / (dx2+dy2)
      error = 0
      do 200,i=2,n-1
         do 100,j=2,m-1
            tmp = u(i,j)
            u(i,j) = ((u(i-1,j) + u(i+1,j))*dy2+
     &           (u(i,j-1) + u(i,j+1))*dx2)*dnr_inv 
            diff = u(i,j) - tmp
            error = error + diff*diff
 100     continue
 200  continue
      error = sqrt(error)
      end

-------------------
# modifications to laplace.py

import flaplace

    # New method to LaplaceSolver class
    def ftimeStep(self,dt = 0.0):
        u,err = flaplace.timestep(self.g.u,self.g.dx,self.g.dy)
        return err

    # Slightly modified solve method:
    def solve(self, n_iter=0, eps=1e-16,fortran = 0):
        if fortran:
            timeStep = self.ftimeStep
        else:
            timeStep = self.timeStep
        err = timeStep()
        #print err
        count = 1
        while err > eps:
            if n_iter and (count > n_iter):
                return err;
            err = timeStep()
            #print err
            count += 1
        return count

    # Slightly modified runner:
    t1 = time.time()
    err = s.solve(100,fortran)
    t2 = time.time()
    print "Iterations took ", t2 - t1, " seconds."
    print "Error: ", err
    t1 = time.time()
    err = s.solve(100,fortran=1)
    t2 = time.time()
    print "Fortran iterations took ", t2 - t1, " seconds."
    print "Error: ", err

---------------
EOF message