/** Within CART_COMM, processes find about their new rank numbers, their cartesian coordinates,
  and their neighbors  */
inline void VCtopology3D::setup_vctopology(MPI_Comm old_comm) {
  // create a matrix with ranks, and neighbours for fields
  MPI_Cart_create(old_comm, 3, divisions, periods, reorder, &CART_COMM);
  // create a matrix with ranks, and neighbours for Particles
  MPI_Cart_create(old_comm, 3, divisions, periods_P, reorder, &CART_COMM_P);
  // field Communicator
  if (CART_COMM != MPI_COMM_NULL) {
    MPI_Comm_rank(CART_COMM, &cartesian_rank);
    MPI_Comm_size(CART_COMM, &nproc);
    MPI_Cart_coords(CART_COMM, cartesian_rank, 3, coordinates);

    MPI_Cart_shift(CART_COMM, XDIR, RIGHT, &xleft_neighbor, &xright_neighbor);
    MPI_Cart_shift(CART_COMM, YDIR, RIGHT, &yleft_neighbor, &yright_neighbor);
    MPI_Cart_shift(CART_COMM, ZDIR, RIGHT, &zleft_neighbor, &zright_neighbor);
  }
  else {
    // EXCEPTION
    cout << "A process is trown away from the new topology for fields. VCtopology3D.h" << endl;
  }
  // Particles Communicator
  if (CART_COMM_P != MPI_COMM_NULL) {
    MPI_Comm_rank(CART_COMM_P, &cartesian_rank);
    MPI_Comm_size(CART_COMM, &nproc);
    MPI_Cart_coords(CART_COMM_P, cartesian_rank, 3, coordinates);

    MPI_Cart_shift(CART_COMM_P, XDIR, RIGHT, &xleft_neighbor_P, &xright_neighbor_P);
    MPI_Cart_shift(CART_COMM_P, YDIR, RIGHT, &yleft_neighbor_P, &yright_neighbor_P);
    MPI_Cart_shift(CART_COMM_P, ZDIR, RIGHT, &zleft_neighbor_P, &zright_neighbor_P);
  }
  else {
    // EXCEPTION
    cout << "A process is trown away from the new topology for Particles. VCtopology3D.h" << endl;
  }

}

/**
 * Creates a Cartesian communicator of size c_dims[0]xc_dims[1] from its input.
 * If c_dims[0]xc_dims[1] would not match the size of in_comm, then the function prints
 * an error and automatically sets c_dims to the correct values.
 *
 * @param in_comm Input MPI communicator handle
 * @param c_dims A 2D integer array, which sets the size of the Cartesian array to c_dims[0]xc_dims[1]
 * @param c_comm A pointer to the Cartesian communicator which will be created
 */
void accfft_create_comm(MPI_Comm in_comm,int * c_dims,MPI_Comm *c_comm){

  int nprocs, procid;
  MPI_Comm_rank(in_comm, &procid);
  MPI_Comm_size(in_comm, &nprocs);

  if(c_dims[0]*c_dims[1]!=nprocs){
    PCOUT<<"ERROR c_dims!=nprocs --> "<<c_dims[0]<<"*"<<c_dims[1]<<" !="<<nprocs<<std::endl;
    c_dims[0]=0;c_dims[1]=0;
    MPI_Dims_create(nprocs,2, c_dims);
    //std::swap(c_dims[0],c_dims[1]);
    PCOUT<<"Switching to c_dims_0="<< c_dims[0]<<" , c_dims_1="<<c_dims[1]<<std::endl;
  }

  /* Create Cartesian Communicator */
  int period[2], reorder;
  int coord[2];
  period[0]=0; period[1]=0;
  reorder=1;

  MPI_Cart_create(in_comm, 2, c_dims, period, reorder, c_comm);
  //PCOUT<<"dim[0]= "<<c_dims[0]<<" dim[1]= "<<c_dims[1]<<std::endl;

  //MPI_Cart_coords(c_comm, procid, 2, coord);

  return;

}

