/** Compile and link with:
  /opt/mpich_mx/bin/mpicc -static -o transp transp.c
*/ 

#include <stdio.h>
#include "mpi.h"


int n=4, nb=4, lda=16;

int main(int argc, char *argv[])
{ float A[lda][n],AB[nb][nb],aa; 
  /* The matrix is stored in a column vector of blocks : A'=(A0,A1,A2,A3) */
  int i,j,ib,count=n*n,rank,size,tag;
  MPI_Status status;
  int ip[4]={0,2,1,3};

/** Initialize MPI */
  ...
  ...
  
/** Test the number of processors; print a message and exit if there are not enough */
  ...
  ...
  
/* Processor 0 put the matrix into the array A */ 
  if (rank==0)
  {  for (ib=0; ib<nb; ib++)
       for (j=0; j<n; j++)
         for (i=0; i<nb; i++)  A[i+nb*ib][j]=(ib%2)*nb+j-(i+(ib/2)*nb);
		 
/*  Print the matrix to check it's OK */
      printf("\n ORIGINAL MATRIX \n");
      for (i=0; i<nb; i++)
      {  printf("\n");
         ib = i;
         for (j=0; j<nb; j++) printf("%4.1f ",A[ib][j]);
	 ib=i+nb;
         for (j=0; j<nb; j++) printf("%4.1f ",A[ib][j]);
      }
      for (i=0; i<nb; i++)
      {  printf("\n");
         ib = i+2*nb;
         for (j=0; j<nb; j++) printf("%4.1f ",A[ib][j]);
	 ib=i+3*nb;
         for (j=0; j<nb; j++) printf("%4.1f ",A[ib][j]);
      }
      printf("\n\n");
  } 
  
/* Scatter the blocks to AB on each processor */
  ...
  
/* Transpose the blocks */
  ...
  ...

/* Before gathering the results, be sure every processor is done */
  ...

/** Gather the transpose matrix 
 * Caution: one can not user MPI_Gather because blocks 1 and 2 have to be swapped. 
 * Hence: Process #ib sends bloc #ip[ib] using tag ib 
 * and process #0 receives it into array A */
  ...
  ...
    
/** Process #0 prints the tanspose matrix */  
  if (rank==0)
  {
      printf("\n TRANSPOSED MATRIX \n");
      for (i=0; i<nb; i++)
      {  printf("\n");
         ib = i;
         for (j=0; j<nb; j++) printf("%4.1f ",A[ib][j]);
	 ib=i+nb;
         for (j=0; j<nb; j++) printf("%4.1f ",A[ib][j]);
      }
      for (i=0; i<nb; i++)
      {  printf("\n");
         ib = i+2*nb;
         for (j=0; j<nb; j++) printf("%4.1f ",A[ib][j]);
	 ib=i+3*nb;
         for (j=0; j<nb; j++) printf("%4.1f ",A[ib][j]);
      }
      printf("\n");
  } 
/** Stop MPI */  

  ...
  
  return 0;
}