29 const std::size_t locally_owned_size,
30 const Number *current_vector,
31 const std::vector<const Number *> &orthogonal_vectors,
39 static constexpr unsigned int n_lanes =
43 for (
unsigned int i = 0; i < n_vectors; ++i)
46 correct[delayed_reorthogonalization ? 129 : 1];
47 if (delayed_reorthogonalization)
48 for (
unsigned int i = 0; i < n_vectors + 1; ++i)
53 constexpr unsigned int inner_batch_size =
54 delayed_reorthogonalization ? 6 : 12;
56 for (; c < locally_owned_size / n_lanes / inner_batch_size;
57 ++c, j += n_lanes * inner_batch_size)
60 for (
unsigned int k = 0; k < inner_batch_size; ++k)
61 vvec[k].load(current_vector + j + k * n_lanes);
63 for (
unsigned int k = 0; k < inner_batch_size; ++k)
64 prev_vector[k].load(orthogonal_vectors[n_vectors - 1] + j +
70 prev_vector[0] * prev_vector[0];
72 for (
unsigned int k = 1; k < inner_batch_size; ++k)
74 local_sum_0 += prev_vector[k] * vvec[k];
75 if (delayed_reorthogonalization)
77 local_sum_1 += prev_vector[k] * prev_vector[k];
78 local_sum_2 += vvec[k] * vvec[k];
81 hs[n_vectors - 1] += local_sum_0;
82 if (delayed_reorthogonalization)
84 correct[n_vectors - 1] += local_sum_1;
85 correct[n_vectors] += local_sum_2;
89 for (
unsigned int i = 0; i < n_vectors - 1; ++i)
95 temp.
load(orthogonal_vectors[i] + j);
98 delayed_reorthogonalization ? temp * prev_vector[0] : 0.;
99 for (
unsigned int k = 1; k < inner_batch_size; ++k)
101 temp.
load(orthogonal_vectors[i] + j + k * n_lanes);
102 local_sum_0 += temp * vvec[k];
103 if (delayed_reorthogonalization)
104 local_sum_1 += temp * prev_vector[k];
106 hs[i] += local_sum_0;
107 if (delayed_reorthogonalization)
108 correct[i] += local_sum_1;
112 c *= inner_batch_size;
113 for (; c < locally_owned_size / n_lanes; ++c, j += n_lanes)
116 vvec.
load(current_vector + j);
117 prev_vector.
load(orthogonal_vectors[n_vectors - 1] + j);
118 hs[n_vectors - 1] += prev_vector * vvec;
119 if (delayed_reorthogonalization)
121 correct[n_vectors - 1] += prev_vector * prev_vector;
122 correct[n_vectors] += vvec * vvec;
125 for (
unsigned int i = 0; i < n_vectors - 1; ++i)
128 temp.
load(orthogonal_vectors[i] + j);
129 hs[i] += temp * vvec;
130 if (delayed_reorthogonalization)
131 correct[i] += temp * prev_vector;
135 for (
unsigned int i = 0; i < n_vectors; ++i)
138 if (delayed_reorthogonalization)
139 h(i + n_vectors) += correct[i].
sum();
141 if (delayed_reorthogonalization)
142 h(n_vectors + n_vectors) += correct[n_vectors].
sum();
147 for (; j < locally_owned_size; ++j)
149 const double vvec = current_vector[j];
150 const double prev_vector = orthogonal_vectors[n_vectors - 1][j];
151 h(n_vectors - 1) += prev_vector * vvec;
152 if (delayed_reorthogonalization)
154 h(n_vectors + n_vectors - 1) += prev_vector * prev_vector;
155 h(n_vectors + n_vectors) += vvec * vvec;
157 for (
unsigned int i = 0; i < n_vectors - 1; ++i)
159 const double temp = orthogonal_vectors[i][j];
161 if (delayed_reorthogonalization)
162 h(n_vectors + i) += temp * prev_vector;
172 const std::size_t locally_owned_size,
173 const std::vector<const Number *> &orthogonal_vectors,
175 Number *current_vector)
177 double norm_vv_temp = 0;
179 Number *previous_vector =
180 const_cast<Number *
>(orthogonal_vectors[n_vectors - 1]);
181 const double inverse_norm_previous =
182 delayed_reorthogonalization ? 1. / h(n_vectors + n_vectors - 1) : 0.;
183 const double scaling_factor_vv =
184 delayed_reorthogonalization ?
185 (h(n_vectors + n_vectors) > 0.0 ?
186 inverse_norm_previous / h(n_vectors + n_vectors) :
187 inverse_norm_previous / h(n_vectors + n_vectors - 1)) :
192 constexpr unsigned int inner_batch_size =
193 delayed_reorthogonalization ? 6 : 12;
197 for (; c < locally_owned_size / n_lanes / inner_batch_size;
198 ++c, j += n_lanes * inner_batch_size)
203 const double last_factor = h(n_vectors - 1);
204 for (
unsigned int k = 0; k < inner_batch_size; ++k)
206 temp[k].
load(current_vector + j + k * n_lanes);
207 prev_vector[k].
load(previous_vector + j + k * n_lanes);
208 if (!delayed_reorthogonalization)
209 temp[k] -= last_factor * prev_vector[k];
212 for (
unsigned int i = 0; i < n_vectors - 1; ++i)
214 const double factor = h(i);
215 const double correction_factor =
216 (delayed_reorthogonalization ? h(n_vectors + i) : 0.0);
217 for (
unsigned int k = 0; k < inner_batch_size; ++k)
220 vec.
load(orthogonal_vectors[i] + j + k * n_lanes);
221 temp[k] -= factor * vec;
222 if (delayed_reorthogonalization)
223 prev_vector[k] -= correction_factor * vec;
227 if (delayed_reorthogonalization)
228 for (
unsigned int k = 0; k < inner_batch_size; ++k)
230 prev_vector[k] = prev_vector[k] * inverse_norm_previous;
231 prev_vector[k].
store(previous_vector + j + k * n_lanes);
232 temp[k] -= last_factor * prev_vector[k];
233 temp[k] = temp[k] * scaling_factor_vv;
234 temp[k].
store(current_vector + j + k * n_lanes);
237 for (
unsigned int k = 0; k < inner_batch_size; ++k)
239 temp[k].
store(current_vector + j + k * n_lanes);
240 norm_vv_temp_vectorized += temp[k] * temp[k];
244 c *= inner_batch_size;
245 for (; c < locally_owned_size / n_lanes; ++c, j += n_lanes)
248 temp.
load(current_vector + j);
249 prev_vector.
load(previous_vector + j);
250 if (!delayed_reorthogonalization)
251 temp -= h(n_vectors - 1) * prev_vector;
253 for (
unsigned int i = 0; i < n_vectors - 1; ++i)
256 vec.
load(orthogonal_vectors[i] + j);
258 if (delayed_reorthogonalization)
259 prev_vector -= h(n_vectors + i) * vec;
262 if (delayed_reorthogonalization)
264 prev_vector = prev_vector * inverse_norm_previous;
265 prev_vector.
store(previous_vector + j);
266 temp -= h(n_vectors - 1) * prev_vector;
267 temp = temp * scaling_factor_vv;
268 temp.
store(current_vector + j);
272 temp.
store(current_vector + j);
273 norm_vv_temp_vectorized += temp * temp;
277 if (!delayed_reorthogonalization)
278 norm_vv_temp += norm_vv_temp_vectorized.
sum();
280 for (; j < locally_owned_size; ++j)
282 double temp = current_vector[j];
283 double prev_vector = previous_vector[j];
284 if (delayed_reorthogonalization)
286 for (
unsigned int i = 0; i < n_vectors - 1; ++i)
288 const double vec = orthogonal_vectors[i][j];
290 prev_vector -= h(n_vectors + i) * vec;
292 prev_vector *= inverse_norm_previous;
293 previous_vector[j] = prev_vector;
294 temp -= h(n_vectors - 1) * prev_vector;
295 temp *= scaling_factor_vv;
299 temp -= h(n_vectors - 1) * prev_vector;
300 for (
unsigned int i = 0; i < n_vectors - 1; ++i)
301 temp -= h(i) * orthogonal_vectors[i][j];
302 norm_vv_temp += temp * temp;
304 current_vector[j] = temp;
315 const std::size_t locally_owned_size,
316 const std::vector<const Number *> &tmp_vectors,
322 constexpr unsigned int inner_batch_size = 12;
326 for (; c < locally_owned_size / n_lanes / inner_batch_size;
327 ++c, j += n_lanes * inner_batch_size)
334 for (
unsigned int k = 0; k < inner_batch_size; ++k)
335 temp[k].load(output + j + k * n_lanes);
337 for (
unsigned int i = 0; i < n_vectors; ++i)
339 const double h_i = h(i);
340 for (
unsigned int k = 0; k < inner_batch_size; ++k)
343 v_ij.
load(tmp_vectors[i] + j + k * n_lanes);
344 temp[k] += v_ij * h_i;
348 for (
unsigned int k = 0; k < inner_batch_size; ++k)
349 temp[k].store(output + j + k * n_lanes);
352 c *= inner_batch_size;
353 for (; c < locally_owned_size / n_lanes; ++c, j += n_lanes)
357 temp.
load(output + j);
359 for (
unsigned int i = 0; i < n_vectors; ++i)
362 v_ij.
load(tmp_vectors[i] + j);
366 temp.
store(output + j);
369 for (; j < locally_owned_size; ++j)
371 double temp = zero_out ? 0.0 : output[j];
372 for (
unsigned int i = 0; i < n_vectors; ++i)
373 temp += tmp_vectors[i][j] * h(i);