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

Тема в разделе "C/C++/C#", создана пользователем LudmilaUt, 16 янв 2012.

  1. LudmilaUt

    LudmilaUt New Member

    Регистрация:
    16 янв 2012
    Сообщения:
    1
    Симпатии:
    0
    Задача: обратить матрицу методом Жордана с выбором максимального элемента по всей матрице. (Метод Жордана - это приведение матрицы к диагональному виду; сначала выбирается максимальный элемент и ставится на первое место первой строки, потом удаляются элементы первого столбца со второго до конца. Потом то же самое для подматрицы, и т.д. потом то же самое для верхнего правого угла матрицы). Использовать параллельное программирование. Выбор главного элемента, обращение матрицы, функция 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;
    }
     

    Вложения:

Загрузка...
Похожие Темы - Обратная Матрица
  1. egor811
    Ответов:
    1
    Просмотров:
    1.787
  2. Normann
    Ответов:
    3
    Просмотров:
    3.548
  3. fatpunk
    Ответов:
    0
    Просмотров:
    1.137
  4. нини
    Ответов:
    3
    Просмотров:
    1.745
  5. phobiaxx
    Ответов:
    1
    Просмотров:
    1.278

Поделиться этой страницей