GCC Code Coverage Report

Line	Branch	Exec	Source
1			#include "ilin_a_gaussian_method_horizontal_band_scheme/mpi/include/ops_mpi.hpp"
2
3			#include <mpi.h>
4
5			#include <algorithm>
6			#include <cmath>
7			#include <cstddef>
8			#include <vector>
9
10			#include "ilin_a_gaussian_method_horizontal_band_scheme/common/include/common.hpp"
11
12			namespace ilin_a_gaussian_method_horizontal_band_scheme {
13
14		✗	IlinAGaussianMethodMPI::IlinAGaussianMethodMPI(const InType &in) {
15			SetTypeOfTask(GetStaticTypeOfTask());
16		✗	GetInput() = in;
17		✗	GetOutput() = std::vector<double>();
18		✗	}
19
20		✗	bool IlinAGaussianMethodMPI::ValidationImpl() {
21		✗	MPI_Comm_rank(MPI_COMM_WORLD, &rank_);
22
23		✗	if (rank_ != 0) {
24			return true;
25			}
26
27			const auto &input = GetInput();
28		✗	if (input.size() < 4) {
29			return false;
30			}
31
32		✗	int size = static_cast<int>(input[0]);
33		✗	int band_width = static_cast<int>(input[1]);
34
35		✗	if (size <= 0 || band_width <= 0 || band_width > size) {
36			return false;
37			}
38
39		✗	size_t expected_count = 2 + (size * band_width) + size;
40		✗	return input.size() == expected_count;
41			}
42
43		✗	bool IlinAGaussianMethodMPI::PreProcessingImpl() {
44			InitializeMPI();
45		✗	BroadcastInputData();
46		✗	ScatterLocalData();
47		✗	return true;
48			}
49
50		✗	void IlinAGaussianMethodMPI::InitializeMPI() {
51		✗	MPI_Comm_rank(MPI_COMM_WORLD, &rank_);
52		✗	MPI_Comm_size(MPI_COMM_WORLD, &size_);
53		✗	}
54
55		✗	void IlinAGaussianMethodMPI::BroadcastInputData() {
56		✗	if (rank_ == 0) {
57			const auto &input = GetInput();
58		✗	data_.size = static_cast<int>(input[0]);
59		✗	data_.band_width = static_cast<int>(input[1]);
60			}
61
62		✗	MPI_Bcast(&data_.size, 1, MPI_INT, 0, MPI_COMM_WORLD);
63		✗	MPI_Bcast(&data_.band_width, 1, MPI_INT, 0, MPI_COMM_WORLD);
64
65		✗	n_ = data_.size;
66		✗	band_ = data_.band_width;
67
68		✗	rows_per_proc_ = n_ / size_;
69		✗	remainder_ = n_ % size_;
70
71		✗	if (rank_ < remainder_) {
72		✗	local_rows_ = rows_per_proc_ + 1;
73		✗	row_start_ = rank_ * local_rows_;
74			} else {
75		✗	local_rows_ = rows_per_proc_;
76		✗	row_start_ = (remainder_ * (rows_per_proc_ + 1)) + ((rank_ - remainder_) * rows_per_proc_);
77			}
78		✗	row_end_ = row_start_ + local_rows_;
79
80		✗	if (rank_ == 0) {
81			const auto &input = GetInput();
82		✗	int mat_size = n_ * band_;
83		✗	data_.matrix.resize(static_cast<size_t>(mat_size));
84		✗	data_.vector.resize(static_cast<size_t>(n_));
85
86		✗	std::copy(input.begin() + 2, input.begin() + 2 + mat_size, data_.matrix.begin());
87		✗	std::copy(input.begin() + 2 + mat_size, input.end(), data_.vector.begin());
88			} else {
89		✗	data_.matrix.resize(static_cast<size_t>(n_) * static_cast<size_t>(band_));
90		✗	data_.vector.resize(static_cast<size_t>(n_));
91			}
92
93		✗	MPI_Bcast(data_.matrix.data(), n_ * band_, MPI_DOUBLE, 0, MPI_COMM_WORLD);
94		✗	MPI_Bcast(data_.vector.data(), n_, MPI_DOUBLE, 0, MPI_COMM_WORLD);
95		✗	}
96
97		✗	void IlinAGaussianMethodMPI::ScatterLocalData() {
98		✗	local_matrix_.resize(static_cast<size_t>(local_rows_) * static_cast<size_t>(band_));
99		✗	local_vector_.resize(static_cast<size_t>(local_rows_));
100
101		✗	for (int i = 0; i < local_rows_; ++i) {
102		✗	int global_row = row_start_ + i;
103		✗	const auto global_row_large = static_cast<std::ptrdiff_t>(global_row);
104		✗	std::copy(data_.matrix.begin() + (global_row_large * band_),
105		✗	data_.matrix.begin() + ((global_row_large + 1) * band_),
106		✗	local_matrix_.begin() + static_cast<std::ptrdiff_t>(i) * band_);
107		✗	local_vector_[static_cast<size_t>(i)] = data_.vector[static_cast<size_t>(global_row)];
108			}
109
110		✗	solution_.resize(static_cast<size_t>(n_), 0.0);
111		✗	pivot_row_buf_.resize(static_cast<size_t>(band_));
112		✗	}
113
114		✗	bool IlinAGaussianMethodMPI::RunImpl() {
115		✗	ProcessForwardElimination();
116		✗	ProcessBackwardSubstitution();
117		✗	GetOutput() = solution_;
118		✗	return true;
119			}
120
121		✗	void IlinAGaussianMethodMPI::ProcessForwardElimination() {
122			const double eps = 1e-12;
123		✗	pivot_row_buf_.assign(band_, 0.0);
124		✗	double pivot_b = 0.0;
125
126		✗	for (int k = 0; k < n_; ++k) {
127			double local_max = 0.0;
128			int local_max_row = -1;
129
130		✗	local_max = FindLocalPivotValue(k, local_max_row);
131
132			struct {
133			double max_val;
134			int max_row;
135		✗	} local_data{.max_val = local_max, .max_row = local_max_row}, global_data{.max_val = 0.0, .max_row = -1};
136
137		✗	MPI_Allreduce(&local_data, &global_data, 1, MPI_DOUBLE_INT, MPI_MAXLOC, MPI_COMM_WORLD);
138
139		✗	int global_max_row = global_data.max_row;
140		✗	double global_max_val = global_data.max_val;
141
142		✗	if (global_max_val < eps) {
143		✗	solution_[static_cast<size_t>(k)] = 0.0;
144		✗	continue;
145			}
146
147			int pivot_owner = CalculatePivotOwner(global_max_row);
148		✗	BroadcastPivotData(pivot_owner, pivot_row_buf_, pivot_b, global_max_row);
149
150		✗	if (global_max_row != k) {
151			int k_owner = CalculateRowOwner(k);
152
153			int k_local_idx = FindLocalRowIndex(k);
154		✗	if (k_owner == rank_ && k_local_idx >= 0) {
155		✗	HandleRowSwapLocal(k_local_idx, pivot_owner, global_max_row);
156		✗	} else if (pivot_owner == rank_) {
157		✗	HandleRowSwapRemote(k_owner, global_max_row);
158			}
159			}
160
161		✗	double pivot = pivot_row_buf_[static_cast<size_t>(band_ - 1)];
162		✗	if (std::fabs(pivot) < eps) {
163		✗	continue;
164			}
165
166		✗	for (int i = 0; i < local_rows_; ++i) {
167		✗	EliminateRow(i, k, pivot_row_buf_, pivot_b);
168			}
169			}
170		✗	}
171
172		✗	double IlinAGaussianMethodMPI::FindLocalPivotValue(int k, int &local_max_row) const {
173			double local_max = 0.0;
174		✗	local_max_row = -1;
175
176		✗	for (int i = 0; i < local_rows_; ++i) {
177		✗	int global_row = row_start_ + i;
178		✗	if (global_row < k) {
179		✗	continue;
180			}
181
182		✗	int diag_idx = band_ - 1 - (global_row - k);
183		✗	if (diag_idx < 0 || diag_idx >= band_) {
184		✗	continue;
185			}
186
187		✗	double val = std::fabs(local_matrix_[(static_cast<size_t>(i) * static_cast<size_t>(band_)) + diag_idx]);
188		✗	if (val > local_max) {
189			local_max = val;
190		✗	local_max_row = global_row;
191			}
192			}
193
194		✗	return local_max;
195			}
196
197		✗	void IlinAGaussianMethodMPI::BroadcastPivotData(int pivot_owner, std::vector<double> &pivot_row, double &pivot_b,
198			int global_max_row) const {
199		✗	if (pivot_owner == rank_) {
200			int local_pivot_idx = FindLocalRowIndex(global_max_row);
201		✗	if (local_pivot_idx >= 0) {
202		✗	const double *row_start = &local_matrix_[(static_cast<size_t>(local_pivot_idx) * static_cast<size_t>(band_))];
203		✗	std::ranges::copy(row_start, row_start + band_, pivot_row.begin());
204		✗	pivot_b = local_vector_[static_cast<size_t>(local_pivot_idx)];
205			}
206			}
207
208		✗	MPI_Bcast(pivot_row.data(), band_, MPI_DOUBLE, pivot_owner, MPI_COMM_WORLD);
209		✗	MPI_Bcast(&pivot_b, 1, MPI_DOUBLE, pivot_owner, MPI_COMM_WORLD);
210		✗	}
211
212		✗	int IlinAGaussianMethodMPI::CalculateRowOwner(int row) const {
213		✗	if (row < (remainder_ * (rows_per_proc_ + 1))) {
214		✗	return row / (rows_per_proc_ + 1);
215			}
216		✗	return remainder_ + ((row - remainder_ * (rows_per_proc_ + 1)) / rows_per_proc_);
217			}
218
219		✗	int IlinAGaussianMethodMPI::CalculatePivotOwner(int global_max_row) const {
220		✗	if (global_max_row < 0) {
221			return -1;
222			}
223			return CalculateRowOwner(global_max_row);
224			}
225
226		✗	int IlinAGaussianMethodMPI::FindLocalRowIndex(int global_row) const {
227		✗	for (int i = 0; i < local_rows_; ++i) {
228		✗	if (row_start_ + i == global_row) {
229			return i;
230			}
231			}
232			return -1;
233			}
234
235		✗	void IlinAGaussianMethodMPI::HandleRowSwapLocal(int k_local_idx, int pivot_owner, int global_max_row) {
236		✗	if (k_local_idx < 0) {
237			return;
238			}
239
240		✗	if (pivot_owner == rank_) {
241			int pivot_local_idx = FindLocalRowIndex(global_max_row);
242		✗	if (pivot_local_idx >= 0) {
243		✗	SwapRowsLocally(k_local_idx, pivot_local_idx);
244			}
245			} else {
246		✗	ExchangeRowsWithRemote(k_local_idx, pivot_owner);
247			}
248			}
249
250		✗	void IlinAGaussianMethodMPI::HandleRowSwapRemote(int k_owner, int global_max_row) {
251			int pivot_local_idx = FindLocalRowIndex(global_max_row);
252		✗	if (pivot_local_idx < 0) {
253			return;
254			}
255
256		✗	ReceiveRowFromRemote(pivot_local_idx, k_owner);
257			}
258
259		✗	void IlinAGaussianMethodMPI::SwapRowsLocally(int k_local_idx, int pivot_local_idx) {
260		✗	std::swap_ranges(&local_matrix_[(static_cast<size_t>(k_local_idx) * static_cast<size_t>(band_))],
261		✗	&local_matrix_[(static_cast<size_t>(k_local_idx) * static_cast<size_t>(band_))] + band_,
262		✗	&local_matrix_[(static_cast<size_t>(pivot_local_idx) * static_cast<size_t>(band_))]);
263			std::swap(local_vector_[static_cast<size_t>(k_local_idx)], local_vector_[static_cast<size_t>(pivot_local_idx)]);
264		✗	}
265
266		✗	void IlinAGaussianMethodMPI::ExchangeRowsWithRemote(int k_local_idx, int pivot_owner) {
267		✗	std::vector<double> temp_row(band_);
268		✗	double temp_b = 0.0;
269			MPI_Status status;
270
271		✗	MPI_Recv(temp_row.data(), band_, MPI_DOUBLE, pivot_owner, 0, MPI_COMM_WORLD, &status);
272		✗	MPI_Recv(&temp_b, 1, MPI_DOUBLE, pivot_owner, 1, MPI_COMM_WORLD, &status);
273
274		✗	std::ranges::copy(temp_row, &local_matrix_[(static_cast<size_t>(k_local_idx) * static_cast<size_t>(band_))]);
275		✗	local_vector_[static_cast<size_t>(k_local_idx)] = temp_b;
276
277		✗	MPI_Send(&local_matrix_[(static_cast<size_t>(k_local_idx) * static_cast<size_t>(band_))], band_, MPI_DOUBLE,
278			pivot_owner, 2, MPI_COMM_WORLD);
279		✗	MPI_Send(&local_vector_[static_cast<size_t>(k_local_idx)], 1, MPI_DOUBLE, pivot_owner, 3, MPI_COMM_WORLD);
280		✗	}
281
282		✗	void IlinAGaussianMethodMPI::ReceiveRowFromRemote(int pivot_local_idx, int k_owner) {
283		✗	MPI_Send(&local_matrix_[(static_cast<size_t>(pivot_local_idx) * static_cast<size_t>(band_))], band_, MPI_DOUBLE,
284			k_owner, 0, MPI_COMM_WORLD);
285		✗	MPI_Send(&local_vector_[static_cast<size_t>(pivot_local_idx)], 1, MPI_DOUBLE, k_owner, 1, MPI_COMM_WORLD);
286
287		✗	std::vector<double> temp_row(band_);
288		✗	double temp_b = 0.0;
289			MPI_Status status;
290
291		✗	MPI_Recv(temp_row.data(), band_, MPI_DOUBLE, k_owner, 2, MPI_COMM_WORLD, &status);
292		✗	MPI_Recv(&temp_b, 1, MPI_DOUBLE, k_owner, 3, MPI_COMM_WORLD, &status);
293
294		✗	std::ranges::copy(temp_row, &local_matrix_[(static_cast<size_t>(pivot_local_idx) * static_cast<size_t>(band_))]);
295		✗	local_vector_[static_cast<size_t>(pivot_local_idx)] = temp_b;
296		✗	}
297
298		✗	void IlinAGaussianMethodMPI::EliminateRow(int i, int k, const std::vector<double> &pivot_row, double pivot_b) {
299			const double eps = 1e-12;
300		✗	int global_row = row_start_ + i;
301
302		✗	if (global_row <= k) {
303			return;
304			}
305
306		✗	int factor_idx = band_ - 1 - (global_row - k);
307		✗	if (factor_idx < 0 || factor_idx >= band_) {
308			return;
309			}
310
311		✗	double pivot = pivot_row[static_cast<size_t>(band_ - 1)];
312		✗	double factor = local_matrix_[(static_cast<size_t>(i) * static_cast<size_t>(band_)) + factor_idx] / pivot;
313
314		✗	if (std::fabs(factor) <= eps) {
315			return;
316			}
317
318		✗	UpdateRowValues(i, k, factor, pivot_row);
319		✗	local_vector_[static_cast<size_t>(i)] -= factor * pivot_b;
320			}
321
322		✗	void IlinAGaussianMethodMPI::UpdateRowValues(int i, int k, double factor, const std::vector<double> &pivot_row) {
323		✗	int global_row = row_start_ + i;
324
325		✗	for (int j = 0; j < band_; ++j) {
326		✗	int src_idx = j - (global_row - k);
327		✗	if (src_idx < 0 || src_idx >= band_) {
328		✗	continue;
329			}
330
331		✗	local_matrix_[(static_cast<size_t>(i) * static_cast<size_t>(band_)) + j] -=
332		✗	factor * pivot_row[static_cast<size_t>(src_idx)];
333			}
334		✗	}
335
336		✗	void IlinAGaussianMethodMPI::ProcessBackwardSubstitution() {
337		✗	std::vector<double> recv_matrix;
338		✗	std::vector<double> recv_vector;
339
340		✗	GatherAllData(recv_matrix, recv_vector);
341
342		✗	if (rank_ != 0) {
343		✗	MPI_Bcast(solution_.data(), n_, MPI_DOUBLE, 0, MPI_COMM_WORLD);
344			return;
345			}
346
347		✗	std::vector<double> full_matrix(static_cast<size_t>(n_) * static_cast<size_t>(band_));
348		✗	std::vector<double> full_vector(static_cast<size_t>(n_));
349
350		✗	ReconstructFullMatrix(recv_matrix, recv_vector, full_matrix, full_vector);
351		✗	SolveBackwardSubstitution(full_matrix, full_vector);
352
353		✗	MPI_Bcast(solution_.data(), n_, MPI_DOUBLE, 0, MPI_COMM_WORLD);
354			}
355
356		✗	void IlinAGaussianMethodMPI::GatherAllData(std::vector<double> &recv_matrix, std::vector<double> &recv_vector) {
357		✗	if (rank_ == 0) {
358		✗	recv_matrix.resize(static_cast<size_t>(n_) * static_cast<size_t>(band_));
359		✗	recv_vector.resize(static_cast<size_t>(n_));
360			}
361
362		✗	std::vector<int> recv_counts(static_cast<size_t>(size_));
363		✗	std::vector<int> displs(static_cast<size_t>(size_));
364		✗	std::vector<int> vec_counts(static_cast<size_t>(size_));
365		✗	std::vector<int> vec_displs(static_cast<size_t>(size_));
366
367		✗	for (int i = 0; i < size_; ++i) {
368		✗	int rows_for_i = (i < remainder_) ? (rows_per_proc_ + 1) : rows_per_proc_;
369		✗	recv_counts[static_cast<size_t>(i)] = rows_for_i * band_;
370		✗	vec_counts[static_cast<size_t>(i)] = rows_for_i;
371
372		✗	if (i == 0) {
373		✗	displs[0] = 0;
374		✗	vec_displs[0] = 0;
375			} else {
376		✗	displs[static_cast<size_t>(i)] = displs[static_cast<size_t>(i - 1)] + recv_counts[static_cast<size_t>(i - 1)];
377		✗	vec_displs[static_cast<size_t>(i)] =
378		✗	vec_displs[static_cast<size_t>(i - 1)] + vec_counts[static_cast<size_t>(i - 1)];
379			}
380			}
381
382		✗	MPI_Gatherv(local_matrix_.data(), local_rows_ * band_, MPI_DOUBLE, recv_matrix.data(), recv_counts.data(),
383			displs.data(), MPI_DOUBLE, 0, MPI_COMM_WORLD);
384		✗	MPI_Gatherv(local_vector_.data(), local_rows_, MPI_DOUBLE, recv_vector.data(), vec_counts.data(), vec_displs.data(),
385			MPI_DOUBLE, 0, MPI_COMM_WORLD);
386		✗	}
387
388		✗	void IlinAGaussianMethodMPI::ReconstructFullMatrix(const std::vector<double> &recv_matrix,
389			const std::vector<double> &recv_vector,
390			std::vector<double> &full_matrix,
391			std::vector<double> &full_vector) const {
392		✗	std::vector<int> displs(static_cast<size_t>(size_));
393		✗	std::vector<int> vec_displs(static_cast<size_t>(size_));
394
395		✗	for (int i = 0; i < size_; ++i) {
396		✗	if (i == 0) {
397		✗	displs[0] = 0;
398		✗	vec_displs[0] = 0;
399			} else {
400		✗	int prev_rows = (i - 1 < remainder_) ? (rows_per_proc_ + 1) : rows_per_proc_;
401		✗	displs[static_cast<size_t>(i)] = displs[static_cast<size_t>(i - 1)] + (prev_rows * band_);
402		✗	vec_displs[static_cast<size_t>(i)] = vec_displs[static_cast<size_t>(i - 1)] + prev_rows;
403			}
404			}
405
406		✗	for (int i = 0; i < size_; ++i) {
407		✗	int rows_for_i = (i < remainder_) ? (rows_per_proc_ + 1) : rows_per_proc_;
408		✗	int start_row = vec_displs[static_cast<size_t>(i)];
409
410		✗	for (int j = 0; j < rows_for_i; ++j) {
411		✗	int global_row = start_row + j;
412		✗	const std::ptrdiff_t offset = static_cast<std::ptrdiff_t>(displs[static_cast<size_t>(i)]) +
413		✗	(static_cast<std::ptrdiff_t>(j) * static_cast<std::ptrdiff_t>(band_));
414			std::copy(&recv_matrix[static_cast<size_t>(offset)],
415		✗	&recv_matrix[static_cast<size_t>(offset) + static_cast<size_t>(band_)],
416		✗	&full_matrix[(static_cast<size_t>(global_row) * static_cast<size_t>(band_))]);
417		✗	full_vector[static_cast<size_t>(global_row)] =
418		✗	recv_vector[static_cast<size_t>(vec_displs[static_cast<size_t>(i)] + static_cast<size_t>(j))];
419			}
420			}
421		✗	}
422
423		✗	void IlinAGaussianMethodMPI::SolveBackwardSubstitution(const std::vector<double> &full_matrix,
424			const std::vector<double> &full_vector) {
425			const double eps = 1e-12;
426
427		✗	for (int i = n_ - 1; i >= 0; --i) {
428			double sum = 0.0;
429
430		✗	for (int j = i + 1; j < std::min(n_, i + band_); ++j) {
431		✗	int idx = band_ - 1 + (j - i);
432		✗	if (idx >= band_) {
433		✗	continue;
434			}
435
436		✗	sum +=
437		✗	full_matrix[(static_cast<size_t>(i) * static_cast<size_t>(band_)) + idx] * solution_[static_cast<size_t>(j)];
438			}
439
440		✗	int diag_idx = band_ - 1;
441		✗	double diag = full_matrix[(static_cast<size_t>(i) * static_cast<size_t>(band_)) + diag_idx];
442
443		✗	if (std::fabs(diag) > eps) {
444		✗	solution_[static_cast<size_t>(i)] = (full_vector[static_cast<size_t>(i)] - sum) / diag;
445			} else {
446		✗	solution_[static_cast<size_t>(i)] = 0.0;
447			}
448			}
449
450		✗	for (int i = 0; i < n_; ++i) {
451		✗	if (!std::isfinite(solution_[static_cast<size_t>(i)])) {
452		✗	solution_[static_cast<size_t>(i)] = 0.0;
453			}
454			}
455		✗	}
456
457		✗	bool IlinAGaussianMethodMPI::PostProcessingImpl() {
458		✗	MPI_Barrier(MPI_COMM_WORLD);
459
460		✗	if (rank_ == 0) {
461		✗	GetOutput() = solution_;
462			}
463
464		✗	MPI_Bcast(solution_.data(), n_, MPI_DOUBLE, 0, MPI_COMM_WORLD);
465
466		✗	return true;
467			}
468
469			} // namespace ilin_a_gaussian_method_horizontal_band_scheme
470