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

#include <iostream>
using namespace std;

int k=0;
double temp;
double a=0, b=1, c=2, nondiag, diag;
int i,j,N=500,M=500,nonzero=200;

void createMatrix(double** Matrix){
        for(i=0; i<N; i++) {
                for(j=0; j<nonzero; j++) {
                        int val = rand()%N;
                        while(val==i) {
                                val=rand()%N;
                        }
                        nondiag = a + rand() * (b-a) / RAND_MAX;
                        Matrix[i][val] = nondiag;
                }
        }

        for(i=0; i<N; i++) {
                diag = b + rand() * (c-b) / RAND_MAX;
                Matrix[i][i] = diag;
        }

}

void printMatrix(double** Matrix){
        for(i=0; i<N; i++) {
                for(j=0; j<M; j++) {
                        printf("%f ",Matrix[i][j]);
                }
                printf("\n");
        }
}

void printMatrix_Q(double** Matrix, int rank){
        printf("%d\n",rank);
        for(i=0; i<N/2; i++) {
                for(j=0; j<M/2; j++) {
                        printf("%f ",Matrix[i][j]);
                }
                printf("\n");
        }
}

void multiplyMatrixByVector(double** Matrix, double* Result, double* Vector){
        double sum = 0;
        for(i=0; i<N; i++) {
                for(j=0; j<M; j++) {
                        sum+=Matrix[i][j]*Vector[j];
                }
                Result[i] = sum;
                sum = 0;
        }
}


void generateVector(double* Vector){
        for(i=0; i<N; i++) {
                Vector[i] = 1 + rand() * (5-1) / RAND_MAX;
        }
}

void createCRS(double** Matrix, double** CRS, int* CRS_ROWS){
        int k = 0;
        double temp;
        for(i=0; i<N/2; i++) {
                for(j=0; j<M/2; j++) {
                        temp = Matrix[i][j];
                        if(temp!=0.0) {
                                CRS[k][0] = temp;
                                CRS[k++][1] = j;
                        }
                }
                CRS_ROWS[i] = k;
        }
}

void decompressCRS(double** Matrix_validate_CRS, double** CRS, int* CRS_ROWS){
        j=0;
        for(i=0; i<N; i++) {
                int k = CRS_ROWS[i];
                for(j; j<k; j++) {
                        Matrix_validate_CRS[i][(int)CRS[j][1]] = CRS[j][0];
                }
        }
}

void validateTwoMultipliedMatrixes(double* Matrix_one, double* Matrix_two, double* Matrix_three){
        int errors = 0;
        for(i=0; i<N; i++) {
                if(Matrix_one[i]!=Matrix_two[i] || Matrix_one[i]!=Matrix_three[i]||Matrix_two[i]!=Matrix_three[i])
                        errors++;
        }
        if(errors==0)
                printf("Both multiplied matrixes are the same\n");
}

void printCRS(double** CRS){
        printf("\n\nCRS \n");
        for(i=0; i<(M * (nonzero + 1)); i++) {
                for(j=0; j<2; j++) {
                        printf("%lf ", CRS[i][j]);
                }
                printf("\n");
        }
}

void printResult(double* result){
        for(i=0; i<N; i++)
                printf("%lf ",result[i]);

        printf("\n");
}

void multiplyCRSbyVector(double **CRS, int* CRS_ROWS, double* Vector, double* multiply_result, int rank){
        int vecshift=0,rowshift=0;
        double sum=0.0;
        j=0;

        if(rank==1)
                vecshift = N/2;
        else if(rank==2)
                rowshift = N/2;
        else if(rank==3) {
                rowshift = N/2;
                vecshift = N/2;
        }
        for(i=0; i<N/2; i++) {
                for(j; j<CRS_ROWS[i]; j++) {
                        sum+=CRS[j][0] * Vector[vecshift+(int)CRS[j][1]];
                }
                multiply_result[i+rowshift] = sum;
                sum=0.0;
        }
}

void compareTwoResults(double* resultSeq, double* resultPar){
        int errors=0;
        for(i=0; i<N; i++) {
                if(resultSeq[i]!=resultPar[i])
                        errors++;
        }
        printf("Erros: %d\n", errors);
}