void create_grid(int myrank, int gd,
	MPI_Comm* comm_grid, MPI_Comm* comm_row, MPI_Comm* comm_col)
{
	int dims[2] = {gd, gd};
	int coords[2]; // coords[0] = i, coords[1] = j
	int periods[2];
	int reorder;

	int grid_rank;
	int subdivision[2];

	periods[0] = 0 ; 
	periods[1] = 1 ;
	reorder = 1 ;
	MPI_Cart_create(MPI_COMM_WORLD, 2, dims, periods, reorder, comm_grid);

	MPI_Cart_coords(*comm_grid, myrank, 2, coords); //Outputs the i,j coordinates of the process
	MPI_Cart_rank(*comm_grid, coords, &grid_rank);  //Outputs the rank of the process

	subdivision[0] = 1;
	subdivision[1] = 0;
 	MPI_Cart_sub (*comm_grid,subdivision,comm_col); // Communicator between lines
 	subdivision[0] = 0;
 	subdivision[1] = 1; 
 	MPI_Cart_sub (*comm_grid,subdivision,comm_row); // Communicator between row
}

 /* A two-dimensional torus of 12 processes in a 4x3 grid */
int main(int argc, char *argv[])
{
    int rank, size;
    MPI_Comm comm;
    int dim[2], period[2], reorder;
    int coord[2], id;

    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &size);
    if (size != 12)
    {
        printf("Please run with 12 processes.\n");fflush(stdout);
        MPI_Abort(MPI_COMM_WORLD, 1);
    }

    dim[0]=4; dim[1]=3;
    period[0]=0; period[1]=1;
    reorder=1;
    MPI_Cart_create(MPI_COMM_WORLD, 2, dim, period, reorder, &comm);
    if (rank == 5)
    {
        MPI_Cart_coords(comm, rank, 2, coord);
        printf("Rank %d coordinates are %d %d\n", rank, coord[0], coord[1]);fflush(stdout);
    }
    if(rank==0)
    {
        coord[0]=3; coord[1]=1;
        MPI_Cart_rank(comm, coord, &id);
        printf("The processor at position (%d, %d) has rank %d\n", coord[0], coord[1], id);fflush(stdout);
    }
    MPI_Finalize();
    return 0;
}

void ompi_cart_create_f(MPI_Fint *old_comm, MPI_Fint *ndims, MPI_Fint *dims,
                       ompi_fortran_logical_t *periods, ompi_fortran_logical_t *reorder,
                       MPI_Fint *comm_cart, MPI_Fint *ierr)
{
    MPI_Comm c_comm1, c_comm2;
    int size, c_ierr;
    OMPI_ARRAY_NAME_DECL(dims);
    OMPI_LOGICAL_ARRAY_NAME_DECL(periods);

    c_comm1 = MPI_Comm_f2c(*old_comm);

    size = OMPI_FINT_2_INT(*ndims);
    OMPI_ARRAY_FINT_2_INT(dims, size);
    OMPI_ARRAY_LOGICAL_2_INT(periods, size);

    c_ierr = MPI_Cart_create(c_comm1, OMPI_FINT_2_INT(*ndims),
                             OMPI_ARRAY_NAME_CONVERT(dims),
                             OMPI_LOGICAL_ARRAY_NAME_CONVERT(periods),
                             OMPI_LOGICAL_2_INT(*reorder),
                             &c_comm2);
    if (NULL != ierr) *ierr = OMPI_INT_2_FINT(c_ierr);

    if (MPI_SUCCESS == c_ierr) {
        *comm_cart = MPI_Comm_c2f(c_comm2);
    }

    /*
     * Need to convert back into Fortran, to not surprise the user
     */
    OMPI_ARRAY_FINT_2_INT_CLEANUP(dims);
    OMPI_ARRAY_INT_2_LOGICAL(periods, size);
}

