linear_algebra.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <cmath>
#include <cstdint>
#include <stdexcept>
#include <utility>
#include <vector>
 
namespace pulsatrix {
 
[[nodiscard]] inline std::vector<float> SolveLinearSystem(std::vector<std::vector<float>> a,
                                                           std::vector<float> b) {
    const auto n = static_cast<int64_t>(b.size());
    for (int64_t col = 0; col < n; ++col) {
        int64_t pivot_row = col;
        float pivot_magnitude = std::fabs(a[static_cast<size_t>(col)][static_cast<size_t>(col)]);
        for (int64_t row = col + 1; row < n; ++row) {
            float magnitude = std::fabs(a[static_cast<size_t>(row)][static_cast<size_t>(col)]);
            if (magnitude > pivot_magnitude) {
                pivot_magnitude = magnitude;
                pivot_row = row;
            }
        }
        if (pivot_magnitude < 1e-8f) {
            throw std::runtime_error("SolveLinearSystem: near-singular system");
        }
        std::swap(a[static_cast<size_t>(col)], a[static_cast<size_t>(pivot_row)]);
        std::swap(b[static_cast<size_t>(col)], b[static_cast<size_t>(pivot_row)]);
 
        for (int64_t row = col + 1; row < n; ++row) {
            float factor = a[static_cast<size_t>(row)][static_cast<size_t>(col)] /
                           a[static_cast<size_t>(col)][static_cast<size_t>(col)];
            for (int64_t c = col; c < n; ++c) {
                a[static_cast<size_t>(row)][static_cast<size_t>(c)] -=
                    factor * a[static_cast<size_t>(col)][static_cast<size_t>(c)];
            }
            b[static_cast<size_t>(row)] -= factor * b[static_cast<size_t>(col)];
        }
    }
 
    std::vector<float> x(static_cast<size_t>(n), 0.0f);
    for (int64_t row = n - 1; row >= 0; --row) {
        float sum = b[static_cast<size_t>(row)];
        for (int64_t c = row + 1; c < n; ++c) {
            sum -= a[static_cast<size_t>(row)][static_cast<size_t>(c)] * x[static_cast<size_t>(c)];
        }
        x[static_cast<size_t>(row)] = sum / a[static_cast<size_t>(row)][static_cast<size_t>(row)];
    }
    return x;
}
 
}  // namespace pulsatrix