/* builds a communicator from the first n processes
 * in a communicator containing p = m*n processes.
 *
 * here m is implicitly the number of rows in a virtual
 * process grid and n is the number of columns in a virtual
 * process grid.
 *
 * Input: m (the global total number of rows)
 *        Note: with m given n is then derived via n = p/m
 *        im (the global row to create a communicator from)
 *
 * Output: n -- program tests correct creation of new communicator
 *     by broadcasting the value 1 to its members -- all other 
 *     processes have the value 0 -- global sum computed across
 *     ALL the processes.
 *
 * Note: Assumes that MPI_COMM_WORLD contains p = m*n processes
 *
 * See Chap 7, pp. 132 & ff in PPMPI
 * by John L. Weatherwax
 */
#include <stdio.h>
#include "mpi.h"
#include <math.h>
#include <stdlib.h>

main(int argc, char* argv[]) {
    int        p,m,im,n;
    int        my_rank;
    MPI_Group  group_world;
    MPI_Group  im_row_group;
    MPI_Comm   im_row_comm;
    int*       process_ranks;
    int        proc;
    int        test = 0;
    int        sum;
    int        my_rank_in_im_row;

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

    /* Input arguments */ 
    if( my_rank==0 ){
      printf("Enter m: "); scanf("%d",&m);
      printf("Virtual grid is of size [%d,%d]\n",m,(int) (p/m)); 
      printf("Enter im: "); scanf("%d",&im); 
    }
    MPI_Bcast(&m,1,MPI_INT,0,MPI_COMM_WORLD); 
    MPI_Bcast(&im,1,MPI_INT,0,MPI_COMM_WORLD);

    n = (int) (p/m); 

    /* Make a list of the processors to include in the new communicator */
    process_ranks = (int*) malloc(n*sizeof(int));
    for (proc = 0; proc < n; proc++)
      process_ranks[proc] = n*im + proc;

    /* Get the group underlying MPI_COMM_WORLD */
    MPI_Comm_group(MPI_COMM_WORLD, &group_world);

    /* Create the new group */
    MPI_Group_incl(group_world, n, process_ranks, &im_row_group);

    /* Create the new communicator */
    MPI_Comm_create(MPI_COMM_WORLD, im_row_group, &im_row_comm);

    /* Now check whether we can do collective ops in im_row_comm */
    if ((my_rank >= im*n) && (my_rank < (im+1)*n)) {
        MPI_Comm_rank(im_row_comm, &my_rank_in_im_row);
	/* The first processor in im_row_comm should broadcast */ 
        if (my_rank_in_im_row == 0) test = 1;
        MPI_Bcast(&test, 1, MPI_INT, 0, im_row_comm);
    }
    MPI_Reduce(&test, &sum, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);

    if (my_rank == 0) {
      printf("sum = %d\n", sum);
    }

    MPI_Finalize();
}  /* main */