Основы параллельных вычислений (1119586), страница 3

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 3)

можно умножать...// Умножение матрицы A на U(x)CD.m = a; // ставим указатель на матрицу Aint TotalCols = n - (k+1); // Сколько всего столбцов (зависит от шага!)int Rest = TotalCols % ThreadCount; //Остатокif (!Rest)// делится?FullThreads = ThreadCount; // да, еще как!else// нет, не делится...FullThreads = ThreadCount-1;// вычисляем количество строк на один ’полный’ потокint ColsPerThread = 0;if (FullThreads) ColsPerThread = TotalCols/FullThreads;for (int i = 0; i < FullThreads; i++) // run "full" threads{// заполняем структуруT[i].m=a;T[i].x=x;T[i].k=k;T[i].n=n;T[i].cs = k+1 + i*ColsPerThread;T[i].ce = k+1 + (i+1)*ColsPerThread;pthread_create(&(T[i].id), 0, Thread, &T[i]);}if (Rest) // если есть остаточные столбцы, доделать...{7T[FullThreads].m=a;T[FullThreads].x=x;T[FullThreads].k=k;T[FullThreads].n=n;T[FullThreads].cs = k+1 + FullThreads*ColsPerThread;T[FullThreads].ce = n; // до самого концаpthread_create(&(T[FullThreads].id), 0, Thread, &T[FullThreads]);}// А теперь подождем немного...for (int i = 0; i < ThreadCount; i++) pthread_join(T[i].id, 0);A(k,k) = Na;// Умножение обратной матрицы на U(x)CD.m = ai; // ставим указатель на матрицу AI (обратную к A),// далее все точно так же...for (int i = 0; i < iFullThreads; i++) // run "full" threads{T[i].m=a;T[i].x=x;T[i].k=k;T[i].n=n;T[i].cs = i*iColsPerThread;T[i].ce = (i+1)*iColsPerThread;pthread_create(&(T[i].id), 0, Thread, &T[i]);}if (iRest > 0) // the rest...{T[FullThreads].m = a;T[FullThreads].x = x;T[FullThreads].k = k;T[FullThreads].n = n;T[iFullThreads].cs = iFullThreads*iColsPerThread;T[iFullThreads].ce = n;pthread_create(&(T[iFullThreads].id), 0, Thread, &T[iFullThreads]);}for (int i = 0; i < ThreadCount; i++) pthread_join(T[i].id, 0);}// Ну, все, готово: матрица верхнетреугольная.delete x;delete T;// Обратный ходfor (int k = 0; k < n; k++) GaussReverse(a, ai, k, n);return true;}// Обратный ход метода Гауссаvoid GaussReverse(double * a, double * ai, int col, int n){for (int k = 1; k <= n; k++){double koeff = AI(n-k, col);for (int i = n-k+1; i < n; i++) koeff -= AI(i, col)*A(n-k,i);AI(n-k,col) = koeff / A(n-k,n-k);}}4.2.

