Add new numerical methods to solve linear system:

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/machineEpsilonPrecision.kt

@@ -0,0 +1,22 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.functions
+
+import space.kscience.kmath.UnstableKMathAPI
+
+private var cachedMachineEpsilonPrecision: Double? = null
+
+@UnstableKMathAPI
+public val machineEpsilonPrecision: Double
+    get() = cachedMachineEpsilonPrecision ?: run {
+        var eps = 1.0
+        while (1 + eps > 1) {
+            eps /= 2.0
+        }
+        cachedMachineEpsilonPrecision = eps
+
+        eps
+    }

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/linearsystemsolving/jacobimethod/JacobiMethod.kt

@@ -0,0 +1,110 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.linearsystemsolving.jacobimethod
+
+import space.kscience.kmath.UnstableKMathAPI
+import space.kscience.kmath.functions.machineEpsilonPrecision
+import space.kscience.kmath.linear.Matrix
+import space.kscience.kmath.linear.Point
+import space.kscience.kmath.nd.Floa64FieldOpsND.Companion.mutableStructureND
+import space.kscience.kmath.nd.Floa64FieldOpsND.Companion.structureND
+import space.kscience.kmath.nd.MutableStructure1D
+import space.kscience.kmath.nd.ShapeND
+import space.kscience.kmath.nd.as1D
+import kotlin.math.abs
+
+/**
+ * Jacobi's method implementation.
+ *
+ * In numerical linear algebra, the Jacobi method is an iterative algorithm for determining the solutions of a strictly
+ * diagonally dominant system of linear equations.
+ *
+ * Asymptotic complexity is O(n^3).
+ *
+ * **See Also:** [https://en.wikipedia.org/wiki/Jacobi_method], [https://ru.wikipedia.org/wiki/Метод_Якоби]
+ *
+ * @param [A] is the input matrix of the system.
+ * @param [B] is the input vector of the right side of the system.
+ * @param [initialApproximation] is the input initial approximation vector.
+ * If the user does not pass custom initial approximation vector, then it will be calculated with default algorithm.
+ * @param [epsilonPrecision] is the input precision of the result.
+ * The user can use, for example, 'epsilonPrecision = 0.001' if he need quickly solution with the small solution error.
+ * If the user does not pass [epsilonPrecision],
+ * then will be used default machine precision as the most accurate precision, but with smaller performance.
+ *
+ * @return vector X - solution of the system 'A*X=B'.
+ */
+@UnstableKMathAPI
+public fun solveSystemByJacobiMethod(
+    A: Matrix<Double>,
+    B: Point<Double>,
+    initialApproximation: Point<Double> = getDefaultInitialApproximation(A, B),
+    epsilonPrecision: Double = machineEpsilonPrecision,
+): Point<Double> {
+    val n: Int = A.rowNum
+
+    require(n == A.colNum) {
+        "The number of rows of matrix 'A' must match the number of columns."
+    }
+    require(n == B.size) {
+        "The number of matrix 'A' rows must match the number of vector 'B' elements."
+    }
+    require(B.size == initialApproximation.size) {
+        "The size of vector 'B' must match the size of 'initialApproximation' vector."
+    }
+
+    // Check the sufficient condition of the convergence of the method: must be diagonal dominance in the matrix 'A'
+    for (i in 0 until n) {
+        var sumOfRowWithoutDiagElem = 0.0
+        var diagElemAbs = 0.0
+        for (j in 0 until n) {
+            if (i != j) {
+                sumOfRowWithoutDiagElem += abs(A[i, j])
+            } else {
+                diagElemAbs = abs(A[i, j])
+            }
+        }
+        require(sumOfRowWithoutDiagElem < diagElemAbs) {
+            "The sufficient condition for the convergence of the Jacobi method is not satisfied: there is no diagonal dominance in the matrix 'A'."
+        }
+    }
+
+    var X: Point<Double> = initialApproximation
+    var norm: Double
+
+    do {
+        val xTmp: MutableStructure1D<Double> = mutableStructureND(ShapeND(n)) { 0.0 }.as1D()
+
+        for (i in 0 until n) {
+            var sum = 0.0
+            for (j in 0 until n) {
+                if (j != i) {
+                    sum += A[i, j] * X[j]
+                }
+            }
+            xTmp[i] = (B[i] - sum) / A[i, i]
+        }
+
+        norm = calcNorm(X, xTmp)
+        X = xTmp
+    } while (norm > epsilonPrecision)
+
+    return X
+}
+
+private fun calcNorm(x1: Point<Double>, x2: Point<Double>): Double {
+    var sum = 0.0
+    for (i in 0 until x1.size) {
+        sum += abs(x1[i] - x2[i])
+    }
+
+    return sum
+}
+
+private fun getDefaultInitialApproximation(A: Matrix<Double>, B: Point<Double>): Point<Double> =
+    structureND(ShapeND(A.rowNum)) { (i) ->
+        B[i] / A[i, i]
+    }.as1D()

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/linearsystemsolving/seidelmethod/SeidelMethod.kt

