Date: Tue, 19 Jun 2001 15:30:18 -0400 (EDT)
From: Robert Edwards <edwards@jlab.org>
To: Richard Brower <brower@ctp.mit.edu>
Subject: pargrid.c


/* $Id: pargrid.c,v 12.6 2001/02/15 21:44:53 edwards Exp $
 *
 * PARGRID: generic parallel grid routines
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <parsmp_rts.h>

#define MAX_GEOMETRY  2

#undef DEBUG

/* Return the coordinate of the underlying grid */
void
PARSMP_grid_coord(PARSMP_geometry_t geom, int coord[ND])
{
  int phys_coord[ND];
  unsigned m;

  PARSMP_physical_coord(phys_coord);

  for(m=0; m < ND; ++m)
  {
    coord[m] = phys_coord[geom->ordering[m]];
  }
}
  

/* Return the logical size of the machine */
void
PARSMP_machine_size(PARSMP_geometry_t geom, int coord[ND])
{
  unsigned m;

  for(m=0; m < geom->dimension; ++m)
  {
    coord[m] = geom->machine_size[m];
  }
}
  

/* Check if a is a factor of b */
static unsigned is_factor(unsigned a, unsigned b)
{
  return ((b % a) == 0 ? 1 : 0);
}



/* 
 * PARSMP_create_detailed_geometry
 *
 * Detailed geometry constructor 
 */
static struct PARSMP_geometry  PAR_geometries[MAX_GEOMETRY];
static unsigned current_geometry = 0;