Метод ЖорданаВ предыдущем пункте было замечено, что реализация массива с данными для потоков не очень эффективна.Здесь мы рассмотрим реализацию еще одного метода, но на этот раз поступим более грамотно — с использованием общей для всех потоков структуры данных. Так как метод Жордана крайне прост и компактен, приведём8полный текст программы.#include#include#include#include<stdio.h><math.h><sys/time.h><pthread.h>#define A(i,j) a[(i)*n+(j)]// matrires &#define AS(i,#define XS(i)#define BS(i)vectors with swapped coordsj) a[(i) * n + index[j]]x[index[i]]b[index[i]]void PrintMatrix(double *a, double *b, int n){for (int i = 0; i < n; i++){for (int j = 0; j < n; j++) printf("%1.4lf ", A(i, j));printf(" | %1.4lf\n", b[i]);}}class CCommonData // общие данные{public:double * a; // matrix pointerdouble * b; //int * index;int n; // total size of matrixint i; // current step};class CThreadData // данные для каждого потока{public:CCommonData * CD;pthread_t id; // Thread identifierint rs; // start rowint re; // end row};void * Thread(void * Ptr) // subtracts max row from rows [rs..re) of matrix A{CThreadData * TD = (CThreadData *) Ptr;int rs = TD->rs; // startint re = TD->re; // enddouble * a = TD->CD->a; // matrixdouble * b = TD->CD->b; // vector (right part)int * index = TD->CD->index; // variable swapint i = TD->CD->i; // current stepint n = TD->CD->n; // sizefor (int j = rs; j < re; j++){double p = AS(j, i);for (int k = i; k < n; k++) AS(j, k) -= p * AS(i, k);b[j] -= p * b[i];}return 0;}bool SolveSystem(int n, double *a, double *b, double *x, int *index, int ThreadCount){CCommonData CD;9// заполняем общие данныеCD.a=a;CD.b=b;CD.index=index;CD.n=n;CThreadData * T = new CThreadData[ThreadCount];int i;for (int t = 0; t < ThreadCount; t++) T[t].CD = &CD;for (i = 0; i < n; i++) index[i] = i; // перестановка (исходно тождественная)for (i = 0; i < n; i++){int k, j; // min indexfor (k = i, j = i; j < n; j++)if (fabs(AS(i, k)) < fabs(AS(i, j))) k = j;j = index[k];index[k] = index[i];index[i] = j;double p = AS(i, i);if (fabs(p) < 1e-100) return false;p = 1.0/p;for (k = i; k < n; k++) AS(i, k) *= p;b[i] *= p;CD.i = i; // set current stepint Rows = n - (i+1); // количество строкint Rest = Rows % ThreadCount;int FullThreads = ThreadCount;if (Rest) FullThreads = ThreadCount-1;int RowsPerThread = 0;if (FullThreads > 0) RowsPerThread = Rows/FullThreads;for (int t = 0; t < FullThreads; t++) // run "full" threads{T[t].rs = i+1 + t*RowsPerThread;T[t].re = i+1 + (t+1)*RowsPerThread;pthread_create(&(T[t].id), 0, Thread, &T[t]);}if (Rest) // the rest...{T[FullThreads].rs = i+1 + FullThreads*RowsPerThread;T[FullThreads].re = n;pthread_create(&(T[FullThreads].id), 0, Thread, &T[FullThreads]);}for (int t = 0; t < ThreadCount; t++) pthread_join(T[t].id, 0);}// second stepfor (i = n - 1; i > 0; i--){for (int j = i - 1; j >= 0; j--){double p = AS(j, i);AS(j, i) = 0; b[j] -= p * b[i];}}// обратная перестановка10for (i = 0; i < n; i++) XS(i) = b[i];delete T;return true;}double f(int i, int j){if (i==j) return 1.0+1.0 / double(i+j+1); elsereturn 1.0/double(i+j+1.0);}void FillMatrix(double * a, double * b, int n){for (int i = 0; i < n; i++){double s = 0;for (int j = 0; j < n; j++) { a[i*n+j] = f(i,j); s += f(i,j); }b[i] = s;}}/* ||A*x - b|| */double GetError(double *a, double *b, double *x, int n){double y = 0;for (int i = 0; i < n; i++){double z = -b[i];for (int j = 0; j < n; j++) z += A(i, j) * x[j];y += z * z;}if (y > 0) y = sqrt(y);return y;}void PrintSolution(double *x, int n){for (int i = 0; i < n; i++) printf("%1.4lf ", x[i]);}int main(){int n, TC;printf("Input dimension (n): ");scanf("%d", &n);if (n <= 0) return -1;printf("Input number of the threads: ");scanf("%d", &TC);if (n <= 0) return -2;double * a = new double[n*n];double * b = new double[n];double * x = new double[n];int * index = new int[n];double * ac = new double[n*n];double * bc = new double[n*n];FillMatrix(a, b, n);for (int i = 0; i < n*n; i++) ac[i] = a[i]; // copy of Afor (int i = 0; i < n; i++) bc[i] = b[i]; // copy of BPrintMatrix(a, b, n);timeval tv1;gettimeofday(&tv1, 0);if (!SolveSystem(n, a, b, x, index, TC))11{printf("Bad matrix!\n");return -3;}timeval tv2;gettimeofday(&tv2,0);double Time = double(tv2.tv_sec)*1000000.0+double(tv2.tv_usec)double(tv1.tv_sec)*1000000.0+double(tv1.tv_usec);printf("Result:\n"); PrintSolution(x, n);printf("\n\nThreads: %d,\nError: %1.17lf,\nTime=%1.4lf sec.\n", TC, GetError(ac, bc, x, n), Time/1000000.0);delete a; delete b;delete x; delete index;delete ac;return 0;}4.3.

