LQCD on the CELL BE

Main achievement

	Full implementation of the Dirac Matrix Inversion on the CELL Engine (PPU + SPEs)
	Effective integration into the tmLDCQ package
	Succesful tests on one CELL Blade (Moving to a cluster of CELLs should be «natural»)

Illustrative results

A 32×16³ configuration solved using the CGR algorithm on the PPU in 138 seconds

The same configuration and algorithm on the (PPU + 8 SPEs) double precision in 4.58 seconds

CGR:	57 iterations in	4.58 s (QS20)	3.68 s (QS22)	28.34s s (Intel 2.83Ghz)
CG:	685 iterations in	51.70 s (QS20)	38.80 s (QS22)	362.45 s (Intel 2.83Ghz)

Compare to the PPU, a speedup around 30 was obtained (also valid per iteration)

Material

	Our CELL executable version of the adapted tmLQCD is available here QS20 or QS22
	Instructions for running the code are the following: download the code (cell_invert.exe), look for a data file and prepare a configuration file (more info and data here) set the LD_LIBRARY_PATH to the /lib directory of the CELL SDK run the program with the sequence ./cell_invert -f config_file
	Online documentation of the tmLQCD package and the DD-HMC package
	More information about the problem, project and data could be found on the PetaQCD webpage.

SOLVING THE DIRAC EQUATION ON THE CELL
C. TADONKI (claude.tadonki@u-psud.fr). PetaQCD project, funded by the french agency ANR.

	Cell is a heterogeneous chip multiprocessor consisting of a 64-bit Power core (Power Processing Unit), augmented with 8 specialized SIMD co-processors (Synergistic Processor Units)
	Peak performance: single precision 256 GFlops / double precision 26 GFlops
	Hardware constraints and special programming model (skillful programmers!)

	In LQCD, we need to solve large-scale and numerically intractable linear system of equations
	Having a very powerful node is crucial for an efficient computation system, thus we tried the CELL

SOLVING THE DIRAC EQUATION ON THE CELL
C. TADONKI (claude.tadonki@u-psud.fr). PetaQCD project, funded by the french agency ANR.

What was done

A vector-SPU implementation of the critical operations

The linear algebra operation so far revisited are:

Hopping-Matrix (the biggest cake!)
Scalar product of two spinors
Square norm of a spinor

Multiplying a spinor by a complex
Assigning a spinor to another
Accumulating a weighted serie of spinors

Efficient iterative data communication mechanism and a total DMA/computation overlapping

Operational integration and connexion from the original tmLQD package

Why such a speedup

	Fully optimized implementation on the SPU (mainly SIMD)
	Memory latency hidden by overlapping DMA with computation (intensive DMA lists)
	Keeping SPU threads alive during «hot iterations»
	Good scaling using 8 SPUs (see one example with Hopping Matrix here)

What are the key points

	Optimizing each iteration at the best as we use iterative methods
	Really care about memory overhead as we manipulate huge arrays everytime and everywhere
	Whenever possible, aggregate some operations in order to break the sequentiality

SOLVING THE DIRAC EQUATION ON THE CELL
C. TADONKI (claude.tadonki@u-psud.fr). PetaQCD project, funded by the french agency ANR.

SIMD Illustration

The multiplication of a 9×9 complex matrix by a complex vector (see the scalar version) is written using SPU intrinsics as follows:

