Cuda 1 xp arrays

conftest.py

@@ -1,6 +1,7 @@
 """Root-level pytest configuration."""
 import pytest
 import sys
+import importlib.util
 from pathlib import Path
 
 
@@ -32,13 +33,17 @@ def pytest_collection_modifyitems(config, items):
 
     items_to_remove = []
     skip = pytest.mark.skip(reason="Requires optional dependency (sympde)")
+    petsc_available = importlib.util.find_spec("petsc4py") is not None
 
     for item in items:
         # Skip if module is in skip list
         if item.fspath.basename in skip_modules:
             items_to_remove.append(item)
             continue
 
+        if item.get_closest_marker("petsc") and not petsc_available:
+            item.add_marker(pytest.mark.skip(reason="petsc4py is not installed"))
+
         # If running with xdist, automatically skip mpi and petsc tests
         if config.pluginmanager.has_plugin("xdist"):
             if item.get_closest_marker("mpi") or item.get_closest_marker("petsc"):

feectools/core/bsplines.py

@@ -16,6 +16,7 @@
 
 """
 import cunumpy as xp
+from cunumpy import PyccelKernel
 from cunumpy.xp import array_backend
 import numpy as np
 
@@ -38,6 +39,27 @@
                                           cell_index_p,
                                           basis_ders_on_irregular_grid_p)
 
+# Kernels generated by Pyccel only understand NumPy arrays; wrap them so they
+# can also be called with CuPy arrays (see cunumpy.kernel.PyccelKernel).
+find_span_p = PyccelKernel(find_span_p)
+find_spans_p = PyccelKernel(find_spans_p)
+basis_funs_p = PyccelKernel(basis_funs_p)
+basis_funs_array_p = PyccelKernel(basis_funs_array_p)
+basis_funs_1st_der_p = PyccelKernel(basis_funs_1st_der_p)
+basis_funs_all_ders_p = PyccelKernel(basis_funs_all_ders_p)
+collocation_matrix_p = PyccelKernel(collocation_matrix_p)
+histopolation_matrix_p = PyccelKernel(histopolation_matrix_p)
+greville_p = PyccelKernel(greville_p)
+breakpoints_p = PyccelKernel(breakpoints_p)
+elements_spans_p = PyccelKernel(elements_spans_p)
+make_knots_p = PyccelKernel(make_knots_p)
+elevate_knots_p = PyccelKernel(elevate_knots_p)
+quadrature_grid_p = PyccelKernel(quadrature_grid_p)
+basis_ders_on_quad_grid_p = PyccelKernel(basis_ders_on_quad_grid_p)
+basis_integrals_p = PyccelKernel(basis_integrals_p)
+cell_index_p = PyccelKernel(cell_index_p)
+basis_ders_on_irregular_grid_p = PyccelKernel(basis_ders_on_irregular_grid_p)
+
 __all__ = ('find_span',
            'find_spans',
            'basis_funs',
@@ -84,7 +106,7 @@ def find_span(knots, degree, x):
         Knot span index.
     """
     x = float(x)
-    knots = xp.ascontiguousarray(knots, dtype=float)
+    knots = xp.ascontiguousarray(xp.asarray(knots), dtype=float)
     return find_span_p(knots, degree, x)
 
 #==============================================================================
@@ -116,8 +138,8 @@ def find_spans(knots, degree, x, out=None):
     spans : array of ints
         Knots span indexes.
     """
-    knots = xp.ascontiguousarray(knots, dtype=float)
-    x = xp.ascontiguousarray(x, dtype=float)
+    knots = xp.ascontiguousarray(xp.asarray(knots), dtype=float)
+    x = xp.ascontiguousarray(xp.asarray(x), dtype=float)
     if out is None:
         out = xp.zeros_like(x, dtype=int)
     else:
@@ -155,7 +177,7 @@ def basis_funs(knots, degree, x, span, out=None):
         1D array containing the values of ``degree + 1`` non-zero
         Bsplines at location ``x``.
     """
-    knots = xp.ascontiguousarray(knots, dtype=float)
+    knots = xp.ascontiguousarray(xp.asarray(knots), dtype=float)
     # Get native float
     x = float(x)
     if out is None:
