[nsd/11arraydesign] Deprecate xtensor from the Laplace equation

.github/workflows/nsd.yml

@@ -202,7 +202,7 @@ jobs:
     - name: 11arraydesign_solve_cpp
       run: |
         make -C ${{ env.ARRAYDESIGNROOT }} solve_cpp.so
-        make -C ${{ env.ARRAYDESIGNROOT }} solve_cpp_simplearray.so
+        make -C ${{ env.ARRAYDESIGNROOT }} solve_cpp_xarray.so
 
     - name: 11arraydesign_data_prep
       # macos does not have MKL.

nsd/11arraydesign/arraydesign.rst

@@ -276,10 +276,9 @@ While numpy does speed up the naive Python loops, it has two drawbacks:
    method.
 
 Here we will see how much faster it can get to move the code from numpy to C++.
-To simplify the algorithm implementation, I use a library named `xtensor
-<http://xtensor.readthedocs.io/>`_.  It defines the multi-dimensional array
-data structure suitable for compile-time optimization.  An array library like
-that makes it easy to organize the code and achieve faster runtime.
+To simplify the implementation, I use an array library from `modmesh
+<https://github.com/solvcon/modmesh>`_.  It helps to organize the code and
+achieve faster runtime.
 
 Except the parentheses, the C++ version looks almost the same as the Python
 version.
@@ -1534,33 +1533,25 @@ elsewhere <https://github.com/solvcon/modmesh>`__):
 Exercises
 =========
 
-1. xtensor allows writing array-based code in C++, just like what numpy does
-   for Python.  Use xtensor to write array-based code in C++ by modifying the
-   C++ version of the point-Jacobi solver.  The array-based C++ version should
-   not have the inner loops.
-2. By allowing changing the signature of the function ``fit_poly()``, how can
+1. By allowing changing the signature of the function ``fit_poly()``, how can
    we ensure the shapes of ``xarr`` and ``yarr`` to be the same, without the
    explicit check with ``"xarr and yarr size mismatch"``?  Write code to show.
 
 References
 ==========
 
-.. [1] xtensor; multi-dimensional arrays with broadcasting and lazy computing:
-   https://xtensor.readthedocs.io
+.. [1] modmesh: https://github.com/solvcon/modmesh
 
-.. [2] xtensor-python; Python bindings for the xtensor C++ multi-dimensional
-   array library: https://xtensor-python.readthedocs.io
-
-.. [3] pybind11 Seamless operability between C++11 and Python:
+.. [2] pybind11 Seamless operability between C++11 and Python:
    https://pybind11.readthedocs.io/en/stable/
 
-.. [4] IPython / Jupyter integration for pybind11:
+.. [3] IPython / Jupyter integration for pybind11:
    https://github.com/aldanor/ipybind
 
-.. [5] Memory Layout Transformations:
+.. [4] Memory Layout Transformations:
    https://software.intel.com/content/www/us/en/develop/articles/memory-layout-transformations.html
 
-.. [6] Chapter 31 - Abstraction for AoS and SoA Layout in C++, *GPU Computing
+.. [5] Chapter 31 - Abstraction for AoS and SoA Layout in C++, *GPU Computing
    Gems Jade Edition*, pp. 429-441, Robert Strzodka, January 1, 2012.
    http://www.sciencedirect.com/science/article/pii/B9780123859631000319
 

nsd/11arraydesign/code/03_solve_cpp_simplearray.py → nsd/11arraydesign/code/03_solve_cpp_xarray.py

@@ -16,7 +16,7 @@ def make_grid():
 nx, x, uoriginal = make_grid()
 
 # [begin example]
-import solve_cpp_simplearray
+import solve_cpp_xarray
 # [end example]
 
 import time
@@ -30,15 +30,15 @@ def __exit__(self, *args, **kw):
 
 # [begin pycon]
 with Timer():
