MPI 带状划分矩阵转置

MPI矩阵操作分两部分：

1.读文件

2.操作

读文件的时候可以使用MPI_File_read来直接分配数据。

矩阵存储的格式：连续内存空间，头两个int表示行列数，后续位置使用double表示数据，行优先存储。

文件使用二进制方式直接存储矩阵内存块。

转置非方阵，使用新申请的内存。如果使用O(1)的内存，并没有规则的交换的过程，所以不大现实，逻辑太复杂，太慢。

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

/*
matrix transpostion using mpi 
linear division
author GT
2013.4.18
*/
int rank,size;
/*
file is in binary form
initial two int are size
afterwards is double data
all matrix is in the same form
*/
/*
read the size of @fileName
return in @sizes 
*/
void readSize(char *fileName, int* sizes)
{
	FILE *in;
	in = fopen(fileName, "rb");
	if(in == NULL)
	{
		return ;
	}
	fread(sizes, sizeof(int), 2, in);
	fclose(in);
}
/*
divide the matrix columns averagely among all threads
while reading the @fileName
the matrix is @row*@col
return in @dest 
*/
void readMatrixByColumn (char *fileName, int row, int col, double *dest)
{
	//amazingly simple function to averagely divide array
	//[colStart,colEnd)
	int colStart = rank * col / size;
	int colEnd = (rank + 1) * col / size; 
	int myCol = colEnd - colStart;
	//columns for a thread 
	MPI_Datatype columns; 
	MPI_Type_vector (row, myCol, col, MPI_DOUBLE, &columns);
	MPI_Type_commit(&columns);
	//file operation
	MPI_File file;
	MPI_Status status;
	MPI_File_open (MPI_COMM_WORLD, fileName, MPI_MODE_CREATE | MPI_MODE_RDWR, MPI_INFO_NULL, &file);
	//current start position
	MPI_Offset offset = colStart * sizeof(double) + 2 * sizeof(int);
	//read
	MPI_File_set_view (file, offset, MPI_DOUBLE, columns, "native", MPI_INFO_NULL);
	MPI_File_seek (file, 0, MPI_SEEK_SET);
	double *data = (double*)((int*)dest + 2);
	MPI_File_read (file, data, myCol * row, MPI_DOUBLE, &status);
	//set coordinate
	*((int*) dest) = row;
	*((int*) dest + 1) = myCol;
	//clean up
	MPI_File_close (&file);
	MPI_Type_free (&columns);
}
	

/*
convert 1d index @one into 2d coordinate @xy 
the matrix is row*@col
*/
extern inline void one2two (int one, int col,int *xy)
{
	xy[0] = one / col;
	xy[1] = one % col;
}
/*
convert using two ints
*/
extern inline int two2one3 (int x, int y, int col)
{
	return x * col + y;
}

/*
convert 2d @xy into 1d @return 
matrix is row * @col
*/
extern inline int two2one2 (int *xy, int col)
{
	return two2one3 (xy[0], xy[1], col);
}
/*
transposition 2d @xy to 1d @return
*/
extern inline int transpositionTwo2one (int *xy, int row)
{
	return two2one3 (xy[1], xy[0], row);
}
/*
transposite @matrix and store in @tran
*/

void transposite (double *matrix, double *tran)
{
	int row = *((int*) matrix);
	int col = *((int*) matrix + 1);
	*((int*) tran) = col;
	*((int*) tran + 1) = row;
	int xy[2];
	double *matrixData = (double*) ((int*) matrix + 2);
	double *tranData = (double*) ((int*) tran + 2);
	for ( int i = 0; i < row * col; ++ i )
	{
		one2two ( i, col, xy );
		*(tranData + transpositionTwo2one (xy, row)) = *(matrixData + i);	
	}
}

void show (double *matrix)
{
	double *data =(double*) ((int*) matrix + 2);
	int row = *((int*) matrix);
	int col = *((int*) matrix + 1);
	printf ("%d %d matrix \n", row, col);
	for (int i = 0; i < row; ++i)
	{
		for (int j = 0; j < col; ++j)
		{
			printf("%lf ",*(data + i * col + j));
		}
		printf("\n");
	}
}
void testRead (char *fileName, int *sizes)
{
	int err;
	MPI_File file;
	MPI_Status status;
	MPI_File_open (MPI_COMM_WORLD, fileName, MPI_MODE_CREATE | MPI_MODE_RDWR, MPI_INFO_NULL, &file);
	MPI_File_set_view (file, 0, MPI_INT, MPI_INT, "native", MPI_INFO_NULL);
	if(0 == rank)
		printf("e %d %d %d \n", err, sizes[0], sizes[1] );
	err = MPI_File_read (file, sizes, 2, MPI_INT, &status);
	if(0 == rank)
		printf("ee %d %d %d \n", err, sizes[0], sizes[1] );
	MPI_File_close (&file);
}
int main (int argc, char* argv[])
{
	MPI_Init (&argc, &argv);
	MPI_Comm_rank (MPI_COMM_WORLD, &rank);
	MPI_Comm_size (MPI_COMM_WORLD, &size);

	int sizes[2];
	readSize (argv[1], sizes);
	double *matrix = (double*) malloc(sizes[0] * sizes[1] * sizeof(double) + 2 *sizeof(int));
	double *transp = (double*) malloc(sizes[0] * sizes[1] * sizeof(double) + 2 *sizeof(int));
	readMatrixByColumn (argv[1], sizes[0], sizes[1], matrix);
	transposite(matrix ,transp);
	if (2 == rank)
	{
		printf ("matrix\n");
		show (matrix);
		printf ("transposition \n");
		show (transp);
	}
	MPI_Finalize ();
	free (matrix);
	return 0;
}