MPI - Sending and Receiving a column

Question

I need to send a matrix column from one process and receive it from another. I tried running the following program, but I'm getting a weird result (at least I think so); Only the first element of the matrix is copied, and some matrix elements change unexpectedly.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "mpi.h"

void swap(int* a,int* b){
    int temp;
    temp=*a;
    *a=*b;
    *b=temp;
}
void print_matrix(double** A,int n){
    int i,j;
    for(i=0;i<n;i++){
        for(j=0;j<n;j++){
            printf("%f ",A[i][j]);
        }
        printf("\n");
    }
}

int main(int argc, char *argv[]){
    int i,j,k,l,n,myid,p,maxp;
    double **A;
    MPI_Datatype col_type;
    MPI_Status status;

    n=3;
    A=malloc(n*sizeof(double*)); /*allocating memory */
    for(i=0;i<n;i++)
        A[i]=malloc(n*sizeof(double));

    A[0][0]=-1;
    A[0][1]=2;
    A[0][2]=-1;
    A[1][0]=2;
    A[1][1]=-1;
    A[1][2]=0;
    A[2][0]=1;
    A[2][1]=7;
    A[2][2]=-3;

    MPI_Init(&argc,&argv);

    MPI_Type_vector(n, 1, n, MPI_DOUBLE,&col_type);
    MPI_Type_commit(&col_type);
    MPI_Comm_size(MPI_COMM_WORLD,&p);
    MPI_Comm_rank(MPI_COMM_WORLD,&myid);

    if(myid==0){
        printf("Starting Method with p=%d\n",p);
        print_matrix(A,n);
    }
    if(myid==0){
            maxp=2;
            A[0][0]=-43;
            A[1][0]=-33;
            A[2][0]=-23;
            printf("BEFORE SENDING\n");
            print_matrix(A,n);
            for(l=0;l<p;l++)
                if(l!=myid){ 
                    MPI_Send(&A[0][0], 1, col_type,l,0,MPI_COMM_WORLD);
                    MPI_Send(&maxp,1,MPI_INT,l,1,MPI_COMM_WORLD);
                }
            printf("AFTER SENDING\n");
            print_matrix(A,n);
    }
    else{
            //receive(k)
            printf("BEFORE RECIEVING\n");
            print_matrix(A,n);
            MPI_Recv(&A[0][1],1,col_type,0,0,MPI_COMM_WORLD,&status);
            MPI_Recv(&maxp,1,MPI_INT,0,1,MPI_COMM_WORLD,&status);
            printf("Just Recieved\n");
            print_matrix(A,n);
    }

    MPI_Finalize();
}

Answer 1

The problem is with your allocation:

A=malloc(n*sizeof(double*)); /*allocating memory */
for(i=0;i<n;i++)
    A[i]=malloc(n*sizeof(double));

This is perfectly fine, but it doesn't necessarily allocate a contiguous array of n*n doubles; it allocates n arrays of n doubles, and those could be scattered all over memory relative to each other. Which (aside from potential cache issues) is also fine, except that when you define a column in this way:

MPI_Type_vector(n, 1, n, MPI_DOUBLE,&col_type);

Eg, n doubles, each of which are n doubles away from the previous one, you're assuming that all of the data is laid out in one big block.

The easiest thing to change is your allocation, to make sure that it is all contiguous and in order (this is almost always what you want for scientific computing):

A=malloc(n*sizeof(double*));        /*allocating pointers */
A[0] = malloc(n*n*sizeof(double));  /* allocating data */
for(i=1;i<n;i++)
    A[i]=&(A[0][i*n]);

/* ... */

free(A[0]);
free(A);