Jacobi Iteration Method  https:en.wikipedia.orgwikiJacobimethod Method to find solution of system of linear equations Jacobi Iteration Method: An iterative algorithm to determine the solutions of strictly diagonally dominant system of linear equations 4x1   x2   x3   2 x1  5x2  2x3  6 x1  2x2  4x3  4 xinit  0.5, 0.5 , 0.5 Examples:  coefficient  np.array4, 1, 1, 1, 5, 2, 1, 2, 4  constant  np.array2, 6, 4  initval  0.5, 0.5, 0.5  iterations  3  jacobiiterationmethodcoefficient, constant, initval, iterations 0.909375, 1.14375, 0.7484375  coefficient  np.array4, 1, 1, 1, 5, 2  constant  np.array2, 6, 4  initval  0.5, 0.5, 0.5  iterations  3  jacobiiterationmethodcoefficient, constant, initval, iterations Traceback most recent call last: ... ValueError: Coefficient matrix dimensions must be nxn but received 2x3  coefficient  np.array4, 1, 1, 1, 5, 2, 1, 2, 4  constant  np.array2, 6  initval  0.5, 0.5, 0.5  iterations  3  jacobiiterationmethod ...     coefficient, constant, initval, iterations ...    doctest: NORMALIZEWHITESPACE Traceback most recent call last: ... ValueError: Coefficient and constant matrices dimensions must be nxn and nx1 but received 3x3 and 2x1  coefficient  np.array4, 1, 1, 1, 5, 2, 1, 2, 4  constant  np.array2, 6, 4  initval  0.5, 0.5  iterations  3  jacobiiterationmethod ...     coefficient, constant, initval, iterations ...    doctest: NORMALIZEWHITESPACE Traceback most recent call last: ... ValueError: Number of initial values must be equal to number of rows in coefficient matrix but received 2 and 3  coefficient  np.array4, 1, 1, 1, 5, 2, 1, 2, 4  constant  np.array2, 6, 4  initval  0.5, 0.5, 0.5  iterations  0  jacobiiterationmethodcoefficient, constant, initval, iterations Traceback most recent call last: ... ValueError: Iterations must be at least 1 Iterates the whole matrix for given number of times for  in rangeiterations: newval   for row in rangerows: temp  0 for col in rangecols: if col  row: denom  tablerowcol elif col  cols  1: val  tablerowcol else: temp  1  tablerowcol  initvalcol temp  temp  val  denom newval.appendtemp initval  newval denominator  a list of values along the diagonal vallast  values of the last column of the table array masks  boolean mask of all strings without diagonal elements array coefficientmatrix nodiagonals  coefficientmatrix array values without diagonal elements Here we get 'icol'  these are the column numbers, for each row without diagonal elements, except for the last column. 'icol' is converted to a twodimensional list 'ind', which will be used to make selections from 'initval' 'arr' array see below. Iterates the whole matrix for given number of times Checks if the given matrix is strictly diagonally dominant  table  np.array4, 1, 1, 2, 1, 5, 2, 6, 1, 2, 4, 4  strictlydiagonallydominanttable True  table  np.array4, 1, 1, 2, 1, 5, 2, 6, 1, 2, 3, 4  strictlydiagonallydominanttable Traceback most recent call last: ... ValueError: Coefficient matrix is not strictly diagonally dominant Test Cases
from __future__ import annotations

import numpy as np
from numpy import float64
from numpy.typing import NDArray


