Обратная Матрица

LudmilaUt

New member
16.01.2012
1
0
#1
Задача: обратить матрицу методом Жордана с выбором максимального элемента по всей матрице. (Метод Жордана - это приведение матрицы к диагональному виду; сначала выбирается максимальный элемент и ставится на первое место первой строки, потом удаляются элементы первого столбца со второго до конца. Потом то же самое для подматрицы, и т.д. потом то же самое для верхнего правого угла матрицы). Использовать параллельное программирование. Выбор главного элемента, обращение матрицы, функция TestMatrix ("невязка") должны быть распаралелены.
Моя программа не работает. Почему - не понятно.

/*Обращение матрицы методом Жордана с выбором элемента по матрице*/

#include <sys/resource.h>
#include <sys/time.h>
#include <stdio.h>
#include <math.h>
#include <pthread.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

#define MAX_OUTPUT_SIZE 8

void synchronize(int total_threads){
static pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;
static pthread_cond_t condvar_in = PTHREAD_COND_INITIALIZER;
static pthread_cond_t condvar_out = PTHREAD_COND_INITIALIZER;
static int threads_in = 0;
static int threads_out = 0;

pthread_mutex_lock(&mutex);

threads_in++;
if (threads_in >= total_threads){
threads_out = 0;
pthread_cond_broadcast(&condvar_in);
}
else
while (threads_in < total_threads)
pthread_cond_wait(&condvar_in,&mutex);

threads_out++;
if (threads_out >= total_threads){
threads_in = 0;
pthread_cond_broadcast(&condvar_out);
}
else
while (threads_out < total_threads)
pthread_cond_wait(&condvar_out,&mutex);

pthread_mutex_unlock(&mutex);
}

/*Выбор главного элемента в функции InvMatrix нужно тоже
распараллеливать. (Каждый тред выбирает "локально" максимальный
элемент в своей подматрице, складывает его в массив, затем нулевой
тред из локально максимальных элементов выбирает уже глобальный
максимум).*/

int InvMatrix(int n, double *a, double *x, int *index, int thread_num, int total_threads){
int i, j, k, l;
int j1, j2, k1, k2, s, t;
int first_row;
int last_row;
double tmp;

// if (thread_num == 0){
for (i = 0; i < n; i++)
index = i;

for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
x[i * n + j] = (double)(i == j);
// }

for (i = 0; i < n; i++){
first_row = (n - i - 1) * thread_num;
first_row = first_row/total_threads + i + 1;
last_row = (n - i - 1) * (thread_num + 1);
last_row = last_row/total_threads + i + 1;

k1 = i;
k2 = i;

for(j1=first_row; j1<last_row; j1++){
for (j2=i; j2<n; j2++){
if (fabs(a[k1*n + k2]) < fabs(a[j1*n + j2])){
a[n*n + 2*thread_num] = j1;
a[n*n + (2*thread_num + 1)] = j2;
}
else{
a[n*n + 2*thread_num] = i;
a[n*n + (2*thread_num + 1)] = i;
}
}
}
synchronize(total_threads);

if (thread_num == 0){
k1=i;
k2=i;
for(l=0; l<2*total_threads; l++){
s = a[n*n + 2*thread_num];
t = a[n*n + (2*thread_num + 1)];
if(fabs(a[k1*n + k2]) < fabs(a[s*n + t])){
k1 = s;
k2 = t;
}
l++;
}
}
synchronize(total_threads);

for (i = 0; i < n; i++){
/* if (thread_num == 0){
k1 = i;
k2 = i;
for (j1 = i; j1 < n; j1++)
for (j2 = i; j2 < n; j2++)
if (fabs(a[k1 * n + k2]) < fabs(a[j1 * n + j2])){
k1 = j1;
k2 = j2;
}
*/
j = index;
index = index[k2];
index[k2] = j;

for (j = 0; j < n; j++){
tmp = a[i * n + j];
a[i * n + j] = a[k1 * n + j];
a[k1 * n + j] = tmp;
}

for (j = 0; j < n; j++){
tmp = a[j * n + i];
a[j * n + i] = a[j * n + k2];
a[j * n + k2] = tmp;
}

for (j = 0; j < n; j++){
tmp = x[i * n + j];
x[i * n + j] = x[k1 * n + j];
x[k1 * n + j] = tmp;
}

// tmp = a[i * n + i];
tmp = 1.0/a[i * n + i];
for (j = i; j < n; j++)
a[i * n + j] *= tmp;
for (j = 0; j < n; j++)
x[i * n + j] *= tmp;
}
synchronize(total_threads);

first_row = (n - i - 1) * thread_num;
first_row = first_row/total_threads + i + 1;
last_row = (n - i - 1) * (thread_num + 1);
last_row = last_row/total_threads + i + 1;

for (j = first_row; j < last_row; j++){
tmp = a[j * n + i];
for (k = i; k < n; k++)
a[j * n + k] -= tmp * a[i * n + k];
for (k = 0; k < n; k++)
x[j * n + k] -= tmp * x[i * n + k];
}
synchronize(total_threads);
}

first_row = n * thread_num;
first_row = first_row/total_threads;
last_row = n * (thread_num + 1);
last_row = last_row/total_threads;

for (k = first_row; k < last_row; k++)
for (i = n - 1; i >= 0; i--){
tmp = x[i * n + k];
for (j = i + 1; j < n; j++)
tmp -= a[i * n + j] * x[j * n + k];
x[i * n + k] = tmp;
}
synchronize(total_threads);

if (thread_num == 0){
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
a[index * n + j] = x[i * n + j];

for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
x[i * n + j] = a[i * n + j];
}

return 0;
}