@@ -0,0 +1,122 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.linearsystemsolving.seidelmethod
+
+import space.kscience.kmath.UnstableKMathAPI
+import space.kscience.kmath.functions.machineEpsilonPrecision
+import space.kscience.kmath.linear.Matrix
+import space.kscience.kmath.linear.Point
+import space.kscience.kmath.nd.*
+import space.kscience.kmath.nd.Floa64FieldOpsND.Companion.mutableStructureND
+import space.kscience.kmath.nd.Floa64FieldOpsND.Companion.structureND
+import kotlin.math.abs
+
+/**
+ * Seidel method implementation.
+ *
+ * In numerical linear algebra, the Gauss–Seidel method, also known as the Liebmann method or
+ * the method of successive displacement,
+ * is an iterative method used to solve a system of linear equations and is similar to the 'Jacobi Method'.
+ * The Gauss-Seidel method can be considered as a modification of the Jacobi method, which, as practice shows,
+ * requires approximately half the number of iterations compared to the Jacobi method.
+ * Though it can be applied to any matrix with non-zero elements on the diagonals,
+ * convergence is only guaranteed if the matrix is either strictly diagonally dominant,
+ * or symmetric and positive definite.
+ *
+ * Asymptotic complexity: O(n^3)
+ *
+ * **See Also:** [https://en.wikipedia.org/wiki/Gauss–Seidel_method], [https://ru.wikipedia.org/wiki/Метод_Гаусса_—_Зейделя_решения_системы_линейных_уравнений]
+ *
+ * @param [A] is the input matrix of the system.
+ * @param [B] is the input vector of the right side of the system.
+ * @param [initialApproximation] is the input initial approximation vector.
+ * If the user does not pass custom initial approximation vector, then it will be calculated with default algorithm.
+ * @param [epsilonPrecision] is the input precision of the result.
+ * The user can use, for example, 'epsilonPrecision = 0.001' if he need quickly solution with the small solution error.
+ * If the user does not pass [epsilonPrecision],
+ * then will be used default machine precision as the most accurate precision, but with smaller performance.
+ *
+ * @return vector X - solution of the system 'A*X=B'.
+ */
+@UnstableKMathAPI
+public fun solveSystemBySeidelMethod(
+    A: Matrix<Double>,
+    B: Point<Double>,
+    initialApproximation: Point<Double>? = null,
+    epsilonPrecision: Double = machineEpsilonPrecision,
+): Point<Double> {
+    val n: Int = A.rowNum
+
+    require(n == A.colNum) {
+        "The number of rows of matrix 'A' must match the number of columns."
+    }
+    require(n == B.size) {
+        "The number of matrix 'A' rows must match the number of vector 'B' elements."
+    }
+    initialApproximation?.let {
+        require(B.size == initialApproximation.size) {
+            "The size of vector 'B' must match the size of 'initialApproximation' vector."
+        }
+    }
+
+    // Check the sufficient condition of the convergence of the method: must be diagonal dominance in the matrix 'A'
+    for (i in 0 until n) {
+        var sumOfRowWithoutDiagElem = 0.0
+        var diagElemAbs = 0.0
+        for (j in 0 until n) {
+            if (i != j) {
+                sumOfRowWithoutDiagElem += abs(A[i, j])
+            } else {
+                diagElemAbs = abs(A[i, j])
+            }
+        }
+        require(sumOfRowWithoutDiagElem < diagElemAbs) {
+            "The sufficient condition for the convergence of the Jacobi method is not satisfied: there is no diagonal dominance in the matrix 'A'."
+        }
+    }
+
+    val X: MutableStructure1D<Double> = initialApproximation?.let {
+        mutableStructureND(ShapeND(n)) { (i) ->
+            initialApproximation[i]
+        }.as1D()
+    } ?: getDefaultInitialApproximation(A, B)
+
+    var xTmp: Structure1D<Double> = structureND(ShapeND(n)) { 0.0 }.as1D()
+    var norm: Double
+
+    do {
+        for (i in 0 until n) {
+            var sum = B[i]
+            for (j in 0 until n) {
+                if (j != i) {
+                    sum -= A[i, j] * X[j]
+                }
+            }
+            X[i] = sum / A[i, i]
+        }
+
+        norm = calcNorm(X, xTmp)
+        xTmp = structureND(ShapeND(n)) { (i) ->
+            X[i]
+        }.as1D() // TODO find another way to make copy, for example, support 'clone' method in the library
+    } while (norm > epsilonPrecision)
+
+    return X
+}
+
+private fun calcNorm(x1: Point<Double>, x2: Point<Double>): Double {
+    var sum = 0.0
+    for (i in 0 until x1.size) {
+        sum += abs(x1[i] - x2[i])
+    }
+
+    return sum
+}
+
+private fun getDefaultInitialApproximation(A: Matrix<Double>, B: Point<Double>): MutableStructure1D<Double> =
+    mutableStructureND(ShapeND(A.rowNum)) { (i) ->
+        B[i] / A[i, i]
+    }.as1D()

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/linearsystemsolving/thomasmethod/ThomasMethod.kt

