Skip to content

Commit

Permalink
Enable compilation with language C
Browse files Browse the repository at this point in the history
  • Loading branch information
camillo81 committed Nov 14, 2023
1 parent d393a7c commit 9a016b1
Show file tree
Hide file tree
Showing 5 changed files with 191 additions and 42 deletions.
94 changes: 94 additions & 0 deletions psydac/core/arrays.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,94 @@
from pyccel.decorators import pure
from numpy import shape

@pure
def matmul(a: 'float[:,:]', b: 'float[:,:]', c: 'float[:,:]'):
"""
Performs the matrix-matrix product a*b and writes the result into c.
Parameters
----------
a : array[float]
The first input array (matrix).
b : array[float]
The second input array (matrix).
c : array[float]
The output array (matrix) which is the result of the matrix-matrix product a.dot(b).
"""

c[:, :] = 0.

sh_a = shape(a)
sh_b = shape(b)

for i in range(sh_a[0]):
for j in range(sh_b[1]):
for k in range(sh_a[1]):
c[i, j] += a[i, k] * b[k, j]


@pure
def sum_vec(a: 'float[:]') -> float:
"""
Sum the elements of a 1D vector.
Parameters
----------
a : array[float]
The 1d vector.
"""

out = 0.

sh_a = shape(a)

for i in range(sh_a[0]):
out += a[i]

return out


@pure
def min_vec(a: 'float[:]') -> float:
"""
Compute the minimum a 1D vector.
Parameters
----------
a : array[float]
The 1D vector.
"""

out = a[0]

sh_a = shape(a)

for i in range(sh_a[0]):
if a[i] < out:
out = a[i]

return out


@pure
def max_vec(a: 'float[:]') -> float:
"""
Compute the maximum a 1D vector.
Parameters
----------
a : array[float]
The 1D vector.
"""

out = a[0]

sh_a = shape(a)

for i in range(sh_a[0]):
if a[i] > out:
out = a[i]

return out
15 changes: 8 additions & 7 deletions psydac/core/bsplines_pyccel.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
# like Numpy functions.

import numpy as np
from psydac.core.arrays import matmul, sum_vec, min_vec, max_vec

