MPIonCell/examples/matmult/src/spu/matmult_spu.c

#include <moc.h>
#include <stdlib.h>
#include <math.h>
#include <stdio.h>

#include <sys/time.h>

#define MASTER  0

void initMatrix(double** M, int x, int y, double value)
    {
    int i, j;

    for (i = 0; i < x; i++)
	{
	for (j = 0; j < y; j++)
	    M[i][j] = value;
	}
    }

void initMatrixInc(double** M, int x, int y, double value, int offset)
    {
    int i, j;

    for (i = 0; i < x; i++)
	{
	for (j = 0; j < y; j++)
	    M[i][j] = (i+1) * value +( double)offset;
	    }
	}

int getChunkSize(int processes, int rank, int lines)
    {
    int result = lines/processes;

    if (lines%processes > rank)
	result++;

    return result;
    }

int getOffset(int processes, int rank, int lines)
    {
    int result = 0;
    int i;

    for (i = 0; i < rank; i++)
	result += getChunkSize(processes, i, lines);

    return result;
    }

double vecMult(double* vec1, double* vec2, int size)
    {
    double result = 0.0;
    int i;

    for (i = 0; i < size; i++)
	result += vec1[i] * vec2[i];

    return result;
    }

int main()
    {
    struct timeval t;
    struct timezone z;

    double sec = 0.0;

    gettimeofday(&t, &z);
    sec =( double)t.tv_usec/1e6 + t.tv_sec;

    int testFlag;

    int n, m, p; //dimensions (m x n) * (n x p) = (m x p) 
	    int i, j, k; //indeces 

	    int offset;

	    double** A; // matrices 
	    double** B;
	    double** C;

	    int np, myrank, chunksize;
	    MOC_Status status;

	    double startTime, computeTime, tmpTime, communicationTime;

	    communicationTime = 0.0;
	    computeTime = 0.0;

	    MOC_Init( 0, 0, NULL);

	    startTime = MOC_Wtime();

	    MOC_Comm_size(MOC_COMM_WORLD, &np); /*get #(processors)*/
	    MOC_Comm_rank(MOC_COMM_WORLD, &myrank);/*find rank of this process*/

	    m = n = p = 150;

	    // B allocation, each column is a vector!!!!  
	    B = (double**)malloc( p * sizeof(double*) );

	    for(i = 0; i < p; i++ )
		{
		B[i] = (double*)malloc( n * sizeof(double) );
		}

	    initMatrix(B, p, n, 0.0);

	    chunksize = getChunkSize(np, myrank, m);
	    printf("RANK[%d] chunksize: %d\n", myrank, chunksize);

	    if(myrank > 0)
		{

		tmpTime = MOC_Wtime();

		communicationTime += MOC_Wtime() - tmpTime;

		for(i = 0; i < p;i++) //get Matrix B from MASTER

		    {
		    tmpTime = MOC_Wtime();
		    MOC_Recv(B[i], n ,MOC_DOUBLE,MASTER,MASTER,MOC_COMM_WORLD,&status);
		    communicationTime += MOC_Wtime() - tmpTime;
		    }

		printf("RANK[%d] has Matrix B\n", myrank);

		A = (double**)malloc(chunksize * sizeof(double)); //a chunk of A 
		for(i = 0; i < chunksize; i++)
		A[i] = (double*)malloc(n * sizeof(double));

		C = (double**)malloc(chunksize * sizeof(double)); //a chunk of A 
		for(i = 0; i < chunksize; i++)
		C[i] = (double*)malloc(p * sizeof(double));

		offset = getOffset(np, myrank, m);

		initMatrix(A, chunksize, n, 0.0);
		initMatrix(C, chunksize, p, 0.0);

		printf("RANK[%d] empfange daten\n",myrank);

		for(k=0; k < chunksize; k++) // Recv chunk of A

		    {
		    tmpTime = MOC_Wtime();

		    MOC_Recv(A[k],n ,MOC_DOUBLE,MASTER,myrank,MOC_COMM_WORLD,&status);

		    communicationTime += MOC_Wtime() - tmpTime;
		    }

		computeTime = MOC_Wtime();
		printf("RANK[%d] bearbeite daten\n",myrank);
		for(i = 0; i < chunksize; i++) // compute

		    {
		    for(j = 0; j < p; j++)
			{
			C[i][j] = vecMult(A[i], B[j], n);
			}
		    }

		computeTime = MOC_Wtime() - computeTime;

		printf("RANK[%d] sende daten\n",myrank);
		for(k = 0; k < chunksize;k++)
		    {
		    tmpTime = MOC_Wtime();

		    testFlag = MOC_Send(C[k], p, MOC_DOUBLE, MASTER, myrank, MOC_COMM_WORLD);
		    if(testFlag != 0)
			{
			k--;
			}
		    communicationTime += MOC_Wtime() - tmpTime;
		    }

		for(i = 0; i < chunksize; i++)
		free(A[i]);
		free(A);

		for(i = 0; i < chunksize; i++)
		free(C[i]);
		free(C);

		}

	    for(i = 0; i < chunksize; i++)
	    free(B[i]);

	    free(B);

	    printf("RANK[%d] Communication Time: %f\n", myrank, communicationTime);
	    printf("RANK[%d] Computation Time: %f\n", myrank, computeTime);

	    MOC_Finalize();

	    gettimeofday(&t, &z);
	    sec = (double)t.tv_usec/1e6 + t.tv_sec - sec;

	    printf("RANK[%d] Complete time = %f\n",myrank, sec);

	    return 0;
	    }