PARSMP_geometry_t
PARSMP_create_detailed_geometry(int *length, unsigned rank, unsigned *ordering)
{
  int i;
  int k;
  int phys_coord[ND];
  int phys_size[ND];
  int wrap_around[ND];
  int phys_rank;
  int phys_num_nodes;
  unsigned m;
  int reorder;
  PARSMP_geometry_t geom;

  if (rank < 1 || rank > ND)
  {
    fprintf(stderr,"PARSMP_create_geometry: invalid rank = %d\n", rank);
    exit(1);
  }

  if (current_geometry < MAX_GEOMETRY)
  {
    fprintf(stderr,"current_geometry = %d\n", current_geometry);
    geom = &(PAR_geometries[current_geometry++]);
  }
  else
  {
    fprintf(stderr,"PARSMP_create_geometry: too many allocated geometries = %d\n",
	   current_geometry);
    exit(1);
  }

  /* Geometry should already have been allocated. Fill it in */
  geom->dimension = rank;
  geom->subgrid_vol = 1;

  for(i=0; i < ND; ++i)
  {
    geom->ordering[i] = 1;
    geom->length[i] = 1;
    geom->subgrid_length[i] = 1;
  }

  for(i=0; i < rank; ++i)
  {
    geom->ordering[i] = ordering[i];
  }

  /* Get info about the machine configuration */
  phys_rank      = PARSMP_physical_dimension();
  phys_num_nodes = PARSMP_physical_num_nodes();

  /* If physical rank is zero, this machine allows configuration of size */
  if (phys_rank == 0)
  {
#define NUM_PRIME 4
    unsigned prime_factors[NUM_PRIME] = {2, 3, 5, 7};
    unsigned num_nodes;
    unsigned virt_size[ND];
    unsigned found_prime;
    int k;

    for(i=0; i < rank; ++i)
    {
      virt_size[i] = length[i];
      phys_size[i] = 1;
    }

    /* A simple (primitive) method to determine lattice sizes */
    /* Loop over each dimension repeatedly pulling out prime factors */
    /* No attempt is made to determine if the resulting subgrid is optimal */
    /* in surface area to volume */
    num_nodes = phys_num_nodes;

    do
    {
      /* Reset prime factor count */
      found_prime = 0;

      /* Check and keep track if any prime factor found */
      for(k=0; k < NUM_PRIME; ++k)
      {
	unsigned p = prime_factors[k];

	for(i=rank-1; i >= 0; --i)
	{
	  if (is_factor(p, num_nodes))
	  {
	    if (is_factor(p, virt_size[i]))
	    {
	      phys_size[i] *= p;
	      virt_size[i] /= p;
	      num_nodes    /= p;
	      found_prime = 1;
	    }
	  }
	}
      }
    }
    while (num_nodes > 1 && found_prime);


    /* Sanity checks */
    if (num_nodes != 1)
    {
      fprintf(stderr,"PARSMP_create_geometry: node number not composed of selected prime factors, still %d remaining\n",
	      num_nodes);
      exit(1);
    }

    num_nodes = 1;
    for(i=0; i < rank; ++i)
      num_nodes *= phys_size[i];

    if (num_nodes != phys_num_nodes)
    {
      fprintf(stderr,"PARSMP_create_geometry: some error matching physical problem size to number of nodes\n");
      exit(1);
    }

    /* Create the physical machine geometry */
    PARSMP_create_physical_geometry(phys_size, rank);

    /* Layout should be okay, proceed with subgrid assignment */
    for(i=0; i < rank; ++i)
    {
      unsigned s = phys_size[i];

      geom->length[i] = length[i];

      if (geom->length[i] % s != 0)
      {
	fprintf(stderr,"PARSMP_create_geometry: virtual lattice not divisible by physical lattice\n");
      }
      if (geom->length[i] < s)
      {
	fprintf(stderr,"PARSMP_create_geometry: virtual lattice smaller than physical lattice\n");
      }

      num_nodes *= s;
      geom->machine_size[i]   = s;
      geom->subgrid_length[i] = geom->length[i] / s;
      geom->subgrid_vol *= geom->subgrid_length[i];
    }
  }
  else /* Physical machine size is fixed already */
  {
    if (phys_rank != rank)
    {
      fprintf(stderr,"PARSMP_create_geometry: currently only support equal physical and problem rank\n");
      exit(1);
    }

    /* Get info about the machine configuration */
    PARSMP_physical_size(phys_size);

    for(i=0; i < rank; ++i)
    {
      unsigned a = geom->ordering[i];
      unsigned s = phys_size[a];

      geom->length[i] = length[i];

      if (geom->length[i] % s != 0)
      {
	fprintf(stderr,"PARSMP_create_geometry: virtual lattice not divisible by physical lattice\n");
      }
      if (geom->length[i] < s)
      {
	fprintf(stderr,"PARSMP_create_geometry: virtual lattice smaller than physical lattice\n");
      }

      geom->machine_size[i]   = s;
      geom->subgrid_length[i] = geom->length[i] / s;
      geom->subgrid_vol *= geom->subgrid_length[i];
    }

    for(i=rank; i < ND; ++i)
    {
      unsigned a = geom->ordering[i];
      unsigned s = phys_size[a];

      if (s > 1)
      {
	fprintf(stderr,"PARSMP_create_geometry: unused machine axis is nontrivial: %d\n", a);
      }
    }
  }

  /* Construct virtual node layout and physical neighbours */
  PARSMP_physical_coord(phys_coord);
  geom->nodeid = 0;
  for(k=rank-1;k>=0;k--)
  {
    unsigned a,s,c;
    a = geom->ordering[k];
    s = phys_size[a];
    c = phys_coord[a];
    geom->coord[k] = c;
    geom->nodeid = geom->nodeid * s + c; 
  }


  /* Need some better mechanism here */
  geom->PAR_dim = rank;

  return geom;
}

/* 
 * PARSMP_create_geometry
 *
 * General geometry constructor. Automatically determine sthe 
 * optimal layout -> decides the optimal ordering of axes.
 */
