omp2/solve.cpp - master - NenNil/vmk2025-parallel-computing

omp2/solve.cpp
128 строк3 KB
Нил Ненахов
some not essential changes
03 дек 2025, 16:49
03 дек 2025, 16:4902715ff
#include "solve.h"
#include "matrix.h"
#include <cmath>
#include <iostream>
#include <mpi.h>
#include <string>
#include <vector>
#include <chrono>

using namespace std;

void printVector(const std::vector<double>& vec) {
    for (const auto& element : vec) { // Используем цикл range-based for
        std::cout << element << ' ';
    }
    std::cout << '\n';
}

// Решение системы линейных уравнений методом сопряженных градиентов
void solve(
    vector<double> &w,
    vector<double> &a,
    vector<double> &b,
    vector<double> &f,
    double eps,
    int maxit,
    int &n,
    double h1,
    double h2,
    int x, int y
) {

    double rho = 0.0, beta = 0.0, alpha = 0.0;

    std::vector<double> nw(a.size());
    std::vector<double> r(a.size());
    std::vector<double> nr(a.size());
    std::vector<double> buf(a.size());
    std::vector<double> z(a.size());
    std::vector<double> nz(a.size());
    std::vector<double> p(a.size());
    double check;

    A(r, w, a, b, x, y);
    for (int i=0; i<r.size(); i++)
        r[i] = f[i] - r[i];

    D1(z, r, a, b, x, y);
    p = z;
    
    ++n;

    do {

        ++n;

        A(buf, p, a, b, x, y);
        alpha = dot(nz, nr, h1, h2)/dot(buf, p, h1, h2);
        axpy(nw, w, alpha, p);
        axpy(nr, r, -alpha, buf);

        D1(nz, nr, a, b, x, y);

        beta = dot(nz, nr, h1, h2)/dot(z, r, h1, h2);

        axpy(p, nz, beta, p);

        z = nz;
        r = nr;
        rho = norm(r);
        w = nw;

    } while ((rho > eps) && (n < maxit));

}

void solve2(
    vector<double> &w,
    vector<double> &a,
    vector<double> &b,
    vector<double> &f,
    double eps,
    int max_it,
    int &n,
    double h1,
    double h2,
    int x, int y
) {
    double rho = 0.0, beta = 0.0, alpha = 0.0;
    double old_dot_zr = 0.0; // Для хранения (z^{k-1}, r^{k-1})

    vector<double> r(f.size()), z(f.size()), p(f.size()), Ap(f.size()), temp(f.size());

    vfill(w.size(), w, 0.0); // Начальное приближение w=0

    A(r, w, a, b, x, y);
    for (size_t i = 0; i < r.size(); ++i) {
        r[i] = f[i] - r[i];
    }
    rho = norm(r);
    if (rho < eps) {
        n = 0;
        return;
    }

    D1(z, r, a, b, x, y);
    vcopy(p, z);
    old_dot_zr = dot(z, r, h1, h2);

    n = 0;
    while (n < max_it && rho >= eps) {
        A(Ap, p, a, b, x, y); // Вычисляем Ap для текущего p
        alpha = old_dot_zr / dot(Ap, p, h1, h2);
        axpy(w, w, alpha, p);
        axpy(r, r, -alpha, Ap); // temp = r - alpha * Ap

        rho = norm(r); // Проверяем норму невязки
        if (rho < eps) break;

        D1(z, r, a, b, x, y);
        double new_dot_zr = dot(z, r, h1, h2);
        beta = new_dot_zr / old_dot_zr;
        axpy(p, z, beta, p); // temp = z + beta * p

        old_dot_zr = new_dot_zr;
        ++n;
    }
}