/*  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;
}