final year project：C++手写numpy并移植到RISC-V上——纪念我在中科院实习 ...

  我毕设做的项目是用C++去实现一个Numpy，因为我是大数据专业，Numpy又是跟数据分析有关的工具，所以我计划自己动手去实现一个小型的Numpy，现在代码规模大概在六千多行左右，并且可以乐成移植到OpenEuler RISC-V上面。在这个项目当中，我实现了比较多的数学函数，并且用到了各种高性能有关的技术，如：SIMD，OpenMP，OpenBlas，分别用来做数学运算的加速，向量化循环以及矩阵运算加速，我首先是在x86架构下完成了项目标大部分，后面才移植到了RISC-V上面，现在RISC-V有关的优化只有OpenMP以及OpenBlas，这两者都已经在OpenEuler RISC-V上乐成移植并且可以乐成运行。至于RVV指令集，现在该操作系统似乎还不太支持，所以如果偶尔间后面我再别的想办法。

1. template <typename T>
2. class ndarray {
3. private:
4.     std::vector<T> __data;
5.     std::vector<size_t> __shape;
6.     std::vector<size_t> __strides;
7.     size_t __size;
9.     void compute_strides();
11.     size_t calculate_offset(size_t row, size_t col) const noexcept;
12.     ...   
13. }

复制代码

  项目标数据结构大概如上，__data用于存储实际的数据，不管是一维还是二维的数据，都存在__data内里，至于二维怎么存，可以通过__strides数组结合calculate_offset去映射到一维数组上面；__shape是存储形状的数据结构，如果数组是一维的并且有五个元素，那么__shape就是{5}，如果是二乘三的数组，那么__shape就是{2, 3}；__size则是存储数组实际元素数目，也就是__data字段的巨细。
  

  这是项目标结构，logical部分包罗了位运算的实当代码，math部分包罗了数学函数的实现，shift部分包罗了位移的部分，matrix_operations部分包罗了矩阵算法的部分，sort是排序部分，parallel_for是向量化for循环的部分。而simd_traits.cpp和xsimd_traits.cpp则负责编译期萃取类型，因为SIMD函数有许多类型，我将它封装成了一个结构体，内里各类型的函数都封装成一个统一接口，如许就能减少许多重复代码。

  这里则是单元测试部分，我现在使用了GoogleTest对其进行了单元测试，分别对一维数组和二维数组的各种数学函数，并且也有异常相关的测试。该项目提供了CMake和Meson的构建方式，用户可以一键式构建，只需要提前下载gcc、GoogleTest、OpenBlas以及xsimd即可，SIMD指令集以及OpenMP套件一般在GNU工具链内里有。
  

1. template <typename T>
2.     std::vector<T> add1(const std::vector<T>& A, const std::vector<T>& B) {
3.         static_assert(std::is_arithmetic_v<T>, "Type must be arithmetic");
4.         static_assert(!std::is_same_v<T, char>);
6.         if (A.size() != B.size()) {
7.             throw std::invalid_argument("Vector dimension mismatch");
8.         }
10.         size_t N = A.size();
11.         std::vector<T> C(N);
13.         if constexpr (std::is_same_v<T, float>) {
14.             cblas_scopy(N, B.data(), 1, C.data(), 1);
15.             cblas_saxpy(N, 1.0, A.data(), 1, C.data(), 1);
16.         } else if constexpr (std::is_same_v<T, double>) {
17.             cblas_dcopy(N, B.data(), 1, C.data(), 1);
18.             cblas_daxpy(N, 1.0, A.data(), 1, C.data(), 1);
19.         } else {
20.             std::vector<float> float_A(N), float_B(N), float_C(N);
21.             for (size_t i = 0; i < N; ++i) {
22.                 float_A[i] = static_cast<float>(A[i]);
23.                 float_B[i] = static_cast<float>(B[i]);
24.             }
25.             cblas_scopy(N, float_B.data(), 1, float_C.data(), 1);
26.             cblas_saxpy(N, 1.0, float_A.data(), 1, float_C.data(), 1);
27.             for (size_t i = 0; i < N; ++i) {
28.                 C[i] = static_cast<T>(float_C[i]);
29.             }
30.         }
32.         return C;
33.     }

复制代码

  这里是向量加法部分，因为OpenBlas的函数接口有单精度和双精度类型的，所以我用了编译期条件判断去选择相应的函数实现，对于非浮点数类型，则将其转换为浮点类型再调用Blas的接口（因为Blas实在是快！），但实在这里另有优化的点，后面偶尔间再想想。

