//
// 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;
}