eigen_sor.cpp

#include <iostream>
#include <Eigen/Core>

/**
 * Solve the linear system Ax=b using the Successive Over Relaxation (SOR) method.
 *
 * @tparam Scalar The scalar type of the matrix and vectors (e.g., float, double).
 * @tparam Rows The number of rows of the matrix and vectors.
 * @tparam Cols The number of columns of the matrix (ignored for vectors).
 * @param A The matrix A of the linear system (Matrix<Scalar, Rows, Cols>).
 * @param b The vector b of the linear system (Vector<Scalar, Rows>).
 * @param x0 The initial guess for the solution (Vector<Scalar, Rows>).
 * @param omega The relaxation parameter (Scalar).
 * @param tol The tolerance for convergence (Scalar).
 * @param maxiter The maximum number of iterations (int).
 *
 * @return The solution to the linear system (Vector<Scalar, Rows>).
 */
template<typename Scalar, auto Rows, auto Cols>
Eigen::Vector<Scalar, Rows> sor(Eigen::Matrix<Scalar, Rows, Cols> A, Eigen::Vector<Scalar, Rows> b, Eigen::Vector<Scalar, Rows> x0, Scalar omega, Scalar tol, int maxiter) {
    
    auto n = b.size();
    Eigen::Vector<Scalar, Rows> x = x0;
    auto iter = 0;
    auto err = tol + 1;
    
    while (err > tol && iter < maxiter) {
        Eigen::Vector<Scalar, Rows> x_old = x;
        
        for (Eigen::Index i = 0; i < n; i++) {
            Scalar sigma = A.row(i).dot(x);
            x(i) = (1 - omega) * x(i) + omega * (b(i) - sigma + A(i, i) * x(i)) / A(i, i);
        }
        
        err = (x - x_old).norm() / x.norm();
        iter++;
    }
    
    return x;
}

int main() {
    Eigen::Matrix<double, 2, 2> A;
    A << 2, 1, 1, 2;
    Eigen::Vector<double, 2> b{5,6};
    Eigen::Vector<double, 2> x0{0,0};
    double omega = 1.2;
    double tol = 1e-6;
    int maxiter = 1000;

    Eigen::Vector<double, 2> x = sor(A, b, x0, omega, tol, maxiter);
    std::cout << "Solution:\n" << x << std::endl;
    
    return 0;
}

https://godbolt.org/z/MGxh1ch84