int compute_communicator(int nb_proc_tot, int* nb_proc_row, MPI_Comm* new_comm, int* rank){
  if (nb_proc_row == NULL || new_comm == NULL)
    return EXIT_FAILURE;

  // classical algo to compute quickly integer square root
  // this square root will be the amout of process in each dimension of the matrix
  register unsigned int op, res, one;  
  op = (unsigned int) nb_proc_tot;  
  res = 0;  
  /* "one" starts at the highest power of four <= than the argument. */  
  one = 1 << 30;  /* second-to-top bit set */  
  while (one > op) one >>= 2;  
  while (one != 0) {
    if (op >= res + one) {  
      op -= res + one;
      res += one << 1;  // <-- faster than 2 * one  
    }  
    res >>= 1;  
    one >>= 2;  
  }  
  *nb_proc_row = (int) res;

  // creation of the communicator for the grid
  const int dims[] =  {*nb_proc_row, *nb_proc_row};
  const int periods[] = {1, 1};

  MPI_Cart_create(MPI_COMM_WORLD, 2, dims, periods, 0, new_comm);
  if (*new_comm != MPI_COMM_NULL)
    MPI_Comm_rank(*new_comm, rank);
  else 
    return EXIT_FAILURE;
  // some process may be unnecessary (if nb_proc_tot isn't a integer power 2)
  // so they are stopped
  return EXIT_SUCCESS;
}

int
main(int argc, char *argv[])
{
	int			my_rank;
	int			comm_sz;
	MPI_Init(&argc, &argv);
	MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
	MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);

	MPI_Comm grid;
	int dim[2] = {4, 3};
	int period[2] = {1, 0};
	int reorder = 1;
	int grid_rank = 0;
	int coord[2];

	MPI_Cart_create(MPI_COMM_WORLD, 2, dim, period, reorder, &grid);
	MPI_Comm_rank(MPI_COMM_WORLD, &grid_rank);
	MPI_Cart_coords(grid, grid_rank, 2, coord);

	fprintf(stdout, "my grid rank is %02d, my grid coordinate is (%d, %d) -- my global rank is %d of %d; \n",
			grid_rank, coord[0], coord[1], my_rank, comm_sz);

	MPI_Finalize();
	return 0;
}

/*-----------------------------------------------------------*/
int findNeighbours(){
	int dims[1];
	int periods[1];
	int reorder;

	/* find a good process distribution */
	dims[0] = 0; /* zero so that dims_create tries to rearrange */
	MPI_Dims_create(numMPIprocs, 1, dims);

	/* set periods equal to TURE for periodic boundary conditions ... */
	periods[0] = TRUE;
	/* ...and reorder = FALSE */
	reorder = FALSE;

	/* Create the cartesian topology */
	MPI_Cart_create(comm, 1, dims, periods, reorder, &commCart);

	/* Find the ranks of the left and right neighbour */
	MPI_Cart_shift(commCart, 0, 1, &leftNeighbour, &rightNeighbour);

    MPI_Isend(&myMPIRank, 1, MPI_INT, leftNeighbour, TAG, commCart, &requestArray[0]);

    MPI_Isend(&myMPIRank, 1, MPI_INT, rightNeighbour, TAG, commCart, &requestArray[1]);

    MPI_Irecv(&myGASPIrightNeighbour, 1, MPI_INT, leftNeighbour, TAG, commCart, &requestArray[2]);

    MPI_Irecv(&myGASPIleftNeighbour, 1, MPI_INT, rightNeighbour, TAG, commCart, &requestArray[3]);

    MPI_Waitall(4, requestArray, statusArray);

	return 0;
}

int main(int argc, char *argv[])
{
    int errs = 0, err;
    int dims[2];
    int periods[2];
    int size, rank;
    MPI_Comm comm;

    MTest_Init(&argc, &argv);

    MPI_Comm_size(MPI_COMM_WORLD, &size);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);

    dims[0] = size;
    dims[1] = size;
    periods[0] = 0;
    periods[1] = 0;

    MPI_Comm_set_errhandler(MPI_COMM_WORLD, MPI_ERRORS_RETURN);
    err = MPI_Cart_create(MPI_COMM_WORLD, 2, dims, periods, 0, &comm);
    if (err == MPI_SUCCESS) {
        errs++;
        printf("Cart_create returned success when dims > size\n");
    } else if (comm != MPI_COMM_NULL) {
        errs++;
        printf("Expected a null comm from cart create\n");
    }

    MTest_Finalize(errs);


    return MTestReturnValue(errs);
}

