Re: DAcreate2d process layout order

[Sean Dettrick <sean@xxxxxxxxxxxxxxxxxx>]

Barry Smith wrote:

> On Thu, 18 May 2006, Sean Dettrick wrote:
>
> > Hi Barry,
> > the order is determined by MPI_Cart_create.
>
>
>    Do you mean that MPI_Cart_create() orders across the 2nd (y-axis)
> fastest and then the first (x-axis)? Hmmm, maybe we should change the
> DA? Changing it once and for all (not supporting both) is probably
> not a big deal and shouldn't break much (I hope).

Hi Barry,

it depends, what do you call x and what do you call y?
MPI_Cart_coords returns a vector, coords - I tend to say x is coords[0], 
y is coords[1] and z is coords[2]. 

For what it's worth, there's a short code appended to this email, which 
produces:

rank = 0 has Cartesian coords = { 0, 0 }
rank = 1 has Cartesian coords = { 0, 1 }
rank = 2 has Cartesian coords = { 1, 0 }
rank = 3 has Cartesian coords = { 1, 1 }
rank = 0 has DA range x=[0,50) and y=[0,50)
rank = 1 has DA range x=[50,100) and y=[0,50)
rank = 2 has DA range x=[0,50) and y=[50,100)
rank = 3 has DA range x=[50,100) and y=[50,100)

> > >    I don't completely understand what goes wrong. Is it because YOUR
> > > application orders the processors related to geometry in the 
> > > following way?
> > >
> > >     ^ y direction
> > >     |
> > >        2   5  8
> > >        1   4  7
> > >        0   3  6
> > >
> > >                     -> x direction
> > >
> > > Or is this something inherent in MPI_Cart_create?

For my interpretation of x and y, MPI_Cart_create produces the above 
layout.  But if I said x=coords[1] and y=coords[0], then it would match 
the one below.

> > > PETSc does it so
> > >
> > >     ^ y direction
> > >     |
> > >        6   7  8
> > >        3   4  5
> > >        0   1  2
> > >
> > >                     -> x direction
> > >
> > >   

Code and makefile attached ... hopefully within the message size limit.
Just make cartcommtest.

Sean



## makefile for cartcommtest

include ${PETSC_DIR}/bmake/common/base

cartcommtest: cartcommtest.o  chkopts
	-${CLINKER} -o cartcommtest cartcommtest.o  ${PETSC_KSP_LIB}
	${RM} cartcommtest.o

include ${PETSC_DIR}/bmake/common/test

## end makefile

/*  cartcommtest.c  */

static char help[] = "Look at MPI Cartesian communicator process order";


#define CART_COMM
/*ifdef CART_COMM then supply a cartesian communicator.*/


#include "petscda.h"

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **args)
{
  DA             da;             /* distributed array data structure */
  PetscInt       mx = 100,my = 100;
  PetscErrorCode ierr;
  PetscTruth     QUIET=0;
  PetscMPIInt    rank,size;
  int id;

  
#ifdef CART_COMM
  int dims[2] = {0,0};
  int periods[2] = {0,0};
  int coords[2] = {0,0};
  MPI_Comm CartComm2D;

  MPI_Init( &argc, &args );
  MPI_Comm_rank( MPI_COMM_WORLD, &rank );
  MPI_Comm_size( MPI_COMM_WORLD, &size );
  MPI_Dims_create( size, 2, dims);

  MPI_Cart_create( MPI_COMM_WORLD, 2, dims, periods, 0, &CartComm2D );

  if ( rank == 0 ) /* print Cart coords from root:*/
    {
      for (id=0; id<size; id++)
	{
	  MPI_Cart_coords( CartComm2D, id, 2, coords );
	  printf( "rank = %d has Cartesian coords = { %d, %d } \n", id, coords[0], coords[1] );
	}
    }

  PETSC_COMM_WORLD = CartComm2D;
#endif

  PetscInitialize(&argc,&args,(char *)0,help);

  ierr = PetscOptionsGetInt(PETSC_NULL,"-mx",&mx,PETSC_NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetInt(PETSC_NULL,"-my",&my,PETSC_NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetTruth(PETSC_NULL,"-quiet",&QUIET,PETSC_NULL);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Create distributed array (DA) to manage parallel grid and vectors
  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

#ifdef CART_COMM

  /* manually specify domain decomposition */

  PetscInt *lx = malloc( dims[0] * sizeof(PetscReal));
  PetscInt *ly = malloc( dims[1] * sizeof(PetscReal));

  /* make lx[j] = mx/dims[0], and correct for integer rounding to ensure sum(lx[j])=mx  */

  PetscInt i,isum=0;
  for (i=0; i<dims[0]; i++){
    lx[i]=mx/dims[0];
    isum += lx[i];
  }
  lx[dims[0]-1] += mx-isum;
  isum=0;
  for (i=0; i<dims[1]; i++){
    ly[i]=my/dims[1];
    isum += ly[i];
  }
  ly[dims[1]-1] += my-isum;

  ierr = DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,mx,my,dims[0],dims[1],1,1,lx,ly, &da);CHKERRQ(ierr);
#else
  /* automatic domain decomposition */
  ierr = DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,mx,my,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL, &da);CHKERRQ(ierr);
#endif


  PetscInt xs,ys,xm,ym;

  ierr = DAGetInfo( da, 0, &mx, &my, 0,0,0,0,0,0,0,0);CHKERRQ(ierr);
  ierr = DAGetCorners( da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);

  MPI_Comm_rank( MPI_COMM_WORLD, &rank );
  MPI_Comm_size( MPI_COMM_WORLD, &size );
 
  for (id=0;id<size;id++) // print from each CPU in order:
    {
      if ( id == rank ) printf( "rank = %d has DA range x=[%d,%d) and y=[%d,%d) \n", rank, xs, xs+xm, ys, ys+ym );
      MPI_Barrier( MPI_COMM_WORLD );
    }

  ierr = PetscFinalize();CHKERRQ(ierr);
  return 0;
}