/*
This program demonstrate the use of MPI_Type_indexed. 
Sending the lower triangular portion of a matrix to all processes. */
/****
Note that even though the blocks are uniformly spaced in memory we couldn't use 
MPI_Type_vector because the block lengths are different for each row.
*****/

/* Compile : mpicc -o type_vector type_indexed.c */
/* Run     : mpirun -np 2 type_indexed      */

#include "mpi.h"
#include <stdio.h>
#define SIZE 4

int main(int argc,char *argv[]) {
        int pid, np, tag=1, i,j;
        double a[SIZE][SIZE] = {1.0, 2.0, 3.0, 4.0,  
                               5.0, 6.0, 7.0, 8.0, 
                               9.0, 10.0, 11.0, 12.0,
                               13.0, 14.0, 15.0, 16.0};
        double b[SIZE][SIZE];
        int blocklength[SIZE];
        int displacement[SIZE];
        MPI_Status stat;
        MPI_Datatype index_t;

MPI_Init(&argc,&argv);
MPI_Comm_rank(MPI_COMM_WORLD, &pid);
MPI_Comm_size(MPI_COMM_WORLD, &np);
   
for (i=0;i<SIZE;i++)
   for (j=0;j<SIZE;j++)
       b[i][j]=0;

for (i=0;i<SIZE;i++){
    blocklength[i]=i+1;
    displacement[i]=SIZE*i;
}

MPI_Type_indexed(SIZE, blocklength,displacement, MPI_DOUBLE, &index_t);
MPI_Type_commit(&index_t);

if (np == 2)
   if (pid == 0) {
       MPI_Send(a, 1, index_t, 1, tag, MPI_COMM_WORLD);
       for (i=0;i<SIZE;i++)
           printf("rank= %d:  a[%d,:]= %4.1f %4.1f %4.1f %4.1f\n", 
           pid,i,a[i][0],a[i][1],a[i][2],a[i][3]);
   }
   else {
       MPI_Recv(b, 1, index_t, 0, tag, MPI_COMM_WORLD, &stat); 
       for (i=0;i<SIZE;i++)
           printf("rank= %d:  b[%d,:]= %4.1f %4.1f %4.1f %4.1f\n", 
           pid,i,b[i][0],b[i][1],b[i][2],b[i][3]);
   }
else
  printf("Must specify 2 processors. Terminating.\n");
   
MPI_Finalize();
}