Поиск:

Ответ в темуСоздание новой темы Создание опроса
> Распараллеливание алгоритма Гаусса, нахождения обратной матрицы 
V
    Опции темы
k0c9k
Дата 24.12.2006, 23:18 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 36
Регистрация: 24.12.2006

Репутация: нет
Всего: нет



Товарищи ! нужно реализовать параллельный алгоритм нахождения обратной амтрицы методом Гаусса. Не поделитесь соображениями ? Просто не могу сообразить как его вообще можно распараллелить… Но можно: надыбал какой-то исходник решения СЛАУ методом Гаусса (сам разобраться в нем как-то немогу =( )
Код

void main(int argc, char **argv)
{
    int  id, N, i, j, s, k, l, index, m, n, p, locN, t;
    double starttime, endtime, max, c;
    double **a, *b, *buf;

    MPI_Init (&argc, &argv);
    starttime = MPI_Wtime ();

    //n = atoi (argv[1]);
    n = 16;     // matrix size

    MPI_Comm_size (MPI_COMM_WORLD,  &N);
    MPI_Comm_rank (MPI_COMM_WORLD,  &id);

    if (n % N)
    {
        MPI_Finalize ();
        printf ("Matrix size (%i) can't be divided to the number of processors (%i)\n",  n,  N);
        exit (1);
    }

    locN = n/N;
    if (id == 0)
        printf ("locN = %i\n\n", locN);

    a  =  (double**) malloc (locN * sizeof (double));       // create matrix a (for each process)
    if (id == 0)
        b = (double*) malloc (n * sizeof (double));         // create vector b (for bm only)

    buf = (double*) malloc (n * sizeof (double));

    // fill matrix a[locN x n]
    for (i=0; i<locN; i++)
    { 
        a[i] = (double*) malloc (n * sizeof (double));
        for (j=0; j<n; j++)
        {
            if (j == i + id * locN)
                a[i][j] = j+1;
            else if (j == i + id * locN + 1)
                a[i][j] = 2;
            else
                a[i][j] = 0;
        }
    }

    // fill vector b
    if (id == 0)
        for (i = 0;i<n;i++)
            b[i] = (i+1) * 2;

    if (id == 0)
    {
        printf("Vector:\n");

        for (i=0; i<n; i++)
            printf("% 7.1f", b[i]);

        printf ("\n\n");
    }

    if (id == 0)
    {
        printf ("Matrix:\n");
        for (i=0; i<locN; i++)
        {
            for (j=0; j<n; j++)
            {
                printf ("% 7.1f", a[i][j]);
            }
            printf ("\n");
        }
        printf ("\n");
    }
   
    for (p=0; p<N; p++)
    {
        for (i=0; i<locN; i++)
        {
            // find max element in matrix a for the current process
            if (id == p)
            {
                max = a[i][i + p*locN];
                index = i + p*locN;

                for (j=i + p*locN + 1; j<n; j++)
                {
                    if (a[i][j] > max)
                    {
                        max = a[i][j];
                        index = j;
                    }
                }
            }

            // broadcast to other processes
            MPI_Bcast (&index, 1, MPI_INT, p, MPI_COMM_WORLD);
            MPI_Bcast (&max, 1, MPI_DOUBLE, p, MPI_COMM_WORLD);
           
            /*
            printf("id: %d   max: %f  index: %d\n", id, max, index);   
            printf("Matrix before subst ... \n");

            for (k = 0;k<locN;k++)
            {
            for (s = 0;s<n;s++)
            printf("id: %d  a: %.2f ", id, a[k][s]);
            printf("\n");
            }
            */
           
            if (index != i + p*locN)
            {   
                for (l=0; l<locN; l++)
                {
                    c = a[l][i + p*locN];
                    a[l][i + p*locN] = a[l][index];
                    a[l][index] = c;
                }

                c = b[i + p*locN];
                b[i + p*locN] = b[index];
                b[index] = c;
            }


            /*    printf("Matrix after subst ... \n");
            for (k = 0;k<locN;k++)
            {
            for (s = 0;s<n;s++)
            printf("id: %d  a: %.2f ", id, a[k][s]);
            printf("\n");
            }
            */       
            if (id == p)
            {
                for (k = 0;k<n;k++)
                {
                    buf[k] = a[i][k];
                    /* printf("id: %d a[%d][%d]: %.2f\n", id, i, k, a[i][k]);*/
                }
            }

            MPI_Bcast(buf, n, MPI_DOUBLE, p, MPI_COMM_WORLD);
            /*    for (s = 0;s<n;s++)
            printf("id: %d  buf: %f\n", id, buf[s]);

            printf("New matrix before ... \n");
            for (s = 0;s<locN;s++)
            {
            for (k = 0;k<n;k++)
            printf("id: %d  a: %.2f ", id, a[s][k]);
            printf("\n");
            }
            */       
            for (s = 0;s<locN;s++)
                for (t = 0;t<n;t++)
                    if (t!= i+p*locN)
                        a[s][t] -= a[s][i+p*locN]*(buf[t]/buf[i+p*locN]);     

            if (id == 0)
                for (k = 0;k<n;k++)
                    if (k!= i+p*locN)
                        b[k] -= b[i+p*locN]*(buf[k]/buf[i+p*locN]);

            MPI_Barrier(MPI_COMM_WORLD);
            /*    printf("New matrix after next step... \n", i+1);
            for (k = 0;k<locN;k++)
            {
            for (j = 0;j<n;j++)
            printf("id: %d  a: %.2f ", id, a[k][j]);
            printf("\n");
            }
            */
        }         
    }           

    MPI_Barrier(MPI_COMM_WORLD);
    /* if (id == 0) printf("END MATRIX ... \n");
    for (i = 0;i<locN;i++)
    {
    for (j = 0;j<n;j++)
    printf("id: %d  a: %.2f ", id, a[i][j]);
    printf("\n");
    }
    if (id == 0)
    {     
    printf("END VECTOR ... \n");
    for (k = 0;k<n;k++)
    printf("%f\n", b[k]);
    }
    */
    MPI_Barrier(MPI_COMM_WORLD);
    MPI_Bcast(b, n, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    /*for (i = 0;i<n;i++)
    printf("id: %d   b: %f\n", id, b[i]);
    if (id == 0)
    printf("So,  we find the solution of SLAE ... \n");
    for (p = 0;p<locN;p++)
    printf("id: %d res: %f\n", id, b[p+id*locN]/a[p][p+id*locN]);
    */   
    MPI_Barrier(MPI_COMM_WORLD);
    endtime = MPI_Wtime()-starttime;     
    if (id == 0) 
        printf("Time %f\n", endtime);

    for (i = 0;i<locN;i++)
        free(a[i]);
    free(a);
    free(b);
    free(buf);
    MPI_Finalize();
    exit(1);



Помогите плиз ! 

PM MAIL   Вверх
SoWa
Дата 25.12.2006, 05:03 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Харекришна
****


Профиль
Группа: Комодератор
Сообщений: 2422
Регистрация: 18.10.2004

Репутация: 6
Всего: 74



И ты хочешь, чтобы твой код разобрали? Нет. Лучше сам напиши. Не сложно ведь?

Что касается алгоритма. Что значит параллельный?


--------------------
Всем добра smile
PM MAIL ICQ   Вверх
k0c9k
Дата 25.12.2006, 10:00 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 36
Регистрация: 24.12.2006

Репутация: нет
Всего: нет



Не сложно. Почти написал. Такая шиза получилась. Кому интересно - могу выложить. 
Параллельный в смысле, что в несколько процессов выполняется. Реализован с помощью MPI
Как доделаю расскажу свои домыслы.

Это сообщение отредактировал(а) k0c9k - 25.12.2006, 10:04
PM MAIL   Вверх
V.A.KeRneL
  Дата 25.12.2006, 10:46 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Vadim A. Kazantsev
**


Профиль
Группа: Участник
Сообщений: 291
Регистрация: 3.12.2006
Где: Moscow, Russia

Репутация: 1
Всего: 14



Цитата(k0c9k @  25.12.2006, 10:00 Найти цитируемый пост)

Кому интересно - могу выложить. 

Вообще, имхо, надо по дефолту выкладывать... Наверняка кому-нибудь в будущем пригодится! И не надо будет заново всё с нуля писать, а можно будет просто ссылочку сюда дать. А те, кто умеют пользоваться поиском по форуму, вообще новую тему создавать не станут!.. Так что, выкладывай, конечно!


Это сообщение отредактировал(а) V_A_KeRneL - 25.12.2006, 10:47


--------------------
«C'est un pense-creux d'ici. C'est le meilleur et le plus irascible homme du monde...» © Ф.М. Достоевский, «Бесы»
---/)/)---(\.../)---(\(\
--(':'=)---(=';'=)---(=':')
(")(")..)-(").--.(")-(..(")(")

PM MAIL IM ICQ AOL YIM MSN   Вверх
SoWa
Дата 25.12.2006, 15:03 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Харекришна
****


Профиль
Группа: Комодератор
Сообщений: 2422
Регистрация: 18.10.2004

Репутация: 6
Всего: 74



Нет. Выкладывать не стоит. Кому надо- сами напишут.
Или если уж совсем невтерпеж похвастаться- подредактируй свой первый пост и в нем выложи


--------------------
Всем добра smile
PM MAIL ICQ   Вверх
k0c9k
Дата 26.12.2006, 06:41 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 36
Регистрация: 24.12.2006

Репутация: нет
Всего: нет



Пока исходник не отлажен расскажу что я сам надумал.
Напомню, метод Гаусса нахождения обратной матрицы заключается в приписании к правой части единичной матрицы такой-же размерности. И  с помощью перестановок сведением левой матрицы к единичной. Таким образом в правой части получается обратная.
Сам процесс вычисления состоит в последовательном делении i-ой строки на элемент a[i,i] (имеется ввиду что на этом этапе все элементы a[i,k], k<i уже равны нулю). Элемент a[i,i] станет равным 1. после этого можно нужно будет из всех строк вычесть i-ю строку домноженную на a[k,i],k!=i. таким образом мы сводим левую (исходную) матрицу к единичной, а в правой части появится обратная.
Теперь как я всё это дело решил запараллелить:
В принципе, если у нас матрица больших размеров, самой трудоёмкой частью будет вычитание текущей строки из остальных строк.
Поэтому решил что у всех процессов будет храниться своя часть матрицы. После того как диспетчерский процесс считает саму матрицу он поледовательно будет выдавать всем процессам по строке. В итоге получаем что у каждого процесса будет в среднем (размерность/кол-во потоков). Теперь диспетчерский процесс посылает остальным сообщение о том, что i-строка делится на a[i,i]. Тот процесс который понял что i-строка принадлежит ему производит деление строки, после чего отсылает её остальным процессам. Процессы вычитают из своих строк полученную. Цикл повторяется. По окончании цикла, все процессы отсылают диспетчерскому строки правой части. 


Это сообщение отредактировал(а) k0c9k - 26.12.2006, 07:34
PM MAIL   Вверх
Silver
Дата 26.12.2006, 09:05 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 7
Регистрация: 1.9.2005

Репутация: 1
Всего: 1



PM MAIL   Вверх
k0c9k
Дата 26.12.2006, 22:02 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 36
Регистрация: 24.12.2006

Репутация: нет
Всего: нет



Спасибо за хелп который здесь появился  smile !
Программную  реализацию выкладывать не буду, ИМХО нечитабельные исходники выкладывать бесполезно. Так что если комуто ОЧЕНЬ понадобится: стучите в ПМ.
PM MAIL   Вверх
DENNN
Дата 28.12.2006, 12:49 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Эксперт
****


Профиль
Группа: Участник Клуба
Сообщений: 3878
Регистрация: 27.3.2002
Где: Москва

Репутация: 1
Всего: 43



Цитата(k0c9k @  26.12.2006,  06:41 Найти цитируемый пост)
после чего отсылает её остальным процессам. Процессы вычитают из своих строк полученную. 

И в чем смысл? Параллельность - это когда несколько действий выполняются одновременно. А здесь одновременность только в вычитании, а потом все процессы опять ждут когда один из них поделит строку на число и разошлеть остальным.
Кроме того, у меня большие сомнения, что при продуманной организации данных  вычитание скажем N строк N процессами будет заметно бытсрее, чем последовательная обработка этих строк одним процессом. Вариант с многопроцессорным кластром не берем ...  smile 
PM ICQ   Вверх
maxim1000
Дата 28.12.2006, 12:55 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Эксперт
****


Профиль
Группа: Участник
Сообщений: 3334
Регистрация: 11.1.2003
Где: Киев

Репутация: 33
Всего: 110



честно говоря, я вижу здесь одну вещь, котогрую можно распараллелить - вычитание одной строки из другой (ну и деление тогда же), каждая компонента векторов обрабатывается независимо, так что можно без проблем поручить это дело разным процессам
но с другой стороны, ещё неизвестно какие накладные расходы вызовет такое мелкое деление...


--------------------
qqq
PM WWW   Вверх
k0c9k
Дата 28.12.2006, 13:28 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 36
Регистрация: 24.12.2006

Репутация: нет
Всего: нет



Да. у меня тоже была мысль ещё деление распределить, но из-за такой мелочи как-то не намного скорость я думаю увеличится.
Да вообще лучше для подсчета обратной матрицы в несколько процессов использовать рекурсивный алгоритм (с кватернарными деревьями)! Вот он параллелится замечательно! А гаусс - так, по приколу   smile 

Это сообщение отредактировал(а) k0c9k - 28.12.2006, 13:29
PM MAIL   Вверх
DENNN
Дата 28.12.2006, 13:51 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Эксперт
****


Профиль
Группа: Участник Клуба
Сообщений: 3878
Регистрация: 27.3.2002
Где: Москва

Репутация: 1
Всего: 43



Угу. Мое скромное ИМХО, паралелльность нужна для независимых процессов/потоков. Когда получится сформулировать алгоритм Гаусса как набор нескольки независимых процессов, то решение о реализации будет очевидным ))
PM ICQ   Вверх
  
Ответ в темуСоздание новой темы Создание опроса
Правила форума "Алгоритмы"

maxim1000

Форум "Алгоритмы" предназначен для обсуждения вопросов, связанных только с алгоритмами и структурами данных, без привязки к конкретному языку программирования и/или программному продукту.


Если Вам понравилась атмосфера форума, заходите к нам чаще! С уважением, maxim1000.

 
1 Пользователей читают эту тему (1 Гостей и 0 Скрытых Пользователей)
0 Пользователей:
« Предыдущая тема | Алгоритмы | Следующая тема »


 




[ Время генерации скрипта: 0.0863 ]   [ Использовано запросов: 20 ]   [ GZIP включён ]


Реклама на сайте     Информационное спонсорство

 
По вопросам размещения рекламы пишите на vladimir(sobaka)vingrad.ru
Отказ от ответственности     Powered by Invision Power Board(R) 1.3 © 2003  IPS, Inc.