Распараллеливание метода Жордана с поиском главного элемента по всей матрице(мьютексы)

//
//  main.cpp
//  serg_инфа
//
//  Created by Антон Лоскутов on 07/12/14.
//  Copyright (c) 2014 Антон Лоскутов. All rights reserved.
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <unistd.h>
#include <pthread.h>

typedef struct Parameters {
    double *a, c;
    int i, j, k, n;
} Parameters;

int *IndexRow, *IndexCol;

int indexRow(int);
int indexCol(int);
void changeRow(int, int);
void changeCol(int, int);
void initIndexRow(int);
void initIndexCol(int);

double* initMatrix(int);
double* initMatrix2(double*, int);
int MatrixSubRow(double*, int, int, int, double);
int MatrixDivRow(double*, int, int, double);
double* MatrixJordan(double*, int);

void* run_sub(void* arg);
void* run_div(void* arg);

static pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;

double eps = 1e-10;

int main(int argc, const char * argv[]) {
    double *res;
    int i, j, n;
    printf("n = ");
    scanf("%d",&n);
    initIndexRow(n);
    initIndexCol(n);
    double *a;
    a = initMatrix(n);
    res = MatrixJordan(a, n);
    if (res == NULL)
        return 1;
    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++)
            printf("%5.2lf ", res[i*n+j]);
        printf("\n");
    }
    free(IndexRow);
    free(IndexCol);
    free(a);
    free(res);
    return 0;
}

int indexRow(int a) {
    return IndexRow[a];
}

int indexCol(int a) {
    return IndexCol[a];
}

void changeRow(int a, int b) {
    int dub = IndexRow[a];
    IndexRow[a] = IndexRow[b];
    IndexRow[b] = dub;
}

void changeCol(int a, int b) {
    int dub = IndexCol[a];
    IndexCol[a] = IndexCol[b];
    IndexCol[b] = dub;
}

void initIndexRow(int n) {
    int i;
    IndexRow = (int*)malloc(n*sizeof(int));
    for (i = 0; i < n; i++)
        IndexRow[i] = i;
}

void initIndexCol(int n) {
    int i;
    IndexCol = (int*)malloc(n*sizeof(int));
    for (i = 0; i < n; i++)
        IndexCol[i] = i;
}

double* initMatrix(int n) {
    double *res = (double*)calloc(n*n, sizeof(double));
    int i,j;
    for (i = 0; i < n; i++)
        for (j = 0; j < n; j++)
            res[i*n+j] = 1./(i+j+1)+1;
    //scanf("%lf", &res[i*n+j]);
    return res;
}

double* initMatrix2(double *a, int n) {
    int i;
    double *res = (double*)calloc(n*n, sizeof(double));
    for (i = 0; i < n*n; i++)
        res[i] = a[i];
    return res;
}

int MatrixSubRow(double *a, int n, int i, int j, double c) {
    int k, res;
    if (fabs(c) < eps)
        return -1;
    Parameters *masParameters;
    masParameters = (Parameters*)malloc(n*sizeof(Parameters));
    pthread_t *thread;
    thread = (pthread_t*)malloc(n*sizeof(pthread_t));
    for (k = 0; k < n; k++) {
        masParameters[k].a = a;
        masParameters[k].i = i*n+k;
        masParameters[k].j = j*n+k;
        masParameters[k].c = c;
        res = pthread_create(&thread[k], NULL, &run_sub, &masParameters[k]);
        if (res != 0) {
            printf("Cannot create thread %d", i+1);
            exit(1);
        }
    }
    for (k = 0; k < n; k++)
        pthread_join(thread[k], NULL);
    free(masParameters);
    free(thread);
    return 0;
}

int MatrixDivRow(double *a, int n, int i, double c) {
    int k, res;
    if (fabs(c) < eps)
        return -1;
    Parameters *masParameters;
    masParameters = (Parameters*)malloc(n*sizeof(Parameters));
    pthread_t *thread;
    thread = (pthread_t*)malloc(n*sizeof(pthread_t));
    for (k = 0; k < n; k++) {
        masParameters[k].a = a;
        masParameters[k].i = i*n+k;
        masParameters[k].c = c;
        res = pthread_create(&thread[k], NULL, &run_div, &masParameters[k]);
        if (res != 0) {
            printf("Cannot create thread %d", i+1);
            exit(1);
        }
    }
    for (k = 0; k < n; k++)
        pthread_join(thread[k], NULL);
    pthread_mutex_destroy(&mutex);
    free(masParameters);
    free(thread);
    return 0;
}

double* MatrixJordan(double* mas, int n) {
    double max, buf, *res, *a = initMatrix2(mas, n);
    int i, j, k, indexColMax, indexRowMax;
    res = (double*)calloc(n*n, sizeof(double));
    for (i = 0; i < n; i++)
        res[i*n+i] = 1.;
    for (i = 0; i < n; i++) {
        max = fabs(a[indexCol(i)*n+indexRow(i)]);
        indexColMax = i;
        indexRowMax = i;
        for (j = i; j < n; j++) {
            for (k = i; k < n; k++) {
                if (fabs(a[indexCol(j)*n+indexRow(k)]) > max) {
                    max = fabs(a[indexCol(j)*n+indexRow(k)]);
                    indexRowMax = k;
                    indexColMax = j;
                }
            }
        }
        if (fabs(max) < eps) {
            printf("Determinant is 0\n");
            return NULL;
        }
        if (i != indexRowMax)
            changeRow(i, indexRowMax);
        if (i != indexColMax)
            changeCol(i, indexColMax);
        buf = a[indexCol(i)*n + indexRow(i)];
        if (MatrixDivRow(a, n, indexCol(i), buf) != 0) {
            printf("Panic...\n");
            return NULL;
        }
        if (MatrixDivRow(res, n, indexCol(i), buf) != 0) {
            printf("Panic!\n");
            return NULL;
        }
        for (j = 0; j < i; j++) {
            buf = a[indexCol(j)*n + indexRow(i)];
            if (fabs(buf) < eps) {
                printf("cont\n");
                continue;
            }
            if (MatrixSubRow(a, n, indexCol(j), indexCol(i), buf) != 0) {
                printf("Panic.1\n");
                return NULL;
            }
            if (MatrixSubRow(res, n, indexCol(j), indexCol(i), buf) != 0) {
                printf("Panic.2\n");
                return NULL;
            }
        }
        for (j = i+1; j < n; j++) {
            buf = a[indexCol(j)*n + indexRow(i)];
            if (fabs(buf) < eps) {
                printf("cont2\n");
                continue;
            }
            if (MatrixSubRow(a, n, indexCol(j), indexCol(i), buf) != 0) {
                printf("Panic.3\n");
                return NULL;
            }
            if (MatrixSubRow(res, n, indexCol(j), indexCol(i), buf) != 0) {
                printf("Panic.4\n");
                return NULL;
            }
        }
    }
    free(a);
    return res;
}

void* run_div(void* arg) {
    int i;
    double *a, c;
    Parameters param = *((Parameters *) arg);
    
    c = param.c;
    i = param.i;
    a = param.a;
    
    pthread_mutex_lock(&mutex);
    a[i] /= c;
    pthread_mutex_unlock(&mutex);
    return NULL;
}

void* run_sub(void* arg) {
    int i, j;
    double *a, c;
    Parameters param = *((Parameters *) arg);
    
    a = param.a;
    i = param.i;
    j = param.j;
    c = param.c;
            
    pthread_mutex_lock(&mutex);
    a[i] -= c*a[j];
    pthread_mutex_unlock(&mutex);
    return NULL;
}