From 58931944b9bf5c832b800e0cdafed6e491758c50 Mon Sep 17 00:00:00 2001 From: Max Date: Fri, 13 Sep 2024 11:48:06 +0200 Subject: [PATCH 01/12] add sqrt option for diagonal in stencil and block matrices. Corresponds to commit 121215d in stefan-psydac. --- psydac/linalg/block.py | 8 ++++++-- psydac/linalg/stencil.py | 9 ++++++++- 2 files changed, 14 insertions(+), 3 deletions(-) diff --git a/psydac/linalg/block.py b/psydac/linalg/block.py index d54e63626..8c77ddf79 100644 --- a/psydac/linalg/block.py +++ b/psydac/linalg/block.py @@ -777,7 +777,7 @@ def __sub__(self, M): #-------------------------------------- # New properties/methods #-------------------------------------- - def diagonal(self, *, inverse = False, out = None): + def diagonal(self, *, inverse = False, sqrt = False, out = None): """Get the coefficients on the main diagonal as another BlockLinearOperator object. Parameters @@ -785,6 +785,10 @@ def diagonal(self, *, inverse = False, out = None): inverse : bool If True, get the inverse of the diagonal. (Default: False). + sqrt : bool + If True, get the square root of the diagonal. (Default: False). + Can be combined with inverse to get the inverse square root + out : BlockLinearOperator If provided, write the diagonal entries into this matrix. (Default: None). @@ -815,7 +819,7 @@ def diagonal(self, *, inverse = False, out = None): # Store the diagonal (or its inverse) into `out` for i, j in self.nonzero_block_indices: if i == j: - out[i, i] = self[i, i].diagonal(inverse = inverse, out = out[i, i]) + out[i, i] = self[i, i].diagonal(inverse = inverse, sqrt=sqrt, out = out[i, i]) return out diff --git a/psydac/linalg/stencil.py b/psydac/linalg/stencil.py index 0c818af26..a8b0402fe 100644 --- a/psydac/linalg/stencil.py +++ b/psydac/linalg/stencil.py @@ -1372,7 +1372,7 @@ def _exchange_assembly_data_serial(self): self._data[idx_to] += self._data[idx_from] # ... - def diagonal(self, *, inverse = False, out = None): + def diagonal(self, *, inverse = False, sqrt = False, out = None): """ Get the coefficients on the main diagonal as a StencilDiagonalMatrix object. @@ -1381,6 +1381,10 @@ def diagonal(self, *, inverse = False, out = None): inverse : bool If True, get the inverse of the diagonal. (Default: False). + sqrt : bool + If True, get the square root of the diagonal. (Default: False). + Can be combined with inverse to get the inverse square root + out : StencilDiagonalMatrix If provided, write the diagonal entries into this matrix. (Default: None). @@ -1410,6 +1414,9 @@ def diagonal(self, *, inverse = False, out = None): diag = self._data[diagonal_indices] data = out._data if out else None + if sqrt: + diag = np.sqrt(diag) + # Calculate entries of StencilDiagonalMatrix if inverse: data = np.divide(1, diag, out=data) From b15daf77b8c7fa74d7b50d461ce734281775ef87 Mon Sep 17 00:00:00 2001 From: maxlin-ipp Date: Wed, 18 Sep 2024 12:13:54 +0200 Subject: [PATCH 02/12] Update psydac/linalg/block.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Yaman Güçlü --- psydac/linalg/block.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/psydac/linalg/block.py b/psydac/linalg/block.py index 8c77ddf79..63e5fb445 100644 --- a/psydac/linalg/block.py +++ b/psydac/linalg/block.py @@ -787,7 +787,7 @@ def diagonal(self, *, inverse = False, sqrt = False, out = None): sqrt : bool If True, get the square root of the diagonal. (Default: False). - Can be combined with inverse to get the inverse square root + Can be combined with inverse to get the inverse square root. out : BlockLinearOperator If provided, write the diagonal entries into this matrix. (Default: None). From cd4016e84e3317fc90dcfdee32fa930a7ce721d1 Mon Sep 17 00:00:00 2001 From: maxlin-ipp Date: Wed, 18 Sep 2024 12:15:31 +0200 Subject: [PATCH 03/12] Update psydac/linalg/stencil.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Yaman Güçlü --- psydac/linalg/stencil.py | 1 + 1 file changed, 1 insertion(+) diff --git a/psydac/linalg/stencil.py b/psydac/linalg/stencil.py index a8b0402fe..60d8a948f 100644 --- a/psydac/linalg/stencil.py +++ b/psydac/linalg/stencil.py @@ -1380,6 +1380,7 @@ def diagonal(self, *, inverse = False, sqrt = False, out = None): ---------- inverse : bool If True, get the inverse of the diagonal. (Default: False). + Can be combined with sqrt to get the inverse square root. sqrt : bool If True, get the square root of the diagonal. (Default: False). From 643de9afb74e85546e17daa0d711a39188d4bc44 Mon Sep 17 00:00:00 2001 From: maxlin-ipp Date: Wed, 18 Sep 2024 12:15:38 +0200 Subject: [PATCH 04/12] Update psydac/linalg/stencil.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Yaman Güçlü --- psydac/linalg/stencil.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/psydac/linalg/stencil.py b/psydac/linalg/stencil.py index 60d8a948f..c132aad9e 100644 --- a/psydac/linalg/stencil.py +++ b/psydac/linalg/stencil.py @@ -1384,7 +1384,7 @@ def diagonal(self, *, inverse = False, sqrt = False, out = None): sqrt : bool If True, get the square root of the diagonal. (Default: False). - Can be combined with inverse to get the inverse square root + Can be combined with inverse to get the inverse square root. out : StencilDiagonalMatrix If provided, write the diagonal entries into this matrix. (Default: None). From fdd2730282ea39bc9f72bc0ebd8b9d2facfc7a56 Mon Sep 17 00:00:00 2001 From: maxlin-ipp Date: Wed, 18 Sep 2024 12:15:52 +0200 Subject: [PATCH 05/12] Update psydac/linalg/block.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Yaman Güçlü --- psydac/linalg/block.py | 1 + 1 file changed, 1 insertion(+) diff --git a/psydac/linalg/block.py b/psydac/linalg/block.py index 63e5fb445..4e9cbadc6 100644 --- a/psydac/linalg/block.py +++ b/psydac/linalg/block.py @@ -784,6 +784,7 @@ def diagonal(self, *, inverse = False, sqrt = False, out = None): ---------- inverse : bool If True, get the inverse of the diagonal. (Default: False). + Can be combined with sqrt to get the inverse square root. sqrt : bool If True, get the square root of the diagonal. (Default: False). From 62499e09f017f3ede0a537e6f1c2e2e32bb998ac Mon Sep 17 00:00:00 2001 From: Max Lindqvist Date: Wed, 18 Sep 2024 13:58:26 +0200 Subject: [PATCH 06/12] Use np.sqrt(data, out=data) for in-place sqrt --- psydac/linalg/stencil.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/psydac/linalg/stencil.py b/psydac/linalg/stencil.py index c132aad9e..6941e2832 100644 --- a/psydac/linalg/stencil.py +++ b/psydac/linalg/stencil.py @@ -1415,17 +1415,18 @@ def diagonal(self, *, inverse = False, sqrt = False, out = None): diag = self._data[diagonal_indices] data = out._data if out else None - if sqrt: - diag = np.sqrt(diag) # Calculate entries of StencilDiagonalMatrix if inverse: - data = np.divide(1, diag, out=data) + np.divide(1, diag, out=data) elif out: np.copyto(data, diag) else: data = diag.copy() + if sqrt: + np.sqrt(data, out=data) + # If needed create a new StencilDiagonalMatrix object if out is None: out = StencilDiagonalMatrix(V, W, data) From 1a65a28a73daf517137664870a190f8280b84899 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 18 Sep 2024 17:03:39 +0200 Subject: [PATCH 07/12] Apply minor suggestions from my very last code review --- psydac/linalg/block.py | 2 +- psydac/linalg/stencil.py | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/psydac/linalg/block.py b/psydac/linalg/block.py index 4e9cbadc6..49e29b60d 100644 --- a/psydac/linalg/block.py +++ b/psydac/linalg/block.py @@ -820,7 +820,7 @@ def diagonal(self, *, inverse = False, sqrt = False, out = None): # Store the diagonal (or its inverse) into `out` for i, j in self.nonzero_block_indices: if i == j: - out[i, i] = self[i, i].diagonal(inverse = inverse, sqrt=sqrt, out = out[i, i]) + out[i, i] = self[i, i].diagonal(inverse = inverse, sqrt = sqrt, out = out[i, i]) return out diff --git a/psydac/linalg/stencil.py b/psydac/linalg/stencil.py index 6941e2832..afad8e0d2 100644 --- a/psydac/linalg/stencil.py +++ b/psydac/linalg/stencil.py @@ -1415,7 +1415,6 @@ def diagonal(self, *, inverse = False, sqrt = False, out = None): diag = self._data[diagonal_indices] data = out._data if out else None - # Calculate entries of StencilDiagonalMatrix if inverse: np.divide(1, diag, out=data) From 50530cf1e460bc1ea014ad25831ca11ed740a37f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 18 Sep 2024 17:26:58 +0200 Subject: [PATCH 08/12] Update continuous-integration.yml Try `apt-get update` to avoid download error. --- .github/workflows/continuous-integration.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/continuous-integration.yml b/.github/workflows/continuous-integration.yml index 6181f9861..ae8d0ab61 100644 --- a/.github/workflows/continuous-integration.yml +++ b/.github/workflows/continuous-integration.yml @@ -58,6 +58,7 @@ jobs: - name: Set default MPI and HDF5 C compilers on Ubuntu if: matrix.os == 'ubuntu-latest' run: | + sudo apt-get update sudo apt-get install --reinstall openmpi-bin libhdf5-openmpi-dev - name: Install non-Python dependencies on macOS From 25f3a90e4d35f9a676e7c34f6ce4b7922a7daf0f Mon Sep 17 00:00:00 2001 From: Max Date: Thu, 19 Sep 2024 13:48:21 +0200 Subject: [PATCH 09/12] Changed in-place division operation to assignment to handle cases where data is None. --- psydac/linalg/stencil.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/psydac/linalg/stencil.py b/psydac/linalg/stencil.py index afad8e0d2..882d2c854 100644 --- a/psydac/linalg/stencil.py +++ b/psydac/linalg/stencil.py @@ -1417,7 +1417,7 @@ def diagonal(self, *, inverse = False, sqrt = False, out = None): # Calculate entries of StencilDiagonalMatrix if inverse: - np.divide(1, diag, out=data) + data = np.divide(1, diag, out=data) elif out: np.copyto(data, diag) else: From d770331b6f615f131d6a2291da7a6a3fb62bd681 Mon Sep 17 00:00:00 2001 From: Max Date: Wed, 25 Sep 2024 13:38:11 +0200 Subject: [PATCH 10/12] added a test for diagonal(sqrt=True) --- psydac/linalg/tests/test_linalg.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/psydac/linalg/tests/test_linalg.py b/psydac/linalg/tests/test_linalg.py index b34e222c1..ef7c0e98a 100644 --- a/psydac/linalg/tests/test_linalg.py +++ b/psydac/linalg/tests/test_linalg.py @@ -528,8 +528,9 @@ def test_in_place_operations(n1, n2, p1, p2, P1=False, P2=False): @pytest.mark.parametrize('n2', n2array) @pytest.mark.parametrize('p1', p1array) @pytest.mark.parametrize('p2', p2array) +@pytest.mark.parametrize('sqrt', [False, True]) -def test_inverse_transpose_interaction(n1, n2, p1, p2, P1=False, P2=False): +def test_inverse_transpose_interaction(n1, n2, p1, p2, sqrt, P1=False, P2=False): # 1. Initiate square LOs: S (V->V, StencilMatrix), S1 (W->W, StencilMatrix) #Initiate BlockLO: B (VxW -> VxW) and a BlockVector u element of VxW @@ -618,6 +619,12 @@ def test_inverse_transpose_interaction(n1, n2, p1, p2, P1=False, P2=False): ### -1,T & T,-1 --- -1,T,T --- -1,T,-1 --- T,-1,-1 --- T,-1,T (the combinations I test) ### + # Square root test + scaled_matrix = B * np.random.random() # Ensure the diagonal elements != 1 + diagonal_values = scaled_matrix.diagonal(sqrt=False).toarray() + sqrt_diagonal_values = scaled_matrix.diagonal(sqrt=True).toarray() + assert np.allclose(sqrt_diagonal_values, np.sqrt(diagonal_values)) + tol = 1e-5 C = inverse(B, 'cg', tol=tol) P = B.diagonal(inverse=True) From d890058bb0c76cccf48a6786664a5cb3a0bbbbb5 Mon Sep 17 00:00:00 2001 From: Max Date: Thu, 26 Sep 2024 16:32:46 +0200 Subject: [PATCH 11/12] removed unused sqrt parameter --- psydac/linalg/tests/test_linalg.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/psydac/linalg/tests/test_linalg.py b/psydac/linalg/tests/test_linalg.py index ef7c0e98a..64f0548b4 100644 --- a/psydac/linalg/tests/test_linalg.py +++ b/psydac/linalg/tests/test_linalg.py @@ -528,9 +528,8 @@ def test_in_place_operations(n1, n2, p1, p2, P1=False, P2=False): @pytest.mark.parametrize('n2', n2array) @pytest.mark.parametrize('p1', p1array) @pytest.mark.parametrize('p2', p2array) -@pytest.mark.parametrize('sqrt', [False, True]) -def test_inverse_transpose_interaction(n1, n2, p1, p2, sqrt, P1=False, P2=False): +def test_inverse_transpose_interaction(n1, n2, p1, p2, P1=False, P2=False): # 1. Initiate square LOs: S (V->V, StencilMatrix), S1 (W->W, StencilMatrix) #Initiate BlockLO: B (VxW -> VxW) and a BlockVector u element of VxW From 8faa3eccf324e49308165f9a6c51c3fa1922a941 Mon Sep 17 00:00:00 2001 From: maxlin-ipp Date: Fri, 27 Sep 2024 13:25:38 +0200 Subject: [PATCH 12/12] Update psydac/linalg/tests/test_linalg.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Yaman Güçlü --- psydac/linalg/tests/test_linalg.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/psydac/linalg/tests/test_linalg.py b/psydac/linalg/tests/test_linalg.py index 64f0548b4..c9962c616 100644 --- a/psydac/linalg/tests/test_linalg.py +++ b/psydac/linalg/tests/test_linalg.py @@ -622,7 +622,7 @@ def test_inverse_transpose_interaction(n1, n2, p1, p2, P1=False, P2=False): scaled_matrix = B * np.random.random() # Ensure the diagonal elements != 1 diagonal_values = scaled_matrix.diagonal(sqrt=False).toarray() sqrt_diagonal_values = scaled_matrix.diagonal(sqrt=True).toarray() - assert np.allclose(sqrt_diagonal_values, np.sqrt(diagonal_values)) + assert np.array_equal(sqrt_diagonal_values, np.sqrt(diagonal_values)) tol = 1e-5 C = inverse(B, 'cg', tol=tol)