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

seq/solve.cpp
83 строки2 KB
Нил Ненахов
added seq and omp program
31 окт 2025, 20:17
31 окт 2025, 20:177c054ab
#include "solve.h"
#include "matrix.h"
#include <cmath>
#include <iostream>
#include <mpi.h>
#include <string>
#include <vector>

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(CSRMatrix* matrix, const std::vector<double>& b, double eps, int maxit,
           std::vector<double>& x, int* n, std::vector<double>& z,
           std::vector<double>& p, std::vector<double>& q, CSRMatrix* M) {

    int nc = matrix->rows;
    // Инициализация вектора x нулями
    vfill(nc, x, 0.0);

    std::vector<double> r = b;
    double rho = 0.0, prho = 0.0, beta = 0.0, alpha = 0.0;
    int k = 0;

    std::vector<double> h(r.size());
    double check;

    // Копирование диагонали из matrix в M
    // copy_diag(M, matrix);

    do {
        ++k;

        // Вычисление SpMV: z = M * r
        SpMV(M, r, z);

        // Обновление rho: rho = dot(r, z)
        prho = rho;
        rho = dot(nc, r, z);

        if (k == 1) {
            // Первый шаг: p = z
            vcopy(nc, z, p);
        } else {
            // Вычисление beta: beta = rho / prho
            beta = rho / prho;

            // Обновление p: p = z + beta * p
            axpy(nc, z, p, beta, p);
        }

        // Вычисление SpMV: q = matrix * p
        SpMV(matrix, p, q);

        // Вычисление alpha: alpha = rho / dot(p, q)
        alpha = rho / dot(matrix->IA.size()-1, p, q);

        // Обновление x: x = x + alpha * p
        axpy(nc, x, p, alpha, x);

        // Обновление r: r = r - alpha * q
        axpy(nc, r, q, -alpha, r);

        /*
        //Проверка
        cout << "J(w) = ";
        for (int i = 0; i<r.size(); i++)
        {
            h[i] = r[i] + b[i]; 
        }
        cout << " " << -0.5*dot(r.size(), h, x) << endl;
        */

    } while ((rho > eps * eps) && (k < maxit));

    *n = k;
}