Метод ГауссаНу, и напоследок, еще один метод — самый простой. . . Эх, пусть простит меня читатель за обилие иностранных комментариев в тексте. . .#include#include#include#include<stdio.h><math.h><sys/time.h><pthread.h>#define A(i,j) a[(i)*n+(j)]// matrires &#define AS(i,#define XS(i)#define BS(i)vectors with swapped coordsj) a[(i) * n + index[j]]x[index[i]]b[index[i]]// MATRIX FUNCTION: A=(Hilbert(n)+E)double f(int i, int j){if (i==j) return 1.0+1.0 / double(i+j+1); elsereturn 1.0/double(i+j+1.0);}void PrintMatrix(double *a, double *b, int n){for (int i = 0; i < n; i++){for (int j = 0; j < n; j++) printf("%1.5lf ", A(i, j));printf(" | %1.5lf\n", b[i]);}}// print solution Xvoid PrintX(double*x,int n) { for (int i=0;i<n;i++) printf("%1.4lf ",x[i]); }class CCommonData{public:double * a; // matrix pointerdouble * b; //int * index;int n; // total size of matrixint i; // current step};class CThreadData // Data for each Thread{public:12CCommonData * CD;pthread_t id; // Thread identifierint rs; // start rowint re; // end row};// thread functionvoid * Thread(void * Ptr) // subtracts one row from rows [rs..re) of matrix A{CThreadData * TD = (CThreadData *) Ptr;int rs = TD->rs; // startint re = TD->re; // enddouble * a = TD->CD->a; // matrixdouble * b = TD->CD->b; // vector (right part)int * index = TD->CD->index; // variable swapint i = TD->CD->i; // current stepint n = TD->CD->n; // sizefor (int j = rs; j < re; j++){double p = AS(j, i);for (int k = i; k < n; k++) AS(j, k) -= p * AS(i, k);b[j] -= p * b[i];}return 0;}// main solving functionbool SolveSystem(int n, double *a, double *b, double *x, int *index, int ThreadCount){CCommonData CD; // common thread data: pointers to matrix A and vector B,// matrix size and current stepCD.a=a;CD.b=b;CD.index=index; // rows substtution arrayCD.n=n;CThreadData * T = new CThreadData[ThreadCount]; // array for thread dataint i;for (int t = 0; t < ThreadCount; t++) T[t].CD = &CD; // fill pointers// to common datafor (i = 0; i < n; i++) index[i] = i; // init substitution: s=idfor (i = 0; i < n; i++) // steps (’i’ is inmber of step){int k, j; // min indexfor (k = i, j = i; j < n; j++) // search for MAX elementif (fabs(AS(i, k)) < fabs(AS(i, j))) k = j;j = index[k]; // exchange linesindex[k] = index[i];index[i] = j;double p = AS(i, i);if (fabs(p) < 1e-300) return false; // det A ~= 0 => very bad matrix :(p = 1.0/p;for (k = i; k < n; k++) AS(i, k) *= p;b[i] *= p;CD.i = i; // set current step in common structureint Rows = n - (i+1); // total rowsint Rest = Rows % ThreadCount; // rest number of rows13int FullThreads = ThreadCount;if (Rest) FullThreads = ThreadCount-1;int RowsPerThread = 0;if (FullThreads > 0) RowsPerThread = Rows/FullThreads;// run threadsfor (int t = 0; t < FullThreads; t++) // run "full" threads{T[t].rs = i+1 + t*RowsPerThread;T[t].re = i+1 + (t+1)*RowsPerThread;pthread_create(&(T[t].id), 0, Thread, &T[t]);}if (Rest) // the rest...{T[FullThreads].rs = i+1 + FullThreads*RowsPerThread;T[FullThreads].re = n;pthread_create(&(T[FullThreads].id), 0, Thread, &T[FullThreads]);}// wait while threads finish their workfor (int t = 0; t < ThreadCount; t++) pthread_join(T[t].id, 0);}// Gauss’s "reverse step"for (i = n - 1; i >= 0; i--){double p = BS(i);for (int j = i + 1; j < n; j++) p -= x[j] * AS(i, j);x[i] = p/AS(i, i);}delete T;return true;}// fill matrix and vector b using f().// here we make solution x=(1,1...1)void FillMatrix(double * a, double * b, int n){for (int i = 0; i < n; i++){double s = 0;for (int j = 0; j < n; j++) { a[i*n+j] = f(i,j); s += f(i,j); }b[i] = s;}}/* return |A*x - b| */double GetError(double *a, double *b, double *x, int n){double y = 0;for (int i = 0; i < n; i++){double z = -b[i];for (int j = 0; j < n; j++) z += A(i, j) * x[j];y += z * z;}if (y > 0) y = sqrt(y);return y;}////////////////////////////////////////////////////////////////int main()14{int n, TC;printf("Input dimension (n): ");scanf("%d", &n);if (n <= 0) return -1;printf("Input number of the threads: ");scanf("%d", &TC);if (n <= 0) return -2;double * a = new double[n*n];double * b = new double[n];double * x = new double[n];int * index = new int[n];double * ac = new double[n*n];double * bc = new double[n*n];FillMatrix(a, b, n);for (int i = 0; i < n*n; i++) ac[i] = a[i]; // copy of Afor (int i = 0; i < n; i++) bc[i] = b[i]; // copy of BPrintMatrix(a, b, n);timeval tv1;gettimeofday(&tv1, 0); // get current timeif (!SolveSystem(n, a, b, x, index, TC)){printf("Bad matrix!\n");return -3;}timeval tv2;gettimeofday(&tv2,0); // again get current time// Get time in microseconds, then divide by 1.000.000double Time = (double(tv2.tv_sec)*1000000.0+double(tv2.tv_usec)double(tv1.tv_sec)*1000000.0+double(tv1.tv_usec))/1000000.0;printf("Result:\n"); PrintX(x, n);printf("\n\nThreads: %d,\nError: %1.17lf,\nTime=%1.4lf sec.\n", TC, GetError(ac, bc, x, n), Time);delete a; delete b;delete x; delete index;delete ac;return 0;}5.

Заключение5.1. Ещё более эффективные реализации многопоточных программВсе примеры, приведённые выше, страдают одним недостатком. В них потоки рождаются и умирают накаждом шаге алгоритма. Вообще говоря, это плохо. Дело в том, что создание потока требует некоторых усилий от операционной системы, и на это уходит драгоценное время. Есть некоторый механизм, позволяющийизбежать этой утраты. Он связан с объектами синхронизации типа mutex. Эти объекты позволяют управлятьвыполнением потоков и представляю собой что-то вроде глобальных «выключателей», которые доступны изразных потоков. Можно в самом начале работы алгоритма создать нужное число потоков, а на каждом шагесинхронизировать выполнение очередной порции вычислений с помощью mutex’ов.

Характеристики

Список файлов учебной работы