Cellang est un ensemble logiciel destiné à représenter, simuler et analyser une grande variété d'automates cellulaires. Il dispose de son propre langage, Cellang, de son compilateur, cellc, d'une machine virtuelle, avcam, et d'un programme de visualisation, cellview. Le tout est distribué gratuitement avec la licence GPL sur le site de son créateur, Dana Eckart :

            Malgré les faibles performances et le code rudimentaire, Cellang permet de "prototyper" des automates cellulaires de complexité arbitraire et d'en extraire les caractéristiques statistiques avant de commencer le codage de bas niveau. De plus, Cellang est un mélange de Fortran, de C et de Pascal qui n'est pas difficile à comprendre et à apprendre. Il reste tout de même incertain que cellc puisse générer du code très efficace (il semble générer du C) et donc un investissement très important doit être fourni pour améliorer la vitesse de simulation. Générer du code natif n'est qu'un premier pas, il faut aussi traduire efficacement l'abstraction de Cellang en un code qui soit en phase avec le microprocesseur. Comme en VHDL, le compilateur a plus de possibilités pour choisir les algorithmes et les structures de données, mais rien ne dit que le compilateur n'effectue pas une traduction simple, trandformant un forall en une véritable boucle (alors qu'en VHDL une boucle n'est pas traduite en un système itératif). Il reste donc beaucoup d'efforts à faire dans ce domaine mais cela ne pourra être entrepris que lorsque tous les mécanismes sous-jacents auront été étudié.

#  fhp
#  Copyright (C) 1994, 1995, 1997  J Dana Eckart
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 1, or (at your option)
#  any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  You should have received a copy of the GNU General Public License
#  along with CELLULAR; see the file COPYING.  If not, write to the 
#  Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

# An fhp latice gas automata with Cellang.
#
2 dimensions of 
        const which of 0..1
        count of 0..6

        # The presence (1) of a particle traveling in the given direction.
        #
        particle[] for 6 of 0..1
end

# Hexagonal neighborhood.
#
neighbor[] for 6
        := [-1, 1], [0, 1], [1, 0], [0, -1], [-1, -1], [-1, 0] when cell.which
        := [0, 1], [1, 1], [1, 0], [1, -1], [0, -1], [-1, 0] otherwise

# Find out how many particles are incident on this cell.
#
count := 0
first := -1
forall i
        if neighbor[i].particle[i+%3] then
                count := count + 1
                if first < 0 then
                        first := i
                else
                        last  := i
                end
        end
end
cell.count := count

if count = 2 & first+3 = last then
        if random%2 then
                # Rotate Clockwise
                #
                forall i
                        cell.particle[i+%1] := neighbor[i+%3].particle[i]
                end
        else
                # Rotate Counter-Clockwise
                #
                forall i
                        cell.particle[i-%1] := neighbor[i+%3].particle[i]
                end
        end
else
        # The incoming particles pass through unabated.
        #
        forall i
                cell.particle[i] := neighbor[i+%3].particle[i]
        end
end

            En 1996, Oleh Baran (oleh@physics.mcgill.ca) m'a envoyé ce code en Fortran 77. Il ne m'a pas été utile car à l'époque je travaillais sur la technique de consultation de table et j'avais besoin d'une table, pas d'un calcul. De plus, les règles de collision sont FHPI, donc non saturées et sous-efficaces, sans même une particule immobile.

            En revanche, il est intéressant de constater l'organisation du tableau principal, NGAS: c'est un tableau tridimensionnel, il suffit donc de changer l'ordre de déclaration des dimensions pour passer sans effort d'une organisation à une autre. Il ne faut pas oublier cependant que les ordres de déclaration sont différents en C et en Fortran. Enfin, la "consultation de table" dans ce code se résume à reconnaitre 5 cas (collisions frontales et triangulaires) et effectuer une rotation de l'octet si le cas est reconnu. Cette technique astucieuse n'est pas applicable à une échelle différente, par exemple en SIMD et/ou pour d'autres types de collisions.

type de structure : tableau multidimensionnel
collisions : FHPI, calcul (pas de LUT) des collisions frontales et triangulaires
performance : non mesurée. Selon l'auteur, l'efficacité est décevante, puisque les collisions sont réduites au minimum. Il n'a apparemment pas connaissance de méthodes plus efficace.

reproduit avec permission.

c********** one timestep in life of LGA *********
c************* by Oleh Baran ********************
c*********** Au 8th, 1994 ***********************

		SUBROUTINE timestep(CMP,NGAS)
		IMPLICIT REAL*4 (A-H,O-Z)
		COMMON /const/ LL1,LL2,DI
		complex CMP(0:LL1,0:LL2,1:2)
		integer NGAS(0:LL1+2,0:LL2+2,1:2)
		iidump=1

