#include <stdio.h>
#include <stdlib.h>
#include "mmio.h"
#include "math.h"
#include <sys/time.h>
#include "omp.h"
#include <unistd.h>

double wtime()
{
    struct timeval t;
    gettimeofday(&t, NULL);
    return (double)t.tv_sec + (double)t.tv_usec * 1E-6;
}

void quicksort(double* a, double* vindex, int* rindex, int* cindex, int n)
{
    int i, j, m;
    double p, t, s;
    if (n < 2)
        return;
    p = vindex[n / 2];

    for (i = 0, j = n - 1;; i++, j--) {
        while (vindex[i]<p)
            i++;
        while (p<vindex[j])
            j--;
        if (i >= j)
            break;
        t = a[i];
        a[i] = a[j];
        a[j] = t;

        s = vindex[i];
        vindex[i] = vindex[j];
        vindex[j] = s;

        m = rindex[i];
        rindex[i] = rindex[j];
        rindex[j] = m;

        m = cindex[i];
        cindex[i] = cindex[j];
        cindex[j] = m;
    }
    quicksort(a, vindex, rindex, cindex, i);
    quicksort(a + i, vindex + i, rindex + i, cindex + i, n - i);
}

int main(int argc, char *argv[])
{
    int j;
    int i;
    double t;
    int ret_code;
    MM_typecode matcode;
    FILE *f;
    int M, N, nz, *row, *col;
    double *val, *vector, *result, *result_t, *check_result, *check_result_t;

    if (argc < 2)
    {
        fprintf(stderr, "Usage: %s [martix-market-filename]\n", argv[0]);
        exit(1);
    }
    else
    {
        if ((f = fopen(argv[1], "r")) == NULL)
            exit(1);
    }

    if (mm_read_banner(f, &matcode) != 0)
    {
        printf("Could not process Matrix Market banner.\n");
        exit(1);
    }

    if (mm_is_complex(matcode) && mm_is_matrix(matcode) &&
            mm_is_sparse(matcode) )
    {
        printf("Sorry, this application does not support ");
        printf("Market Market type: [%s]\n", mm_typecode_to_str(matcode));
        exit(1);
    }

    if ((ret_code = mm_read_mtx_crd_size(f, &M, &N, &nz)) !=0)
        exit(1);

    row = (int *) malloc(nz * sizeof(int));
    col = (int *) malloc(nz * sizeof(int));
    val = (double *) malloc(nz * sizeof(double));

    for (int i = 0; i < nz; i++)
    {
        fscanf(f, "%d %d %lg\n", &row[i], &col[i], &val[i]);
        row[i]--;
        col[i]--;
    }

    if (f !=stdin) fclose(f);

//->> Сортировка строк
    double *vIndex;
    vIndex = (double*)malloc(nz*sizeof(double));
    memset(vIndex, 0, nz*sizeof(double));
    for (int i = 0; i < nz; i++)
    {
        vIndex[i] = (double)row[i] * N + col[i];
        if (vIndex[i] < 0)
        {
            printf("Error!\n");
            exit(1);
        }
    }

    quicksort(val, vIndex, row, col, nz);

//->> Конвертер
    int *RowIndex;
    RowIndex = (int*)malloc((N + 1)*sizeof(int));
    RowIndex[0] = 0;
    int row1 = 0;
    int num, num2;
    i = 0;
    for (i = 1; i < nz; i++){
        if (row[i] != row[i-1]){
            row1++;
            RowIndex[row1] = i;
            num = row[i] - row[i-1];
            if (num > 1){
                num2 = row1 + num - 1;
                for (int j = row1; j < num2; j++){
                    row1++;
                    RowIndex[row1] = i;
                }
            }
        }
    }
    row1++;
    RowIndex[row1] = i;

    t = wtime();
    vector = (double*)malloc(sizeof(double)*nz);
    result = (double*)malloc(sizeof(double)*nz);
    result_t = (double*)malloc(sizeof(double)*nz);
    for (int i = 0; i < N; i++)
    {
        vector[i] = i + 2.00;
    }

//************************************************************* умножение *************************************************



            #pragma omp parallel for private(i, j)
            for(i = 0; i < N; i++ )
            {
                result[i] = 0;
                for (j = RowIndex[i]; j < RowIndex[i + 1]; j++)
                    result[i] += val[j] * vector[col[j]];

            }




//----------------------------------------------------------------------------------------------------------------------------------------
//************************************************************* ТРАНСПОНИРОВАНИЕ *************************************************

        for (int i = 0; i < N; ++i)
        {
            result_t[i] = 0.0;
        }

        #pragma omp parallel for private(i, j)
        for (i = 0; i < N; i++)
        {

            for(j = RowIndex[i]; j < RowIndex[i+1]; j++)
            {
                if (i != col[j])
                {

                    result_t[col[j]] += val[j] * vector[i];
                }
            }
        }


//----------------------------------------------------------------------------------------------------------------------------------------

//************************************************************* проверка *************************************************
    check_result = (double*)malloc(sizeof(double)*nz);
    check_result_t = (double*)malloc(sizeof(double)*nz);
    for (i = 0; i < N; i++)
    {
        check_result[i] = 0;
        check_result_t[i] = 0;
        for ( j=RowIndex[i]; j<RowIndex[i + 1]; j++)
        {
            check_result[i] += val[j] * vector[col[j]];
            if (i != col[j])
                check_result_t[col[j]] += val[j] * vector[i];

        }
    }

    ///check multiplication
    for (int i = 0; i < N; ++i)
    {
        if (result[i] == check_result[i])
            continue;
        else
        {
            printf("error-->>%d\n", i);
            printf("check_result_t[%d] = %0.12f\t  result[%d] = %0.12f\n",i, check_result[i], i, result[i]);
            printf("Vectors not match\n");
            break;
        }
    }

    ///check transposition
    for (int i = 0; i < N; ++i)
    {
        if (result_t[i] == check_result_t[i])
            continue;
        else
        {
            printf("error-->>%d\n", i);
            printf("check_result_t[%d] = %0.12f\t  result_t[%d] = %0.12f\n",i, check_result_t[i], i, result_t[i]);
            printf("Vectors not match\n");
            break;
        }
    }
    return 0;
}