1. template <typename T>
2. struct round_simd_traits;
4. template <>
5. struct round_simd_traits<float> {
6.     using scalar_type = float;
7.     using simd_type = __m256;
8.     static constexpr size_t step = 8;
10.     static simd_type load(const scalar_type *ptr) noexcept {
11.         return _mm256_loadu_ps(ptr);
12.     }
14.     static void store(scalar_type *ptr, simd_type val) noexcept {
15.         _mm256_storeu_ps(ptr, val);
16.     }
18.     static simd_type op(simd_type a) noexcept {
19.         return _mm256_round_ps(a, _MM_FROUND_TO_NEAREST_INT);
20.     }
21. };
23. template <>
24. struct round_simd_traits<double> {
25.     using scalar_type = double;
26.     using simd_type = __m256d;
27.     static constexpr size_t step = 4;
29.     static simd_type load(const scalar_type *ptr) noexcept {
30.         return _mm256_loadu_pd(ptr);
31.     }
33.     static void store(scalar_type *ptr, simd_type val) noexcept {
34.         _mm256_storeu_pd(ptr, val);
35.     }
37.     static simd_type op(simd_type a) noexcept {
38.         return _mm256_round_pd(a, _MM_FROUND_TO_NEAREST_INT);
39.     }
40. };

复制代码

  这里是SIMD萃取的过程，如许子外部调用的时候就可以直接统一用OP去调用，省去了函数重载的过程。

1. template <typename T>
2.     std::vector<T> acos1_simd(const std::vector<T>& A) {
3.         if (A.size() < 32)
4.             return apply_unary_op_plain(A, [](const T& a) {
5.                 return std::acos(a);
6.             });
8.         #ifdef __riscv
9.             return apply_unary_op_plain(A, [](const T& a) {
10.                 return std::acos(a);
11.             });
12.         #endif
14.         #ifdef __AVX2__
15.             return apply_unary_op_simd<T, acos_simd_traits<T>>(A, [](const T& a) {
16.                 return std::acos(a);
17.             });
18.         #endif
19.     }

复制代码

  这里是math函数部分，对于数据规模较小的，直接调用朴素循环，对于RISC-V架构，现在我的做法也是直接调用朴素循环（还可以优化），如果是支持AVX2的环境，则可以调用SIMD版本的函数。

1. if(CMAKE_SYSTEM_PROCESSOR MATCHES "x86_64|AMD64|i386|i686")
2.     add_definitions(-mavx2 -fopenmp -O3)
3. elseif(CMAKE_SYSTEM_PROCESSOR MATCHES "riscv64")
4.     add_definitions(-fopenmp -march=rv64gcv -O3)
5. endif()

复制代码

  像CMake部分，经过我的测试，在O3优化的情况下性能最高，所以我开了O3优化。另有一点区别就是RISC-V需要开V扩展，而x86开-mavx2用来支持AVX2指令集。

1. name: CI
3. on:
4.   push:
5.     branches:
6.       - main
7.   pull_request:
8.     branches:
9.       - main
10.   workflow_dispatch:
13. jobs:
14.   build-and-test:
15.     runs-on: ${{ matrix.os }}
16.     strategy:
17.       matrix:
18.         os: [ubuntu-latest]
20.     steps:
21.       - name: Checkout code
22.         uses: actions/checkout@v2
24.       - name: Install system dependencies
25.         run: |
26.           if [[ "${{ matrix.os }}" == "ubuntu-latest" ]]; then
27.             sudo apt-get update
28.             sudo apt-get install -y cmake build-essential libopenblas-dev libxsimd-dev libboost-all-dev
29.           elif [[ "${{ matrix.os }}" == "fedora-latest" ]]; then
30.             sudo dnf update -y
31.             sudo dnf install -y cmake make gcc-c++ openblas-devel xsimd-devel boost-devel
32.           fi
33.       
34.       - name: Download and install Google Test
35.         run: |
36.           git clone https://github.com/google/googletest.git
37.           cd googletest
38.           mkdir build
39.           cd build
40.           cmake ..
41.           make
42.           sudo make install
44.       - name: Create build directory
45.         run: mkdir -p build
47.       - name: Configure CMake
48.         working-directory: build
49.         run: cmake ..
51.       - name: Build project
52.         working-directory: build
53.         run: make -j$(nproc)
55.       - name: Run tests
56.         working-directory: build/test
57.         run: ./run_all_tests

复制代码

  别的项目当中还提供了Ubuntu最新版的CI/CD，也就是Ubuntu 24.04，项目已开源，偶尔间我会多多commit，欢迎关注~
  项目地点：https://github.com/Thomas134/numpy_project.git

免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。