From f89e92053a7c15e5f682a8f283f9070ecbb1daf0 Mon Sep 17 00:00:00 2001 From: zagriadskov maksim Date: Thu, 16 Jul 2026 21:16:21 +0300 Subject: [PATCH] Initial commit --- include/HCSR.h | 104 ++++++++++++++ src/rvv/spmv_mtx_hcsr.cpp | 253 +++++++++++++++++++++++++++++++++++ src/scalar/spmv_mtx_hcsr.cpp | 123 +++++++++++++++++ 3 files changed, 480 insertions(+) create mode 100644 include/HCSR.h create mode 100644 src/rvv/spmv_mtx_hcsr.cpp create mode 100644 src/scalar/spmv_mtx_hcsr.cpp diff --git a/include/HCSR.h b/include/HCSR.h new file mode 100644 index 0000000..2600024 --- /dev/null +++ b/include/HCSR.h @@ -0,0 +1,104 @@ +#pragma once +#include +#include +#include +#include "CRS.h" +#include +#include "convert.h" +#include + +namespace SparseMatrixLib +{ + +template +struct inCRS { + size_t m = 0; + size_t n = 0; + size_t nz = 0; + std::vector row_ptr; // int16_t + std::vector col_idx; // int16_t + std::vector values; + + inCRS() = default; + inCRS(const inCRS& crs) = default; + inCRS& operator=(const inCRS& crs) = default; +}; + +template +class spMtxHCSR { +public: + size_t R; + size_t C; + size_t m = 0; // row + size_t n = 0; // column + size_t nz = 0; // non-zeros + size_t b_m; + size_t b_n; + std::vector block_ptr; + std::vector row_offset; + std::vector row_ptr; + std::vector col_idx; + std::vector values; + // my_container row_ptr; + // my_container col_idx; + // my_container> values; // can improve by adding bit-masks to elements + + spMtxHCSR() : R(0), C(0), m(0), n(0), nz(0), b_m(0), b_n(0), row_ptr(), col_idx(), values(), block_ptr() {} + + spMtxHCSR(size_t _m, size_t _n, size_t _nz, size_t _R, size_t _C) : m(_m), n(_n), nz(_nz), R(_R), C(_C) { + b_m = (m + R - 1) / R; + b_n = (n + C - 1) / C; + } + + spMtxHCSR(const spMtxHCSR& copy) : m(copy.m), n(copy.n), nz(copy.nz), b_m(copy.b_m), b_n(copy.b_n), + R(copy.R), C(copy.C), row_ptr(copy.row_ptr), col_idx(copy.col_idx), values(copy.values), block_ptr(copy.block_ptr) + {} + + spMtxHCSR(spMtxHCSR&& mov) : m(mov.m), n(mov.n), nz(mov.nz), b_m(mov.b_m), b_n(mov.b_n), + R(mov.R), C(mov.C), row_ptr(std::move(mov.row_ptr)), col_idx(std::move(mov.col_idx)), + values(std::move(mov.values)), block_ptr(std::move(mov.block_ptr)) { + + } + + ~spMtxHCSR() {} + + spMtxHCSR& operator=(const spMtxHCSR& copy) { + if (this == ©) + return *this; + + m = copy.m; + n = copy.n; + nz = copy.nz; + b_m = copy.b_m; + b_n = copy.b_n; + R = copy.R; + C = copy.C; + row_ptr = copy.row_ptr; + col_idx = copy.col_idx; + values = copy.values; + block_ptr = copy.block_ptr; + + return *this; + } + + spMtxHCSR& operator=(spMtxHCSR&& mov) { + if (this == &mov) + return *this; + + m = mov.m; + n = mov.n; + nz = mov.nz; + b_m = mov.b_m; + b_n = mov.b_n; + R = mov.R; + C = mov.C; + row_ptr = std::move(mov.row_ptr); + col_idx = std::move(mov.col_idx); + values = std::move(mov.values); + block_ptr = std::move(mov.block_ptr); + + return *this; + } +}; + +} diff --git a/src/rvv/spmv_mtx_hcsr.cpp b/src/rvv/spmv_mtx_hcsr.cpp new file mode 100644 index 0000000..207025b --- /dev/null +++ b/src/rvv/spmv_mtx_hcsr.cpp @@ -0,0 +1,253 @@ +#include "sparse_matrix.h" +#include "spmv_mtx.h" + +#include + +namespace SparseMatrixLib +{ + +inline void HCSR_spmv_ker_crs_double(double alpha, const uint32_t* row_ptr, const uint16_t* col_idx, const double* values, const double* b, double* y, int m) { + int vlmax = __riscv_vsetvlmax_e64m1(); + for (int i = 0; i < m; ++i) { + double tmp = 0.0; + int vl = vlmax; + vfloat64m1_t res = __riscv_vfmv_v_f_f64m1(0.0, vl); + + for (int j = row_ptr[i]; j < row_ptr[i + 1]; j += vlmax) { + vl = __riscv_vsetvl_e64m1(row_ptr[i + 1] - j); + vfloat64m1_t val = __riscv_vle64_v_f64m1(values + j, vl); + vuint16mf4_t index = __riscv_vle16_v_u16mf4(col_idx + j, vl); // reinterpret cast + vuint16mf4_t index_shift = __riscv_vsll_vx_u16mf4(index, 3, vl); + vfloat64m1_t b_ = __riscv_vloxei16_v_f64m1(b, index_shift, vl); + res = __riscv_vfmacc_vv_f64m1_tu(res, val, b_, vl); + } + vl = __riscv_vsetvl_e64m1(vlmax); + vfloat64m1_t sum = __riscv_vfmv_v_f_f64m1(0.0, 1); + sum = __riscv_vfredosum_vs_f64m1_f64m1(res, sum, vl); + tmp = __riscv_vfmv_f_s_f64m1_f64(sum); + + y[i] = alpha * tmp + y[i]; + } +} + +inline void HCSR_spmv_ker_crs_float(float alpha, const uint32_t* row_ptr, const uint16_t* col_idx, const float* values, const float* b, float* y, int m) { + int vlmax = __riscv_vsetvlmax_e32m1(); + for (int i = 0; i < m; ++i) { + float tmp = 0.0f; + int vl = vlmax; + vfloat32m1_t res = __riscv_vfmv_v_f_f32m1(0.0, vl); + + for (int j = row_ptr[i]; j < row_ptr[i + 1]; j += vlmax) { + vl = __riscv_vsetvl_e32m1(row_ptr[i + 1] - j); + vfloat32m1_t val = __riscv_vle32_v_f32m1(values + j, vl); + vuint16mf2_t index = __riscv_vle16_v_u16mf2(col_idx + j, vl); // reinterpret cast + vuint16mf2_t index_shift = __riscv_vsll_vx_u16mf2(index, 2, vl); + vfloat32m1_t b_ = __riscv_vloxei16_v_f32m1(b, index_shift, vl); + res = __riscv_vfmacc_vv_f32m1_tu(res, val, b_, vl); + } + vl = __riscv_vsetvl_e32m1(vlmax); + vfloat32m1_t sum = __riscv_vfmv_v_f_f32m1(0.0, 1); + sum = __riscv_vfredosum_vs_f32m1_f32m1(res, sum, vl); + tmp = __riscv_vfmv_f_s_f32m1_f32(sum); + + y[i] = alpha * tmp + y[i]; + } +} + +inline void HCSR_spmv_ker_coo_double(double alpha, const uint32_t* row_ptr, const uint16_t* col_idx, const double* values, const double* b, double* y, int m, int nz) { + for (int i = 0; i < nz; ++i) { + y[row_ptr[i]] += alpha * values[i] * b[col_idx[i]]; + } +} + +inline void HCSR_spmv_ker_coo_float(float alpha, const uint32_t* row_ptr, const uint16_t* col_idx, const float* values, const float* b, float* y, int m, int nz) { + for (int i = 0; i < nz; ++i) { + y[row_ptr[i]] += alpha * values[i] * b[col_idx[i]]; + } +} + +inline void HCSR_spmv_double(double alpha, const spMtxHCSR& mat, const std::vector& b, double beta, std::vector& y) { +#pragma omp parallel for schedule(dynamic) + for (int i = 0; i < mat.m; i += mat.R) { + int maxr = std::min(mat.R, mat.m - i); + for (int r = 0; r < maxr; ++r) + y[i + r] *= beta; + + for (int j = 0; j < mat.n; j += mat.C) { // for every block in row: + int block = j / mat.C + (i / mat.R) * mat.b_n; + int bptr = mat.block_ptr[block]; + uint32_t row_ptr_offset = mat.row_offset[block]; + + uint32_t mask = 0x3FFF'FFFFu; + uint32_t key = row_ptr_offset >> 30; + row_ptr_offset &= mask; + + int nz = mat.block_ptr[block + 1] - bptr; + if (key == 1) + HCSR_spmv_ker_crs_double(alpha, mat.row_ptr.data() + row_ptr_offset, mat.col_idx.data() + bptr, mat.values.data() + bptr, b.data() + j, y.data() + i, maxr); + else + HCSR_spmv_ker_coo_double(alpha, mat.row_ptr.data() + row_ptr_offset, mat.col_idx.data() + bptr, mat.values.data() + bptr, b.data() + j, y.data() + i, maxr, nz); + } + + } +} + +inline void HCSR_spmv_float(float alpha, const spMtxHCSR& mat, const std::vector& b, float beta, std::vector& y) { +#pragma omp parallel for schedule(dynamic) + for (int i = 0; i < mat.m; i += mat.R) { + int maxr = std::min(mat.R, mat.m - i); + for (int r = 0; r < maxr; ++r) + y[i + r] *= beta; + + for (int j = 0; j < mat.n; j += mat.C) { // for every block in row: + int block = j / mat.C + (i / mat.R) * mat.b_n; + int bptr = mat.block_ptr[block]; + uint32_t row_ptr_offset = mat.row_offset[block]; + + uint32_t mask = 0x3FFF'FFFFu; + uint32_t key = row_ptr_offset >> 30; + row_ptr_offset &= mask; + + int nz = mat.block_ptr[block + 1] - bptr; + if (key == 1) + HCSR_spmv_ker_crs_float(alpha, mat.row_ptr.data() + row_ptr_offset, mat.col_idx.data() + bptr, mat.values.data() + bptr, b.data() + j, y.data() + i, maxr); + else + HCSR_spmv_ker_coo_float(alpha, mat.row_ptr.data() + row_ptr_offset, mat.col_idx.data() + bptr, mat.values.data() + bptr, b.data() + j, y.data() + i, maxr, nz); + } + + } +} + +template<> +sparse_matrix_status sparse_mv( + sparse_operation_t type_op, + double alpha, + const spMtxHCSR &mat, + sparse_matrix_descr descr, + const std::vector &b, + double beta, + std::vector &y){ + sparse_matrix_status status; + + return status; +} + + +template<> +sparse_matrix_status sparse_mv( + sparse_operation_t type_op, + double alpha, + const spMtxHCSR &mat, + sparse_matrix_descr descr, + const std::vector &b, + double beta, + std::vector &y){ + sparse_matrix_status status; + + HCSR_spmv_double(alpha, mat, b, beta, y); + + return status; +} + +template<> +sparse_matrix_status sparse_mv( + sparse_operation_t type_op, + double alpha, + const spMtxHCSR &mat, + sparse_matrix_descr descr, + const std::vector &b, + double beta, + std::vector &y){ + sparse_matrix_status status; + + return status; +} + +template<> +sparse_matrix_status sparse_mv( + sparse_operation_t type_op, + double alpha, + const spMtxHCSR &mat, + sparse_matrix_descr descr, + const std::vector &b, + double beta, + std::vector &y){ + sparse_matrix_status status; + + sparse_mv + (type_op, alpha, mat, descr, b, beta, y); + sparse_mv + (type_op, alpha, mat, descr, b, beta, y); + sparse_mv + (type_op, alpha, mat, descr, b, beta, y); + + return status; +} + +template<> +sparse_matrix_status sparse_mv( + sparse_operation_t type_op, + float alpha, + const spMtxHCSR &mat, + sparse_matrix_descr descr, + const std::vector &b, + float beta, + std::vector &y){ + sparse_matrix_status status; + + return status; +} + +template<> +sparse_matrix_status sparse_mv( + sparse_operation_t type_op, + float alpha, + const spMtxHCSR &mat, + sparse_matrix_descr descr, + const std::vector &b, + float beta, + std::vector &y){ + sparse_matrix_status status; + + HCSR_spmv_float(alpha, mat, b, beta, y); + + return status; +} + +template<> +sparse_matrix_status sparse_mv( + sparse_operation_t type_op, + float alpha, + const spMtxHCSR &mat, + sparse_matrix_descr descr, + const std::vector &b, + float beta, + std::vector &y){ + sparse_matrix_status status; + + return status; +} + + +template<> +sparse_matrix_status sparse_mv( + sparse_operation_t type_op, + float alpha, + const spMtxHCSR &mat, + sparse_matrix_descr descr, + const std::vector &b, + float beta, + std::vector &y){ + sparse_matrix_status status; + + sparse_mv + (type_op, alpha, mat, descr, b, beta, y); + sparse_mv + (type_op, alpha, mat, descr, b, beta, y); + sparse_mv + (type_op, alpha, mat, descr, b, beta, y); + + return status; +} + +} \ No newline at end of file diff --git a/src/scalar/spmv_mtx_hcsr.cpp b/src/scalar/spmv_mtx_hcsr.cpp new file mode 100644 index 0000000..106267c --- /dev/null +++ b/src/scalar/spmv_mtx_hcsr.cpp @@ -0,0 +1,123 @@ +#include "sparse_matrix.h" + +#include +#include + +namespace SparseMatrixLib +{ + +inline void HCSR_spmv_ker_crs_double(double alpha, const uint32_t* row_ptr, const uint16_t* col_idx, const double* values, const double* b, double* y, int m) { + for (int i = 0; i < m; ++i) { + double tmp = 0.0; + for (int j = row_ptr[i]; j < row_ptr[i + 1]; ++j) + tmp += values[j] * b[col_idx[j]]; + y[i] += alpha * tmp; + } +} + +inline void HCSR_spmv_ker_crs_float(float alpha, const uint32_t* row_ptr, const uint16_t* col_idx, const float* values, const float* b, float* y, int m) { + for (int i = 0; i < m; ++i) { + float tmp = 0.0f; + for (int j = row_ptr[i]; j < row_ptr[i + 1]; ++j) + tmp += values[j] * b[col_idx[j]]; + y[i] += alpha * tmp; + } +} + +inline void HCSR_spmv_ker_coo_double(double alpha, const uint32_t* row_ptr, const uint16_t* col_idx, const double* values, const double* b, double* y, int m, int nz) { + for (int i = 0; i < nz; ++i) { + y[row_ptr[i]] += alpha * values[i] * b[col_idx[i]]; + } +} + +inline void HCSR_spmv_ker_coo_float(float alpha, const uint32_t* row_ptr, const uint16_t* col_idx, const float* values, const float* b, float* y, int m, int nz) { + for (int i = 0; i < nz; ++i) { + y[row_ptr[i]] += alpha * values[i] * b[col_idx[i]]; + } +} + +inline void HCSR_spmv_double(double alpha, const spMtxHCSR& mat, const std::vector& b, double beta, std::vector& y) { +#pragma omp parallel for + for (int i = 0; i < mat.m; i += mat.R) { + int maxr = std::min(mat.R, mat.m - i); + for (int r = 0; r < maxr; ++r) + y[i + r] *= beta; + + for (int j = 0; j < mat.n; j += mat.C) { // for every block in row: + int block = j / mat.C + (i / mat.R) * mat.b_n; + int bptr = mat.block_ptr[block]; + uint32_t row_ptr_offset = mat.row_offset[block]; + + uint32_t mask = ~(0xC000'0000u); + uint32_t key = row_ptr_offset >> 30; + row_ptr_offset &= mask; + + int nz = mat.block_ptr[block + 1] - bptr; + if (key == 1) + HCSR_spmv_ker_crs_double(alpha, mat.row_ptr.data() + row_ptr_offset, mat.col_idx.data() + bptr, mat.values.data() + bptr, b.data() + j, y.data() + i, maxr); + else + HCSR_spmv_ker_coo_double(alpha, mat.row_ptr.data() + row_ptr_offset, mat.col_idx.data() + bptr, mat.values.data() + bptr, b.data() + j, y.data() + i, maxr, nz); + } + + } +} + +inline void HCSR_spmv_float(float alpha, const spMtxHCSR& mat, const std::vector& b, float beta, std::vector& y) { +#pragma omp parallel for + for (int i = 0; i < mat.m; i += mat.R) { + int maxr = std::min(mat.R, mat.m - i); + for (int r = 0; r < maxr; ++r) + y[i + r] *= beta; + + for (int j = 0; j < mat.n; j += mat.C) { // for every block in row: + int block = j / mat.C + (i / mat.R) * mat.b_n; + int bptr = mat.block_ptr[block]; + uint32_t row_ptr_offset = mat.row_offset[block]; + + uint32_t mask = ~(0xC000'0000u); + uint32_t key = row_ptr_offset >> 30; + row_ptr_offset &= mask; + + int nz = mat.block_ptr[block + 1] - bptr; + if (key == 1) + HCSR_spmv_ker_crs_float(alpha, mat.row_ptr.data() + row_ptr_offset, mat.col_idx.data() + bptr, mat.values.data() + bptr, b.data() + j, y.data() + i, maxr); + else + HCSR_spmv_ker_coo_float(alpha, mat.row_ptr.data() + row_ptr_offset, mat.col_idx.data() + bptr, mat.values.data() + bptr, b.data() + j, y.data() + i, maxr, nz); + } + + } +} + +template<> +sparse_matrix_status sparse_mv( + sparse_operation_t type_op, + double alpha, + const spMtxHCSR &mat, + sparse_matrix_descr descr, + const std::vector &b, + double beta, + std::vector &y){ + sparse_matrix_status status; + +HCSR_spmv_double(alpha, mat, b, beta, y); + + return status; +} + +template<> +sparse_matrix_status sparse_mv( + sparse_operation_t type_op, + float alpha, + const spMtxHCSR &mat, + sparse_matrix_descr descr, + const std::vector &b, + float beta, + std::vector &y){ + sparse_matrix_status status; + +HCSR_spmv_float(alpha, mat, b, beta, y); + + return status; +} + +} \ No newline at end of file