# Method to find solution of system of linear equations
def jacobi_iteration_method(
    coefficient_matrix: NDArray[float64],
    constant_matrix: NDArray[float64],
    init_val: list[float],
    iterations: int,
) -> list[float]:
    """
    Jacobi Iteration Method:
    An iterative algorithm to determine the solutions of strictly diagonally dominant
    system of linear equations

    4x1 +  x2 +  x3 =  2
     x1 + 5x2 + 2x3 = -6
     x1 + 2x2 + 4x3 = -4

    x_init = [0.5, -0.5 , -0.5]

    Examples:

    >>> coefficient = np.array([[4, 1, 1], [1, 5, 2], [1, 2, 4]])
    >>> constant = np.array([[2], [-6], [-4]])
    >>> init_val = [0.5, -0.5, -0.5]
    >>> iterations = 3
    >>> jacobi_iteration_method(coefficient, constant, init_val, iterations)
    [0.909375, -1.14375, -0.7484375]


    >>> coefficient = np.array([[4, 1, 1], [1, 5, 2]])
    >>> constant = np.array([[2], [-6], [-4]])
    >>> init_val = [0.5, -0.5, -0.5]
    >>> iterations = 3
    >>> jacobi_iteration_method(coefficient, constant, init_val, iterations)
    Traceback (most recent call last):
        ...
    ValueError: Coefficient matrix dimensions must be nxn but received 2x3

    >>> coefficient = np.array([[4, 1, 1], [1, 5, 2], [1, 2, 4]])
    >>> constant = np.array([[2], [-6]])
    >>> init_val = [0.5, -0.5, -0.5]
    >>> iterations = 3
    >>> jacobi_iteration_method(
    ...     coefficient, constant, init_val, iterations
    ... )  # doctest: +NORMALIZE_WHITESPACE
    Traceback (most recent call last):
        ...
    ValueError: Coefficient and constant matrices dimensions must be nxn and nx1 but
                received 3x3 and 2x1

    >>> coefficient = np.array([[4, 1, 1], [1, 5, 2], [1, 2, 4]])
    >>> constant = np.array([[2], [-6], [-4]])
    >>> init_val = [0.5, -0.5]
    >>> iterations = 3
    >>> jacobi_iteration_method(
    ...     coefficient, constant, init_val, iterations
    ... )  # doctest: +NORMALIZE_WHITESPACE
    Traceback (most recent call last):
        ...
    ValueError: Number of initial values must be equal to number of rows in coefficient
                matrix but received 2 and 3

    >>> coefficient = np.array([[4, 1, 1], [1, 5, 2], [1, 2, 4]])
    >>> constant = np.array([[2], [-6], [-4]])
    >>> init_val = [0.5, -0.5, -0.5]
    >>> iterations = 0
    >>> jacobi_iteration_method(coefficient, constant, init_val, iterations)
    Traceback (most recent call last):
        ...
    ValueError: Iterations must be at least 1
    """

    rows1, cols1 = coefficient_matrix.shape
    rows2, cols2 = constant_matrix.shape

    if rows1 != cols1:
        msg = f"Coefficient matrix dimensions must be nxn but received {rows1}x{cols1}"
        raise ValueError(msg)

    if cols2 != 1:
        msg = f"Constant matrix must be nx1 but received {rows2}x{cols2}"
        raise ValueError(msg)

    if rows1 != rows2:
        msg = (
            "Coefficient and constant matrices dimensions must be nxn and nx1 but "
            f"received {rows1}x{cols1} and {rows2}x{cols2}"
        )
        raise ValueError(msg)

    if len(init_val) != rows1:
        msg = (
            "Number of initial values must be equal to number of rows in coefficient "
            f"matrix but received {len(init_val)} and {rows1}"
        )
        raise ValueError(msg)

    if iterations <= 0:
        raise ValueError("Iterations must be at least 1")

    table: NDArray[float64] = np.concatenate(
        (coefficient_matrix, constant_matrix), axis=1
    )

    rows, cols = table.shape

    strictly_diagonally_dominant(table)

    """
    # Iterates the whole matrix for given number of times
    for _ in range(iterations):
        new_val = []
        for row in range(rows):
            temp = 0
            for col in range(cols):
                if col == row:
                    denom = table[row][col]
                elif col == cols - 1:
                    val = table[row][col]
                else:
                    temp += (-1) * table[row][col] * init_val[col]
            temp = (temp + val) / denom
            new_val.append(temp)
        init_val = new_val
    """

    # denominator - a list of values along the diagonal
    denominator = np.diag(coefficient_matrix)

    # val_last - values of the last column of the table array
    val_last = table[:, -1]

    # masks - boolean mask of all strings without diagonal
    # elements array coefficient_matrix
    masks = ~np.eye(coefficient_matrix.shape[0], dtype=bool)

    # no_diagonals - coefficient_matrix array values without diagonal elements
    no_diagonals = coefficient_matrix[masks].reshape(-1, rows - 1)

    # Here we get 'i_col' - these are the column numbers, for each row
    # without diagonal elements, except for the last column.
    i_row, i_col = np.where(masks)
    ind = i_col.reshape(-1, rows - 1)

    #'i_col' is converted to a two-dimensional list 'ind', which will be
    # used to make selections from 'init_val' ('arr' array see below).

    # Iterates the whole matrix for given number of times
    for _ in range(iterations):
        arr = np.take(init_val, ind)
        sum_product_rows = np.sum((-1) * no_diagonals * arr, axis=1)
        new_val = (sum_product_rows + val_last) / denominator
        init_val = new_val

    return new_val.tolist()


# Checks if the given matrix is strictly diagonally dominant
def strictly_diagonally_dominant(table: NDArray[float64]) -> bool:
    """
    >>> table = np.array([[4, 1, 1, 2], [1, 5, 2, -6], [1, 2, 4, -4]])
    >>> strictly_diagonally_dominant(table)
    True

    >>> table = np.array([[4, 1, 1, 2], [1, 5, 2, -6], [1, 2, 3, -4]])
    >>> strictly_diagonally_dominant(table)
    Traceback (most recent call last):
        ...
    ValueError: Coefficient matrix is not strictly diagonally dominant
    """

    rows, cols = table.shape

    is_diagonally_dominant = True

    for i in range(rows):
        total = 0
        for j in range(cols - 1):
            if i == j:
                continue
            else:
                total += table[i][j]

        if table[i][i] <= total:
            raise ValueError("Coefficient matrix is not strictly diagonally dominant")

    return is_diagonally_dominant


# Test Cases
if __name__ == "__main__":
    import doctest

    doctest.testmod()