c-----------periodic boundary conditions---------
            do 38 ivr = 1,2
            do 35 iv2 = 1,LL2+1
            NGAS(0,iv2,ivr) = NGAS(LL1+1,iv2,ivr)
            NGAS(LL1+2,iv2,ivr) = NGAS(1,iv2,ivr)
35          continue
            do 37 iv1 = 0,LL1+2
            NGAS(iv1,0,ivr) = NGAS(iv1,LL2+1,ivr)
            NGAS(iv1,LL2+2,ivr) = NGAS(iv1,1,ivr)
37          continue
38          continue
c-----------end of boundary conditions-----------

c-------PROPAGATION------------------------------
		do 41 iv1 = 1,LL1+1
		do 43 iv2 = 1,LL2+1

		kr = NGAS(iv1,iv2+1,1)
		kk0 = and(1,kr)
		kr = NGAS(iv1-1,iv2,2)
		kk1 = and(2,kr)
		kr = NGAS(iv1-1,iv2-1,2)
		kk2 = and(4,kr)
		kr = NGAS(iv1,iv2-1,1)
		kk3 = and(8,kr) Au 8th, 1994
		kr = NGAS(iv1,iv2-1,2)
		kk4 = and(16,kr)
		kr = NGAS(iv1,iv2,2)
		kk5 = and(32,kr)

		rkk= real(kk0+kk1+kk2+kk3+kk4+kk5)
		CMP(iv1-1,iv2-1,1) = cmplx(rkk,0.0)
43		continue
41		continue
                do 45 iv1 = 1,LL1+1
				do 47 iv2 = 1,LL2+1
				kk0 = and(1,NGAS(iv1,iv2+1,2))
				kk1 = and(2,NGAS(iv1,iv2+1,1))
				kk2 = and(4,NGAS(iv1,iv2,1))
				kk3 = and(8,NGAS(iv1,iv2-1,2))
				kk4 = and(16,NGAS(iv1+1,iv2,1))
				kk5 = and(32,NGAS(iv1+1,iv2+1,1))
				rkk= real(kk0+kk1+kk2+kk3+kk4+kk5)
				CMP(iv1-1,iv2-1,2) = cmplx(rkk,0.0)
47              continue
45              continue


c-------end of propagation-----------------------

c-------COLLISIONS-------------------------------
		do 77 iv = 1,2
		do 77 iv1 = 1,LL1+1
		do 77 iv2 = 1,LL2+1
		kk = int(real(CMP(iv1-1,iv2-1,iv)))
c - - - - - table of collisions - - - - - - - - -

			if (kk.EQ.18.OR.kk.EQ.36.OR.kk.EQ.9) THEN
			kk = rotr(kk)
			r3=ran3(iidump)
			if (r3.LT.0.5) THEN
			kk = rotr(kk)
			end if

			elseif (kk.EQ.42.OR.kk.EQ.21) THEN
			kk = rotr(kk)
			end if

c - - - - - end of table of collisions  - - - - -
		NGAS(iv1,iv2,iv) = kk
		bkk = real(bsum(kk)) - 6.*DI
c		bkk = real(bsum(kk))*(-1.)**(iv1+iv2) - 6.*DI
		CMP(iv1-1,iv2-1,iv) = cmplx(bkk,0.0)
77		continue
c-------end of collisions------------------------

		end

            LGAPACK est presque une implémentation de référence dans le domaine, il est écrit par Daniel H. Rothman (MIT) et Stéphane Zaleski (Paris 6).

            Ce code est curieux pour plusieurs raisons, surtout par la manière dont il traite la structure des données, par exemple en utilisant une syntaxe du langage C qui permet une utilisation indépendante de la représentation en mémoire. Il reste à vérifier l'impact de ce choix sur les performances (mesurer avec une configuration et une autre). Ce code montre qu'on peut programmer quelque chose de classique tout en étant étrange.

type de structure : selon #define, mais orienté octets, avec lignes paires et impaires.
collisions : FHPI (5 cas), 6 bits
buffer temporaire : 2 plans
performance : Non mesurée.
 

/*  fhplattice.c

    downloaded from
    http://www.media.mit.edu/people/wad/mas864/src/fhplattice.txt

    This application solves for the behavior of a lattice gas,
    using cellular automata in a hexagonal lattice (FHP)

    J. Watlington, 10/25/95

    Usage:  fhplattice [ block | random | test ] [ -s <array_size> ]

    Modified:
*/

#include <stdio.h>
#include <stdlib.h>
#include "display.h" /* YG : non inclus */

#define SCATTER_ALTERNATELY
#define DEFAULT_MAX_ITERATIONS  2000
#define MAX_ARRAY_SIZE          1024
#define DEFAULT_ARRAY_SIZE      512

/*  The following define how the array is seeded  */