@@ -193,8 +215,8 @@ def basis_funs_array(knots, degree, span, x, out=None):
         2D array of shape ``(len(x), degree + 1)`` containing the values of ``degree + 1`` non-zero
         Bsplines at each location in ``x``.
     """
-    knots = xp.ascontiguousarray(knots, dtype=float)
-    x = xp.ascontiguousarray(x, dtype=float)
+    knots = xp.ascontiguousarray(xp.asarray(knots), dtype=float)
+    x = xp.ascontiguousarray(xp.asarray(x), dtype=float)
     if out is None:
         out = xp.zeros(x.shape + (degree + 1,), dtype=float)
     else:
@@ -240,7 +262,7 @@ def basis_funs_1st_der(knots, degree, x, span, out=None):
     ----------
     .. [2] SELALIB, Semi-Lagrangian Library. http://selalib.gforge.inria.fr
     """
-    knots = xp.ascontiguousarray(knots, dtype=float)
+    knots = xp.ascontiguousarray(xp.asarray(knots), dtype=float)
     # Get native float to work on windows
     x = float(x)
     if out is None:
@@ -291,7 +313,7 @@ def basis_funs_all_ders(knots, degree, x, span, n, normalization='B', out=None):
         ders[i,j] = (d/dx)^i B_k(x) with k=(span-degree+j),
         for 0 <= i <= n and 0 <= j <= degree+1.
     """
-    knots = xp.ascontiguousarray(knots, dtype=float)
+    knots = xp.ascontiguousarray(xp.asarray(knots), dtype=float)
     # Get native float to work on windows
     x = float(x)
     if out is None:
@@ -346,8 +368,8 @@ def collocation_matrix(knots, degree, periodic, normalization, xgrid, out=None,
     if xgrid.size == 1:
         return xp.ones((1, 1), dtype=float)
 
-    knots = xp.ascontiguousarray(knots, dtype=float)
-    xgrid = xp.ascontiguousarray(xgrid, dtype=float)
+    knots = xp.ascontiguousarray(xp.asarray(knots), dtype=float)
+    xgrid = xp.ascontiguousarray(xp.asarray(xgrid), dtype=float)
     if out is None:
         nb = len(knots) - degree - 1
         if periodic:
@@ -430,8 +452,8 @@ def histopolation_matrix(knots, degree, periodic, normalization, xgrid, multipli
     if not xp.all(xp.diff(xgrid) > 0):
         raise ValueError("Grid points must be ordered, with no repetitions: {}".format(xgrid))
 
-    knots = xp.ascontiguousarray(knots, dtype=float)
-    xgrid = xp.ascontiguousarray(xgrid, dtype=float)
+    knots = xp.ascontiguousarray(xp.asarray(knots), dtype=float)
+    xgrid = xp.ascontiguousarray(xp.asarray(xgrid), dtype=float)
     elevated_knots = elevate_knots(knots, degree, periodic, multiplicity=multiplicity)
 
     normalization = normalization == "M"
@@ -477,7 +499,7 @@ def breakpoints(knots, degree, tol=1e-15, out=None):
     breaks : numpy.ndarray (1D)
         Abscissas of all breakpoints.
     """
-    knots = xp.ascontiguousarray(knots, dtype=float)
+    knots = xp.ascontiguousarray(xp.asarray(knots), dtype=float)
     if out is None:
         out = xp.zeros(len(knots), dtype=float)
     else:
@@ -518,8 +540,7 @@ def greville(knots, degree, periodic, out=None, multiplicity=1):
     # Greville points are index arrays, keep on NumPy
     if isinstance(knots, (list, tuple)):
         knots = np.asarray(knots, dtype=float)
-    if hasattr(knots, 'get'):
-        knots = knots.get()  # Convert CuPy to NumPy
+    knots = xp.to_numpy(knots)
     knots = np.ascontiguousarray(knots, dtype=float)
     if out is None:
         n = len(knots) - 2 * degree - 2 + multiplicity if periodic else len(knots) - degree - 1