#define _su3_multiply(r,u,s) \
   (r).c0.re= (u).c00.re*(s).c0.re-(u).c00.im*(s).c0.im  \
             +(u).c01.re*(s).c1.re-(u).c01.im*(s).c1.im  \
             +(u).c02.re*(s).c2.re-(u).c02.im*(s).c2.im; \
   (r).c0.im= (u).c00.re*(s).c0.im+(u).c00.im*(s).c0.re  \
             +(u).c01.re*(s).c1.im+(u).c01.im*(s).c1.re  \
             +(u).c02.re*(s).c2.im+(u).c02.im*(s).c2.re; \
   (r).c1.re= (u).c10.re*(s).c0.re-(u).c10.im*(s).c0.im  \
             +(u).c11.re*(s).c1.re-(u).c11.im*(s).c1.im  \
             +(u).c12.re*(s).c2.re-(u).c12.im*(s).c2.im; \
   (r).c1.im= (u).c10.re*(s).c0.im+(u).c10.im*(s).c0.re  \
             +(u).c11.re*(s).c1.im+(u).c11.im*(s).c1.re  \
             +(u).c12.re*(s).c2.im+(u).c12.im*(s).c2.re; \
   (r).c2.re= (u).c20.re*(s).c0.re-(u).c20.im*(s).c0.im  \
             +(u).c21.re*(s).c1.re-(u).c21.im*(s).c1.im  \
             +(u).c22.re*(s).c2.re-(u).c22.im*(s).c2.im; \
   (r).c2.im= (u).c20.re*(s).c0.im+(u).c20.im*(s).c0.re  \
             +(u).c21.re*(s).c1.im+(u).c21.im*(s).c1.re  \
             +(u).c22.re*(s).c2.im+(u).c22.im*(s).c2.re;

#define _SPU_su3_multiply(r,u,s) \
   v_r = (vector double *)&(r); \
   v_u = (vector double *)&(u);  \
   v_s = (vector double *)&(s); \	         
   v_a = v_s[0];     \
   v_b = v_s[1];     \
   v_c = v_s[2];         \
   vv_a = spu_rlqwbyte(v_a, 8);     \
   vv_b = spu_rlqwbyte(v_b, 8);     \
   vv_c = spu_rlqwbyte(v_c, 8);         \
   v_a = spu_mul(v_a, v_i);     \
   v_b = spu_mul(v_b, v_i);     \
   v_c = spu_mul(v_c, v_i);         \
   v_d   = spu_add(spu_add(spu_mul(v_u[0], v_a),  spu_mul(v_u[1], v_b)),  spu_mul(v_u[2], v_c)); \
   vv_d  = spu_add(spu_add(spu_mul(v_u[0], vv_a), spu_mul(v_u[1], vv_b)), spu_mul(v_u[2], vv_c)); \
   v_e   = (vector double){spu_extract(v_d,0),spu_extract(vv_d,1)}; \
   v_f   = (vector double){spu_extract(v_d,1),spu_extract(vv_d,0)}; \
   v_r[0] = spu_add(v_e, v_f);\
   v_d   = spu_add(spu_add(spu_mul(v_u[3], v_a),  spu_mul(v_u[4], v_b)),  spu_mul(v_u[5], v_c)); \
   vv_d  = spu_add(spu_add(spu_mul(v_u[3], vv_a), spu_mul(v_u[4], vv_b)), spu_mul(v_u[5], vv_c)); \
   v_e   = (vector double){spu_extract(v_d,0),spu_extract(vv_d,1)}; \
   v_f   = (vector double){spu_extract(v_d,1),spu_extract(vv_d,0)}; \
   v_r[1] = spu_add(v_e, v_f);\
   v_d   = spu_add(spu_add(spu_mul(v_u[6], v_a),  spu_mul(v_u[7], v_b)),  spu_mul(v_u[8], v_c)); \
   vv_d  = spu_add(spu_add(spu_mul(v_u[6], vv_a), spu_mul(v_u[7], vv_b)), spu_mul(v_u[8], vv_c)); \
   v_e   = (vector double){spu_extract(v_d,0),spu_extract(vv_d,1)}; \
   v_f   = (vector double){spu_extract(v_d,1),spu_extract(vv_d,0)}; \
   v_r[2] = spu_add(v_e, v_f);

SOLVING THE DIRAC EQUATION ON THE CELL
C. TADONKI (claude.tadonki@u-psud.fr). PetaQCD project, funded by the french agency ANR.