216{
217 const size_t maxIter = nAllVars + 1;
218 std::vector<double> x(nAllVars, 0.0);
219 std::vector<double> y(nAllVars, 0.0);
220 std::vector<double> residual(nAllVars, 0.0);
221 std::vector<double> r(nAllVars, 0.0);
222 std::vector<double> w(nAllVars, 0.0);
223 std::vector<double> wOld(nAllVars, 0.0);
224 std::vector<double> g(nAllVars, 0.0);
225 std::vector<double> c(nAllVars, 0.0);
226 std::vector<double> s(nAllVars, 0.0);
227 std::vector<double> hisGs(nAllVars + 1, 0.0);
228 std::vector<double> v(nAllVars * (maxIter + 1), 0.0);
229
230 std::vector<double> vecFromV(nAllVars, 0.0);
231
232 std::vector<double> h((nAllVars + 1) * maxIter, 0.0);
233 double beta = 0.0, gs, normb = 0.0;
234 size_t m = 0;
235
236 for (size_t i = 0; i < nAllVars; ++i)
237 normb += rhsDir[i] * rhsDir[i];
238 normb = sqrt(normb);
239
240 for (size_t i = 0; i < nAllVars; ++i)
241 residual[i] = rhsDir[i];
242
243#pragma omp parallel for
244 for (int aflI = 0; aflI < nafl; ++aflI)
245 {
246
247 r = residual;
248 }
249
250 for (size_t i = 0; i < nAllVars; ++i)
251 beta += r[i] * r[i];
252 beta = sqrt(beta);
253 gs = beta;
254 g[0] = beta;
255
256
257 if (beta == 0)
258 {
259 for (size_t aflI = 0; aflI < nafl; ++aflI)
260 {
261 for (size_t i = 0; i < vsize[aflI]; ++i)
262 gam[aflI][i] = 0.0;
263 R[aflI] = 0.0;
264 }
265 return;
266 }
267
268 for (size_t aflI = 0; aflI < nafl; ++aflI)
269 R[aflI] = x[x.size() - nafl + aflI];
270
271 for (size_t i = 0; i < nAllVars; ++i)
272 v[i * (maxIter + 1) + 0] = r[i] / beta;
273
274 for (size_t j = 0; j < nAllVars; ++j)
275 {
276 double timer1 = omp_get_wtime();
277
278 for (size_t p = 0; p < nAllVars; ++p)
279 vecFromV[p] = v[p * (maxIter + 1) + j];
280
281
282
283 for (int i = 0; i < nAllVars; ++i)
284 {
285 wOld[i] = 0.0;
286 for (int jj = 0; jj < nAllVars; ++jj)
287 {
288 wOld[i] += mtrDir[i * nAllVars + jj] * vecFromV[jj];
289 }
290 }
291
292 double timer2 = omp_get_wtime();
293
294#pragma omp parallel for
295 for (int aflI = 0; aflI < nafl; ++aflI)
296 {
297
298 w = wOld;
299 }
300
301 for (size_t i = 0; i <= j; ++i)
302 {
303 h[i * maxIter + j] = 0.0;
304 for (size_t q = 0; q < nAllVars; ++q)
305 h[i * maxIter + j] += w[q] * v[q * (maxIter + 1) + i];
306
307 for (size_t q = 0; q < nAllVars; ++q)
308 w[q] -= h[i * maxIter + j] * v[q * (maxIter + 1) + i];
309 }
310
311 h[(j + 1) * maxIter + j] = 0.0;
312 for (size_t p = 0; p < nAllVars; ++p)
313 h[(j + 1) * maxIter + j] += w[p] * w[p];
314 h[(j + 1) * maxIter + j] = sqrt(h[(j + 1) * maxIter + j]);
315
316 for (size_t p = 0; p < nAllVars; ++p)
317 v[p * (maxIter + 1) + (j + 1)] = w[p] / h[(j + 1) * maxIter + j];
318
319 for (size_t i = 0; i < j; ++i)
320 {
321 double buf = h[i * maxIter + j];
322 h[i * maxIter + j] = c[i] * buf + s[i] * h[(i + 1) * maxIter + j];
323 h[(i + 1) * maxIter + j] = -s[i] * buf + c[i] * h[(i + 1) * maxIter + j];
324 }
325
326 double zn = sqrt(h[j * maxIter + j] * h[j * maxIter + j] + h[(j + 1) * maxIter + j] * h[(j + 1) * maxIter + j]);
327 c[j] = fabs(h[j * maxIter + j] / zn);
328 s[j] = c[j] * h[(j + 1) * maxIter + j] / h[j * maxIter + j];
329
330 gs *= -s[j];
331
332 h[j * maxIter + j] = c[j] * h[j * maxIter + j] + s[j] * h[(j + 1) * maxIter + j];
333 h[(j + 1) * maxIter + j] = 0.0;
334 hisGs[j] = fabs(gs) / normb;
335
337 {
338 m = j + 1;
339 break;
340 }
341 }
342
343
344 for (size_t i = 0; i < m; ++i)
345 {
346 double buf = g[i];
347 g[i] *= c[i];
348 if (i < (m - 1))
349 g[i + 1] = -s[i] * buf;
350
351 }
352 y[m - 1] = g[m - 1] / h[(m - 1) * maxIter + (m - 1)];
353
354 for (size_t i = m - 2; i + 1 > 0; --i)
355 {
356 double sum = 0;
357 for (size_t j = i + 1; j < m; ++j)
358 sum += h[i * maxIter + j] * y[j];
359
360 y[i] = (g[i] - sum) / h[i * maxIter + i];
361 }
362
363 for (size_t p = 0; p < nAllVars; ++p)
364 {
365 for (size_t q = 0; q < m; ++q)
366 x[p] += v[p * (maxIter + 1) + q] * y[q];
367 }
368
369 std::cout << "#iterations: " << m << std::endl;
370
371 for (size_t aflI = 0; aflI < nafl; ++aflI)
372 {
373 for (size_t i = 0; i < vsize[aflI]; ++i)
374 gam[aflI][i] = x[pos[aflI] + i];
375 }
376
377 for (size_t aflI = 0; aflI < nafl; ++aflI)
378 R[aflI] = x[x.size() - nafl + aflI];
379}
NumericalSchemes numericalSchemes
Структура с используемыми численными схемами
const Passport & getPassport() const
Возврат константной ссылки на паспорт
1.9.8