#define SEED_BLOCK_SIZE         ((array_size) / 4)
#define SEED_BLOCK_START        ((3 *(array_size)) / 8)
#define SEED_RAND_SEED          0xFF67533
#define SEED_RAND_THRESHOLD     22000
#define SEED_BLOCK              0
#define SEED_BLOCK_AND_RND      1
#define SEED_TEST               2
#define DEFAULT_ARRAY_SEED      SEED_BLOCK_AND_RND

/*  The following global variables are used for command line parameters */

int   array_size = DEFAULT_ARRAY_SIZE;
int   array_seed = DEFAULT_ARRAY_SEED;
int   max_iterations = DEFAULT_MAX_ITERATIONS;

/*  This is the structure that will represent the lattice of nodes  */

#define DEGREE_OF_INTERCONNECT   6
#define NUM_RULES                64   /* 2^DEGREE_OF_INTERCONNECT */
#define MAX_PARTICLES            6    /* used only for display  */

#ifndef USE_BIT_STORAGE
typedef struct {
  char    northeast;
  char    northwest;
  char    east;
  char    west;
  char    southeast;
  char    southwest;
} lattice_node;
#else
typedef struct {
  char    northeast:1,
          northwest:1,
          east:1,
          west:1,
          southeast:1,
	  southwest:1,
          reserved:2;
} lattice_node;
#endif

/*  One of the few peculiarities of this program is the manner in
    which I represent the lattice arrays.  For ease of access, I
    store an array of pointers (d) to the rows of the array.   */

typedef struct {
  int           size;     /*  Assume square array !  */
  lattice_node  *d[ MAX_ARRAY_SIZE ];
} lattice;

/*  These are the rules the determine the behavior of this lattice gas */

lattice_node  rules[ NUM_RULES ] = {  /* { NE, NW, E, W, SE, SW }   Col. Scatter */
  { 0,0,0,0,0,0 }, { 1,0,0,0,0,0 }, { 0,1,0,0,0,0 }, { 1,1,0,0,0,0 },
  { 0,0,1,0,0,0 }, { 1,0,1,0,0,0 }, { 0,1,1,0,0,0 }, { 1,1,1,0,0,0 },
  { 0,0,0,1,0,0 }, { 1,0,0,1,0,0 }, { 0,1,0,1,0,0 }, { 1,1,0,1,0,0 },
  { 0,1,0,0,1,0 }, { 1,0,1,1,0,0 }, { 0,1,1,1,0,0 }, { 1,1,1,1,0,0 },  /* 0 2b */

  { 0,0,0,0,1,0 }, { 1,0,0,0,1,0 }, { 0,0,1,1,0,0 }, { 1,1,0,0,1,0 },  /* 2 2b */
  { 0,0,1,0,1,0 }, { 1,0,1,0,1,0 }, { 0,1,1,0,1,0 }, { 1,1,1,0,1,0 },
  { 0,0,0,1,1,0 }, { 0,1,1,0,0,1 }, { 0,1,0,1,1,0 }, { 1,1,0,1,1,0 },  /* 1 3b */
  { 0,0,1,1,1,0 }, { 1,0,1,1,1,0 }, { 0,1,1,1,1,0 }, { 1,1,1,1,1,0 },

  { 0,0,0,0,0,1 }, { 0,0,1,1,0,0 }, { 0,1,0,0,0,1 }, { 1,1,0,0,0,1 },  /* 1 2b */
  { 0,0,1,0,0,1 }, { 1,0,1,0,0,1 }, { 1,0,0,1,1,0 }, { 1,1,1,0,0,1 },  /* 2 3b */
  { 0,0,0,1,0,1 }, { 1,0,0,1,0,1 }, { 0,1,0,1,0,1 }, { 1,1,0,1,0,1 },
  { 0,0,1,1,0,1 }, { 1,0,1,1,0,1 }, { 0,1,1,1,0,1 }, { 1,1,1,1,0,1 },

  { 0,0,0,0,1,1 }, { 1,0,0,0,1,1 }, { 0,1,0,0,1,1 }, { 1,1,0,0,1,1 },
  { 0,0,1,0,1,1 }, { 1,0,1,0,1,1 }, { 0,1,1,0,1,1 }, { 1,1,1,0,1,1 },
  { 0,0,0,1,1,1 }, { 1,0,0,1,1,1 }, { 0,1,0,1,1,1 }, { 1,1,0,1,1,1 },
  { 0,0,1,1,1,1 }, { 1,0,1,1,1,1 }, { 0,1,1,1,1,1 }, { 1,1,1,1,1,1 },
};

#ifdef SCATTER_ALTERNATELY
static int           scatter_dir = 0;
static lattice_node  scatter[ 4 ] = {
  { 0,0,1,1,0,0 }, { 0,1,0,0,1,0 }, { 1,0,0,0,0,1 }, { 0,0,1,1,0,0 }
};
#endif