void grid_setup(struct grid_info *grid) {
    /* obtener datos globales */
    MPI_Comm_size(MPI_COMM_WORLD, &(grid->nr_world_processes));
    MPI_Comm_rank(MPI_COMM_WORLD, &(grid->my_world_rank));

    /* calcular cuantos procesos por lado tendra la grilla */
    grid->ppside = intsqrt(grid->nr_world_processes);

    /* crear comunicador para topologia de grilla */
    int dimensions[2]  = {grid->ppside, grid->ppside};
    int wrap_around[2] = {TRUE, TRUE};
    int reorder = TRUE;
    MPI_Cart_create(MPI_COMM_WORLD, 2, dimensions, wrap_around, reorder, &(grid->comm));
    MPI_Comm_rank(grid->comm, &(grid->my_rank));

    /* obtener coordenadas grillisticas del proceso */
    int coordinates[2];
    MPI_Cart_coords(grid->comm, grid->my_rank, 2, coordinates);
    grid->my_row = coordinates[0];
    grid->my_col = coordinates[1];

    /* obtener comunicadores para la fila y la columna del proceso */
    int free_coords_for_rows[] = {FALSE, TRUE};
    int free_coords_for_cols[] = {TRUE, FALSE};
    MPI_Cart_sub(grid->comm, free_coords_for_rows, &(grid->row_comm));
    MPI_Cart_sub(grid->comm, free_coords_for_cols, &(grid->col_comm));
}

void grid_changed_n_nodes()
{
  int per[3] = { 1, 1, 1 };
  GRID_TRACE(fprintf(stderr,"%d: grid_changed_n_nodes:\n",this_node));

  MPI_Comm_free(&comm_cart);
  
  MPI_Cart_create(MPI_COMM_WORLD, 3, node_grid, per, 0, &comm_cart);

  MPI_Comm_rank(comm_cart, &this_node);

  MPI_Cart_coords(comm_cart, this_node, 3, node_pos);

  calc_node_neighbors(this_node);

#ifdef GRID_DEBUG
  fprintf(stderr,"%d: node_pos=(%d,%d,%d)\n",this_node,node_pos[0],node_pos[1],node_pos[2]);
  fprintf(stderr,"%d: node_neighbors=(%d,%d,%d,%d,%d,%d)\n",this_node,
	  node_neighbors[0],node_neighbors[1],node_neighbors[2],
	  node_neighbors[3],node_neighbors[4],node_neighbors[5]);
  fprintf(stderr,"%d: boundary=(%d,%d,%d,%d,%d,%d)\n",this_node,
	  boundary[0],boundary[1],boundary[2],boundary[3],boundary[4],boundary[5]);
#endif

  grid_changed_box_l();
}

static void init_mpi(void)
{
    MPI_Comm_size(MPI_COMM_WORLD, &nproc); //プロセス数の取得
    int dim = 2;          //number of dimension
    int procs[2] = {0,0}; //[0]: x方向の分割数, [1]:y方向の分割数 がはいる
    int period[2] = {0,0};//境界条件, 0は固定境界
    MPI_Comm grid_comm;
    int reorder = 1;   //re-distribute rank flag

    MPI_Dims_create(nproc, dim, procs); //縦横を何分割にするか自動的に計算
    MPI_Cart_create(MPI_COMM_WORLD, 2, procs, period, reorder, &grid_comm); //領域を自動分割 => procs, grid_commは変更される
    MPI_Cart_shift(grid_comm, 0, 1, &ltRank, &rtRank);
    MPI_Cart_shift(grid_comm, 1, 1, &bmRank, &tpRank);

    //プロセス座標において自分がどの位置に居るのか求める(何行何列に居るか)
    int coordinates[2];
    MPI_Comm_rank(grid_comm, &rank);
    MPI_Cart_coords(grid_comm, rank, 2, coordinates);

    SUB_N_X = N_PX / procs[0];
    SUB_N_Y = N_PY / procs[1];
    SUB_N_PX = SUB_N_X + 2; //のりしろ(となりの領域の値が入る部分)の分2大きい
    SUB_N_PY = SUB_N_Y + 2;
    SUB_N_CELL = SUB_N_PX*SUB_N_PY;
    offsetX = coordinates[0] * SUB_N_X; //ランクのインデックスではなく, セル単位のオフセットなのでSUB_N_Xずれる
    offsetY = coordinates[1] * SUB_N_Y;

    /*これだと, 1個のデータをSUB_N_PY跳び(次のデータまでSUB_N_PY-1個隙間がある),SUB_N_X行ぶん取ってくる事になる */
    MPI_Type_vector(SUB_N_X, 1, SUB_N_PY, MPI_C_DOUBLE_COMPLEX, &X_DIRECTION_DOUBLE_COMPLEX);
    MPI_Type_commit(&X_DIRECTION_DOUBLE_COMPLEX);
}