-    u, step, norm = solve_cpp_simplearray.solve_cpp(uoriginal)
+    u, step, norm = solve_cpp_xarray.solve_cpp(uoriginal)
 # [end pycon]
 
 from matplotlib import pyplot as plt
 
 def show_result(u, step, norm, size=7):
     print("step", step, "norm", norm)
     fig, ax = plt.subplots(figsize=(size,size))
-    cs = ax.contour(x, x, u.ndarray.T)
+    cs = ax.contour(x, x, u.T)
     ax.clabel(cs, inline=1, fontsize=10)
 
     ax.set_xticks(np.linspace(0,1,6))

nsd/11arraydesign/code/Makefile

@@ -67,25 +67,22 @@ testimport: libst
 data_prep.so: data_prep.cpp Makefile
 	g++ $< -o $@ -O3 -fPIC -shared -std=c++17 -lpython$(PYTHON_VERSION) -L$(PYTHON_LIB) ${INC} -I${MKLINC} ${MKLLINKLINE}
 
-solve_cpp.so: solve_cpp.cpp Makefile
-	g++ solve_cpp.cpp -o $@ -O3 -fPIC -shared -std=c++17 -lpython$(PYTHON_VERSION) -L$(PYTHON_LIB) ${INC}
+solve_cpp_xarray.so: solve_cpp_xarray.cpp Makefile
+	g++ $< -o $@ -O3 -fPIC -shared -std=c++17 -lpython$(PYTHON_VERSION) -L$(PYTHON_LIB) ${INC}
 
-../image/03_solve_cpp.png: 03_solve_cpp.py solve_cpp.so
+../image/03_solve_cpp_xarray.png: 03_solve_cpp_xarray.py solve_cpp_xarray.so
 	./$<
 
-solve_cpp_simplearray.so: \
-	solve_cpp_simplearray.cpp \
-	$(MODMESH_PYMOD)/wrap_SimpleArray.cpp \
-	$(MODMESH_PYMOD)/wrap_ConcreteBuffer.cpp \
-	$(MODMESH_PYMOD)/buffer_pymod.cpp \
-	Makefile
-	g++ -c solve_cpp_simplearray.cpp -o solve_cpp_simplearray.o -O3 -fPIC -std=c++17 ${INC}
-	g++ -c $(MODMESH_PYMOD)/wrap_ConcreteBuffer.cpp -o wrap_ConcreteBuffer.o -O3 -fPIC -std=c++17 ${INC}
-	g++ -c $(MODMESH_PYMOD)/wrap_SimpleArray.cpp -o wrap_SimpleArray.o -O3 -fPIC -std=c++17 ${INC}
-	g++ -c $(MODMESH_PYMOD)/buffer_pymod.cpp -o buffer_pymod.o -O3 -fPIC -std=c++17 ${INC}
-	g++ solve_cpp_simplearray.o wrap_ConcreteBuffer.o wrap_SimpleArray.o buffer_pymod.o -o $@ -shared -std=c++17 -lpython$(PYTHON_VERSION) -L$(PYTHON_LIB)
-
-../image/03_solve_cpp_simplearray.png: 03_solve_cpp_simplearray.py solve_cpp_simplearray.so
+%.o: $(MODMESH_PYMOD)/%.cpp Makefile
+	g++ -c $< -o $@ -O3 -fPIC -std=c++17 ${INC}
+
+solve_cpp.o: solve_cpp.cpp Makefile
+	g++ -c $< -o $@ -O3 -fPIC -std=c++17 ${INC}
+
+solve_cpp.so: solve_cpp.o wrap_SimpleArray.o wrap_ConcreteBuffer.o buffer_pymod.o Makefile
+	g++ solve_cpp.o wrap_ConcreteBuffer.o wrap_SimpleArray.o buffer_pymod.o -o $@ -shared -std=c++17 -lpython$(PYTHON_VERSION) -L$(PYTHON_LIB)
+
+../image/03_solve_cpp.png: 03_solve_cpp.py solve_cpp.so
 	./$<
 
 ../image/04_fit_poly.png: 04_fit_poly.py data_prep.so

