Subject: [FFC-dev] Benchmark results for new BLAS mode

From: Anders Logg <[email-masked]>
Date: Mon, 10 Oct 2005 10:52:29 -0500

Here are the benchmark results/evaluation of the new BLAS mode as
promised. Full log files attached.

The test compares FFC default mode with FFC BLAS mode for Poisson
in 2D and 3D for q = 1, 2, ..., 8 and different levels of optimization
during the compilation with g++: -O0, -O1, -O2.

Here are some conclusions:

- Optimization (-O1, -O2) has little effect on run-time performance
  in FFC default mode. These are numbers for Poisson degree 5 in 3D:
 
          | Run-time  | Compile-time (g++)
      ------------------------------------
      -O0 | 3.417e-05 | 4.596e+00
      -O1 | 3.116e-05 | 5.725e+01
      -O2 | 3.788e-05 | 6.929e+01

  With -O1, the code runs about 10% faster and with -02, the code
  is actually slower, possibly due to increase in size of object code.
    
  On the other hand, the level of optimization has quite a big
  influence on the compile-time (g++). Compile-time increases with
  more than a factor 10 when going from -O0 to -O1 and even more when
  going to -O2.

  Compiling with optimization also increases the memory usage for g++
  and I can't go beyond q = 5 for Poisson on my 2GB machine with -O2
  and g++ 4.0.2.
  
  Conclusion: -O1 may improve the run-time performance, but only
  little and with a big penalty in compile-time.
 
- Optimization (-O1, -O2) has no effect on run-time performance in
  BLAS mode, at least not for Poisson, where the computation of the
  geometry tensor is only a small contribution to the overall work.
  These are numbers for Poisson degree 5 in 3D:

          | Run-time  | Compile-time (g++)
      ------------------------------------
      -O0 | 6.787e-05 | 1.070e+00
      -O1 | 6.817e-05 | 1.305e+00
      -O2 | 6.805e-05 | 1.448e+00

  Also note that compile-time (g++) is much the same with different
  levels of optimizations in BLAS mode.

  Conclusion: the level of optimization has little effect in BLAS
  mode.

- Comparing now FFC default mode with FFC BLAS mode, we have the
  following timings for Poisson degree 8 in 3D:

               | BLAS mode | default MODE
      -----------------------------------
      FFC      | 6.039e+01 | 6.807e+01
      g++      | 1.383e+00 | 9.053e+01
      Run-time | 5.907e-04 | 1.333e-03

  Compile-time with FFC is about 1 minute in both modes.
  Compile-time with g++ is a factor 6 better in BLAS mode.
  Run-time is more than a factor 2 better in BLAS mode.

  When examining the influence of the order q, one finds that BLAS
  mode is always equal or faster than default mode in terms of g++
  compile-time. In 3D, the difference is significant for q > 3.

  Concerning run-time performance, FFC default mode is faster than
  BLAS mode for small q, and BLAS is faster for large q. For Poisson
  in 3D, the break-even point is at q = 6.

  Conclusion: the FFC blas mode will in some cases generate faster
  code, but it's first and foremost an option that should be
  considered to decrease the compile-time (with g++).

- Overall conclusions:

  It's unecessary to compile with -O1 or -O2. We have -O2 by
  default in DOLFIN, but should probably consider to change to
  -O0, at least for the compilation of forms. There is probably
  a huge improvement in compile-time for Navier-Stokes if we switch
  to -O0.

  BLAS can be an option to reduce compile-time and possibly to
  improve the run-time for high-order forms.

- Finally, note that these benchmarks are for Poisson, where the main
  part of the work is the computation of the element tensor (doing the
  tensor product). For other forms, computing the geometry tensor can
  dominate, and then BLAS mode won't help, but compiling with -O0 may.

/Anders

---

From: Johan Jansson <[email-masked]>
Date: Mon, 10 Oct 2005 19:59:08 +0200

On Mon, Oct 10, 2005 at 10:52:29AM -0500, Anders Logg wrote:

...

>   dominate, and then BLAS mode won't help, but compiling with -O0 may.

Ok, these results look reasonable. BLAS does blocking (reordering
computations to utilize the cache better), and FFC doesn't do that
(yet at least), so it's reasonable that BLAS will be faster the larger
the tensors get (at least that's what I think is happening). BLAS is
compiled by a Fortran compiler though, so it could gain a factor 2-3
there.

I still think code generation is the way to go though. The BLAS mode
probably won't ever be able to exploit the sparsity of the reference
tensor fully (perhaps on a coarse level), and won't be able to use
other Ferari-style optimizations. And if it's really necessary, FFC
could generate Fortran code. It's probably also a significant benefit
to generate code for the mappings as well. In the long term, I don't
see how you're ever going to be able to beat code generation for
runtime speed.

Still, why not have both modes. Evidently the code generation mode
suffers from excessive compile time, so clearly the BLAS mode has an
advantage there.

  Johan

---

From: "Robert C. Kirby" <[email-masked]>
Date: Mon, 10 Oct 2005 16:54:04 -0500

> computations to utilize the cache better), and FFC doesn't do that
BLAS really only gets you blocking advantages for level 3.  I suspend  
all judgments on these benchmarks until level 3 is in place.  The  
only big difference either way seems to be compile time.

