1 #pragma once
2 // jacobi.hpp: Jacobi iterative method to solve Ax = b
3 //
4 // Copyright (C) 2017-2021 Stillwater Supercomputing, Inc.
5 //
6 // This file is part of the universal numbers project, which is released under an MIT Open Source license.
7 #include <cmath>
8 #include <universal/number/posit/posit_fwd.hpp>
9 #include <universal/blas/matrix.hpp>
10
11 namespace sw::universal::blas {
12
13 // Jacobi: Solution of x in Ax=b using Jacobi Method
14 template<typename Matrix, typename Vector, size_t MAX_ITERATIONS = 100>
Jacobi(const Matrix & A,const Vector & b,Vector & x,typename Matrix::value_type tolerance=typename Matrix::value_type (0.00001))15 size_t Jacobi(const Matrix& A, const Vector& b, Vector& x, typename Matrix::value_type tolerance = typename Matrix::value_type(0.00001)) {
16 using Scalar = typename Matrix::value_type;
17 Scalar residual = Scalar(std::numeric_limits<Scalar>::max());
18 size_t m = num_rows(A);
19 size_t n = num_cols(A);
20 size_t itr = 0;
21 while (residual > tolerance && itr < MAX_ITERATIONS) {
22 Vector x_old = x;
23 for (size_t i = 0; i < m; ++i) {
24 Scalar sigma = 0;
25 for (size_t j = 0; j < n; ++j) {
26 if (i != j) sigma += A(i, j) * x(j);
27 }
28 x(i) = (b(i) - sigma) / A(i, i);
29 }
30 residual = normL1(x_old - x);
31 std::cout << '[' << itr << "] " << std::setw(10) << x << " residual " << residual << std::endl;
32 ++itr;
33 }
34
35 return itr;
36 }
37
38 } // namespace sw::universal::blas
39