/* Nhan ma tran dung MPI */


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

#define ROW_LOW(id,p,n)         ((id)*(n)/(p))
#define ROW_HIGH(id,p,n)        (ROW_LOW((id)+1,p,n)-1)
#define ROW_SIZE(id,p,n)        (ROW_HIGH(id,p,n)-ROW_LOW(id,p,n)+1)
#define ROW_OWNER(j,p,n)        (((p)*((j)+1)-1)/(n))
#define THRESHOLD 65535

static double *A,*B,*C;

void MM(int L, int M, int N,
                double *A, double *B, double *C,
                int crow, int ccol , /*corner of C block*/
                int arow, int acol , /* corner of A block*/
                int brow, int bcol , /* corner of B block*/
                int l,
                int m,          /* block B is m x n*/
                int n)          /* block C is l x n*/
{
        int lhalf[3], mhalf[3], nhalf[3];
        int i,j,k;
        double *aptr,*bptr,*cptr;

        if (m*n > THRESHOLD) {
        lhalf[0] = 0; lhalf[1] = l/2; lhalf[2] = l - l/2;
        mhalf[0] = 0; mhalf[1] = m/2; mhalf[2] = m - m/2;
        nhalf[0] = 0; nhalf[1] = n/2; nhalf[2] = n - n/2;

        for (i=0; i<2; i++)
                for (j=0;j<2;j++)
                for (k=0;k<2;k++)
                        MM(L,M,N, A, B, C,
                        crow + lhalf[i], ccol + mhalf[j],
                        arow + lhalf[i], acol + mhalf[k],
                        brow + mhalf[k], bcol + nhalf[j],
                        lhalf[i+1], mhalf[k+1], nhalf[j+1] );
        }else {


                for (i=0; i<l ; i++)
                        for (j=0; j<n; j++) {
                                cptr = C +(crow+i)*N + ccol + j;
                                aptr = A + (arow+i)*M + acol;
                                bptr = B + brow*N + bcol + j;
                                for (k=0; k<m; k++) {
                                        *cptr += *(aptr++) * *bptr;
                                        bptr += N;
                                }
                        }
        }
}

void print_usage(char *s) {
        printf("Usage: mpirun -np <N> %s <l> <m> <n> \n",s);
        printf("\t <N>: number of processes\n");
        printf("\t Matrix size: A=lxm, B=mxn, C=lxn\n");
}

void print_to_file(double *mt, int row, int col, char *name) {
        FILE *f;
        f= fopen(name,"w");
        if (f == NULL){
                printf("Cannot write to file %s\n",name);
                exit(1);
        }
        int i,j;
        for (i=0;i<row;i++) {
                for (j=0;j<col;j++)
                        fprintf(f,"%6.0f ", *(mt + i*col +j));
                fprintf(f,"\n");
        }
        fclose(f);
}
void print(double *mt, int row, int col) {
        int i,j;
        for (i=0;i<row;i++) {
                for (j=0;j<col;j++)
                        printf("%6.0f ", *(mt + i*col +j));
                        printf("\n");
        }
}

void init_matrix(double *MT, int row, int col,int seed) {
        int i,j;
        for (i=0;i<row;i++) {
                for (j=0;j<col;j++) {
                        MT[i*col+j] = 1;
                }
        }
}

int main(int argc, char **argv) {
        int i,j;
        int id; /* process id */
        int p; /* number of processes */

        int localL; /* number of rows each node does the computation*/
        MPI_Status status;
        MPI_Init(&argc, &argv);
        MPI_Barrier(MPI_COMM_WORLD);
        double elapsed_time = -MPI_Wtime();

        MPI_Comm_rank(MPI_COMM_WORLD, &id);
        MPI_Comm_size(MPI_COMM_WORLD, &p);

        if (argc !=4) {
                print_usage(argv[0]);
                exit(1);
        }

        int L,M,N;
        L=atoi(argv[1]);
        M=atoi(argv[2]);
        N=atoi(argv[3]);

        /*compute the number of rows each one has*/
        localL = ROW_SIZE(id,p,L);
        B=(double *) malloc(sizeof(double)*M*N);
        if (id==0)
                init_matrix(B,M,N,0);

        MPI_Barrier(MPI_COMM_WORLD);
        /* send B */
        MPI_Bcast(B,
                        M * N,
                        MPI_DOUBLE,
                        0,
                        MPI_COMM_WORLD);
        /* on root node */
        if (id==0) {

                A=(double *) malloc(sizeof(double)*L*M);
                C=(double *) malloc(sizeof(double)*L*N);

                /* initialize A,B */
                init_matrix(A,L,M,1);
                /*
                printf("A=\n");
                print (A,L,M);
                printf("B=\n");
                print (B, M, N);        
                */
                /* send A, B to other nodes */
                /* each node will receive:
                 * - entire matrix B
                 * - L/p rows of matrix A
                 * and generate:
                 * - L/p rows of C
                 */

                 /* SEND DATA TO OTHER NODES */
                 for (i=1; i<p;i++) {
                        int size = ROW_SIZE(i,p,L);
                        int low = ROW_LOW(i,p,L);
                        //printf(" send %d to %d \n",size,i);
                        fflush(stdout);
                        /* send rows of A */
                        //printf("first data should be %2.0f\n", *(A+(i)*M));
                        MPI_Send( A + low* M, /*pointer to the beginning of the buffer to send */
                                                size * M, /* size of the buffer */
                                                MPI_DOUBLE , /*data type*/
                                                i, /* send to process number i */
                                                100,
                                                MPI_COMM_WORLD);

                 }
                 /* DOES ITS OWN COMPUTATION */
                 MM(L,M,N,A,B,C,0,0,0,0,0,0,localL,M,N);
                 /* RECEIVE DATA FROM OTHER NODES */

                 for (i=1; i<p;i++) {
                         int coffset = ROW_LOW(i,p,L);
                         int size = ROW_SIZE(i,p,L);
                         MPI_Recv(      C + coffset*N ,       /* pointer to receiving buffer */
                                        size * N, /*buffer size */
                                        MPI_DOUBLE,
                                        i, /*send to root node */
                                        200, /* message tag */
                                        MPI_COMM_WORLD,
                                        &status);
                }

                /* PRINT RESULTS */
                printf("C=\n");
                print(C, L, N);    
				//print_to_file(double *mt, int row, int col, char *name) {	
				print_to_file(A, L, M, "inA");
				print_to_file(B, M, N, "inB");
				print_to_file(C, L, N, "outC");


        }

        if (id !=0) {
                A=(double *) malloc(sizeof(double)*localL*M);
                C=(double *) malloc(sizeof(double)*localL*N);

                memset( (double *)C, 0, localL*N);

                MPI_Recv(       A,       /* pointer to receiving buffer */
                                        localL * M, /*buffer size */
                                        MPI_DOUBLE,
                                        0, /*receive from root node */
                                        100,
                                        MPI_COMM_WORLD,
                                        &status);

                MM(L,M,N,A,B,C,0,0,0,0,0,0,localL,M,N);

                //printf("finish computing  on node %d\n", id);
                //fflush(stdout);
                /* send C back to root */
                MPI_Send(       C,       /* pointer to receiving buffer */
                                        localL * N, /*buffer size */
                                        MPI_DOUBLE,
                                        0, /*send to root node */
                                        200, /* message tag */
                                        MPI_COMM_WORLD);

        }

        free(A);
        free(B);
        free(C);

        elapsed_time += MPI_Wtime();
        if (id==0)
                printf("Elapsed time: %10.6f\n",elapsed_time);
        MPI_Finalize();
        return EXIT_SUCCESS;
}