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

[eric]

Hey Pearu,

<snip Fortran compared to Blitz>
> > > Iterations took  6.12530398369  seconds.
> > > Fortran iterations took  3.15447306633  seconds.

I've rewritten the code in C++ as quickly as I know how and added a few
optimizations.  It has the same errors as the previous example, but, again,
is useful for benchmarking. (maybe Prabhu will fix my bugs in this code
also. :) The results are as follows:

Blitz indexing:
Iterations took  3.54500007629  seconds.
After re-write:
iterations took  1.77300000191  seconds.

With the re-write, C++ is in the same performance range as the fortran code.

My timestep loop is included below.  It uses pointer math instead of blitz
arrays to reduce the effort of array element lookups down to incrementing
pointers in the inner loop.  It isn't exactly pretty, but we were looking
for speed, not beauty.  I expect there are some ways to clean it up and
still get about the same speed.

I used the following optimizations:

extra_compile_args = ['-O3','-malign-double','-funroll-loops']

The difference between using and not using -funroll-loops was significant
(maybe 25%).

> Anyway, I don't think that people would use Fortran if
> the gain in speed would be only few percent. It must be higher than that.

I think it used to be, but the recent years have been kinder to C than
Fortran.  Among other things, market forces have pushed C/C++ compilers to
evolve more quickly than Fortran compilers simply because more people are
interested in fast C/C++ compilers.  The result is that C/C++ can be made to
execute algorithms at close to the same speed as Fortran in most cases -- at
least that has been my experience.  As for gcc and g77, I think only the
front ends of the compilers are different.  The same back-end is shared by
both.

Still, there are lots of good reasons besides speed for scientist to use
Fortran:

1. It is the language you know. (probably the main reason many still use it)
2. There are tons of scientific libraries available.
3. All your labs legacy code is in Fortran.
4. You don't have to mess with memory management. (debugging dangling
pointers stinks)
5. You don't have to resort to tricks for speed.
6. You just like the language better.

Of course, there are reasons not to use it also ( number 4 above can go on
both the "for" and "against" lists), but all the above reasons are valid.
Scientist are interested in the fastest way to results -- not in writing
elegant, full featured, re-usable, scriptable programs.  For one or more of
the above reasons, some still consider Fortran their best choice to optimize
this interest.  Other scientist choose Matlab or IDL or even C.  Hopefully
we can increase the share that see Python+Numeric+SciPy+f2py+weave as a good
choice also.

This thread definitely needs to get summarized and put into a web page.

see ya,
eric

----------------------------------------------------------------------------
------------------------------
# just replace the timeStep algorithm from the previous post with this one.

    def timeStep(self, dt = 0.0):
        dx2 = self.g.dx**2
        dy2 = self.g.dy**2
        dnr_inv = .5/(dx2+dy2)
        nx, ny = self.g.shape
        u = self.g.u
        code = """
               #line 44 "laplace.py"
               double tmp, err, diff;
               double *uc, *uu, *ud, *ul, *ur;
               for (int i=1; i<nx-1; ++i) {
                   uc = u_data+i*ny+1;
                   ur = u_data+i*ny+2;     ul = u_data+i*ny;
                   ud = u_data+(i+1)*ny+1; uu = u_data+(i-1)*ny+1;
                   for (int j=1; j<ny-1; ++j) {
                       tmp = *uc;
                       *uc = ((*ul + *ur)*dy2 +
                              (*uu + *ud)*dx2)*dnr_inv;
                       diff = *uc - tmp;
                       err += diff*diff;
                       uc++;ur++;ul++;ud++;uu++;
                   }
               }
               return_val = Py::new_reference_to(Py::Float(sqrt(err)));
               """
        #extra_compile_args =
['-O3','-funroll-loops','-march=i686','-malign-double']
        # compiler keyword only needed on windows with MSVC installed
        err = weave.inline(code, ['u', 'dx2', 'dy2', 'dnr_inv', 'nx','ny'],
                           compiler = 'gcc',
                           extra_compile_args =
['-O3','-malign-double','-funroll-loops'])
        return err