int main(int argc, char *argv[]) {
    int SIZE = atoi(argv[1]);
    MPI_Status status;
    MPI_Comm comm_3d;
    int rank, p;
    int dims[3], periods[3], coords[3];
    dims[0] = dims[1] = dims[2] = periods[0] = periods[1] = periods[2] = 0;
    int i;
    int *b;
    int *my_b = (int *) malloc(SIZE * sizeof(int));

    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &p);

    // Start the vector that is broadcast.
    if (rank == 0){
        b = (int *) malloc(SIZE * sizeof(int));
        for (i = 0; i < SIZE; i++) my_b[i] = b[i] = i;
    }

    // Create the mesh.
    MPI_Dims_create(p, 3, dims);
    MPI_Cart_create(MPI_COMM_WORLD, 3, dims, periods, 0, &comm_3d);

    // Load the coordinates.
    MPI_Cart_coords(comm_3d, rank, 3, coords);

    // The first column will start the broadcast along the rows.
    double start = MPI_Wtime();
    int dim_0_succ, dim_0_pred, dim_1_succ, dim_1_pred, dim_2_succ, dim_2_pred;
    dim_0_succ = dim_0_pred = dim_1_succ = dim_1_pred = dim_2_succ = dim_2_pred = 0;
    MPI_Cart_shift(comm_3d, 0, 1, &dim_0_pred, &dim_0_succ);
    MPI_Cart_shift(comm_3d, 1, 1, &dim_1_pred, &dim_1_succ);
    MPI_Cart_shift(comm_3d, 2, 1, &dim_2_pred, &dim_2_succ);

    if (coords[0] == 0 && coords[1] == 0 && coords[2]) {
        MPI_Send(b, SIZE, MPI_INT, dim_0_succ, 0, MPI_COMM_WORLD);
        MPI_Send(b, SIZE, MPI_INT, dim_1_succ, 0, MPI_COMM_WORLD);
        MPI_Send(b, SIZE, MPI_INT, dim_2_succ, 0, MPI_COMM_WORLD);
    } else if (coords[1] == 0 && coords[2] == 0){
        MPI_Recv(my_b, SIZE, MPI_INT, dim_0_pred, 0, MPI_COMM_WORLD, &status);
        MPI_Send(my_b, SIZE, MPI_INT, dim_0_succ, 0, MPI_COMM_WORLD);
        MPI_Send(my_b, SIZE, MPI_INT, dim_1_succ, 0, MPI_COMM_WORLD);
        MPI_Send(my_b, SIZE, MPI_INT, dim_2_succ, 0, MPI_COMM_WORLD);
    } else if  (coords[3] == 0){
        MPI_Recv(my_b, SIZE, MPI_INT, dim_1_pred, 0, MPI_COMM_WORLD, &status);
        MPI_Send(my_b, SIZE, MPI_INT, dim_1_succ, 0, MPI_COMM_WORLD);
        MPI_Send(my_b, SIZE, MPI_INT, dim_2_succ, 0, MPI_COMM_WORLD);
    } else {
        MPI_Recv(my_b, SIZE, MPI_INT, dim_2_pred, 0, MPI_COMM_WORLD, &status);
        MPI_Send(my_b, SIZE, MPI_INT, dim_2_succ, 0, MPI_COMM_WORLD);
    }

    double end = MPI_Wtime();
 
    if (rank == 0) printf("%d %f\n", SIZE, end - start);
     
    MPI_Finalize();
}