double f(int i, int j){
return 1.0/(sqrt(sqrt(i) + sqrt(j)) + 1.0);
// return 1.0/(i+j+1);
// return fabs(i-j) + 1;
}

void InputMatrix(int n, double *a, int mode, FILE *fp){
int i, j;

if (mode == 1){
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
fscanf(fp, "%lf", &a[i * n + j]);
fclose(fp);
}
else{
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
a[i * n + j] = f(i, j);
a[n*n + i] = 0.0;
}
}

void OutputMatrix(int n, double *a){
int i, j;
int m = (n > MAX_OUTPUT_SIZE) ? MAX_OUTPUT_SIZE : n;

for (i = 0; i < m; i++){
printf("| ");
for (j = 0; j < m; j++)
printf("%10.3g ", a[i * n +j]);
printf(" |\n");
}
}

double TestMatrix(int n, double *a, double *x){
int i, j, k;
double tmp;
double rezult;

rezult = 0.0;
for (i = 0; i < n; i++)
for (j = 0; j < n; j++){
tmp = 0.0;
for (k = 0; k < n; k++)
tmp += a[i * n + k] * x[k * n + j];
if (i == j)
tmp -= 1.0;
rezult += tmp * tmp;
}

return sqrt(rezult);
}

long int get_time(void){
struct rusage buf;

getrusage(RUSAGE_SELF, &buf);

return buf.ru_utime.tv_sec * 100 + buf.ru_utime.tv_usec/10000;
}

long int get_full_time(void){
struct timeval buf;

gettimeofday(&buf, 0);

return buf.tv_sec * 100 + buf.tv_usec/10000;
}

typedef struct{
int n;
double *a;
double *x;
int *index;
int thread_num;
int total_threads;
} ARGS;

long int thread_time = 0;
pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;

void *Inversion(void *p_arg){
ARGS *arg = (ARGS*)p_arg;
long int t1;

t1 = get_time();
InvMatrix(arg->n, arg->a, arg->x, arg->index, arg->thread_num, arg->total_threads);
t1 = get_time() - t1;

pthread_mutex_lock(&mutex);
thread_time += t1;
pthread_mutex_unlock(&mutex);

return NULL;
}

int main(int argc, char **argv){
int i;
int n;
double *a;
double *x;
int *index;
int mode;
FILE *fp;
long int t_full;
int total_threads;
pthread_t *threads;
ARGS *args;

if (argc > 1)
total_threads = atoi(argv[1]);
else{
printf("Incorrect usage!\n");
return -1;
}

input = NULL;

if (argc == 3)
n = atoi(argv[2]);
else{
printf("nado vvodite pravilno \n");
return -3;
}
if (n <= 0){
printf("Incorrect n!\n");
return -4;
}

a = (double*)malloc(n*n * sizeof(double));
x = (double*)malloc(n*n * sizeof(double));
index = (int*)malloc(n*sizeof(int));
threads = (pthread_t*)malloc(total_threads * sizeof(pthread_t));
args = (ARGS*)malloc(total_threads * sizeof(ARGS));

if (!(a && x && index && threads && args)){
printf("Not enough memory!\n");

if (a) free(a);
if (x) free(x);
if (index) free(index);
if (threads) free(threads);
if (args) free(args);

return -4;
}

InputMatrix(n, a, mode, fp);
for (i = 0; i < total_threads; i++){
args.n = n;
args.a = a;
args.x = x;
args.index = index;
args.thread_num = i;
args.total_threads = total_threads;
}

t_full = get_full_time();

for (i = 0; i < total_threads; i++)
if (pthread_create(threads + i, 0, Inversion, args + i)){
printf("Cannot create thread %d!\n", i);

if (a) free(a);
if (x) free(x);
if (index) free(index);
if (threads) free(threads);
if (args) free(args);

return -7;
}

for (i = 0; i < total_threads; i++)
if (pthread_join(threads, 0)){
printf("Cannot wait thread %d!\n", i);

if (a) free(a);
if (x) free(x);
if (index) free(index);
if (threads) free(threads);
if (args) free(args);

return -8;
}

t_full = get_full_time() - t_full;

if (t_full == 0)
t_full = 1;
printf("\n\nInversion time = %ld\nTotal threads time = %ld"\
" (%ld%%)\nPer thread = %ld\n",
t_full, thread_time, thread_time * 100/t_full,
thread_time/total_threads);

InputMatrix(n, a, mode, fp);
printf("\n||A*A^{-1} - E|| = %e\n", TestMatrix(n, a, x));

free(a);
free(x);
free(index);
free(threads);
free(args);

return 0;
}
 

Вложения