> the tensors get (at least that's what I think is happening). BLAS is
This will be very pronounced for level 3 blas.

> compiled by a Fortran compiler though, so it could gain a factor 2-3

Atlas BLAS may be C code.  Don't worry about the language issues.   
It's a black box.

> other Ferari-style optimizations. And if it's really necessary, FFC
FErari only gets you so much.  The factor of 3 on floating-point  
operations I gain on the weighted Laplacian will be more than  
compensated for by level 3 BLAS (factor of 10 or so compared to level  
2).

The best way to go is
i.) Figure out block structure first.
ii.) See whether FErari or level 3 wins after you do the coarse-level  
block structure.  This will depend on the form, the polynomial  
degree, how well ferari does, etc.

> see how you're ever going to be able to beat code generation for
And the build system gets even more complicated :)
Actually, good C code shouldn't lose by a factor of 2-3.  More like  
5-10% atworst.

Rob Kirby

"Mathematical software should be mathematical."

---

From: Johan Jansson <[email-masked]>
Date: Tue, 11 Oct 2005 12:00:10 +0200

On Mon, Oct 10, 2005 at 04:54:04PM -0500, Robert C. Kirby wrote:

...

> "Mathematical software should be mathematical."

All I'm saying is that the theoretical top speed for the generated
code mode is faster than for the BLAS mode. Do you disagree with this
statement?

The reason this is so is because the generated code mode can
potentially use all the optimizations that BLAS does (blocking or
clever cache usage, Fortran) but BLAS cannot use all the optimizations
that the generated code mode uses (Ferari-style optimizations such as
skipping known zeros).

In the short term though, I don't think anybody is going to want to
spend time on these types of optimizations for the generated code mode
(premature optimization is the root of all evil etc.). But in the long
term someone might consider it worthwhile, especially if this type of
assembly becomes an established tool.

So here's my suggested plan: let's not lock ourselves into BLAS, keep
both methods for the foreseeable future. BLAS is a great tool, and may
be faster right now in runtime for certain forms and elements, but the
generated code mode will be faster than BLAS for all forms and all
elements if the BLAS-style optimizations are included (which can
certainly be done).

On a side note, there can be a significant difference in top speed
performance between code produced by Fortran and C compilers. We had a
quite trivial example in a course to demonstrate why Fortran was so
popular for computing, and Fortran turned out to be more than 4 times
faster than C, given full optimization for both the Fortran and C
compilers (Sun compilers on an UltraSPARC). This was due to aliasing,
a Fortran compiler can assume that all variables are distinct in
memory, while a C compiler cannot.

Here's the example in C:

void horner(double* px, double* x, double* coeff, int n)
{
  int           j;
  double        xj;

  for(j = 0; j < n; j++)
    {
      xj = x[j];
      px[j] = coeff[0] + xj * (coeff[1] + xj * (coeff[2] + xj * (coeff[3] + xj * coeff[4])));
    }
}

and Fortran 90:

subroutine horner( px, x, coeff, n )
  implicit none
  integer           n, j
  double precision  px(n), x(n), coeff(5), xj;
  
  do j = 1, n
     xj = x(j)
     px(j) = coeff(1) + xj * (coeff(2) + xj * (coeff(3) + xj * (coeff(4) + xj * coeff(5))))
  end do
  
end subroutine horner

I don't think Fenics should focus on these types of issues though, it
should be enough to compare abstractions of methods, we shouldn't all
have to become Fortran or assembler hackers to prove our points (I
certainly am not, and have no ambition of becoming one).

  Johan

---

From: Robert C Kirby <[email-masked]>
Date: Tue, 11 Oct 2005 07:23:06 -0500

> that the generated code mode uses (Ferari-style optimizations such as
I see you what you mean, and agree somewhat, but  I am skeptical. The 
irregularity introduced by FErari optimizations destroys 
spatial/temporal locality.  The reduced flop count comes at a price 
that may not be able to use the machine as efficiently as the very 
regular computation patterns in BLAS -- you can do fewer flops, but 
they are slower.  Second, as you indicate below, nobody wants to 
recreate something like ATLAS or FLAME.  Besides, these can take a very 
long time just to generate a single matrix-matrix multiply with no 
particular structure (so if you tile one part you know how to tile all 
of it).

> elements if the BLAS-style optimizations are included (which can

I completely agree with the following:
i.) There exist forms for which fine-grained FErari optimizations are 
the way to go.  These are ones with "short G" -- the geometric tensor 
doesn't have much data in it.
ii.) There exist forms for which level 3 BLAS are the way to go.  They 
are the other end of the spectrum.
iii.) For vector-valued Lagrange elements, it will (almost) always pay 
off to utilize block structure; even if there are no zero blocks, the 
nonzero blocks could have lots of structure and repetition that we 
should exploit.

> memory, while a C compiler cannot.

I'm not sure how much this propogates to other compilers, but aliasing 
can be a problem.  On the other hand, since we're generating code and 
know what n is up front, we could generate a very specific struct with 
n fields rather than use an array.  This may explain Anders' 
observation that forming G slowed down when he started putting it in an 
array.  Generating structs is not incompatible with BLAS -- you pass in 
a pointer to the head of the struct, cast to be a double.

Also, most FORTRAN compilers don't prevent you from doing aliasing -- 
you can pass in the same pointer to two arguments, and it won't 
complain at you.  It optimizes in an unsafe way, so we have to be 
careful in the code generated (not that you don't in C).