void PX(split_cart_procmesh_for_3dto2d_remap_q1)(
    const INT *n, MPI_Comm comm_cart_3d,
    MPI_Comm *comm_1d
    )
{
  int p0, p1, q0=0, q1=0;
  int ndims, coords_3d[3];
  int dim_1d, period_1d, reorder=0;
  int color, key;
  MPI_Comm comm;

  if( !PX(is_cart_procmesh)(comm_cart_3d) )
    return;

  MPI_Cartdim_get(comm_cart_3d, &ndims);
  if(ndims != 3)
    return;

  PX(get_mpi_cart_coords)(comm_cart_3d, ndims, coords_3d);
  PX(get_procmesh_dims_2d)(n, comm_cart_3d, &p0, &p1, &q0, &q1);

  /* split into p0*p1*q0 comms of size q1 */
  color = coords_3d[0]*p1*q0 + coords_3d[1]*q0 + coords_3d[2]/q1;
  key = coords_3d[2]%q1;
//   key = coords_3d[2]/q0; /* TODO: delete this line after several tests */
  MPI_Comm_split(comm_cart_3d, color, key, &comm);

  dim_1d = q1; period_1d = 1;
  MPI_Cart_create(comm, ndims=1, &dim_1d, &period_1d, reorder,
      comm_1d);

  MPI_Comm_free(&comm);
}

static MPI_Comm cart_comm_1d(MPI_Comm communicator) {
  /* test whether the communicator is cartesian and correct dimensionality */
  fcs_int status;
  MPI_Topo_test(communicator, &status);
  if (status == MPI_CART) {
    /* Communicator is cartesian, so test dimensionality */
    fcs_int ndims;
    MPI_Cartdim_get(communicator, &ndims);
    if (ndims == 3) {
      /* Correct dimensionality, so get grid and test periodicity */
      fcs_int node_grid[3], periodicity[3], node_pos[3];
      MPI_Cart_get(communicator, 3, node_grid, periodicity, node_pos);
      if (!periodicity[0] && !periodicity[1] && periodicity[2]) {
        return communicator;
      }
    }
  }
  /* Otherwise, we have to set up the cartesian communicator */
  fcs_int comm_size;
  MPI_Comm_size(communicator, &comm_size);
  fcs_int node_grid[3] = {0, 0, 0};
  MPI_Dims_create(comm_size, 3, node_grid);
  fcs_int periodicity[3] = {1, 1, 0};

  MPI_Comm new_communicator;
  MPI_Cart_create(communicator, 3, node_grid, periodicity, 1, &new_communicator);

  return new_communicator;
}

int main( int argc, char *argv[] )
{
    int errs = 0;
    int topo_type, size, dims[1], periods[1];
    MPI_Comm comm;

    MTest_Init( &argc, &argv );

    /* Check that topo test returns the correct type, including 
       MPI_UNDEFINED */

    MPI_Topo_test( MPI_COMM_WORLD, &topo_type );
    if (topo_type != MPI_UNDEFINED) {
	errs++;
	printf( "Topo type of comm world is not UNDEFINED\n" );
    }

    MPI_Comm_size( MPI_COMM_WORLD, &size );
    dims[0] = size;
    periods[0] = 0;
    MPI_Cart_create( MPI_COMM_WORLD, 1, dims, periods, 0, &comm );
    MPI_Topo_test( comm, &topo_type );
    if (topo_type != MPI_CART) {
	errs++;
	printf( "Topo type of cart comm is not CART\n" );
    }

    MPI_Comm_free( &comm );
    /* FIXME: still need graph example */

    MTest_Finalize( errs );
    MPI_Finalize();
    return 0;
  
}

void topo_cartesian(int rank, int *dims, int *coords, int *neigh)
{
  int periods[3];
  
  MPI_Comm commcart;
  
  periods[0]=1;
  periods[1]=1;
  periods[2]=1;

  // creation of cartesian communicator
  MPI_Cart_create(MPI_COMM_WORLD,3,dims,periods,0,&commcart);

  // getting the cartesian position
  MPI_Cart_coords(commcart,rank,3,coords);
  
  // getting the neighbors
  MPI_Cart_shift(commcart,0,1,neigh+0,neigh+1); //X
  MPI_Cart_shift(commcart,1,1,neigh+2,neigh+3); //Y
  MPI_Cart_shift(commcart,2,1,neigh+4,neigh+5); //Z

  printf(" proc #%d has coordinates %d %d %d and neighbors Xm=%d Xp=%d Ym=%d Yp=%d Zm=%d Zp=%d \n",rank,coords[0],coords[1],coords[2],neigh[0],neigh[1],neigh[2],neigh[3],neigh[4],neigh[5]);
  //  printf(" dims = %d %d %d\n",dims[0],dims[1],dims[2]);

}

