GCC Code Coverage Report

Directory: ./
File: tasks/kichanova_k_lin_system_by_conjug_grad/all/src/ops_all.cpp
Date: 2026-06-04 20:25:32
Exec Total Coverage
Lines: 83 94 88.3%
Functions: 13 13 100.0%
Branches: 31 56 55.4%
Line Branch Exec Source
1 #include "kichanova_k_lin_system_by_conjug_grad/all/include/ops_all.hpp"
2
3 #include <omp.h>
4 #include <tbb/tbb.h>
5
6 #include <cmath>
7 #include <cstddef>
8 #include <vector>
9
10 #include "kichanova_k_lin_system_by_conjug_grad/common/include/common.hpp"
11
12 namespace kichanova_k_lin_system_by_conjug_grad {
13
14 namespace {
15
16 298 double ComputeRowProduct(const double *row, const double *v_ptr, int n) {
17 double sum = 0.0;
18 int j = 0;
19
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 298 times.
 298 for (; j <= n - 8; j += 8) {
20 ✗ sum += row[j] * v_ptr[j];
21 ✗ sum += row[j + 1] * v_ptr[j + 1];
22 ✗ sum += row[j + 2] * v_ptr[j + 2];
23 ✗ sum += row[j + 3] * v_ptr[j + 3];
24 ✗ sum += row[j + 4] * v_ptr[j + 4];
25 ✗ sum += row[j + 5] * v_ptr[j + 5];
26 ✗ sum += row[j + 6] * v_ptr[j + 6];
27 ✗ sum += row[j + 7] * v_ptr[j + 7];
28 }
29
2/2
✓ Branch 0 taken 1642 times.
✓ Branch 1 taken 298 times.
 1940 for (; j < n; ++j) {
30 1642 sum += row[j] * v_ptr[j];
31 }
32 298 return sum;
33 }
34
35 138 double DotProductHybrid(const std::vector<double> &a, const std::vector<double> &b, int n) {
36
1/2
✓ Branch 0 taken 138 times.
✗ Branch 1 not taken.
 138 if (n <= 0) {
37 return 0.0;
38 }
39
40 double result = 0.0;
41 138 int num_threads = omp_get_max_threads();
42
43 138 #pragma omp parallel for reduction(+ : result) schedule(static) default(none) shared(a, b, n, num_threads)
44 for (int block = 0; block < num_threads; ++block) {
45 int start = block * (n / num_threads);
46 int end = (block == num_threads - 1) ? n : (block + 1) * (n / num_threads);
47 if (start >= end) {
48 continue;
49 }
50
51 double local_sum = tbb::parallel_reduce(tbb::blocked_range<int>(start, end, 256), 0.0,
52 [&](const tbb::blocked_range<int> &range, double sum) {
53
2/2
✓ Branch 0 taken 670 times.
✓ Branch 1 taken 276 times.
 946 for (int i = range.begin(); i < range.end(); ++i) {
54 670 sum += a[i] * b[i];
55 }
56 return sum;
57 ✗ }, [](double x, double y) { return x + y; });
58 result += local_sum;
59 }
60 138 return result;
61 }
62
63 60 void MatrixVectorProductHybrid(const std::vector<double> &a, const std::vector<double> &v, std::vector<double> &result,
64 int n) {
65
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 60 times.
 60 if (n <= 0) {
66 ✗ return;
67 }
68
69 60 const auto stride = static_cast<size_t>(n);
70 60 const double *v_ptr = v.data();
71 60 int num_threads = omp_get_max_threads();
72
73 60 #pragma omp parallel for schedule(dynamic, 1) default(none) shared(a, result, n, stride, v_ptr, num_threads)
74 for (int block = 0; block < num_threads; ++block) {
75 int start = block * (n / num_threads);
76 int end = (block == num_threads - 1) ? n : (block + 1) * (n / num_threads);
77 if (start >= end) {
78 continue;
79 }
80
81 120 tbb::parallel_for(tbb::blocked_range<int>(start, end, 32), [&](const tbb::blocked_range<int> &range) {
82
2/2
✓ Branch 0 taken 298 times.
✓ Branch 1 taken 120 times.
 418 for (int i = range.begin(); i < range.end(); ++i) {
83 298 const double *row = &a[i * stride];
84 298 result[i] = ComputeRowProduct(row, v_ptr, n);
85 }
86 120 });
87 }
88 }
89
90 60 void UpdateSolutionAndResidualHybrid(std::vector<double> &x, std::vector<double> &r, const std::vector<double> &p,
91 const std::vector<double> &ap, double alpha, int n) {
92
1/2
✓ Branch 0 taken 60 times.
✗ Branch 1 not taken.
 60 if (n <= 0) {
93 return;
94 }
95
96 60 int num_threads = omp_get_max_threads();
97
98 60 #pragma omp parallel for schedule(static) default(none) shared(x, r, p, ap, alpha, n, num_threads)
99 for (int block = 0; block < num_threads; ++block) {
100 int start = block * (n / num_threads);
101 int end = (block == num_threads - 1) ? n : (block + 1) * (n / num_threads);
102 if (start >= end) {
103 continue;
104 }
105
106 120 tbb::parallel_for(tbb::blocked_range<int>(start, end, 512), [&](const tbb::blocked_range<int> &range) {
107
2/2
✓ Branch 0 taken 298 times.
✓ Branch 1 taken 120 times.
 418 for (int i = range.begin(); i < range.end(); ++i) {
108 298 x[i] += alpha * p[i];
109 298 r[i] -= alpha * ap[i];
110 }
111 120 });
112 }
113 }
114
115 42 void UpdateDirectionHybrid(std::vector<double> &p, const std::vector<double> &r, double beta, int n) {
116
1/2
✓ Branch 0 taken 42 times.
✗ Branch 1 not taken.
 42 if (n <= 0) {
117 return;
118 }
119
120 42 int num_threads = omp_get_max_threads();
121
122 42 #pragma omp parallel for schedule(static) default(none) shared(p, r, beta, n, num_threads)
123 for (int block = 0; block < num_threads; ++block) {
124 int start = block * (n / num_threads);
125 int end = (block == num_threads - 1) ? n : (block + 1) * (n / num_threads);
126 if (start >= end) {
127 continue;
128 }
129
130 tbb::parallel_for(tbb::blocked_range<int>(start, end, 512), [&](const tbb::blocked_range<int> &range) {
131
2/4
✓ Branch 0 taken 224 times.
✓ Branch 1 taken 84 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 308 for (int i = range.begin(); i < range.end(); ++i) {
132 224 p[i] = r[i] + (beta * p[i]);
133 }
134 });
135 }
136 }
137
138 18 void InitializeVectorsHybrid(std::vector<double> &r, std::vector<double> &p, const std::vector<double> &b, int n) {
139 18 int num_threads = omp_get_max_threads();
140
141 18 #pragma omp parallel for schedule(static) default(none) shared(r, p, b, n, num_threads)
142 for (int block = 0; block < num_threads; ++block) {
143 int start = block * (n / num_threads);
144 int end = (block == num_threads - 1) ? n : (block + 1) * (n / num_threads);
145 if (start >= end) {
146 continue;
147 }
148
149 tbb::parallel_for(tbb::blocked_range<int>(start, end, 1024), [&](const tbb::blocked_range<int> &range) {
150
2/4
✓ Branch 0 taken 74 times.
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 110 for (int i = range.begin(); i < range.end(); ++i) {
151 74 r[i] = b[i];
152 74 p[i] = r[i];
153 }
154 });
155 }
156 18 }
157
158 } // namespace
159
160
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 18 KichanovaKLinSystemByConjugGradALL::KichanovaKLinSystemByConjugGradALL(const InType &in) {
161 SetTypeOfTask(GetStaticTypeOfTask());
162
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 18 GetInput() = in;
163 18 GetOutput() = OutType();
164 18 }
165
166 18 bool KichanovaKLinSystemByConjugGradALL::ValidationImpl() {
167 const InType &input_data = GetInput();
168
1/2
✓ Branch 0 taken 18 times.
✗ Branch 1 not taken.
 18 if (input_data.n <= 0) {
169 return false;
170 }
171
1/2
✓ Branch 0 taken 18 times.
✗ Branch 1 not taken.
 18 if (input_data.A.size() != static_cast<size_t>(input_data.n) * static_cast<size_t>(input_data.n)) {
172 return false;
173 }
174
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18 times.
 18 if (input_data.b.size() != static_cast<size_t>(input_data.n)) {
175 ✗ return false;
176 }
177 return true;
178 }
179
180 18 bool KichanovaKLinSystemByConjugGradALL::PreProcessingImpl() {
181 18 GetOutput().assign(GetInput().n, 0.0);
182 18 return true;
183 }
184
185
1/2
✓ Branch 0 taken 18 times.
✗ Branch 1 not taken.
 18 bool KichanovaKLinSystemByConjugGradALL::RunImpl() {
186 const InType &input_data = GetInput();
187 OutType &x = GetOutput();
188
189 18 const int n = input_data.n;
190
1/2
✓ Branch 0 taken 18 times.
✗ Branch 1 not taken.
 18 if (n == 0) {
191 return false;
192 }
193
194 18 const std::vector<double> &a = input_data.A;
195 18 const std::vector<double> &b = input_data.b;
196 18 const double epsilon = input_data.epsilon;
197
198 18 std::vector<double> r(n);
199
1/4
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 18 std::vector<double> p(n);
200
1/4
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 18 std::vector<double> ap(n);
201
202 18 InitializeVectorsHybrid(r, p, b, n);
203
204 18 double rr_old = DotProductHybrid(r, r, n);
205 18 double residual_norm = std::sqrt(rr_old);
206
207
1/2
✓ Branch 0 taken 18 times.
✗ Branch 1 not taken.
 18 if (residual_norm < epsilon) {
208 return true;
209 }
210
211 18 const int max_iter = n * 1000;
212
213
1/2
✓ Branch 0 taken 60 times.
✗ Branch 1 not taken.
 60 for (int iter = 0; iter < max_iter; ++iter) {
214 60 MatrixVectorProductHybrid(a, p, ap, n);
215
216 60 const double p_ap = DotProductHybrid(p, ap, n);
217
1/2
✓ Branch 0 taken 60 times.
✗ Branch 1 not taken.
 60 if (std::abs(p_ap) < 1e-30) {
218 break;
219 }
220
221 60 const double alpha = rr_old / p_ap;
222 60 UpdateSolutionAndResidualHybrid(x, r, p, ap, alpha, n);
223
224 60 const double rr_new = DotProductHybrid(r, r, n);
225 60 residual_norm = std::sqrt(rr_new);
226
227
2/2
✓ Branch 0 taken 42 times.
✓ Branch 1 taken 18 times.
 60 if (residual_norm < epsilon) {
228 break;
229 }
230
231 42 const double beta = rr_new / rr_old;
232 42 UpdateDirectionHybrid(p, r, beta, n);
233
234 rr_old = rr_new;
235 }
236
237 return true;
238 }
239
240 18 bool KichanovaKLinSystemByConjugGradALL::PostProcessingImpl() {
241 18 return true;
242 }
243
244 } // namespace kichanova_k_lin_system_by_conjug_grad
245