# =============================================================================
def find_span_p(knots: 'float[:]', degree: int, x: float):
Expand Down Expand Up @@ -322,7 +323,7 @@ def basis_funs_all_ders_p(knots: 'float[:]', degree: int, x: float, span: int, n
j2 = k-1 if (r-1 <= pk) else degree-r

a[s2, j1:j2 + 1] = (a[s1, j1:j2 + 1] - a[s1, j1 - 1:j2]) * ndu[pk + 1, rk + j1:rk + j2 + 1]
temp_d[:, :] = np.matmul(a[s2:s2 + 1, j1:j2 + 1], ndu[rk + j1:rk + j2 + 1, pk: pk + 1])
matmul(a[s2:s2 + 1, j1:j2 + 1], ndu[rk + j1:rk + j2 + 1, pk: pk + 1], temp_d)
d+= temp_d[0, 0]
if r <= pk:
a[s2, k] = - a[s1, k - 1] * ndu[pk + 1, r]
Expand Down Expand Up @@ -574,7 +575,7 @@ def histopolation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, norma
jend = min(spans[i + 1], n)
# Compute non-zero values of histopolation matrix
for j in range(1 + jstart, jend + 1):
s = np.sum(colloc[i, 0:j]) - np.sum(colloc[i + 1, 0:j])
s = sum_vec(colloc[i, 0:j]) - sum_vec(colloc[i + 1, 0:j])
H[i, j - 1] = s

else:
Expand All @@ -586,7 +587,7 @@ def histopolation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, norma
jend = min(spans[i + 1], n)
# Compute non-zero values of histopolation matrix
for j in range(1 + jstart, jend + 1):
s = np.sum(colloc[i, 0:j]) - np.sum(colloc[i + 1, 0:j])
s = sum_vec(colloc[i, 0:j]) - sum_vec(colloc[i + 1, 0:j])
H[i, j - 1] = s * integrals[j - 1]

# Mitigate round-off errors
Expand Down Expand Up @@ -714,10 +715,10 @@ def greville_p(knots: 'float[:]', degree: int, periodic: bool, out:'float[:]'):
# Compute greville abscissas as average of p consecutive knot values
if p == 0:
for i in range(n):
out[i] = sum(T[i:i + 2]) / 2
out[i] = sum_vec(T[i:i + 2]) / 2
else:
for i in range(1, 1+n):
out[i - 1] = sum(T[i:i + p]) / p
out[i - 1] = sum_vec(T[i:i + p]) / p

# Domain boundaries
a = T[p]
Expand Down Expand Up @@ -1056,8 +1057,8 @@ def cell_index_p(breaks: 'float[:]', i_grid: 'float[:]', tol: float, out: 'int[:
nbk = len(breaks)

# Check if there are points outside the domain
if np.min(i_grid) < breaks[0] - tol: return 1
if np.max(i_grid) > breaks[nbk - 1] + tol: return 1
if min_vec(i_grid) < breaks[0] - tol: return 1
if max_vec(i_grid) > breaks[nbk - 1] + tol: return 1

current_index = 0
while current_index < nx:
Expand Down
71 changes: 50 additions & 21 deletions psydac/core/kernels.py
Original file line number Diff line number Diff line change
Expand Up @@ -3267,6 +3267,8 @@ def eval_jacobians_inv_3d(nc1: int, nc2: int, nc3: int, f_p1: int, f_p2: int,
arr_z_x2 = np.zeros_like(global_arr_coeff_z, shape=(k1, k2, k3))
arr_z_x3 = np.zeros_like(global_arr_coeff_z, shape=(k1, k2, k3))

tmp_a = np.zeros((3, 3), dtype=float)

for i_cell_1 in range(nc1):
span_1 = global_spans_1[i_cell_1]

Expand Down Expand Up @@ -3361,12 +3363,14 @@ def eval_jacobians_inv_3d(nc1: int, nc2: int, nc3: int, f_p1: int, f_p2: int,
a_32 = - x_x1 * y_x3 + x_x3 * y_x1
a_33 = x_x1 * y_x2 - x_x2 * y_x1

tmp_a[:] = np.array([[a_11, a_21, a_31],
[a_12, a_22, a_32],
[a_13, a_23, a_33]])

jacobians_inv[i_cell_1 * k1 + i_quad_1,
i_cell_2 * k2 + i_quad_2,
i_cell_3 * k3 + i_quad_3,
:, :] = np.array([[a_11, a_21, a_31],
[a_12, a_22, a_32],
[a_13, a_23, a_33]]) / det
:, :] = tmp_a / det


@template(name='T', types=['float[:,:]', 'complex[:,:]'])
Expand Down Expand Up @@ -3420,6 +3424,8 @@ def eval_jacobians_inv_2d(nc1: int, nc2: int, f_p1: int, f_p2: int, k1: int, k2

arr_y_x1 = np.zeros_like(global_arr_coeff_y, shape=(k1, k2))
arr_y_x2 = np.zeros_like(global_arr_coeff_y, shape=(k1, k2))

tmp_a = np.zeros((2, 2), dtype=float)

for i_cell_1 in range(nc1):
span_1 = global_spans_1[i_cell_1]
Expand Down Expand Up @@ -3468,11 +3474,13 @@ def eval_jacobians_inv_2d(nc1: int, nc2: int, f_p1: int, f_p2: int, k1: int, k2
y_x2 = arr_y_x2[i_quad_1, i_quad_2]

det = x_x1 * y_x2 - x_x2 * y_x1

tmp_a[:] = np.array([[y_x2, - x_x2],
- [- y_x1, x_x1]])

jacobians_inv[i_cell_1 * k1 + i_quad_1,
i_cell_2 * k2 + i_quad_2,
:, :] = np.array([[y_x2, - x_x2],
[- y_x1, x_x1]]) / det
:, :] = tmp_a / det


# -----------------------------------------------------------------------------
Expand Down Expand Up @@ -3551,6 +3559,8 @@ def eval_jacobians_inv_irregular_3d(np1: int, np2: int, np3: int, f_p1: int, f_p
temp_z_x1 = arr_coeffs_x[0,0,0]-arr_coeffs_x[0,0,0]
temp_z_x2 = arr_coeffs_y[0,0,0]-arr_coeffs_y[0,0,0]
temp_z_x3 = arr_coeffs_z[0,0,0]-arr_coeffs_z[0,0,0]

tmp_a = np.zeros((3, 3), dtype=float)

for i_p_1 in range(np1):
i_cell_1 = cell_index_1[i_p_1]
Expand Down Expand Up @@ -3639,12 +3649,14 @@ def eval_jacobians_inv_irregular_3d(np1: int, np2: int, np3: int, f_p1: int, f_p
a_32 = - temp_x_x1 * temp_y_x3 + temp_x_x3 * temp_y_x1
a_33 = temp_x_x1 * temp_y_x2 - temp_x_x2 * temp_y_x1

tmp_a[:] = np.array([[a_11, a_21, a_31],
[a_12, a_22, a_32],
[a_13, a_23, a_33]])

jacobians_inv[i_p_1,
i_p_2,
i_p_3,
:, :] = np.array([[a_11, a_21, a_31],
[a_12, a_22, a_32],
[a_13, a_23, a_33]]) / det
:, :] = tmp_a / det


@template(name='T', types=['float[:,:]', 'complex[:,:]'])
Expand Down Expand Up @@ -3698,7 +3710,8 @@ def eval_jacobians_inv_irregular_2d(np1: int, np2: int, f_p1: int, f_p2: int, ce

temp_y_x1 = arr_coeffs_x[0,0]-arr_coeffs_x[0,0]
temp_y_x2 = arr_coeffs_y[0,0]-arr_coeffs_y[0,0]


tmp_a = np.zeros((2, 2), dtype=float)

for i_p_1 in range(np1):
i_cell_1 = cell_index_1[i_p_1]
Expand Down Expand Up @@ -3742,8 +3755,9 @@ def eval_jacobians_inv_irregular_2d(np1: int, np2: int, f_p1: int, f_p2: int, ce

det = temp_x_x1 * temp_y_x2 - temp_y_x1 * temp_x_x2

jacobians_inv[i_p_1, i_p_2, :, :] = np.array([[temp_y_x2, - temp_x_x2],
[- temp_y_x1, temp_x_x1]]) / det
tmp_a[:] = np.array([[temp_y_x2, - temp_x_x2], [- temp_y_x1, temp_x_x1]])

jacobians_inv[i_p_1, i_p_2, :, :] = tmp_a / det

# -----------------------------------------------------------------------------
# 3: Regular tensor grid with weights
Expand Down Expand Up @@ -3837,6 +3851,8 @@ def eval_jacobians_inv_3d_weights(nc1: int, nc2: int, nc3: int, f_p1: int, f_p2
arr_weights_x2 = np.zeros((k1, k2, k3))
arr_weights_x3 = np.zeros((k1, k2, k3))

tmp_a = np.zeros((3, 3), dtype=float)

for i_cell_1 in range(nc1):
span_1 = global_spans_1[i_cell_1]

Expand Down Expand Up @@ -3983,12 +3999,14 @@ def eval_jacobians_inv_3d_weights(nc1: int, nc2: int, nc3: int, f_p1: int, f_p2
a_32 = - x_x1 * y_x3 + x_x3 * y_x1
a_33 = x_x1 * y_x2 - x_x2 * y_x1

tmp_a[:] = np.array([[a_11, a_21, a_31],
[a_12, a_22, a_32],
[a_13, a_23, a_33]])

jacobians_inv[i_cell_1 * k1 + i_quad_1,
i_cell_2 * k2 + i_quad_2,
i_cell_3 * k3 + i_quad_3,
:, :] = np.array([[a_11, a_21, a_31],
[a_12, a_22, a_32],
[a_13, a_23, a_33]]) / det
:, :] = tmp_a / det


@template(name='T', types=['float[:,:]', 'complex[:,:]'])
Expand Down Expand Up @@ -4056,6 +4074,8 @@ def eval_jacobians_inv_2d_weights(nc1: int, nc2: int, f_p1: int, f_p2: int, k1:
arr_weights_x1 = np.zeros((k1, k2))
arr_weights_x2 = np.zeros((k1, k2))

tmp_a = np.zeros((2, 2), dtype=float)

for i_cell_1 in range(nc1):
span_1 = global_spans_1[i_cell_1]

Expand Down Expand Up @@ -4142,10 +4162,12 @@ def eval_jacobians_inv_2d_weights(nc1: int, nc2: int, f_p1: int, f_p2: int, k1:

det = x_x1 * y_x2 - x_x2 * y_x1

tmp_a[:] = np.array([[y_x2, - x_x2],
[- y_x1, x_x1]])

jacobians_inv[i_cell_1 * k1 + i_quad_1,
i_cell_2 * k2 + i_quad_2,
:, :] = np.array([[y_x2, - x_x2],
[- y_x1, x_x1]]) / det
:, :] = tmp_a / det


# -----------------------------------------------------------------------------
Expand Down Expand Up @@ -4234,6 +4256,8 @@ def eval_jacobians_inv_irregular_3d_weights(np1: int, np2: int, np3: int, f_p1:
temp_z_x2 = arr_coeffs_y[0,0,0]-arr_coeffs_y[0,0,0]
temp_z_x3 = arr_coeffs_z[0,0,0]-arr_coeffs_z[0,0,0]

tmp_a = np.zeros((3, 3), dtype=float)

for i_p_1 in range(np1):
i_cell_1 = cell_index_1[i_p_1]
span_1 = global_spans_1[i_cell_1]
Expand Down Expand Up @@ -4359,12 +4383,14 @@ def eval_jacobians_inv_irregular_3d_weights(np1: int, np2: int, np3: int, f_p1:
a_32 = - x_x1 * y_x3 + x_x3 * y_x1
a_33 = x_x1 * y_x2 - x_x2 * y_x1

tmp_a[:] = np.array([[a_11, a_21, a_31],
[a_12, a_22, a_32],
[a_13, a_23, a_33]])

jacobians_inv[i_p_1,
i_p_2,
i_p_3,
:, :] = np.array([[a_11, a_21, a_31],
[a_12, a_22, a_32],
[a_13, a_23, a_33]]) / det
:, :] = tmp_a / det


@template(name='T', types=['float[:,:]', 'complex[:,:]'])
Expand Down Expand Up @@ -4429,6 +4455,8 @@ def eval_jacobians_inv_irregular_2d_weights(np1: int, np2: int, f_p1: int, f_p2:

temp_y_x1 = arr_coeffs_x[0,0]-arr_coeffs_x[0,0]
temp_y_x2 = arr_coeffs_y[0,0]-arr_coeffs_y[0,0]

tmp_a = np.zeros((2, 2), dtype=float)

for i_p_1 in range(np1):
i_cell_1 = cell_index_1[i_p_1]
Expand Down Expand Up @@ -4501,8 +4529,9 @@ def eval_jacobians_inv_irregular_2d_weights(np1: int, np2: int, f_p1: int, f_p2:

det = x_x1 * y_x2 - x_x2 * y_x1

jacobians_inv[i_p_1, i_p_2, :, :] = np.array([[y_x2, - x_x2],
[- y_x1, x_x1]]) / det
tmp_a[:] = np.array([[y_x2, - x_x2], [- y_x1, x_x1]])

jacobians_inv[i_p_1, i_p_2, :, :] = tmp_a / det


# ==========================================================================
Expand Down
25 changes: 25 additions & 0 deletions psydac/core/tests/test_arrays.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,25 @@
import numpy as np
import pytest

from psydac.core.arrays import matmul, sum_vec, min_vec, max_vec

@pytest.mark.parametrize('seed', [123, 6, 49573])
def test_arrays(seed):

np.random.seed(seed)

a = np.random.rand(30, 20)
b = np.random.rand(20, 40)
c = np.zeros((30, 40))

matmul(a, b, c)
assert np.allclose(np.matmul(a, b), c)

v = np.random.rand(32)
assert np.isclose(np.sum(v), sum_vec(v))
assert np.isclose(np.min(v), min_vec(v))
assert np.isclose(np.max(v), max_vec(v))


if __name__ == '__main__':
test_arrays(459)
Loading

1 comment on commit 9a016b1

@spossann
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Definitely needed in Struphy. Martin wrote an issue about it, since Psydac also at some point wants to enable C compilation.

Please sign in to comment.