void mpla_init_instance_block_rows(struct mpla_instance* instance, MPI_Comm comm)
{
        instance->comm = comm;

        // get number of process
        MPI_Comm_size(comm, &(instance->proc_count));

        // find number of current process
        MPI_Comm_rank(comm, &(instance->cur_proc_rank));
        if (instance->cur_proc_rank==0)
                instance->is_parent = true;
        else
                instance->is_parent = false;

        // compute the process grid, enforcing only a parallelization over rows
        int dims[2];
        dims[0]=instance->proc_count;
	dims[1]=1;
        MPI_Dims_create(instance->proc_count, 2, dims);
        instance->proc_rows = dims[0];
        instance->proc_cols = dims[1];

        // create cartesian communicator and retrieve cartesian coordinates
        int periods[2];
        periods[0]=periods[1]=0;
        MPI_Cart_create(comm, 2, dims, periods, 0, &(instance->comm));
        int cur_proc_coord[2];
        MPI_Cart_get(instance->comm, 2, dims, periods, cur_proc_coord);
        instance->cur_proc_row = cur_proc_coord[0];
        instance->cur_proc_col = cur_proc_coord[1];;

        cublasCreate(&(instance->cublas_handle));

}

FC_FUNC( mpi_cart_create , MPI_CART_CREATE )
	 ( int *comm_old, int *ndims, int *dims, int *periods,
           int *reorder, int *comm_cart, int *ierr)
{
  *ierr = MPI_Cart_create( *comm_old, *ndims, dims, periods, *reorder,
  	                   comm_cart);
}

void initialiseMonde(int argc, char** argv)
{
	int remain[2], periods[2], reorder, i, nb_thread;

	MPI_Init (&argc, &argv);      				  /* starts MPI */
	MPI_Comm_rank (MPI_COMM_WORLD, &rank);        /* get current process id */
	MPI_Comm_size (MPI_COMM_WORLD, &size);        /* get number of processes */

	// On vérifie qu'il existe au moins un deuxieme argument 
	if(argc < 3)
	{
		MPI_Finalize();
		if(rank == 0)
			printf("Le nombre d'arguments n'est pas suffisant\nAssurez vous de préciser en premier argument le nombre de thread puis la taille de matrice");
		exit(1);
	}
	nb_thread = atoi(argv[1]);
	omp_set_num_threads(nb_thread);

	taille_matrice = atoi(argv[2]);

	racine = sqrt(size);
	if (racine * racine != size || taille_matrice % racine != 0)
	{
		MPI_Finalize();
		if(rank == 0)
			printf("Le nombre de processus MPI n'est pas cohérent avec la taille de la matrice\n");
		exit(1);
	}

	taille_block = taille_matrice / racine;

	tab_dim[LIGNE]   = racine;
	tab_dim[COLONNE] = racine;
	periods[LIGNE]   = 0;
	periods[COLONNE] = 0;
	reorder = 0;

	MPI_Cart_create(MPI_COMM_WORLD, 2, tab_dim, periods, reorder, &COMM_CART);
	// Récupération du rang dans le communicateur cartésien
	MPI_Comm_rank(COMM_CART, &rankCart);

	// Récupération des coordonnées dans le communicateur cartésien
	MPI_Cart_coords(COMM_CART, rankCart, 2, coords);

	// Création du communicateur en ligne
	remain[LIGNE]   = 1;
	remain[COLONNE] = 0;
	MPI_Cart_sub(COMM_CART, remain, &COMM_ROWS);
	// Récupération du rang dans le communicateur en ligne
	MPI_Comm_rank(COMM_ROWS, &rankRow);

	// Création du communicateur en colonne
	remain[LIGNE]   = 0;
	remain[COLONNE] = 1;
	MPI_Cart_sub(COMM_CART, remain, &COMM_COLS);
	// Récupération du rang dans le communicateur en colonne
	MPI_Comm_rank(COMM_COLS, &rankCol);
}