int main(int argc, char **argv)
{

        int rank,ranksent,size,source,dest,tag,i,len;
        char name[20];
        len=20;

        MPI_Status status;
        MPI_Init(&argc,&argv);
        MPI_Comm_rank(MPI_COMM_WORLD,&rank);
        MPI_Comm_size(MPI_COMM_WORLD,&size);

        double **Matrix, **CRS;
        int *CRS_ROWS;
        double *Vector;

        Matrix = (double**)malloc(N * sizeof(double*));
        CRS = (double**)malloc(N * (nonzero+1) * sizeof(double*));
        CRS_ROWS = (int*)malloc(N * sizeof(int));
        Vector = (double*)malloc(N * sizeof(double));

        for(i=0; i<N; i++) {
                Matrix[i] = (double*)malloc(M*sizeof(double));
        }

        for(i=0; i<(N*(nonzero+1)); i++) {
                CRS[i] = (double*)malloc(2*sizeof(double));
        }

        double **p0,**p1,**p2,**p3;
        double **dest_matrix;
        double **crs_q;
        int *crs_rows_q;
        double *multiply_result;
        double *result_all;
        double *result_seq;



        p0 = (double**)malloc(N/2*sizeof(double*));
        p1 = (double**)malloc(N/2*sizeof(double*));
        p2 = (double**)malloc(N/2*sizeof(double*));
        p3 = (double**)malloc(N/2*sizeof(double*));
        dest_matrix = (double**)malloc(N/2*sizeof(double*));

        crs_q = (double**)malloc(N/2 * (nonzero+1) * sizeof(double*));
        crs_rows_q = (int*)malloc(N/2 * sizeof(int));
        multiply_result = (double*)malloc(N * sizeof(double));
        result_all = (double*)malloc(N*sizeof(double));
        result_seq = (double*)malloc(N*sizeof(double));

        for(i=0; i<N/2; i++) {
                p0[i] = (double*)malloc(M/2*sizeof(double));
                p1[i] = (double*)malloc(M/2*sizeof(double));
                p2[i] = (double*)malloc(M/2*sizeof(double));
                p3[i] = (double*)malloc(M/2*sizeof(double));
                dest_matrix[i] = (double*)malloc(M/2*sizeof(double));
        }

        for(i=0; i<(N/2 * (nonzero+1)); i++) {
                crs_q[i] = (double*)malloc(2*sizeof(double));
        }

        if(rank==0)
        {
                createMatrix(Matrix);
                // printMatrix(Matrix);
                generateVector(Vector);
                multiplyMatrixByVector(Matrix,result_seq,Vector);

                dest=0;
                tag=0;
                for(i=0; i<N/2; i++) {
                        for(j=0; j<M/2; j++) {
                                dest_matrix[i][j] = Matrix[i][j];
                                p1[i][j] = Matrix[i][j+(N/2)];
                                p2[i][j] = Matrix[i+(M/2)][j];
                                p3[i][j] = Matrix[i+(N/2)][j+(M/2)];
                        }
                }
                // MPI_Send(&(p1[0][0]),(N/2*M/2)*sizeof(double),MPI_DOUBLE,1,10,MPI_COMM_WORLD);
                // MPI_Send(&(p2[0][0]),(N/2*M/2)*sizeof(double),MPI_DOUBLE,2,10,MPI_COMM_WORLD);
                // MPI_Send(&(p3[0][0]),(N/2*M/2)*sizeof(double),MPI_DOUBLE,3,10,MPI_COMM_WORLD);
                // for(int d = 0; d < N/2; d++) {
                //         MPI_Send(&(p1[d][0]),M/2,MPI_DOUBLE,1,10,MPI_COMM_WORLD);
                //         MPI_Send(&(p2[d][0]),M/2,MPI_DOUBLE,2,10,MPI_COMM_WORLD);
                //         MPI_Send(&(p3[d][0]),M/2,MPI_DOUBLE,3,10,MPI_COMM_WORLD);
                // }

        }
        MPI_Barrier(MPI_COMM_WORLD);

        if(rank==0) {
                MPI_Send(&(p1[0][0]),(N/2*M/2)+sizeof(double),MPI_DOUBLE,1,10,MPI_COMM_WORLD);
                MPI_Send(&(p2[0][0]),(N/2*M/2)+sizeof(double),MPI_DOUBLE,2,10,MPI_COMM_WORLD);
                MPI_Send(&(p3[0][0]),(N/2*M/2)+sizeof(double),MPI_DOUBLE,3,10,MPI_COMM_WORLD);
        }
        else{
                MPI_Recv(&(dest_matrix[0][0]),(N/2*M/2)+sizeof(double),MPI_DOUBLE,0,10,MPI_COMM_WORLD,&status);
                // for(int d = 0; d < N/2; d++) {
                //         MPI_Recv(&(dest_matrix[d][0]),M/2,MPI_DOUBLE,0,10,MPI_COMM_WORLD,&status);
                // }
        }

        MPI_Bcast(Vector,N,MPI_DOUBLE,0,MPI_COMM_WORLD);

        // printMatrix_Q(dest_matrix,rank);

        createCRS(dest_matrix,crs_q,crs_rows_q);
        MPI_Barrier(MPI_COMM_WORLD);

        multiplyCRSbyVector(crs_q,crs_rows_q,Vector, multiply_result,rank);
        MPI_Barrier(MPI_COMM_WORLD);


        MPI_Reduce(multiply_result,result_all,N,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
        MPI_Barrier(MPI_COMM_WORLD);

        if(rank==0) {
                compareTwoResults(result_seq,result_all);
        }

        MPI_Finalize();
        return(0);
}