1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
| #include <fstream>
#include <iostream>
#include <mpi.h>
#include <cmath>
#include <malloc.h>
#include <memory.h>
#include <iomanip>
/* myid – номер процесса, numprocs – количество процессов, Root – номер
корневого процесса, sendcounts, displs – массивы для организации обменов */
int myid, numprocs, Root = 0, *sendcounts, *displs;
/* AB – исходная матрица метода, может быть задана любым способом
X – вектор, содержит начальное приближение */
double *AB, *A, *X;
/* задана матрица А, начальное приближение вектора Х_old, размерность
матрицы MATR_SIZE, количество вычисляемых элементов вектора
в данном процессе size, вычисляем новые значение вектора Х с номера first */
void Iter_Jacoby(double *X_old, int size, int MATR_SIZE, int first);
void SolveSLAE(size_t MATR_SIZE, int size, double Error);
int ind(int i, int j, int SIZE); /* процедура обращения к элементу матрицы */
int main(int argc, char *argv[]) {
/* MATR_SIZE – размерность системы, size – количество строк матрицы в
каждом процессе, SIZE – количество элементов матрицы в каждом процессе */
int i, size, MATR_SIZE, SIZE;
double Error; /* точность вычислений */
double startTime, finishTime;
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
MPI_Comm_rank(MPI_COMM_WORLD, &myid);
if (myid == Root) {
/* получаем размерность системы MATR_SIZE */
/* выделяем память для матрицы AB */
std::ifstream fin("Source.txt");
fin >> MATR_SIZE;
fin >> Error;
//AB = new double[1, 0, 0, 2, 0, 1, 0, 2, 0, 0, 1, 2];
AB = (double *)malloc(sizeof(double) * MATR_SIZE * (MATR_SIZE + 1));
for (int i = 0; i < MATR_SIZE * (MATR_SIZE + 1); i++) {
fin >> AB[i];
}
/* получение матрицы АB размерности MATR_SIZE+1(матрица+столбец
свободных членов), задание точности вычислений Error */
}
startTime = MPI_Wtime();
/* рассылка значений всем процессам */
std::cout << &MATR_SIZE << "\n";
MPI_Bcast(&MATR_SIZE, 1, MPI_INT, Root, MPI_COMM_WORLD);
MPI_Bcast(&Error, 1, MPI_DOUBLE, Root, MPI_COMM_WORLD);
/* выделяем память для вектора X */
X = (double *)malloc(sizeof(double) * MATR_SIZE);
if (myid == Root) {
for (int i = 0; i < MATR_SIZE; i++) {
X[i] = 1;
}
/*получаем начальное значение для вектора Х */
}
std::cout << MATR_SIZE << " ; " << X;
/* рассылка вектора Х всем процессам */
MPI_Bcast(X, MATR_SIZE, MPI_DOUBLE, Root, MPI_COMM_WORLD);////////////////////////////////////////////////////////////////////////////////////////
/* определение количества элементов вектора, которые будут
вычисляться в каждом процессе (равно количеству строк
матрицы, находящихся в данном процессе) */
size = (MATR_SIZE / numprocs) + ((MATR_SIZE % numprocs) > myid ? 1 : 0);
/* выделяем память для матрицы А в каждом процессе*/
A = (double *)malloc(sizeof(double) * (MATR_SIZE + 1)*size);
displs = (int *)malloc(numprocs * sizeof(int));
sendcounts = (int *)malloc(numprocs * sizeof(int));
/* рассылка частей матрицы по процессам */
SIZE = (MATR_SIZE + 1) * size;
MPI_Gather(&SIZE, 1, MPI_INT, sendcounts, 1, MPI_INT, Root, MPI_COMM_WORLD);
displs[0] = 0;
for (i = 1; i < numprocs; ++i)
displs[i] = displs[i - 1] + sendcounts[i - 1];
MPI_Scatterv(AB, sendcounts, displs, MPI_DOUBLE, A,
(MATR_SIZE + 1) * size, MPI_DOUBLE, Root, MPI_COMM_WORLD);
/* решение СЛАУ методом простой итерации */
SolveSLAE(MATR_SIZE, size, Error);
/* освобождение памяти */
finishTime = MPI_Wtime();
if (myid == Root) {
std::cout << "algorithm took " << finishTime - startTime << " seconds with eps = " << Error << " on " << numprocs << " processes" << std::endl;
std::ofstream fout("output.txt");
fout << std::fixed << std::setprecision(4);
for (int i = 0; i < MATR_SIZE; i++) {
fout << X[i] << std::endl;
}
}
free(sendcounts); free(displs);
free(AB);
// free(AB);
free(A); free(X);
MPI_Finalize();
//system("pause");
return 0;
}
void SolveSLAE(size_t MATR_SIZE, int size, double Error) {
/* задана матрица А, размерность матрицы MATR_SIZE, количество
элементов, вычисляемых в данном процессе size, погрешность
вычислений Error*/
double *X_old;
int Iter = 0, i, Result, first;
double dNorm = 0, dVal;
/* определение номера элемента first вектора Х, с которого
будут вычисляться новые значения в данном процессе*/
MPI_Scan(&size, &first, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
first -= size;
/* заполнение массива sendcounts значениями size от каждого процесса*/
MPI_Allgather(&size, 1, MPI_INT, sendcounts, 1, MPI_INT, MPI_COMM_WORLD);
/* заполнение массива displs – расстояний между элементами вектора,
распределенных по процессам */
displs[0] = 0;
for (i = 1; i < numprocs; ++i)
displs[i] = displs[i - 1] + sendcounts[i - 1];
/* выделяем память для X_old */
X_old = (double *)malloc(sizeof(double) * MATR_SIZE);
do {
++Iter;
/*сохраняем предыдущее значение вектора Х */
memcpy(X_old, X, sizeof(double)*MATR_SIZE);
/* вычисляем новое значение вектора Х */
//double *R = new double (MATR_SIZE);
Iter_Jacoby(X_old, size, MATR_SIZE, first);
/* Рассылаем вектор Х всем процессам */
MPI_Allgatherv(&X[first], size, MPI_DOUBLE, X, sendcounts, displs,
MPI_DOUBLE, MPI_COMM_WORLD);
/* Считаем норму */
if (myid == Root) {
dNorm = 0;
for (i = 0; i < MATR_SIZE; ++i) {
dVal = fabs(X[i] - X_old[i]);
if (dNorm < dVal) dNorm = dVal;
}
Result = Error < dNorm;
//std::cout << X[0] << std::endl;
}
/* рассылаем результат всем процессам */
MPI_Bcast(&Result, 1, MPI_INT, Root, MPI_COMM_WORLD);
} while (Result);
free(X_old);
return;
}
void Iter_Jacoby(double *X_old, int size, int MATR_SIZE, int first) {
/* задана матрица А, начальное приближение вектора Х_old, размерность
матрицы MATR_SIZE, количество вычисляемых элементов вектора
в данном процессе size, вычисляем новые значение вектора Х с номера first */
int i, j;
double Sum;
for (i = 0; i < size; ++i)
{
Sum = 0;
for (j = 0; j < i + first; ++j)
Sum += A[ind(i, j, MATR_SIZE)] * X_old[j];
for (j = i + 1 + first; j < MATR_SIZE; ++j)
Sum += A[ind(i, j, MATR_SIZE)] * X_old[j];
X[i + first] = (A[ind(i, MATR_SIZE, MATR_SIZE)] - Sum) /
A[ind(i, i + first, MATR_SIZE)];
}
}
int ind(int i, int j, int SIZE) {/* процедура обращения к элементу матрицы */
return (i*(SIZE + 1) + j);
} |
(