PARSMP_geometry_t
PARSMP_create_geometry(int *length, unsigned rank)
{
  unsigned ordering[ND];
  int phys_rank;

  if (rank < 1 || rank > ND)
  {
    fprintf(stderr,"PARSMP_create_geometry: invalid rank = %d\n", rank);
    exit(1);
  }

  /* Get info about the machine configuration */
  phys_rank = PARSMP_physical_dimension();

  /* If physical rank is zero, this machine allows configuration of size */
  if (phys_rank == 0)
  {
    unsigned m;

    for(m=0; m < ND; ++m)
      ordering[m] = m;
  }
  else /* Physical machine is fixed already */
  {
    unsigned m, nd, k;
    int phys_size[ND];
    unsigned min_size, min_len, used_size[ND], used_len[ND];

    /* Get info about the machine configuration */
    PARSMP_physical_size(phys_size);

    for(m=0; m < ND; ++m)
    {
      ordering[m] = 9999;
      used_len[m] = 0;
      used_size[m] = 0;
    }

    nd = rank;

    while(nd > 0)
    {
      unsigned ind;

      /* Find the current smallest physical machine size */
      for(m=0; m < ND; ++m)
      {
	if (! used_size[m])
	{
	  min_size = m;
	  break;
	}
      }

      for(m=0; m < ND; ++m)
      {
	if ((! used_size[m]) && (phys_size[m] < phys_size[min_size]))
	  min_size = m;
      }

      /* Find the first divisible lattice length */
      for(m=0; m < rank; ++m)
      {
	ind = m;
	if ((! used_len[m]) && (length[m] % phys_size[min_size]) == 0)
	  break;
      }

      if (m == rank)
      {
	fprintf(stderr,"PARSMP_create_geometry: lattice size does not fit in the machine\n");
	fprintf(stderr,"physical machine size\n");
	for(m=0; m < ND; ++m)
	  fprintf(stderr," size[%d] = %d\n",m,phys_size[m]);
	fprintf(stderr,"requested lattice size\n");
	for(m=0; m < rank; ++m)
	  fprintf(stderr," len[%d] = %d\n",m,length[m]);
	exit(1);
      }

      for(m=ind; m < rank; ++m)
      {
	if ((! used_len[m]) && 
	    (length[m] % phys_size[min_size] == 0) && 
	    (length[m] < length[ind]))
	  ind = m;
      }

      ordering[ind] = min_size;
      used_len[ind] = 1;
      used_size[min_size] = 1;
      --nd;
    }

    /* Sanity check */
    /* Reset all remaining unsed  */
    {
      for(m=0; m < ND; ++m)
	printf("phys_size[%d] = %d\n",m,phys_size[m]);

      for(m=0; m < ND; ++m)
	printf("ordering[%d] = %d\n",m,ordering[m]);
      
      for(m=0; m < ND; ++m)
	printf("used_len[%d] = %d\n",m,used_len[m]);

      for(m=0; m < ND; ++m)
	printf("used_size[%d] = %d\n",m,used_size[m]);
    }

    for(m=0; m < rank; ++m)
    {
      if (ordering[m] == 9999)
      {
	fprintf(stderr,"PARSMP_create_geometry: ordering not correct\n");
	exit(1);
      }
    }

    for(m=rank; m < ND; ++m)
    {
      if (ordering[m] != 9999)
      {
	fprintf(stderr,"PARSMP_create_geometry: ordering not correct\n");
	exit(1);
      }

      ordering[m] = m;
    }
  }

  return PARSMP_create_detailed_geometry(length, rank, ordering);
}

/* 
 * PARSMP_create_geometry
 *
 * General geometry constructor. Automatically determine sthe 
 * optimal layout -> decides the optimal ordering of axes.
 */
void 
PARSMP_free_geometry(PARSMP_geometry_t geom)
{
  if (geom != &(PAR_geometries[--current_geometry]))
  {
    fprintf(stderr,"PARSMP_free_geometry: error in removing geometry from stack %d\n",
	   current_geometry);
    exit(1);
  }

  fprintf(stderr,"free_geometry: current_geometry = %d\n", current_geometry);
}