/*  These are prototypes of functions internal to this code module */

void parse_arguments( int argc, char **argv );
lattice *alloc_lattice( int size );
void seed_lattice( lattice *l, int seed );
void calc_collisions( lattice *l, lattice *temp );
void calc_collision_row( lattice *in, lattice *out,
   int prev_row_index, int cur_row_index, int next_row_index );
void display_lattice( window_tag wt, lattice *l );

/*  This is the top level of the program  */

void main( int argc, char **argv )
{
  lattice    *nodes[ 2 ];
  int        next_buf = 1;
  int        cur_buf = 0;
  int        num_iters = 0;
  window_tag w;

  parse_arguments( argc, argv );
  nodes[ cur_buf ] = alloc_lattice( array_size );
  nodes[ next_buf ] = alloc_lattice( array_size );
  seed_lattice( nodes[ cur_buf ], array_seed );
  
  if( (w = display_init_window( "FHP", array_size, array_size,
       MAX_PARTICLES + 1)) < 0 )
  {
    printf( "Unable to open display window.\n" );
    exit( -1 );
  }
  display_lattice( w, nodes[ cur_buf ] );

  while( num_iters++ < max_iterations ) {
    calc_collisions( nodes[ cur_buf ], nodes[ next_buf ] );
    display_lattice( w, nodes[ next_buf ] );

    cur_buf = next_buf;
    next_buf = 1 - next_buf;
  }
}

void calc_collisions( lattice *in, lattice *out )
{
  int           size = in->size - 1;
  int           row;

  calc_collision_row( in, out, size, 0, 1 );

  for( row = 1; row < size; row++ )
    calc_collision_row( in, out, row - 1, row, row + 1 );

  calc_collision_row( in, out, size - 1, size, 0 );
}

void calc_collision_row( lattice *in, lattice *out,
   int prev_row_index, int cur_row_index, int next_row_index )
{
  int      adr;
  int      col;
  int      size = in->size - 1;

  lattice_node   *outcome;
  lattice_node   *cur_row_in = in->d[ cur_row_index ];
  lattice_node   *prev_row_out = out->d[ prev_row_index ];
  lattice_node   *cur_row_out = out->d[ cur_row_index ];
  lattice_node   *next_row_out = out->d[ next_row_index ];

  /*  Special case the first node in a row  */

  adr = 0;
  if( cur_row_in->northeast != 0 )  adr |= 32;
  if( cur_row_in->northwest != 0 )  adr |= 16;
  if( cur_row_in->east  != 0 )  adr |= 8;
  if( cur_row_in->west  != 0 )  adr |= 4;
  if( cur_row_in->southeast != 0 )  adr |= 2;
  if( cur_row_in->southwest != 0 )  adr |= 1;
  
  outcome = rules + adr;
  
  if( cur_row_index & 1 )
  {  /* odd row */
    prev_row_out[ 0 ].southwest = outcome->northeast;
    prev_row_out[ size ].southeast = outcome->northwest;  /* boundary */
    cur_row_out[ 1 ].west = outcome->east;
    cur_row_out[ size ].east  = outcome->west;   /* periodic boundary ! */
    next_row_out[ 0 ].northwest = outcome->southeast;
    next_row_out[ size ].northeast = outcome->southwest;  /* boundary */
  }
  else
  {  /* even row - offset in slightly */
    prev_row_out[ 1 ].southwest = outcome->northeast;
    prev_row_out[ 0 ].southeast = outcome->northwest;
    cur_row_out[ 1 ].west = outcome->east;
    cur_row_out[ size ].east  = outcome->west;   /* periodic boundary ! */
    next_row_out[ 1 ].northwest = outcome->southeast;
    next_row_out[ 0 ].northeast = outcome->southwest;
  }
  
  /*  Now process the inner nodes  */
  
  for( col = 1; col < size; col++ ) {
    cur_row_in++;
    adr = 0;
    if( cur_row_in->northeast != 0 )  adr |= 32;
    if( cur_row_in->northwest != 0 )  adr |= 16;
    if( cur_row_in->east  != 0 )  adr |= 8;
    if( cur_row_in->west  != 0 )  adr |= 4;
    if( cur_row_in->southeast != 0 )  adr |= 2;
    if( cur_row_in->southwest != 0 )  adr |= 1;
  
    outcome = rules + adr;

#ifdef SCATTER_ALTERNATELY
    switch( adr ) {
      case 0x21:
	outcome = &scatter[ scatter_dir ];
	scatter_dir = 1 - scatter_dir;
	break;

      case 0x0C:
	outcome = &scatter[ scatter_dir + 1 ];
	scatter_dir = 1 - scatter_dir;
	break;

      case 0x12:
	outcome = &scatter[ scatter_dir + 2 ];
	scatter_dir = 1 - scatter_dir;
	break;
      }
#endif

    if( cur_row_index & 1 )
    {  /* odd row */
      prev_row_out[ col ].southwest = outcome->northeast;
      prev_row_out[ col - 1 ].southeast = outcome->northwest;
      cur_row_out[ col + 1 ].west = outcome->east;
      cur_row_out[ col - 1 ].east  = outcome->west;
      next_row_out[ col ].northwest = outcome->southeast;
      next_row_out[ col - 1 ].northeast = outcome->southwest;
    }
    else
    {  /* even row - offset in slightly */
      prev_row_out[ col + 1 ].southwest = outcome->northeast;
      prev_row_out[ col ].southeast = outcome->northwest;
      cur_row_out[ col + 1 ].west = outcome->east;
      cur_row_out[ col - 1 ].east  = outcome->west;
      next_row_out[ col + 1 ].northwest = outcome->southeast;
      next_row_out[ col ].northeast = outcome->southwest;
    }
  }
  
  /* Special case the last node in a row  */
  
  cur_row_in++;
  adr = 0;
  if( cur_row_in->northeast != 0 )  adr |= 32;
  if( cur_row_in->northwest != 0 )  adr |= 16;
  if( cur_row_in->east  != 0 )  adr |= 8;
  if( cur_row_in->west  != 0 )  adr |= 4;
  if( cur_row_in->southeast != 0 )  adr |= 2;
  if( cur_row_in->southwest != 0 )  adr |= 1;
  
  outcome = rules + adr;
   
  if( cur_row_index & 1 )
  {  /* odd row */
    prev_row_out[ size ].southwest = outcome->northeast;
    prev_row_out[ size - 1 ].southeast = outcome->northwest;
    cur_row_out[ 0 ].west = outcome->east;    /* periodic boundary ! */
    cur_row_out[ size - 1 ].east  = outcome->west;
    next_row_out[ size ].northwest = outcome->southeast;
    next_row_out[ size - 1 ].northeast = outcome->southwest;
  }
  else
  {  /* even row - offset in slightly */
    prev_row_out[ 0 ].southwest = outcome->northeast;  /* boundary */
    prev_row_out[ size ].southeast = outcome->northwest;
    cur_row_out[ 0 ].west = outcome->east;   /* periodic boundary ! */
    cur_row_out[ size - 1 ].east  = outcome->west;
    next_row_out[ 0 ].northwest = outcome->southeast;  /* boundary */
    next_row_out[ size ].northeast = outcome->southwest;
  }
}

