2024-10-20 18:27:51 +02:00
|
|
|
import pytest
|
2024-10-20 16:30:42 +02:00
|
|
|
import numpy as np
|
|
|
|
|
from scipy.sparse.linalg import cg
|
|
|
|
|
from matrix_generator import MatrixGenerator
|
|
|
|
|
from richardson_method import RichardsonMethod
|
2024-10-25 17:01:19 +02:00
|
|
|
from processing_type import ProcessingType
|
2024-11-11 17:25:36 +01:00
|
|
|
from time_measurement import time_measurement, tests_time
|
2024-10-20 16:30:42 +02:00
|
|
|
|
2024-10-20 19:06:15 +02:00
|
|
|
def calculate_norm_numpy(I, w, A):
|
|
|
|
|
# Calculate the difference between I and w * A
|
|
|
|
|
difference = I - w * A
|
|
|
|
|
|
|
|
|
|
# Calculate the Euclidean norm of the difference
|
|
|
|
|
norm = np.linalg.norm(difference)
|
|
|
|
|
|
|
|
|
|
return norm
|
|
|
|
|
|
|
|
|
|
def calculate_eigenvalues(A):
|
|
|
|
|
# Calculate the eigenvalues of matrix A
|
|
|
|
|
eigenvalues = np.linalg.eigvals(A)
|
|
|
|
|
|
|
|
|
|
# Find the minimum and maximum eigenvalues
|
|
|
|
|
lambda_min = np.min(eigenvalues)
|
|
|
|
|
lambda_max = np.max(eigenvalues)
|
|
|
|
|
|
|
|
|
|
return lambda_min, lambda_max
|
|
|
|
|
|
2024-10-20 19:34:12 +02:00
|
|
|
def calcualte_norm_from_matrix_numpy(A, n, max_iterations):
|
2024-10-20 19:06:15 +02:00
|
|
|
lambda_min, lambda_max = calculate_eigenvalues(A)
|
|
|
|
|
omega = 2 / (lambda_min + lambda_max)
|
|
|
|
|
I = np.eye(n)
|
2024-10-20 19:34:12 +02:00
|
|
|
return calculate_norm_numpy(I, omega, A)
|
|
|
|
|
|
2024-10-20 19:20:38 +02:00
|
|
|
@pytest.mark.parametrize("n", [2, 3, 4, 5, 10, 20, 50, 100])
|
2024-11-11 17:25:36 +01:00
|
|
|
@pytest.mark.parametrize("processing_type", [ProcessingType.SEQUENTIAL, ProcessingType.THREADS, ProcessingType.PROCESSES])
|
|
|
|
|
@pytest.mark.parametrize("matrix_type", ["spd", "nemeth12", "poli3"])
|
|
|
|
|
@time_measurement(tests_time)
|
|
|
|
|
def test_richardson_vs_cg(n: int, processing_type: ProcessingType, matrix_type: str, capsys):
|
|
|
|
|
print("matrix type: ", matrix_type)
|
|
|
|
|
print("matrix size: ", n if matrix_type == "spd" else "fixed")
|
2024-10-20 19:20:38 +02:00
|
|
|
tolerance = 1e-5
|
|
|
|
|
max_iterations=1000
|
2024-11-11 17:25:36 +01:00
|
|
|
if matrix_type in ["nemeth12", "poli3"] and n != 2:
|
|
|
|
|
pytest.skip("Fixed matrix size for nemeth12 and poli3, skipping redundant runs.")
|
|
|
|
|
|
|
|
|
|
if matrix_type == "spd":
|
|
|
|
|
A, b = MatrixGenerator.generate_matrix_and_vector('spd', size=n)
|
|
|
|
|
elif matrix_type == "poli3":
|
|
|
|
|
A, b = MatrixGenerator.generate_matrix_and_vector('poli3')
|
|
|
|
|
elif matrix_type == "nemeth12":
|
|
|
|
|
A, b = MatrixGenerator.generate_matrix_and_vector('nemeth12')
|
|
|
|
|
else:
|
|
|
|
|
raise ValueError("Invalid matrix type specified. Choose 'spd', 'poli3', or 'nemeth12'.")
|
|
|
|
|
|
|
|
|
|
richardson_solver = RichardsonMethod(processing_type, A, b, max_iterations, size=A.shape[0], tol=1e-7)
|
2024-10-31 17:50:21 +01:00
|
|
|
# solution_richardson, info_richardson = richardson_solver.solve()
|
|
|
|
|
|
|
|
|
|
solution_richardson, info_richardson = None, None
|
|
|
|
|
with capsys.disabled():
|
|
|
|
|
solution_richardson, info_richardson = richardson_solver.solve()
|
|
|
|
|
|
|
|
|
|
# Przechwytywanie wyjścia po solve
|
|
|
|
|
captured = capsys.readouterr()
|
|
|
|
|
print("Captured output:", captured.out)
|
2024-10-20 18:27:51 +02:00
|
|
|
|
2024-11-11 17:25:36 +01:00
|
|
|
solution_cg, info = cg(A, b, atol=0.)
|
2024-10-20 18:27:51 +02:00
|
|
|
|
|
|
|
|
if info == 0: # SciPy CG converged
|
2024-10-20 19:34:12 +02:00
|
|
|
assert_scipy_converged(solution_richardson, info_richardson, solution_cg, tolerance, A, b, max_iterations, n)
|
2024-10-20 18:27:51 +02:00
|
|
|
else: # SciPy CG did not converge
|
2024-10-20 19:34:12 +02:00
|
|
|
assert_scipy_not_converged(solution_richardson, info_richardson, A, b)
|
2024-10-20 16:30:42 +02:00
|
|
|
|
2024-10-20 19:34:12 +02:00
|
|
|
def assert_scipy_converged(solution_richardson, info_richardson, solution_cg, tolerance, A, b, max_iterations, n):
|
|
|
|
|
if info_richardson == "Richardson method for those values will NOT converge":
|
2024-10-20 18:27:51 +02:00
|
|
|
print("Richardson did not converge, while SciPy did")
|
2024-10-20 19:34:12 +02:00
|
|
|
numpy_norm = calcualte_norm_from_matrix_numpy(A, n, max_iterations)
|
|
|
|
|
print("Numpy norm: ", numpy_norm, " Richardson norm: ", solution_richardson)
|
2024-10-20 18:27:51 +02:00
|
|
|
assert False, "Richardson did not converge, while SciPy did"
|
|
|
|
|
else:
|
2024-10-20 16:30:42 +02:00
|
|
|
difference = np.linalg.norm(solution_richardson - solution_cg)
|
|
|
|
|
print(f"Difference between Richardson and CG solutions: {difference:.8f}")
|
2024-10-20 18:39:01 +02:00
|
|
|
if difference < tolerance:
|
|
|
|
|
print("Both Richardson and CG converged and calculated correct values.")
|
|
|
|
|
print("Solution CG:\n", solution_cg)
|
|
|
|
|
print("Solution Richardson:\n", solution_richardson)
|
|
|
|
|
else:
|
2024-10-20 18:32:12 +02:00
|
|
|
print("Matrix A:\n", A)
|
|
|
|
|
print("Vector b:\n", b)
|
2024-10-20 18:27:51 +02:00
|
|
|
assert difference < tolerance, f"The solutions are different! Difference: {difference:.8f}"
|
|
|
|
|
|
2024-10-20 19:34:12 +02:00
|
|
|
def assert_scipy_not_converged(solution_richardson, info_richardson, A, b):
|
|
|
|
|
if info_richardson == "Richardson method for those values will NOT converge":
|
2024-10-20 18:27:51 +02:00
|
|
|
print("Richardson and SciPy did not converge")
|
|
|
|
|
else:
|
|
|
|
|
print("Richardson converged while SciPy did not:", solution_richardson)
|
2024-10-20 18:32:12 +02:00
|
|
|
print("Matrix A:\n", A)
|
|
|
|
|
print("Vector b:\n", b)
|
2024-10-20 18:27:51 +02:00
|
|
|
assert False, "Richardson converged while SciPy did not"
|
2024-10-20 18:39:01 +02:00
|
|
|
|
|
|
|
|
if __name__ == "__main__":
|
2024-10-27 20:13:43 +01:00
|
|
|
pytest.main()
|