@@ -572,7 +593,7 @@ def elements_spans(knots, degree, out=None):
        spans  = xp.searchsorted( knots, breaks[:-1], side='right' ) - 1
 
     """
-    knots = xp.ascontiguousarray(knots, dtype=float)
+    knots = xp.ascontiguousarray(xp.asarray(knots), dtype=float)
     if out is None:
         out = np.zeros(len(knots), dtype=xp.int64)
     else:
@@ -624,7 +645,7 @@ def make_knots(breaks, degree, periodic, multiplicity=1, out=None):
     # Consistency checks
     assert len(breaks) > 1
     # Convert to numpy for comparison since assertion needs Python bool
-    breaks_np = breaks.get() if hasattr(breaks, 'get') else breaks
+    breaks_np = xp.to_numpy(breaks)
     if isinstance(breaks_np, (list, tuple)):
         breaks_np = np.asarray(breaks_np)
     assert all( np.diff(breaks_np) > 0 )
@@ -638,8 +659,7 @@ def make_knots(breaks, degree, periodic, multiplicity=1, out=None):
 
     # Keep breaks on NumPy for initialization - knots are index arrays needed for CPU operations
     breaks = np.asarray(breaks, dtype=float) if isinstance(breaks, (list, tuple)) else breaks
-    if hasattr(breaks, 'get'):
-        breaks = breaks.get()  # Convert CuPy to NumPy
+    breaks = xp.to_numpy(breaks)
     breaks = np.ascontiguousarray(breaks, dtype=float)
     if out is None:
         # Knots are index arrays, keep them on NumPy
@@ -693,8 +713,7 @@ def elevate_knots(knots, degree, periodic, multiplicity=1, tol=1e-15, out=None):
     multiplicity = int(multiplicity)
     if isinstance(knots, (list, tuple)):
         knots = np.asarray(knots, dtype=float)
-    if hasattr(knots, 'get'):
-        knots = knots.get()  # Convert CuPy to NumPy
+    knots = xp.to_numpy(knots)
     knots = np.ascontiguousarray(knots, dtype=float)
     if out is None:
         if periodic:
@@ -771,14 +790,13 @@ def quadrature_grid(breaks, quad_rule_x, quad_rule_w):
     assert max(quad_rule_x) <= +1
 
     # Convert breaks to numpy if CuPy (breaks/grids should stay on CPU)
-    if hasattr(breaks, 'get'):
-        breaks = breaks.get()
+    breaks = xp.to_numpy(breaks)
     breaks = np.ascontiguousarray(breaks, dtype=float)
 
     if array_backend.backend == "cupy":
         # Convert CuPy arrays to NumPy
-        quad_rule_x = quad_rule_x.get() if hasattr(quad_rule_x, 'get') else quad_rule_x
-        quad_rule_w = quad_rule_w.get() if hasattr(quad_rule_w, 'get') else quad_rule_w
+        quad_rule_x = xp.to_numpy(quad_rule_x)
+        quad_rule_w = xp.to_numpy(quad_rule_w)
 
     quad_rule_x = np.ascontiguousarray(quad_rule_x, dtype=float)
     quad_rule_w = np.ascontiguousarray(quad_rule_w, dtype=float)
@@ -848,8 +866,8 @@ def basis_ders_on_quad_grid(knots, degree, quad_grid, nders, normalization, offs
     """
     offset = int(offset)
     ne, nq = quad_grid.shape
-    knots = xp.ascontiguousarray(knots, dtype=float)
-    quad_grid = xp.ascontiguousarray(quad_grid, dtype=float)
+    knots = xp.ascontiguousarray(xp.asarray(knots), dtype=float)
+    quad_grid = xp.ascontiguousarray(xp.asarray(quad_grid), dtype=float)
     if out is None:
         out = xp.zeros((ne, degree + 1, nders + 1, nq), dtype=float)
     else:
@@ -892,7 +910,7 @@ def basis_integrals(knots, degree, out=None):
     to (len(knots)-degree-1). In the periodic case the last (degree) values in
     the array are redundant, as they are a copy of the first (degree) values.
     """
-    knots = xp.ascontiguousarray(knots, dtype=float)
+    knots = xp.ascontiguousarray(xp.asarray(knots), dtype=float)
     if out is None:
         out = xp.zeros(len(knots) - degree - 1, dtype=float)
     else:
@@ -934,8 +952,8 @@ def cell_index(breaks, i_grid, tol=1e-15, out=None):
         ``cell_index[i]`` is the index of the cell in which
         ``i_grid[i]`` belong.
     """
-    breaks = xp.ascontiguousarray(breaks, dtype=float)
-    i_grid = xp.ascontiguousarray(i_grid, dtype=float)
+    breaks = xp.ascontiguousarray(xp.asarray(breaks), dtype=float)
+    i_grid = xp.ascontiguousarray(xp.asarray(i_grid), dtype=float)
     if out is None:
         out = np.zeros_like(i_grid, dtype=xp.int64)
     else:
@@ -990,8 +1008,8 @@ def basis_ders_on_irregular_grid(knots, degree, i_grid, cell_index, nders, norma
         . il: local basis function   (0 <= il <= degree)
         . id: derivative             (0 <= id <= nders )
     """
-    knots = xp.ascontiguousarray(knots, dtype=float)
-    i_grid = xp.ascontiguousarray(i_grid, dtype=float)
+    knots = xp.ascontiguousarray(xp.asarray(knots), dtype=float)
+    i_grid = xp.ascontiguousarray(xp.asarray(i_grid), dtype=float)
     if out is None:
         nx = i_grid.shape[0]
         out = xp.zeros((nx, degree + 1, nders + 1), dtype=float)

feectools/core/tests/test_bsplines.py

@@ -159,7 +159,7 @@ def test_histopolation_matrix(lims, nc, p, periodic, tol=1e-13):
 def test_cell_index(i_grid, expected):
     breaks = xp.array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.])
     out = cell_index(breaks, xp.asarray(i_grid))
-    assert xp.array_equal(expected, out)
+    assert xp.array_equal(xp.asarray(expected), out)
 
 #==============================================================================
 # SCRIPT FUNCTIONALITY: PLOT BASIS FUNCTIONS

feectools/core/tests/test_bsplines_kernel.py

@@ -2,14 +2,16 @@
 
 import pytest
 import cunumpy as xp
+import numpy as np
 
 
 from feectools.core.bsplines_kernels import cell_index_p
 
 def test_cell_index_p():
-    breaks = xp.array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.])
-    breaks = xp.ascontiguousarray(breaks, dtype=float)
-    out = xp.zeros_like(breaks, dtype=xp.int64)
+    # This directly tests the raw Pyccel kernel, which intentionally accepts
+    # NumPy host arrays only; CuPy coverage belongs to the public wrapper.
+    breaks = np.ascontiguousarray([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.], dtype=float)
+    out = np.zeros_like(breaks, dtype=np.int64)
     tol = 1e-15
 
     # limit case: code should decide wether point is in or out, not fall in infinite loop
@@ -26,15 +28,14 @@ def test_cell_index_p():
         assert status == expected_status
 
     # checking that the values match those of searchsorted (-1) for arbitrary grid points
-    i_grid = xp.array([0.14320482, 0.86569833, 0.77775327, 0.00895956, 0.074629  ,
+    i_grid = np.array([0.14320482, 0.86569833, 0.77775327, 0.00895956, 0.074629  ,
        0.45682646, 0.5384352 , 0.20915311, 0.73121977, 0.01057414,
        0.33756086, 0.17839759, 0.14023414, 0.09846206, 0.79970392,
        0.65330406, 0.82716552, 0.24185731, 0.24054685, 0.72466651,
        0.69125033, 0.3136558 , 0.64794089, 0.47975527, 0.99802844,
        0.64402598, 0.41263526, 0.28178414, 0.57274384, 0.73218562])
-    out = xp.zeros_like(i_grid, dtype=xp.int64)
+    out = np.zeros_like(i_grid, dtype=np.int64)
     status = cell_index_p(breaks, i_grid, tol, out)
     assert status == 0
-    nps = xp.searchsorted(breaks, i_grid)-1
-    assert xp.allclose(out, nps)
-
+    nps = np.searchsorted(breaks, i_grid)-1
+    assert np.allclose(out, nps)

feectools/core/tests/test_bsplines_pyccel.py

@@ -187,8 +187,9 @@ def basis_funs_all_ders_true(knots, degree, x, span, n, normalization='B'):
 
     # Normalization to get M-Splines
     if normalization == 'M':
-        ders *= [(degree + 1) / (knots[i + degree + 1] - knots[i]) \
-                 for i in range(span - degree, span + 1)]
+        scaling = xp.asarray([(degree + 1) / (knots[i + degree + 1] - knots[i])
+                              for i in range(span - degree, span + 1)])
+        ders *= scaling
     return ders
 
 #==============================================================================
@@ -221,7 +222,15 @@ def collocation_matrix_true(knots, degree, periodic, normalization, xgrid):
     for i,x in enumerate( xgrid ):
         span  =  find_span_true( knots, degree, x )
         basis = basis_funs_true( knots, degree, x, span )
-        mat[i,js(span)] = normalize(basis, span)
+        values = normalize(basis, span)
+        if periodic:
+            # NumPy and CuPy differ for indexed assignment with repeated
+            # indices (which occurs when nb <= degree).  The production
+            # kernel assigns in loop order, so make the reference explicit.
+            for j, value in zip(js(span), values):
+                mat[i, j] = value
+        else:
+            mat[i, js(span)] = values
 
     # Mitigate round-off errors
     mat[abs(mat) < 1e-14] = 0.0
@@ -293,7 +302,7 @@ def histopolation_matrix_true(knots, degree, periodic, normalization, xgrid):
     # Compute span for each row (index of last non-zero basis function)
     # TODO: would be better to have this ready beforehand
     # TODO: use tolerance instead of comparing against zero
-    spans = [(row != 0).argmax() + (degree+1) for row in C]
+    spans = [int((row != 0).argmax()) + (degree+1) for row in C]
 
     # Compute histopolation matrix from collocation matrix of higher degree
     m = C.shape[0] - 1