PARSMP_m_site(PARSMP_geometry_t geom,int site,int *s_p,int *s_v)
{
  int i,k;  
  int coord[4];

  /* The lattice coordinates corresponding to this site */
  for(i=0; i<geom->dimension; ++i)
  {
    coord[i] = site % geom->length[i];
    site = site / geom->length[i];
  }

  /* The logical site and logical node */
  k = 0;
  for(i=geom->dimension-1;i>=0;i--)
    k = k*(geom->machine_size[i]) + coord[i]/geom->subgrid_length[i];
  *s_p = k;
  k = 0;
  for(i=geom->dimension-1;i>=0;i--)
    k = k*(geom->subgrid_length[i]) + (coord[i] % geom->subgrid_length[i]);
  *s_v = k;
}


int 
PARSMP_shift(PARSMP_geometry_t geom, int site, unsigned char *data,
	     int prim_size, int sn)
{
  static int save_old_node;
  static int cur_src[4];
  static int cur_lsite[4];
  int nd;
  
  int old_node,cur_node,i,j,pd,shift_size,dir;
  int local_site;
  unsigned int b1,l1,l2;
  int send_dir = 0;
  
  nd = geom->dimension;
  
  if(site == sn)
  {
    save_old_node = 0;
    for(i=0;i<nd;i++)
      cur_lsite[i] = cur_src[i] = 0;
  }  

  PARSMP_m_site(geom,site,&cur_node,&local_site);
  old_node = save_old_node;
  save_old_node = cur_node;

#ifdef DEBUG
  fprintf(stderr," logical site %d node %d local_site %d old_node %d\n",site,cur_node,local_site,old_node);
#endif

  shift_size = prim_size;
  for(dir = 0;dir < nd;dir++)
  {
    shift_size *= geom->subgrid_length[dir];
    b1 =  geom->machine_size[dir];
    i = cur_node % b1;
    j = old_node % b1;
    cur_node = cur_node / b1;
    old_node = old_node / b1;

    if(j != i)
    { 
      /* the coordinates are different, i.e. we have to do a shift */
#ifdef DEBUG
      fprintf(stderr," diff coordinates for dir = %d new %d old %d\n",dir,i,j);
#endif

      if(((i-j+b1) % b1) != 1)
      {
	fprintf(stderr," par_shift: internal error 1\n");
	exit(1);
      }
      if(i != 0)
	cur_lsite[dir] = local_site;
    
#ifdef DEBUG
      fprintf(stderr,"shift in dir %d, offset = %d size = %d\n",
	     dir,cur_src[dir]*geom->subgrid_vol + cur_lsite[dir],shift_size);
#endif

      if (b1 > 1)
      { /* machine_size[dir] > 1 */
	PARSMP_sendrecv_wait(
	  (void *)(data + (cur_src[dir]*geom->subgrid_vol
			   + cur_lsite[dir])*prim_size),
	  (void *)(data + ((1-cur_src[dir])*geom->subgrid_vol
			   + cur_lsite[dir])*prim_size),
	  shift_size, send_dir, geom->ordering[dir]);
      }
      else /* machine_size[dir] = 1 */
      {
	/* Has arguments memcpy(dest, src, count) */
	memcpy((void *)(data + ((1-cur_src[dir])*geom->subgrid_vol
				+ cur_lsite[dir])*prim_size),
	       (void *)(data + (cur_src[dir]*geom->subgrid_vol
				+ cur_lsite[dir])*prim_size),
	       shift_size);
      }

      cur_src[dir] = 1 - cur_src[dir];
      for(j=0;j<dir;j++)
	cur_src[j] = cur_src[dir];
    }
  }

#ifdef DEBUG
  fprintf(stderr,"new offset for read %d\n",local_site + cur_src[0]*geom->subgrid_vol);
#endif

  return local_site + cur_src[0]*geom->subgrid_vol;
}