void parse_arguments( int argc, char **argv )
{
  int   index;
  
  if( argc > 1 )
  {
    index = 1;
	while( index < argc )
	{
	  if( argv[ index ][0] == '-' )
	  {   /*  - option, only -size defined right now  */
	    if( index + 1 >= argc )
		{
		  printf( "%s option requires an argument !\n", argv[ index ] );
		  exit( -1 );
		}
		array_size = atoi( argv[ ++index ] );
		printf( "Simulating a %d x %d array\n", array_size, array_size );
	  }
	  else
	  {   /*  must be array seed - block, rand, or test  */
        switch( argv[ index ][0] ) {
		case 'b':
			array_seed = SEED_BLOCK;
			break;
		case 'r':
			array_seed = SEED_BLOCK_AND_RND;
			break;
		case 't':
		    array_seed = SEED_TEST;
			break;
		default:
		  printf( "I don't understand %s as a seed pattern !\n",
		    argv[ index ] );
		  exit( -1 );
		}
	  }
	  index++;
	}
  }
  srand( SEED_RAND_SEED );
}

void seed_lattice( lattice *l, int seed )
{
  int   row, col;
  int   lattice_size = l->size;
  lattice_node  *n;

  /* Now seed the lattice structure.  It has already
     been cleared to all zeroes.                      */

  if( seed == SEED_BLOCK_AND_RND )
    for( row = 0; row < lattice_size; row++ ) {
      for( n = l->d[ row ], col = 0; col < lattice_size; col++, n++ ) {
	if( rand() > SEED_RAND_THRESHOLD) n->northeast = 1;
	if((rand() > SEED_RAND_THRESHOLD) n->northwest = 1;
	if((rand() > SEED_RAND_THRESHOLD) n->east = 1;
	if((rand() > SEED_RAND_THRESHOLD) n->west = 1;
	if((rand() > SEED_RAND_THRESHOLD) n->southeast = 1;
	if((rand() > SEED_RAND_THRESHOLD) n->southwest = 1;
      }
    }

  if( (seed == SEED_BLOCK) || (seed == SEED_BLOCK_AND_RND) )
    for( row = SEED_BLOCK_START;
         row < (SEED_BLOCK_START + SEED_BLOCK_SIZE); row++ ) {
      n = l->d[ row ] + SEED_BLOCK_START;
      for( col = SEED_BLOCK_START;
	  col < (SEED_BLOCK_START + SEED_BLOCK_SIZE); col++, n++ ) {
	n->northeast = 1;
	n->northwest = 1;
	n->east  = 1;
	n->west  = 1;
	n->southeast = 1;
	n->southwest = 1;
      }
    }
  
  if( seed == SEED_TEST ) {
    for( n = l->d[ 0 ], col = 0; col < l->size; col++, n++ ) {
      n->east  = 1;
      n->southwest = 1;
    }

    for( row = 2; row < (array_size - 1); row++ ) {
      n = l->d[ row ] + 4 + (row / 2);
      n->northeast = 1;
      n->west  = 1;
    }
  }
}

/*  alloc_lattice
    This lattice is stored as an array of lattice_nodes, where each
    row is pointed to independently (although in this case the lines
    are actually stored consecutively.)
*/
lattice *alloc_lattice( int lattice_size )
{
  int           row, col;
  long          total_array_size;
  lattice_node  *n;
  lattice_node  *row_start;
  lattice       *l;

  /* malloc the memory */

  total_array_size = (lattice_size * lattice_size * sizeof(lattice_node))
    + sizeof( lattice );

  if( (l = (lattice *)malloc( total_array_size )) == NULL ) {
    printf( "unable to alloc %d bytes of memory for lattice\n",
    	total_array_size );
    exit( -1 );
  }

  /* Now init the lattice structure.  */

  l->size = lattice_size;
  row_start = (lattice_node *)((long)l + sizeof( lattice ));
  for( row = 0; row < lattice_size; row++ ) {
    l->d[ row ] = row_start;

    for( n = row_start, col = 0; col < lattice_size; col++, n++ ) {
      n->northeast = 0;
      n->northwest = 0;
      n->east = 0;
      n->west = 0;
      n->southeast = 0;
      n->southwest = 0;
    }
    row_start += lattice_size;
  }
  return( l );
}

void display_lattice( window_tag wt, lattice *l )
{
  int           x, y;
  int           size = l->size;
  int           num_particles;
  lattice_node  *n;
	
  for( y = 0; y < size; y++ )
  {
    n = l->d[ y ];
    for( x = 0; x < size; x++, n++ )
    {
      num_particles = n->northeast; 
      num_particles += n->northwest; 
      num_particles += n->east; 
      num_particles += n->west; 
      num_particles += n->southeast; 
      num_particles += n->southwest; 
      display_set_pixel( wt, x, y, num_particles );
    }
  }
  display_update( wt );
}

            En 1987, Fung F. Lee (Stanford University) et George B. Adams (Purdue) ont conduit des recherches (NCC 2-387) sur la parallélisation des codes style FHP dans un hypercube à 32 processeurs au NASA Ames Research Center. Ils observent dans leur rapport que le speedup est quasi linéaire (x28 / 32 CPU). Ils expliquent aussi leur amélioration de la technique de passage de messages (piggy-back) : comme un seul bit doit être envoyé vers le nord-est (ou sud-est) le message est envoyé en deux hops, grâce à un ordonnancement des envois et un tag sur les messages.

            Le modèle à 6 bits est traité par calcul booléen. Ici, les deux phases, translation et collision, sont indépendents pour des raisons de communication intra-CPU. Le mélange entre les deux phase est plus difficile qu'avec un ordinateur à mémoire partagée et mérite qu'on s'y attarde, mais il dépend étroitement de l'architecture cible. Le problème de la complexité du calcul booléen est aussi abordé mais la méthode de résolution n'est pas exposée. Il semble qu'ils aient employé une méthode "directe", descriptive, vu la simplicité des règles de collision (comme en FHP1 mais plus complexe).

type de structure : distribuée, 7 plans de bits par noeud, lignes paires/impaires (ce qui pose un problème de topologie dans un hypercube...)
collisions : FHPII 6 bits (note : mêmes règles que celles utilisées par S. Wolfram sur sa CM1)
buffer temporaire :
performance : 50Kc/s (donc 1,4Mc/s avec 32 noeuds)

#define WSIZE	32
#define MSB	31
typedef unsigned long Word;

/* There are YN by XN*WSIZE bits in each of the state-bit planes */
/* YN must be even */
Word s0[YN][XN],s1[YN][XN],s2[YN][XN],s3[YN][XN],s4[YN][XN],s5[YN][XN];
Word obstacle[YN][XN];

/* This is the main computation loop which updates the cells. */
void fluid_flow()
{
	for (t = 0; t < TMAX; t++) {
		advance();
		collide();
	}
}

/* YG : translation : */

#define AllY  for (y=0; y<YN; y++)
#define EvenY for (y=0; y<YN; y+=2)
#define OddY  for (y=1; y<YN; y+=2)
#define AllX  for (x=0; x<XN; x++)

#define RowLS(s) for (x=0; x<XN-1; x++) s[y][x] = s[y][x]<<1 | s[y][x+1]>>MSB; \
	s[y][XN-1]<<=1;

#define RowRS(s) for (x=XN-1; x>0; x--) s[y][x] = s[y][x]>>1 | s[y][x-1]<<MSB; \
	s[y][0]>>=1;

#define ColDS(s) for (y=0; y<YN-1; y++) s[y][x] = s[y+1][x];

#define ColUS(s) for (y=YN-1; y>0; y--) s[y][x] = s[y-1][x];

/* advance() moves state bits by working on the bit planes one at a 
|| time. This is a over-simplified version ignoring boundary cells and
|| communication.
*/

void advance()
{
	int x, y;

	AllY  {RowRS(s0)}
	EvenY {RowRS(s1)}	AllX {ColUS(s1)}
	OddY  {RowLS(s2)}	AllX {ColUS(s2)}
	AllY  {RowLS(s3)}
	OddY  {RowLS(s4)}	AllX {ColDS(s4)}
	EvenY {RowRS(s5)}	AllX {ColDS(s5)}
}

/* Algorithm for advance with piggyback communication :
|| This is a simplified version of advance() with pseudo code.
|| Conceptually, the following action happens for each of the six
|| non-obstacle bit planes.
|| The boundary bits on a specific side of a bit plane are gathered 
|| into a message buffer, and sent to a specific neighbor.
|| The message is then received by the neighbor, and the bits of the
|| message are scattered to a specific boundary of the same bit plane
|| of the neighbor.
*/
void advance()
{
	int x, y;
	
	/* R-bits, L-bits, U-bits and D-bits are the bits on the
	|| right,left, up (top) and down (bottom) boundaries respectively.
	*/

	{Gather R-bits from s0,s1,s5 into a message, and Send to RIGHT.}
	AllY  {RowRS(s0)}
	EvenY {RowRS(s1)}
	EvenY {RowRS(s5)}
	{Receive a message from LEFT, and scatter bits to s0,s1,s5.}

	{Gather L-bits from s3,s4,s2 into a message, and Send to LEFT.}
	AllY  {RowLS(s3)}
	OddY  {RowLS(s4)}
	OddY  {RowLS(s2)}
	{Receive a message from RIGHT, and scatter bits to s3,s4,s2.}

	{Gather U-bits from s1,s2 into a message, and Send to UP.}
	AllX  {ColUS(s1)}
	AllX  {ColUS(s2)}
	{Receive a message from DOWN, and scatter bits to s1,s2.}

	{Gather D-bits from s4,s5 into a message, and Send to DOWN.}
	AllX  {ColDS(s4)}
	AllX  {ColDS(s5)}
	{Receive a message from UP, and scatter bits to s4,s5.}
}

/* collide() effectively transforms the old states into new states
|| for WSIZE (32) cells at a time in its inner most loop by evaluating
|| the boolean expressions which encode the set of interaction rules.
*/
void collide()
{
	int x, y;
	Word obs, a, b, c, d, e, f;
	Word Ca, Cb, Cc, abcdef;

	for (y = 0; y < YN; y++) {
		for (x = 0; x < XN; x++) {
			obs = obstacle[y][x];
			a = s0[y][x]; b = s1[y][x]; c = s2[y][x];
			d = s3[y][x]; e = s4[y][x]; f = s5[y][x];

			abcdef = a & ~b & c & ~d & e & ~f |
			        ~a & b & ~c & d & ~e & f;

			Ca = ~(a|d) & b & e | a & d & ((c&f) & ~(b|e) |
					~(c|f) & ~(b&e)) | abcdef;
			Cb = ~(b|e) & c & f | b & e & ((d&a) & ~(c|f) |
					~(d|a) & ~(c&f)) | abcdef;
			Cc = ~(c|f) & d & a | c & f & ((e&b) & ~(d|a) |
					~(e|b) & ~(d&a)) | abcdef;

			s0[y][x] = ~obs & (Ca ^ a) | obs & d;
			s1[y][x] = ~obs & (Cb ^ b) | obs & e;
			s2[y][x] = ~obs & (Cc ^ c) | obs & f;
			s3[y][x] = ~obs & (Ca ^ d) | obs & a;
			s4[y][x] = ~obs & (Cb ^ e) | obs & b;
			s5[y][x] = ~obs & (Cc ^ f) | obs & c;
		}
	}
}

            Les codes pour calculer FHP sont très variés et il serait difficile d'en faire une analyse comparée poussée. Chacun peut s'inspirer des exemples de cette annexe pour concevoir et explorer l'algorithme qui lui plaît, en gardant à l'esprit que "Le code le plus rapide n'existe pas" (M. Abrash).

            Toutefois, la lecture de ces exemples permet de dégager quelques axes et règles de conception pour les programmes FHP ou similaires :

• Utiliser les règles de codage saturées (au moins) pour que le calcul soit efficace. L'exemple d'Oleh Baran montre que des règles incomplètes font travailler l'ordinateur vite mais pour presque rien, puisque le mean free path est très long. Les collisions saturées sont plus compliquées mais il faut moins de sites et de pas de calcul pour arriver au résultat. Le caractère memory bound du modèle fait pencher la balance encore plus vers FHP3+.

• Veiller à ce que le code reste indépendent d'une géométrie ou d'une taille de mots particulière. Par exemple en 1), le tunnel est défini par plusieurs #define (ce que l'auteur a modifié plus tard), et la plupart des codes utilisent des tailles de mots de 32 bits. Bien que cela soit portable (puisque c'est le souci de nombreuses personnes), il reste encore des plateformes 16 bits et les plateformes 64 et 128 bits se démocratisent. Il faut donc pouvoir définir la taille des mots et des géométries à un seul endroit. Comme cela a été fait dans le programme du mémoire, il est recommandé de rendre la taille du tunnel variable, afin de pouvoir "accélérer" les calculs en redimensionnant le tunnel en cours de route.

• Les techniques de strip-mining et d'affichage progressif n'ont pas été utilisées : il est recommandé d'en tenir compte dès le début de la conception du programme car elles l'influencent profondément. Les programmes présentés ne visent pas des performances exceptionnelles et ont été écrits pour des stations de travail, ce qui explique que le strip-mining n'ait pas été étudié. Cependant, l'affichage des particules, pas à pas, permet de vérifier formellement le comportement du fluide comme des règles de collision. Une interface graphique, même rudimentaire, aurait permis de trouver rapidement l'erreur commise dans LGAPACK. Une erreur minime arrive tellement vite qu'il est parfois trop tard pour réparer l'erreur, et l'interface graphique du projet est un élément essentiel pour le développement et le débogage du programme, en absence de débogueur.

• Quelle méthode de déplacement choisir ? Celle qui effectue le moins d'opérations possible et qui favorise la localité des données. En particulier, il faut réduire le nombre d'accès désordonnés et redondants à la mémoire, en faisant attention à la taille des lignes de caches. Ainsi, sur une station de travail classique, il faut mélanger la phase de calcul et de déplacement car les données se trouvent au même endroit au même moment. Or, dans la plupart des cas rencontrés, les phases sont séparées : tout le tableau est calculé, puis le déplacement est effectué (voir 1) et LGAPACK). Cette approche ralentit le programme considérablement car il y a deux fois plus de cache miss. En concentrant le calcul et le déplacement, on peut en profiter pour traiter et afficher les données.

• Il est inutile, pour un programme qui veut être rapide, de recourir à de nombreux appels de sous-fonctions : dans l'exemple 1) (void nextstate()) et dans LGAPACK, la boucle centrale contient de nombreux appels de routines qui bénéficieraient de l'inlining. Dans LGAPACK, la fonction de calcul d'un site est appelée pour chaque site, inutilement : il serait possible avec un simple copier/coller de gagner quelques pourcents sur la durée d'exécution du programme...

• L'orientation et la projection du réseau hexagonal sur le réseau carré sont des choix qu'il faut mûrir patiemment. Dans le corps du mémoire, nous avons vu que les lignes paires et impaires ont d'abord été orientées horizontalement, puis verticalement pour le modèle à 32 bits et enfin horizontalement pour le travail du mémoire, à cause de propriétés différentes du traitement. Le réseau rhombohédrique de LGAPACK est une alternative intéressante mais son impact sur les algorithmes annexes est encore inconnu.

• La chiralité est un sujet délicat :
  - certains codes choisissent la chiralité au hasard au début du pas de temps (LGAPACK). Cela peut rendre le développement, les vérifications et le débogage difficiles. La simplicité est assurée mais l'impact sur la chiralité générale du modèle n'a pas été étudiée à ma connaissance.
  - D'autres choisissent la chiralité en fonction de la parité du numéro de la ligne (H. Stockman entre autres). Cette technique est déconseillée.
  - Enfin la technique choisie dans le mémoire est d'utiliser le "bruit brownien" du domaine d'étude pour choisir, site par site, la chiralité. Le compromis est intéressant car la chiralité n'est pas corrélée avec la présence ou l'absence des particules. Dans le cas où il faudrait gagner encore quelques cycles lors du calcul, la première technique (LGAPACK) serait utilisée mais la taille du code de détection risque de doubler et de ralentir le développement du programme.