nsd/11arraydesign/code/solve_cpp.cpp

@@ -1,49 +1,107 @@
 #include <pybind11/pybind11.h>
-#define FORCE_IMPORT_ARRAY
-#include <xtensor-python/pyarray.hpp>
+#include <modmesh/buffer/buffer.hpp>
+#include <modmesh/buffer/pymod/buffer_pymod.hpp>
 
 #include <vector>
 #include <algorithm>
 #include <tuple>
 #include <iostream>
 
-#include <xtensor/xarray.hpp>
-#include <xtensor/xadapt.hpp>
-#include <xtensor/xview.hpp>
+#include <pybind11/numpy.h>
+#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
+#include <numpy/arrayobject.h>
+
+namespace modmesh
+{
+
+namespace python
+{
+
+void import_numpy()
+{
+    auto local_import_numpy = []()
+    {
+        import_array2("cannot import numpy", false); // or numpy c api segfault.
+        return true;
+    };
+    if (!local_import_numpy())
+    {
+        throw pybind11::error_already_set();
+    }
+}
+
+} /* end namespace python */
+
+} /* end namespace modmesh */
+
+template <typename T>
+T calc_norm_amax(modmesh::SimpleArray<T> const & arr0, modmesh::SimpleArray<T> const & arr1)
+{
+    size_t const nelm = arr0.size();
+    T ret = 0;
+    for (size_t it = 0; it < nelm; ++it)
+    {
+        T const val = std::abs(arr0[it] - arr1[it]);
+        if (val > ret)
+        {
+            ret = val;
+        }
+    }
+    return ret;
+}
 
 // [begin example]
-std::tuple<xt::xarray<double>, size_t, double>
-solve1(xt::xarray<double> u)
+std::tuple<modmesh::SimpleArray<double>, size_t, double>
+solve1(modmesh::SimpleArray<double> u)
 {
     const size_t nx = u.shape(0);
-    xt::xarray<double> un = u;
+    modmesh::SimpleArray<double> un = u;
     bool converged = false;
     size_t step = 0;
     double norm;
     while (!converged)
     {
+        norm = 0.0;
         ++step;
         for (size_t it=1; it<nx-1; ++it)
         {
             for (size_t jt=1; jt<nx-1; ++jt)
             {
                 un(it,jt) = (u(it+1,jt) + u(it-1,jt) + u(it,jt+1) + u(it,jt-1)) / 4;
+                double const v = std::abs(un(it,jt) - u(it,jt));
+                if (v > norm)
+                {
+                    norm = v;
+                }
             }
         }
-        norm = xt::amax(xt::abs(un-u))();
+        // additional loop slows down: norm = calc_norm_amax(u, un);
         if (norm < 1.e-5) { converged = true; }
         u = un;
     }
     return std::make_tuple(u, step, norm);
 }
 
+// Expose the C++ helper to Python
 PYBIND11_MODULE(solve_cpp, m)
 {
-    xt::import_numpy();
+    // Boilerplate for using the SimpleArray with Python
+    {
+        modmesh::python::import_numpy();
+        modmesh::python::wrap_ConcreteBuffer(m);
+        modmesh::python::wrap_SimpleArray(m);
+    }
     m.def
     (
         "solve_cpp"
-      , [](xt::pyarray<double> & uin) { return solve1(xt::xarray<double>(uin)); }
+      , [](pybind11::array_t<double> & uin)
+        {
+            auto ret = solve1(modmesh::python::makeSimpleArray(uin));
+            return std::make_tuple(
+                modmesh::python::to_ndarray(std::get<0>(ret)),
+                std::get<1>(ret),
+                std::get<2>(ret));
+        }
     );
 }
 // [end example]