From 9bd7b57cb5ca9e5cf21842b0d8c396660f3f77c3 Mon Sep 17 00:00:00 2001 From: VarenytsiaMykhailo Date: Sun, 14 Apr 2024 19:27:31 +0300 Subject: [PATCH] =?UTF-8?q?Add=20new=20numerical=20methods=20to=20solve=20?= =?UTF-8?q?linear=20system:=201.=20Thomas=20method=20(=D0=BF=D1=80=D0=BE?= =?UTF-8?q?=D0=B3=D0=BE=D0=BD=D0=BA=D0=B8)=202.=20Jacobi=20method=20(itera?= =?UTF-8?q?tive)=203.=20Seidel=20method=20(iterative,=20improvement=20vers?= =?UTF-8?q?ion=20of=20the=20Jacobi=20method)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../functions/machineEpsilonPrecision.kt | 22 ++ .../jacobimethod/JacobiMethod.kt | 110 ++++++ .../seidelmethod/SeidelMethod.kt | 122 ++++++ .../thomasmethod/ThomasMethod.kt | 115 ++++++ .../jacobimethod/JacobiMethodTest.kt | 350 ++++++++++++++++++ .../seidelmethod/SeidelMethodTest.kt | 350 ++++++++++++++++++ .../thomasmethod/ThomasMethodTest.kt | 178 +++++++++ 7 files changed, 1247 insertions(+) create mode 100644 kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/machineEpsilonPrecision.kt create mode 100644 kmath-functions/src/commonMain/kotlin/space/kscience/kmath/linearsystemsolving/jacobimethod/JacobiMethod.kt create mode 100644 kmath-functions/src/commonMain/kotlin/space/kscience/kmath/linearsystemsolving/seidelmethod/SeidelMethod.kt create mode 100644 kmath-functions/src/commonMain/kotlin/space/kscience/kmath/linearsystemsolving/thomasmethod/ThomasMethod.kt create mode 100644 kmath-functions/src/commonTest/kotlin/space/kscience/kmath/linearsystemsolving/jacobimethod/JacobiMethodTest.kt create mode 100644 kmath-functions/src/commonTest/kotlin/space/kscience/kmath/linearsystemsolving/seidelmethod/SeidelMethodTest.kt create mode 100644 kmath-functions/src/commonTest/kotlin/space/kscience/kmath/linearsystemsolving/thomasmethod/ThomasMethodTest.kt diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/machineEpsilonPrecision.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/machineEpsilonPrecision.kt new file mode 100644 index 000000000..e39d4eddc --- /dev/null +++ b/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 + } diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/linearsystemsolving/jacobimethod/JacobiMethod.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/linearsystemsolving/jacobimethod/JacobiMethod.kt new file mode 100644 index 000000000..6c2638959 --- /dev/null +++ b/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, + B: Point, + initialApproximation: Point = getDefaultInitialApproximation(A, B), + epsilonPrecision: Double = machineEpsilonPrecision, +): Point { + 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 = initialApproximation + var norm: Double + + do { + val xTmp: MutableStructure1D = 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, x2: Point): 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, B: Point): Point = + structureND(ShapeND(A.rowNum)) { (i) -> + B[i] / A[i, i] + }.as1D() diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/linearsystemsolving/seidelmethod/SeidelMethod.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/linearsystemsolving/seidelmethod/SeidelMethod.kt new file mode 100644 index 000000000..e790ace29 --- /dev/null +++ b/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, + B: Point, + initialApproximation: Point? = null, + epsilonPrecision: Double = machineEpsilonPrecision, +): Point { + 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 = initialApproximation?.let { + mutableStructureND(ShapeND(n)) { (i) -> + initialApproximation[i] + }.as1D() + } ?: getDefaultInitialApproximation(A, B) + + var xTmp: Structure1D = 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, x2: Point): 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, B: Point): MutableStructure1D = + mutableStructureND(ShapeND(A.rowNum)) { (i) -> + B[i] / A[i, i] + }.as1D() diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/linearsystemsolving/thomasmethod/ThomasMethod.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/linearsystemsolving/thomasmethod/ThomasMethod.kt new file mode 100644 index 000000000..41b8a0544 --- /dev/null +++ b/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, + B: Point, + isStrictMode: Boolean = false, +): Point { + 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 = + 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 +} diff --git a/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/linearsystemsolving/jacobimethod/JacobiMethodTest.kt b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/linearsystemsolving/jacobimethod/JacobiMethodTest.kt new file mode 100644 index 000000000..28e2ac1f2 --- /dev/null +++ b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/linearsystemsolving/jacobimethod/JacobiMethodTest.kt @@ -0,0 +1,350 @@ +/* + * 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.linear.Point +import space.kscience.kmath.linear.linearSpace +import space.kscience.kmath.linear.matrix +import space.kscience.kmath.linear.vector +import space.kscience.kmath.operations.algebra +import kotlin.test.Test +import kotlin.test.assertEquals +import kotlin.test.assertFails +import kotlin.test.assertTrue + +@UnstableKMathAPI +internal class JacobiMethodTest { + + /** + * Positive test, + * WHEN matrix 'A' is 3x3 with sufficient condition of the convergence, + * without the custom initial approximation input, + * with the custom epsilonPrecision input (machineEpsilonPrecision will be used). + */ + @Test + fun positiveTest1() = + Double.algebra.linearSpace.run { + val matrixA = matrix(3, 3)( + 115.0, -20.0, -75.0, + 15.0, -50.0, -5.0, + 6.0, 2.0, 20.0, + ) + + val vectorB = vector( + 20.0, -40.0, 28.0, + ) + + val result: Point = solveSystemByJacobiMethod( + A = matrixA, + B = vectorB, + epsilonPrecision = 0.001, + ) + + assertEquals(3, result.size) + val absoluteTolerance = 0.001 + assertEquals(1.0, result[0], absoluteTolerance) + assertEquals(1.0, result[1], absoluteTolerance) + assertEquals(1.0, result[2], absoluteTolerance) + + val result2 = matrixA.dot(result) + + for (i in 0 until result2.size) { + assertEquals(vectorB[i], result2[i], 0.05) + } + } + + + /** + * Positive test, + * WHEN matrix 'A' is 4x4 with sufficient condition of the convergence, + * without the custom initial approximation input, + * without the custom epsilonPrecision input (machineEpsilonPrecision will be used). + */ + @Test + fun positiveTest2() = + Double.algebra.linearSpace.run { + val matrixA = matrix(4, 4)( + 20.9, 1.2, 2.1, 0.9, + 1.2, 21.2, 1.5, 2.5, + 2.1, 1.5, 19.8, 1.3, + 0.9, 2.5, 1.3, 32.1, + ) + + val vectorB = vector( + 21.70, 27.46, 28.76, 49.72, + ) + + val result: Point = solveSystemByJacobiMethod( + A = matrixA, + B = vectorB, + ) + + assertEquals(4, result.size) + val absoluteTolerance = 0.00000000000001 + assertEquals(0.8, result[0], absoluteTolerance) + assertEquals(1.0, result[1], absoluteTolerance) + assertEquals(1.2, result[2], absoluteTolerance) + assertEquals(1.4, result[3], absoluteTolerance) + + val result2 = matrixA.dot(result) + for (i in 0 until result2.size) { + assertEquals(vectorB[i], result2[i], absoluteTolerance) + } + } + + /** + * Positive test, + * WHEN matrix 'A' is 4x4 with sufficient condition of the convergence, + * with the custom initial approximation input, + * without the custom epsilonPrecision input (machineEpsilonPrecision will be used). + */ + @Test + fun positiveTest3() = + Double.algebra.linearSpace.run { + val matrixA = matrix(4, 4)( + 20.9, 1.2, 2.1, 0.9, + 1.2, 21.2, 1.5, 2.5, + 2.1, 1.5, 19.8, 1.3, + 0.9, 2.5, 1.3, 32.1, + ) + + val vectorB = vector( + 21.70, 27.46, 28.76, 49.72, + ) + + val result: Point = solveSystemByJacobiMethod( + A = matrixA, + B = vectorB, + initialApproximation = vector(1.0, 1.0, 1.0, 1.0), + ) + + assertEquals(4, result.size) + val absoluteTolerance = 0.00000000000001 + assertEquals(0.8, result[0], absoluteTolerance) + assertEquals(1.0, result[1], absoluteTolerance) + assertEquals(1.2, result[2], absoluteTolerance) + assertEquals(1.4, result[3], absoluteTolerance) + + val result2 = matrixA.dot(result) + for (i in 0 until result2.size) { + assertEquals(vectorB[i], result2[i], absoluteTolerance) + } + } + + /** + * Positive test, + * WHEN matrix 'A' is 4x4 with sufficient condition of the convergence, + * without the custom initial approximation input, + * with the custom epsilonPrecision input (machineEpsilonPrecision will be used). + */ + @Test + fun positiveTest4() = + Double.algebra.linearSpace.run { + val matrixA = matrix(4, 4)( + 20.9, 1.2, 2.1, 0.9, + 1.2, 21.2, 1.5, 2.5, + 2.1, 1.5, 19.8, 1.3, + 0.9, 2.5, 1.3, 32.1, + ) + + val vectorB = vector( + 21.70, 27.46, 28.76, 49.72, + ) + + val result: Point = solveSystemByJacobiMethod( + A = matrixA, + B = vectorB, + epsilonPrecision = 0.001, + ) + + assertEquals(4, result.size) + val absoluteTolerance = 0.001 + assertEquals(0.8, result[0], absoluteTolerance) + assertEquals(1.0, result[1], absoluteTolerance) + assertEquals(1.2, result[2], absoluteTolerance) + assertEquals(1.4, result[3], absoluteTolerance) + + val result2 = matrixA.dot(result) + for (i in 0 until result2.size) { + assertEquals(vectorB[i], result2[i], absoluteTolerance) + } + } + + /** + * Positive test, + * WHEN matrix 'A' is 4x4 with sufficient condition of the convergence, + * with the custom initial approximation input, + * with the custom epsilonPrecision input (machineEpsilonPrecision will be used). + */ + @Test + fun positiveTest5() = + Double.algebra.linearSpace.run { + val matrixA = matrix(4, 4)( + 20.9, 1.2, 2.1, 0.9, + 1.2, 21.2, 1.5, 2.5, + 2.1, 1.5, 19.8, 1.3, + 0.9, 2.5, 1.3, 32.1, + ) + + val vectorB = vector( + 21.70, 27.46, 28.76, 49.72, + ) + + val result: Point = solveSystemByJacobiMethod( + A = matrixA, + B = vectorB, + initialApproximation = vector(1.0, 1.0, 1.0, 1.0), + epsilonPrecision = 0.001, + ) + + assertEquals(4, result.size) + val absoluteTolerance = 0.001 + assertEquals(0.8, result[0], absoluteTolerance) + assertEquals(1.0, result[1], absoluteTolerance) + assertEquals(1.2, result[2], absoluteTolerance) + assertEquals(1.4, result[3], absoluteTolerance) + + val result2 = matrixA.dot(result) + for (i in 0 until result2.size) { + assertEquals(vectorB[i], result2[i], absoluteTolerance) + } + } + + /** + * Exception test, + * WHEN matrix 'A' is not square matrix. + */ + @Test + fun exceptionTest1() = + Double.algebra.linearSpace.run { + val matrixA = matrix(3, 4)( + 4.0, 1.0, 0.0, 0.0, + 1.0, 4.0, 1.0, 0.0, + 0.0, 1.0, 4.0, 1.0, + ) + + val vectorB = vector( + 5.0, 6.0, 6.0, 5.0, + ) + + assertFails { + solveSystemByJacobiMethod(matrixA, vectorB) + } + assertTrue(true) // Needs, because single assertFails is not supporting: + // Execution failed for task ':kmath-functions:jvmTest'. + // > No tests found for given includes: [space.kscience.kmath.linearsystemsolving.thomasmethod.ThomasMethodTest.exceptionTestWhenAisNotSquareMatrix](--tests filter) + } + + /** + * Exception test, + * WHEN the dimension of the matrix 'A' is bigger than the dimension of the vector 'B'. + */ + @Test + fun exceptionTest2() = + Double.algebra.linearSpace.run { + val matrixA = matrix(4, 4)( + 4.0, 1.0, 0.0, 0.0, + 1.0, 4.0, 1.0, 0.0, + 0.0, 1.0, 4.0, 1.0, + 0.0, 0.0, 1.0, 4.0, + ) + + val vectorB = vector( + 5.0, 6.0, 6.0, + ) + + assertFails { + solveSystemByJacobiMethod(matrixA, vectorB) + } + assertTrue(true) // Needs, because single assertFails is not supporting: + // Execution failed for task ':kmath-functions:jvmTest'. + // > No tests found for given includes: [space.kscience.kmath.linearsystemsolving.thomasmethod.ThomasMethodTest.exceptionTestWhenAisNotSquareMatrix](--tests filter) + } + + /** + * Exception test, + * WHEN the dimension of the matrix 'A' is smaller than the dimension of the vector 'B'. + */ + @Test + fun exceptionTest3() = + Double.algebra.linearSpace.run { + val matrixA = matrix(4, 4)( + 4.0, 1.0, 0.0, 0.0, + 1.0, 4.0, 1.0, 0.0, + 0.0, 1.0, 4.0, 1.0, + 0.0, 0.0, 1.0, 4.0, + ) + + val vectorB = vector( + 5.0, 6.0, 6.0, 5.0, 7.0, + ) + + assertFails { + solveSystemByJacobiMethod(matrixA, vectorB) + } + assertTrue(true) // Needs, because single assertFails is not supporting: + // Execution failed for task ':kmath-functions:jvmTest'. + // > No tests found for given includes: [space.kscience.kmath.linearsystemsolving.thomasmethod.ThomasMethodTest.exceptionTestWhenAisNotSquareMatrix](--tests filter) + } + + /** + * Exception test, + * WHEN the size of the vector 'B' does not match the size of the vector 'initialApproximation'. + */ + @Test + fun exceptionTest4() = + Double.algebra.linearSpace.run { + val matrixA = matrix(4, 4)( + 20.9, 1.2, 2.1, 0.9, + 1.2, 21.2, 1.5, 2.5, + 2.1, 1.5, 19.8, 1.3, + 0.9, 2.5, 1.3, 32.1, + ) + + val vectorB = vector( + 21.70, 27.46, 28.76, 49.72, + ) + + assertFails { + solveSystemByJacobiMethod( + A = matrixA, + B = vectorB, + initialApproximation = vector(1.0, 1.0, 1.0) + ) + } + assertTrue(true) // Needs, because single assertFails is not supporting: + // Execution failed for task ':kmath-functions:jvmTest'. + // > No tests found for given includes: [space.kscience.kmath.linearsystemsolving.thomasmethod.ThomasMethodTest.exceptionTestWhenAisNotSquareMatrix](--tests filter) + } + + /** + * Exception test, + * WHEN matrix 'A' is 4x4 without the sufficient condition for the convergence of the Jacobi method: + * there is no diagonal dominance of the matrix. + */ + @Test + fun exceptionTest5() = + Double.algebra.linearSpace.run { + val matrixA = matrix(4, 4)( + 20.9, 1.2, 2.1, 0.9, + 1.2, 21.2, 1.5, 2.5, + 2.1, 1.5, 0.00000000000000000, 1.3, + 0.9, 2.5, 1.3, 32.1, + ) + + val vectorB = vector( + 21.70, 27.46, 28.76, 49.72, + ) + + assertFails { + solveSystemByJacobiMethod(matrixA, vectorB) + } + assertTrue(true) // Needs, because single assertFails is not supporting: + // Execution failed for task ':kmath-functions:jvmTest'. + // > No tests found for given includes: [space.kscience.kmath.linearsystemsolving.thomasmethod.ThomasMethodTest.exceptionTestWhenAisNotSquareMatrix](--tests filter) + } +} \ No newline at end of file diff --git a/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/linearsystemsolving/seidelmethod/SeidelMethodTest.kt b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/linearsystemsolving/seidelmethod/SeidelMethodTest.kt new file mode 100644 index 000000000..6c59a62fd --- /dev/null +++ b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/linearsystemsolving/seidelmethod/SeidelMethodTest.kt @@ -0,0 +1,350 @@ +/* + * 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.linear.Point +import space.kscience.kmath.linear.linearSpace +import space.kscience.kmath.linear.matrix +import space.kscience.kmath.linear.vector +import space.kscience.kmath.operations.algebra +import kotlin.test.Test +import kotlin.test.assertEquals +import kotlin.test.assertFails +import kotlin.test.assertTrue + +@UnstableKMathAPI +internal class SeidelMethodTest { + + /** + * Positive test, + * WHEN matrix 'A' is 3x3 with sufficient condition of the convergence, + * without the custom initial approximation input, + * with the custom epsilonPrecision input (machineEpsilonPrecision will be used). + */ + @Test + fun positiveTest1() = + Double.algebra.linearSpace.run { + val matrixA = matrix(3, 3)( + 115.0, -20.0, -75.0, + 15.0, -50.0, -5.0, + 6.0, 2.0, 20.0, + ) + + val vectorB = vector( + 20.0, -40.0, 28.0, + ) + + val result: Point = solveSystemBySeidelMethod( + A = matrixA, + B = vectorB, + epsilonPrecision = 0.001, + ) + + assertEquals(3, result.size) + val absoluteTolerance = 0.0001 + assertEquals(1.0, result[0], absoluteTolerance) + assertEquals(1.0, result[1], absoluteTolerance) + assertEquals(1.0, result[2], absoluteTolerance) + + val result2 = matrixA.dot(result) + + for (i in 0 until result2.size) { + assertEquals(vectorB[i], result2[i], 0.05) + } + } + + + /** + * Positive test, + * WHEN matrix 'A' is 4x4 with sufficient condition of the convergence, + * without the custom initial approximation input, + * without the custom epsilonPrecision input (machineEpsilonPrecision will be used). + */ + @Test + fun positiveTest2() = + Double.algebra.linearSpace.run { + val matrixA = matrix(4, 4)( + 20.9, 1.2, 2.1, 0.9, + 1.2, 21.2, 1.5, 2.5, + 2.1, 1.5, 19.8, 1.3, + 0.9, 2.5, 1.3, 32.1, + ) + + val vectorB = vector( + 21.70, 27.46, 28.76, 49.72, + ) + + val result: Point = solveSystemBySeidelMethod( + A = matrixA, + B = vectorB, + ) + + assertEquals(4, result.size) + val absoluteTolerance = 0.00000000000001 + assertEquals(0.8, result[0], absoluteTolerance) + assertEquals(1.0, result[1], absoluteTolerance) + assertEquals(1.2, result[2], absoluteTolerance) + assertEquals(1.4, result[3], absoluteTolerance) + + val result2 = matrixA.dot(result) + for (i in 0 until result2.size) { + assertEquals(vectorB[i], result2[i], absoluteTolerance) + } + } + + /** + * Positive test, + * WHEN matrix 'A' is 4x4 with sufficient condition of the convergence, + * with the custom initial approximation input, + * without the custom epsilonPrecision input (machineEpsilonPrecision will be used). + */ + @Test + fun positiveTest3() = + Double.algebra.linearSpace.run { + val matrixA = matrix(4, 4)( + 20.9, 1.2, 2.1, 0.9, + 1.2, 21.2, 1.5, 2.5, + 2.1, 1.5, 19.8, 1.3, + 0.9, 2.5, 1.3, 32.1, + ) + + val vectorB = vector( + 21.70, 27.46, 28.76, 49.72, + ) + + val result: Point = solveSystemBySeidelMethod( + A = matrixA, + B = vectorB, + initialApproximation = vector(1.0, 1.0, 1.0, 1.0), + ) + + assertEquals(4, result.size) + val absoluteTolerance = 0.00000000000001 + assertEquals(0.8, result[0], absoluteTolerance) + assertEquals(1.0, result[1], absoluteTolerance) + assertEquals(1.2, result[2], absoluteTolerance) + assertEquals(1.4, result[3], absoluteTolerance) + + val result2 = matrixA.dot(result) + for (i in 0 until result2.size) { + assertEquals(vectorB[i], result2[i], absoluteTolerance) + } + } + + /** + * Positive test, + * WHEN matrix 'A' is 4x4 with sufficient condition of the convergence, + * without the custom initial approximation input, + * with the custom epsilonPrecision input (machineEpsilonPrecision will be used). + */ + @Test + fun positiveTest4() = + Double.algebra.linearSpace.run { + val matrixA = matrix(4, 4)( + 20.9, 1.2, 2.1, 0.9, + 1.2, 21.2, 1.5, 2.5, + 2.1, 1.5, 19.8, 1.3, + 0.9, 2.5, 1.3, 32.1, + ) + + val vectorB = vector( + 21.70, 27.46, 28.76, 49.72, + ) + + val result: Point = solveSystemBySeidelMethod( + A = matrixA, + B = vectorB, + epsilonPrecision = 0.001, + ) + + assertEquals(4, result.size) + val absoluteTolerance = 0.0001 + assertEquals(0.8, result[0], absoluteTolerance) + assertEquals(1.0, result[1], absoluteTolerance) + assertEquals(1.2, result[2], absoluteTolerance) + assertEquals(1.4, result[3], absoluteTolerance) + + val result2 = matrixA.dot(result) + for (i in 0 until result2.size) { + assertEquals(vectorB[i], result2[i], absoluteTolerance) + } + } + + /** + * Positive test, + * WHEN matrix 'A' is 4x4 with sufficient condition of the convergence, + * with the custom initial approximation input, + * with the custom epsilonPrecision input (machineEpsilonPrecision will be used). + */ + @Test + fun positiveTest5() = + Double.algebra.linearSpace.run { + val matrixA = matrix(4, 4)( + 20.9, 1.2, 2.1, 0.9, + 1.2, 21.2, 1.5, 2.5, + 2.1, 1.5, 19.8, 1.3, + 0.9, 2.5, 1.3, 32.1, + ) + + val vectorB = vector( + 21.70, 27.46, 28.76, 49.72, + ) + + val result: Point = solveSystemBySeidelMethod( + A = matrixA, + B = vectorB, + initialApproximation = vector(1.0, 1.0, 1.0, 1.0), + epsilonPrecision = 0.001, + ) + + assertEquals(4, result.size) + val absoluteTolerance = 0.0001 + assertEquals(0.8, result[0], absoluteTolerance) + assertEquals(1.0, result[1], absoluteTolerance) + assertEquals(1.2, result[2], absoluteTolerance) + assertEquals(1.4, result[3], absoluteTolerance) + + val result2 = matrixA.dot(result) + for (i in 0 until result2.size) { + assertEquals(vectorB[i], result2[i], absoluteTolerance) + } + } + + /** + * Exception test, + * WHEN matrix 'A' is not square matrix. + */ + @Test + fun exceptionTest1() = + Double.algebra.linearSpace.run { + val matrixA = matrix(3, 4)( + 4.0, 1.0, 0.0, 0.0, + 1.0, 4.0, 1.0, 0.0, + 0.0, 1.0, 4.0, 1.0, + ) + + val vectorB = vector( + 5.0, 6.0, 6.0, 5.0, + ) + + assertFails { + solveSystemBySeidelMethod(matrixA, vectorB) + } + assertTrue(true) // Needs, because single assertFails is not supporting: + // Execution failed for task ':kmath-functions:jvmTest'. + // > No tests found for given includes: [space.kscience.kmath.linearsystemsolving.thomasmethod.ThomasMethodTest.exceptionTestWhenAisNotSquareMatrix](--tests filter) + } + + /** + * Exception test, + * WHEN the dimension of the matrix 'A' is bigger than the dimension of the vector 'B'. + */ + @Test + fun exceptionTest2() = + Double.algebra.linearSpace.run { + val matrixA = matrix(4, 4)( + 4.0, 1.0, 0.0, 0.0, + 1.0, 4.0, 1.0, 0.0, + 0.0, 1.0, 4.0, 1.0, + 0.0, 0.0, 1.0, 4.0, + ) + + val vectorB = vector( + 5.0, 6.0, 6.0, + ) + + assertFails { + solveSystemBySeidelMethod(matrixA, vectorB) + } + assertTrue(true) // Needs, because single assertFails is not supporting: + // Execution failed for task ':kmath-functions:jvmTest'. + // > No tests found for given includes: [space.kscience.kmath.linearsystemsolving.thomasmethod.ThomasMethodTest.exceptionTestWhenAisNotSquareMatrix](--tests filter) + } + + /** + * Exception test, + * WHEN the dimension of the matrix 'A' is smaller than the dimension of the vector 'B'. + */ + @Test + fun exceptionTest3() = + Double.algebra.linearSpace.run { + val matrixA = matrix(4, 4)( + 4.0, 1.0, 0.0, 0.0, + 1.0, 4.0, 1.0, 0.0, + 0.0, 1.0, 4.0, 1.0, + 0.0, 0.0, 1.0, 4.0, + ) + + val vectorB = vector( + 5.0, 6.0, 6.0, 5.0, 7.0, + ) + + assertFails { + solveSystemBySeidelMethod(matrixA, vectorB) + } + assertTrue(true) // Needs, because single assertFails is not supporting: + // Execution failed for task ':kmath-functions:jvmTest'. + // > No tests found for given includes: [space.kscience.kmath.linearsystemsolving.thomasmethod.ThomasMethodTest.exceptionTestWhenAisNotSquareMatrix](--tests filter) + } + + /** + * Exception test, + * WHEN the size of the vector 'B' does not match the size of the vector 'initialApproximation'. + */ + @Test + fun exceptionTest4() = + Double.algebra.linearSpace.run { + val matrixA = matrix(4, 4)( + 20.9, 1.2, 2.1, 0.9, + 1.2, 21.2, 1.5, 2.5, + 2.1, 1.5, 19.8, 1.3, + 0.9, 2.5, 1.3, 32.1, + ) + + val vectorB = vector( + 21.70, 27.46, 28.76, 49.72, + ) + + assertFails { + solveSystemBySeidelMethod( + A = matrixA, + B = vectorB, + initialApproximation = vector(1.0, 1.0, 1.0) + ) + } + assertTrue(true) // Needs, because single assertFails is not supporting: + // Execution failed for task ':kmath-functions:jvmTest'. + // > No tests found for given includes: [space.kscience.kmath.linearsystemsolving.thomasmethod.ThomasMethodTest.exceptionTestWhenAisNotSquareMatrix](--tests filter) + } + + /** + * Exception test, + * WHEN matrix 'A' is 4x4 without the sufficient condition for the convergence of the Seidel method: + * there is no diagonal dominance of the matrix. + */ + @Test + fun exceptionTest5() = + Double.algebra.linearSpace.run { + val matrixA = matrix(4, 4)( + 20.9, 1.2, 2.1, 0.9, + 1.2, 21.2, 1.5, 2.5, + 2.1, 1.5, 0.00000000000000000, 1.3, + 0.9, 2.5, 1.3, 32.1, + ) + + val vectorB = vector( + 21.70, 27.46, 28.76, 49.72, + ) + + assertFails { + solveSystemBySeidelMethod(matrixA, vectorB) + } + assertTrue(true) // Needs, because single assertFails is not supporting: + // Execution failed for task ':kmath-functions:jvmTest'. + // > No tests found for given includes: [space.kscience.kmath.linearsystemsolving.thomasmethod.ThomasMethodTest.exceptionTestWhenAisNotSquareMatrix](--tests filter) + } +} diff --git a/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/linearsystemsolving/thomasmethod/ThomasMethodTest.kt b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/linearsystemsolving/thomasmethod/ThomasMethodTest.kt new file mode 100644 index 000000000..5c408fa3a --- /dev/null +++ b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/linearsystemsolving/thomasmethod/ThomasMethodTest.kt @@ -0,0 +1,178 @@ +/* + * 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.Point +import space.kscience.kmath.linear.linearSpace +import space.kscience.kmath.linear.matrix +import space.kscience.kmath.linear.vector +import space.kscience.kmath.operations.algebra +import space.kscience.kmath.structures.toDoubleArray +import kotlin.test.* + +@UnstableKMathAPI +internal class ThomasMethodTest { + + /** + * Positive test, + * WHEN matrix 'A' is 3x3 tridiagonal matrix. + */ + @Test + fun positiveTest1() = + Double.algebra.linearSpace.run { + val matrixA = matrix(3, 3)( + 2.0, -1.0, 0.0, + 5.0, 4.0, 2.0, + 0.0, 1.0, -3.0, + ) + + val vectorB = vector( + 3.0, 6.0, 2.0, + ) + + val result: Point = solveSystemByThomasMethod(matrixA, vectorB) + + assertEquals(3, result.size) + assertTrue(result[0].toString().startsWith("1.4883")) + assertTrue(result[1].toString().startsWith("-0.0232")) + assertTrue(result[2].toString().startsWith("-0.6744")) + + assertContentEquals(vectorB.toDoubleArray(), matrixA.dot(result).array) + } + + /** + * Positive test, + * WHEN matrix 'A' is 4x4 tridiagonal matrix. + */ + @Test + fun positiveTest2() = + Double.algebra.linearSpace.run { + val matrixA = matrix(4, 4)( + 4.0, 1.0, 0.0, 0.0, + 1.0, 4.0, 1.0, 0.0, + 0.0, 1.0, 4.0, 1.0, + 0.0, 0.0, 1.0, 4.0, + ) + + val vectorB = vector( + 5.0, 6.0, 6.0, 5.0, + ) + + val result: Point = solveSystemByThomasMethod(matrixA, vectorB) + + assertEquals(4, result.size) + assertEquals(1.0, result[0]) + assertEquals(1.0, result[1]) + assertEquals(1.0, result[2]) + assertEquals(1.0, result[3]) + + assertContentEquals(vectorB.toDoubleArray(), matrixA.dot(result).array) + } + + /** + * Exception test, + * WHEN matrix 'A' is not square matrix. + */ + @Test + fun exceptionTest1() = + Double.algebra.linearSpace.run { + val matrixA = matrix(3, 4)( + 4.0, 1.0, 0.0, 0.0, + 1.0, 4.0, 1.0, 0.0, + 0.0, 1.0, 4.0, 1.0, + ) + + val vectorB = vector( + 5.0, 6.0, 6.0, 5.0, + ) + + assertFails { + solveSystemByThomasMethod(matrixA, vectorB) + } + assertTrue(true) // Needs, because single assertFails is not supporting: + // Execution failed for task ':kmath-functions:jvmTest'. + // > No tests found for given includes: [space.kscience.kmath.linearsystemsolving.thomasmethod.ThomasMethodTest.exceptionTestWhenAisNotSquareMatrix](--tests filter) + } + + /** + * Exception test, + * WHEN the dimension of the matrix 'A' is bigger than the dimension of the vector 'B'. + */ + @Test + fun exceptionTest2() = + Double.algebra.linearSpace.run { + val matrixA = matrix(4, 4)( + 4.0, 1.0, 0.0, 0.0, + 1.0, 4.0, 1.0, 0.0, + 0.0, 1.0, 4.0, 1.0, + 0.0, 0.0, 1.0, 4.0, + ) + + val vectorB = vector( + 5.0, 6.0, 6.0, + ) + + assertFails { + solveSystemByThomasMethod(matrixA, vectorB) + } + assertTrue(true) // Needs, because single assertFails is not supporting: + // Execution failed for task ':kmath-functions:jvmTest'. + // > No tests found for given includes: [space.kscience.kmath.linearsystemsolving.thomasmethod.ThomasMethodTest.exceptionTestWhenAisNotSquareMatrix](--tests filter) + } + + /** + * Exception test, + * WHEN the dimension of the matrix 'A' is smaller than the dimension of the vector 'B'. + */ + @Test + fun exceptionTest3() = + Double.algebra.linearSpace.run { + val matrixA = matrix(4, 4)( + 4.0, 1.0, 0.0, 0.0, + 1.0, 4.0, 1.0, 0.0, + 0.0, 1.0, 4.0, 1.0, + 0.0, 0.0, 1.0, 4.0, + ) + + val vectorB = vector( + 5.0, 6.0, 6.0, 5.0, 7.0, + ) + + assertFails { + solveSystemByThomasMethod(matrixA, vectorB) + } + assertTrue(true) // Needs, because single assertFails is not supporting: + // Execution failed for task ':kmath-functions:jvmTest'. + // > No tests found for given includes: [space.kscience.kmath.linearsystemsolving.thomasmethod.ThomasMethodTest.exceptionTestWhenAisNotSquareMatrix](--tests filter) + } + + /** + * Exception test, + * WHEN the matrix 'A' is not tridiagonal matrix. + */ + @Test + fun exceptionTest4() = + Double.algebra.linearSpace.run { + val matrixA = matrix(4, 4)( + 4.0, 1.0, 0.0, 0.0, + 1.0, 4.0, 1.0, 0.0, + 0.0, 1.0, 4.0, 1.0, + 99999999.9999, 0.0, 1.0, 4.0, + ) + + val vectorB = vector( + 5.0, 6.0, 6.0, 5.0, + ) + + assertFails { + solveSystemByThomasMethod(matrixA, vectorB, isStrictMode = true) + } + assertTrue(true) // Needs, because single assertFails is not supporting: + // Execution failed for task ':kmath-functions:jvmTest'. + // > No tests found for given includes: [space.kscience.kmath.linearsystemsolving.thomasmethod.ThomasMethodTest.exceptionTestWhenAisNotSquareMatrix](--tests filter) + } +} \ No newline at end of file