#include <stdlib.h>
#include <stdio.h>
#include <mpi.h>
#include <mpix.h>

// BLACS
#ifdef GPAW_NO_UNDERSCORE_CBLACS
#define Cblacs_gridexit_   Cblacs_gridexit
#define Cblacs_gridinit_   Cblacs_gridinit
#define Csys2blacs_handle_ Csys2blacs_handle
#define Cblacs2sys_handle_ Cblacs2sys_handle
#endif

void Cblacs_gridexit_(int ConTxt);

void Cblacs_gridinit_(int* ConTxt, char* order, int nprow, int npcol);

int Csys2blacs_handle_(MPI_Comm comm);

MPI_Comm Cblacs2sys_handle_(int ConTxt);

// End of BLACS  

// Daniel Faraj
int isRect(int comm_world_rank, MPI_Comm subcomm)
{
  unsigned int x_size, y_size, z_size, t_size;
  unsigned int coords[8];
  unsigned int min_max_coords[8];
  int subcomm_size;
  MPI_Comm_size(subcomm, &subcomm_size);

  MPIX_rank2torus(comm_world_rank, &coords[0], &coords[1], &coords[2], &coords[3]);

  coords[4] = ~(coords[0]);
  coords[5] = ~(coords[1]);
  coords[6] = ~(coords[2]);
  coords[7] = ~(coords[3]);

  MPIR_Allreduce(coords, min_max_coords, 8, MPI_UNSIGNED, MPI_MAX, subcomm);

  /* find if the communicator is a rectangle */
  t_size = (unsigned) (min_max_coords[3] - ~min_max_coords[7] + 1);
  z_size = (unsigned) (min_max_coords[2] - ~min_max_coords[6] + 1);
  y_size = (unsigned) (min_max_coords[1] - ~min_max_coords[5] + 1);
  x_size = (unsigned) (min_max_coords[0] - ~min_max_coords[4] + 1);

  if (x_size * y_size * z_size * t_size == subcomm_size)
    return 1;
  else
    return 0;
}

int main(int argc, char *argv[])
{
  int nprow = 2; 
  int npcol = 2;

  int iam, nprocs;
  int ConTxt;
  int result = 0;

  char order = 'R';

  // Create blacs grid on this communicator                                        
  MPI_Comm bg_comm;

  MPI_Init(&argc, &argv);
  MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
  MPI_Comm_rank(MPI_COMM_WORLD, &iam);

  if (argc > 1) {
    nprow = strtod(argv[1],NULL);
    npcol = strtod(argv[2],NULL);
  }
    
  if (iam == 0) {
    printf("nprow %d", nprow);
    printf("npcol %d", npcol);
  }

  // Sanity check on MPI_COMM_WORLD, then reset result to zero.
  MPIX_Get_property(MPI_COMM_WORLD, MPIDO_RECT_COMM, &result);
  if(result) fprintf(stderr, "MPI_COMM_WORLD is a rect. comm\n");
  result = 0;

  // Create ConTxt out out MPI_COMM_WORLD, then create the BlacsGrid
  ConTxt = Csys2blacs_handle_(MPI_COMM_WORLD);
  Cblacs_gridinit_(&ConTxt, &order, nprow, npcol);

  // This returns the communicator for the BlacsGrid
  if (ConTxt != -1) {
    bg_comm = Cblacs2sys_handle_(ConTxt);
    MPIX_Get_property(bg_comm, MPIDO_RECT_COMM, &result);
    if(result) fprintf(stderr, "Get_property: bg_comm is a rect. subcomm\n");

    result = 0;
    
    result = isRect(iam, bg_comm);
    if(result) fprintf(stderr, "isRect: bg_comm is a rect. subcomm\n");
  }

  MPI_Barrier(MPI_COMM_WORLD);
  MPI_Finalize();
}