From e3b15d3c121bb4ca2bdfaf681ca79bf5231a6384 Mon Sep 17 00:00:00 2001
From: Simon Bjurek <simbj106@student.liu.se>
Date: Wed, 12 Feb 2025 11:17:21 +0100
Subject: [PATCH] ldlt inverse now verified for matrices up to and including
size 3x3
---
.gitignore | 1 +
b_asic/core_operations.py | 17 ++++++++---
b_asic/sfg_generators.py | 29 +++++++++----------
test/test_sfg_generators.py | 56 +++++++++++++++++++++++++++++++++++++
4 files changed, 85 insertions(+), 18 deletions(-)
diff --git a/.gitignore b/.gitignore
index c5b2148a..d240bbd1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -116,6 +116,7 @@ TODO.txt
b_asic/_version.py
docs_sphinx/_build/
docs_sphinx/examples
+docs_sphinx/sg_execution_times.rst
result_images/
.coverage
Digraph.gv
diff --git a/b_asic/core_operations.py b/b_asic/core_operations.py
index 70d43370..67d6b07f 100644
--- a/b_asic/core_operations.py
+++ b/b_asic/core_operations.py
@@ -1123,7 +1123,7 @@ class MADS(AbstractOperation):
def __init__(
self,
- is_add: Optional[bool] = False,
+ is_add: Optional[bool] = True,
src0: Optional[SignalSourceProvider] = None,
src1: Optional[SignalSourceProvider] = None,
src2: Optional[SignalSourceProvider] = None,
@@ -1142,15 +1142,24 @@ class MADS(AbstractOperation):
latency_offsets=latency_offsets,
execution_time=execution_time,
)
- self._is_add = is_add
- # self.set_param("is_add", is_add)
+ self.set_param("is_add", is_add)
@classmethod
def type_name(cls) -> TypeName:
return TypeName("mads")
def evaluate(self, a, b, c):
- return a + b * c if self._is_add else a - b * c
+ return a + b * c if self.is_add else a - b * c
+
+ @property
+ def is_add(self) -> bool:
+ """Get if operation is an addition."""
+ return self.param("is_add")
+
+ @is_add.setter
+ def is_add(self, is_add: bool) -> None:
+ """Set if operation is an addition."""
+ self.set_param("is_add", is_add)
@property
def is_linear(self) -> bool:
diff --git a/b_asic/sfg_generators.py b/b_asic/sfg_generators.py
index 0ceeb4f6..a0847de4 100644
--- a/b_asic/sfg_generators.py
+++ b/b_asic/sfg_generators.py
@@ -456,40 +456,41 @@ def ldlt_matrix_inverse(N: int, is_complex: bool) -> SFG:
# R*di*R^T factorization
for i in range(N):
- for k in range(1, i):
- D[i] = MADS(False, D[i], M[k][i], R[k][i])
+ for k in range(i):
+ D[i] = MADS(False, D[i], M[k][i], R[k][i], name="M1")
D_inv[i] = Reciprocal(D[i])
- for j in range(i, N):
+ for j in range(i + 1, N):
R[i][j] = A[i][j]
- for k in range(1, i):
- R[i][j] = MADS(False, R[i][j], M[k][i], R[k][j])
+ for k in range(i):
+ R[i][j] = MADS(False, R[i][j], M[k][i], R[k][j], name="M2")
if is_complex:
M[i][j] = ComplexConjugate(R[i][j])
else:
M[i][j] = R[i][j]
- R[i][j] = MADS(True, Constant(0, name="0"), R[i][j], D_inv[i])
+ R[i][j] = MADS(True, Constant(0, name="0"), R[i][j], D_inv[i], name="M3")
# back substitution
A_inv = [[None for _ in range(N)] for _ in range(N)]
for i in reversed(range(N)):
A_inv[i][i] = D_inv[i]
- for j in reversed(range(i)):
- for k in reversed(range(j, N)):
+ for j in reversed(range(i + 1)):
+ for k in reversed(range(j + 1, N)):
if k == N - 1 and i != j:
- # my custom
- if A_inv[k][i]:
+ A_inv[j][i] = MADS(
+ False, Constant(0, name="0"), R[j][k], A_inv[i][k], name="M4"
+ )
+ else:
+ if A_inv[i][k]:
A_inv[j][i] = MADS(
- False, Constant(0, name="0"), R[j][k], A_inv[k][i]
+ False, A_inv[j][i], R[j][k], A_inv[i][k], name="M5"
)
else:
- A_inv[j][i] = Constant(0, name="0")
- else:
- A_inv[j][i] = MADS(False, A_inv[j][i], R[j][k], A_inv[k][i])
+ A_inv[j][i] = MADS(False, A_inv[j][i], R[j][k], A_inv[k][i])
outputs = []
for i in range(N):
diff --git a/test/test_sfg_generators.py b/test/test_sfg_generators.py
index 7064658c..e3f56ec0 100644
--- a/test/test_sfg_generators.py
+++ b/test/test_sfg_generators.py
@@ -12,6 +12,7 @@ from b_asic.sfg_generators import (
direct_form_1_iir,
direct_form_2_iir,
direct_form_fir,
+ ldlt_matrix_inverse,
radix_2_dif_fft,
transposed_direct_form_fir,
wdf_allpass,
@@ -637,3 +638,58 @@ class TestRadix2FFT:
POINTS = 5
with pytest.raises(ValueError, match="Points must be a power of two."):
radix_2_dif_fft(points=POINTS)
+
+
+class TestLdltMatrixInverse:
+ def test_1x1(self):
+ sfg = ldlt_matrix_inverse(N=1, is_complex=False)
+
+ A_input = [Constant(5)]
+
+ sim = Simulation(sfg, A_input)
+ sim.run_for(1)
+
+ res = sim.results
+ assert np.isclose(res["0"], 0.2)
+
+ def test_2x2_simple_spd(self):
+ sfg = ldlt_matrix_inverse(N=2, is_complex=False)
+
+ A = np.array([[1, 2], [2, 1]])
+ A_input = [Constant(1), Constant(2), Constant(1)]
+
+ A_inv = np.linalg.inv(A)
+
+ sim = Simulation(sfg, A_input)
+ sim.run_for(1)
+
+ res = sim.results
+ assert np.isclose(res["0"], A_inv[0, 0])
+ assert np.isclose(res["1"], A_inv[0, 1])
+ assert np.isclose(res["2"], A_inv[1, 1])
+
+ def test_3x3_simple_spd(self):
+ sfg = ldlt_matrix_inverse(N=3, is_complex=False)
+
+ A = np.array([[2, -1, 0], [-1, 3, -1], [0, -1, 2]])
+ A_input = [
+ Constant(2),
+ Constant(-1),
+ Constant(0),
+ Constant(3),
+ Constant(-1),
+ Constant(2),
+ ]
+
+ A_inv = np.linalg.inv(A)
+
+ sim = Simulation(sfg, A_input)
+ sim.run_for(1)
+
+ res = sim.results
+ assert np.isclose(res["0"], A_inv[0, 0])
+ assert np.isclose(res["1"], A_inv[0, 1])
+ assert np.isclose(res["2"], A_inv[0, 2])
+ assert np.isclose(res["3"], A_inv[1, 1])
+ assert np.isclose(res["4"], A_inv[1, 2])
+ assert np.isclose(res["5"], A_inv[2, 2])
--
GitLab