@@ -0,0 +1,115 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.linearsystemsolving.thomasmethod
+
+import space.kscience.kmath.UnstableKMathAPI
+import space.kscience.kmath.linear.Matrix
+import space.kscience.kmath.linear.MutableMatrix
+import space.kscience.kmath.linear.Point
+import space.kscience.kmath.nd.Floa64FieldOpsND.Companion.mutableStructureND
+import space.kscience.kmath.nd.ShapeND
+import space.kscience.kmath.nd.as1D
+import space.kscience.kmath.nd.as2D
+
+/**
+ * Thoma's method (tridiagonal matrix algorithm) implementation.
+ *
+ * In numerical linear algebra, the tridiagonal matrix algorithm, also known as the Thomas algorithm,
+ * is a simplified form of Gaussian elimination that can be used to solve tridiagonal systems of equations.
+ * This method is based on forward and backward sweeps.
+ *
+ * Tridiagonal matrix is a matrix that has nonzero elements on the main diagonal,
+ * the first diagonal below this and the first diagonal above the main diagonal only.
+ *
+ * **See Also:** [https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm], [https://ru.wikipedia.org/wiki/Метод_прогонки]
+ *
+ * @param A is the input matrix of the system.
+ * @param B is the input vector of the right side of the system.
+ * @param isStrictMode is the flag, that says that the method will make validation of input matrix 'A',
+ * which must be tridiagonal.
+ * If [isStrictMode] is true, asymptotic complexity will be O(n^2) and
+ * method will throw an exception if the matrix 'A' is not tridiagonal.
+ * If [isStrictMode] is false (is default value), asymptotic complexity will be O(n) and
+ * solution will be incorrect if the matrix 'A' is not tridiagonal (the responsibility lies with you).
+ *
+ * @return vector X - solution of the system 'A*X=B'.
+ */
+@UnstableKMathAPI
+public fun solveSystemByThomasMethod(
+    A: Matrix<Double>,
+    B: Point<Double>,
+    isStrictMode: Boolean = false,
+): Point<Double> {
+    val n: Int = A.rowNum
+
+    require(n == A.colNum) {
+        "The number of rows of matrix 'A' must match the number of columns."
+    }
+    require(n == B.size) {
+        "The number of matrix 'A' rows must match the number of vector 'B' elements."
+    }
+
+    // Check the condition of the convergence of the Thomas method.
+    // The matrix 'A' must be tridiagonal matrix:
+    // all elements of the matrix must be zero, except for the elements of the main diagonal and adjacent to it.
+    if (isStrictMode) {
+        var matrixAisTridiagonal = true
+        loop@ for (i in 0 until n) {
+            for (j in 0 until n) {
+                if (i == 0 && j != 0 && j != 1 && A[i, j] != 0.0) {
+                    matrixAisTridiagonal = false
+                    break@loop
+                } else if (i == n - 1 && j != n - 2 && j != n - 1 && A[i, j] != 0.0) {
+                    matrixAisTridiagonal = false
+                    break@loop
+                } else if (j != i - 1 && j != i && j != i + 1 && A[i, j] != 0.0) {
+                    matrixAisTridiagonal = false
+                    break@loop
+                }
+            }
+        }
+        require(matrixAisTridiagonal) {
+            "The matrix 'A' must be tridiagonal: all elements must be zero, except elements of the main diagonal and adjacent to it."
+        }
+    }
+
+    // Forward sweep
+    val alphaBettaCoefficients: MutableMatrix<Double> =
+        mutableStructureND(ShapeND(n, 2)) { 0.0 }.as2D()
+    for (i in 0 until n) {
+        val alpha: Double
+        val betta: Double
+        if (i == 0) {
+            val y: Double = A[i, i]
+            alpha = -A[i, i + 1] / y
+            betta = B[i] / y
+        } else {
+            val y: Double = A[i, i] + A[i, i - 1] * alphaBettaCoefficients[i - 1, 0]
+            alpha =
+                if (i != n - 1) { // For the last row, alpha is not needed and should be 0.0
+                    -A[i, i + 1] / y
+                } else {
+                    0.0
+                }
+            betta = (B[i] - A[i, i - 1] * alphaBettaCoefficients[i - 1, 1]) / y
+        }
+        alphaBettaCoefficients[i, 0] = alpha
+        alphaBettaCoefficients[i, 1] = betta
+    }
+
+    // Backward sweep
+    val result = mutableStructureND(ShapeND(n)) { 0.0 }.as1D()
+    for (i in n - 1 downTo 0) {
+        result[i] =
+            if (i == n - 1) {
+                alphaBettaCoefficients[i, 1]
+            } else {
+                alphaBettaCoefficients[i, 0] * result[i + 1] + alphaBettaCoefficients[i, 1]
+            }
+    }
+
+    return result
+}