Bausager/Flux
1
0
Fork 0
Code Issues 36 Pull Requests Actions Packages Projects 1 Releases Wiki Activity

refactor
Sync public mirror / sync (push) Failing after 29s
Details

Browse Source 
almost complete. Need to doublecheck names for functions in *_serial.h
This commit is contained in:
Michelle
2026-01-18 17:51:05 +01:00
parent 9709949a5b
commit 6b8a7ab582
15 changed files with 383 additions and 523 deletions
Download Patch File Download Diff File Expand all files Collapse all files
include/numerics/detail/dot_serial.h
+42
View File
@@ -0,0 +1,42 @@
#pragma once
#include <cstdint> //uint64_t
#include <stdexcept> // std::runtime_error
#include "../utils/vector.h"
#include "../utils/matrix.h"
namespace numerics::detail{
// ---------------- Vector * Vector ----------------
template <typename T>
inline T dot_serial(const utils::Vector<T>& v, const utils::Vector<T>& p) {
const uint64_t N = v.size();
if (N != p.size()) {
throw std::runtime_error("dot_serial: dimension mismatch");
}
T dot_product = T{0};
for (uint64_t i = 0; i < N; ++i){
dot_product += v[i]*p[i];
}
return dot_product;
}
// ---------------- Matrix * Vector ----------------
template <typename T>
inline utils::Vector<T> dot_serial(const utils::Matrix<T>& A, const utils::Vector<T>& v) {
const uint64_t rows = A.rows();
const uint64_t cols = A.cols();
if (cols != v.size()) {
throw std::runtime_error("dot_serial: dimension mismatch");
}
utils::Vector<T> dot_product(rows, T{0});
for (uint64_t j = 0; j < cols; ++j){
for (uint64_t i = 0; i < rows; ++i){
dot_product[i] += A(i,j)*v[j];
}
}
return dot_product;
}
} // namespace numerics
include/numerics/detail/equal_serial.h
+47
View File
@@ -0,0 +1,47 @@
#pragma once
#include <cstdint> //uint64_t
//#include <stdexcept> // std::runtime_error
#include "../utils/vector.h"
#include "../utils/matrix.h"
namespace numerics::detail{
// ---------------- Matrix ----------------
template <typename T>
inline bool equal_serial(const utils::Matrix<T>& A, const utils::Matrix<T> & B) {
const uint64_t rows = A.rows();
const uint64_t cols = A.cols();
if ((rows != B.rows()) || (cols != B.cols())){
return false;
}
for (uint64_t i = 0; i < rows; ++i){
for (uint64_t j = 0; j < cols; ++j){
if (A(i,j) != B(i,j)){
return false;
}
}
}
return true;
}
// ---------------- Vector ----------------
template <typename T>
inline bool equal_serial(const utils::Vector<T>& v, const utils::Vector<T>& p) {
const uint64_t N = v.size();
if (N != p.size()){
return false;
}
for (uint64_t i = 0; i < N; ++i){
if ((v[i] != p[i])){
return false;
}
}
return true;
}
} // namespace numerics
include/numerics/detail/isclose_serial.h
+49
View File
@@ -0,0 +1,49 @@
#pragma once
#include <cstdint> //uint64_t
//#include <stdexcept> // std::runtime_error
#include "../utils/vector.h"
#include "../utils/matrix.h"
#include <cmath> // std::abs
namespace numerics::detail{
// ---------------- Matrix ----------------
template <typename T>
inline bool isclose_serial(const utils::Matrix<T>& A, const utils::Matrix<T> & B, const T tol = T{1e-5}) {
const uint64_t rows = A.rows();
const uint64_t cols = A.cols();
if ((rows != B.rows()) || (cols != B.cols())){
return false;
}
for (uint64_t i = 0; i < rows; ++i){
for (uint64_t j = 0; j < cols; ++j){
if (static_cast<T>(std::abs(A(i,j)-B(i,j))) > tol){
return false;
}
}
}
return true;
}
// ---------------- Vector ----------------
template <typename T>
inline bool isclose_serial(const utils::Vector<T>& v, const utils::Vector<T>& p, const T tol = T{1e-5}) {
const uint64_t N = v.size();
if (N != p.size()){
return false;
}
for (uint64_t i = 0; i < N; ++i){
if (static_cast<T>(std::abs((v[i]-p[i]))) > tol){
return false;
}
}
return true;
}
} // namespace numerics
include/numerics/detail/matmul_serial.h
+35
View File
@@ -0,0 +1,35 @@
#pragma once
#include <cstdint> //uint64_t
#include <stdexcept> // std::runtime_error
#include "../utils/vector.h"
#include "../utils/matrix.h"
namespace numerics::detail{
// ---------------- Matrix * Matrix ----------------
template <typename T>
inline utils::Matrix<T> matmul_serial(const utils::Matrix<T>& A, const utils::Matrix<T>& B){
const uint64_t m = A.rows();
const uint64_t n = A.cols(); // also B.rows()
const uint64_t p = B.cols();
if(n != B.rows()){
throw std::runtime_error("matmul: dimension mismatch");
}
T tmp;
utils::Matrix<T> C(m, p, T{0});
for (uint64_t i = 0; i < m; ++i){
for (uint64_t j = 0; j < n; ++j){
tmp = A(i,j);
for (uint64_t k = 0; k < p; ++k){
C(i,k) += tmp * B(j,k);
}
}
}
return C;
}
} // namespace numerics
include/numerics/detail/random_serial.h
+63
View File
@@ -0,0 +1,63 @@
#pragma once
#include <cstdint> //uint64_t
//#include <stdexcept> // std::runtime_error
#include <random>
#include <type_traits>
#include "../utils/vector.h"
#include "../utils/matrix.h"
namespace numerics::detail{
// Shared engine
inline std::mt19937& rng() {
static std::random_device rd;
static std::mt19937 gen(rd());
return gen;
}
// Integral overload
template <
typename T,
typename std::enable_if<std::is_integral<T>::value, int>::type = 0
>
inline T random(const T low, const T high) {
std::uniform_int_distribution<T> dist(low, high);
return dist(rng());
}
// Floating-point overload
template <
typename T,
typename std::enable_if<std::is_floating_point<T>::value, int>::type = 0
>
inline T random(const T low, const T high) {
std::uniform_real_distribution<T> dist(low, high);
return dist(rng());
}
// ---------------- Matrix ----------------
template <typename T>
inline void inplace_random_serial(utils::Matrix<T>& A, const T low, const T high) {
const uint64_t rows = A.rows();
const uint64_t cols = A.cols();
for (uint64_t i = 0; i < rows; ++i){
for (uint64_t j = 0; j < cols; ++j){
A(i,j) = random(low, high);
}
}
}
// ---------------- Vector ----------------
template <typename T>
inline void inplace_random_serial(utils::Vector<T>& v, const T low, const T high) {
const uint64_t N = v.size();
for (uint64_t i = 0; i < N; ++i){
v[i] = random(low, high);
}
}
} // namespace numerics
include/numerics/dot.h
+21
View File
@@ -0,0 +1,21 @@
#pragma once
#include "./core/omp_config.h"
#include "detail/dot_serial.h"
namespace numerics{
// ---------------- Vector * vector ----------------
template <typename T>
inline T dot(const utils::Vector<T>& v, const utils::Vector<T>& p) {
return detail::dot_serial(v,p);
}
// ---------------- Matrix * vector ----------------
template <typename T>
inline utils::Vector<T> dot(const utils::Matrix<T>& A, const utils::Vector<T>& v) {
return detail::dot_serial(A,v);
}
}
include/numerics/equal.h
+22
View File
@@ -0,0 +1,22 @@
#pragma once
#include "./core/omp_config.h"
#include "detail/equal_serial.h"
namespace numerics{
// ---------------- equal ----------------
template <typename T>
inline bool equal(const utils::Vector<T>& v, const utils::Vector<T>& p) {
return detail::equal_serial(v, p);
}
template <typename T>
inline bool equal(const utils::Matrix<T>& A, const utils::Matrix<T>& B) {
return detail::equal_serial(A, B);
}
}
include/numerics/isclose.h
View File
include/numerics/matdot.h
-63
View File
@@ -1,63 +0,0 @@
#pragma once
#include "./core/omp_config.h"
#include "./utils/matrix.h"
#include "./utils/vector.h"
namespace numerics{
template <typename T>
utils::Vector<T> matdot_row(const utils::Matrix<T>& A, const utils::Matrix<T>& B){
if (A.rows() != B.rows() || A.cols() != B.cols()){
throw std::runtime_error("matmul: dimension mismatch");
}
const uint64_t rows = A.rows();
const uint64_t cols = A.cols();
utils::Vector<T> c(rows, T{0});
for (uint64_t i = 0; i < rows; ++i){
T sum = T{0};
for (uint64_t j = 0; j < cols; ++j){
sum += A(i,j) * A(i,j);
}
c[i] = sum;
}
return c;
}
template <typename T>
utils::Vector<T> matdot_col(const utils::Matrix<T>& A, const utils::Matrix<T>& B){
if (A.rows() != B.rows() || A.cols() != B.cols()){
throw std::runtime_error("matmul: dimension mismatch");
}
const uint64_t rows = A.rows();
const uint64_t cols = A.cols();
utils::Vector<T> c(cols, T{0});
for (uint64_t j = 0; j < cols; ++j){
T sum = T{0};
for (uint64_t i = 0; i < rows; ++i){
sum += A(i,j) * A(i,j);
}
c[j] = sum;
}
return c;
}
} // namespace numerics
include/numerics/matequal.h
-85
View File
@@ -1,85 +0,0 @@
#pragma once
#include "./core/omp_config.h"
#include "./utils/matrix.h"
#include "./numerics/abs.h"
namespace numerics{
// -------------- Serial ----------------
template <typename T>
bool matequal(const utils::Matrix<T>& A, const utils::Matrix<T>& B, double tol = 1e-9) {
if (A.rows() != B.rows() || A.cols() != B.cols()) {
return false;
}
bool decimal = std::is_floating_point<T>::value;
const uint64_t rows=A.rows(), cols=A.cols();
for (uint64_t i = 0; i < rows; ++i)
for (uint64_t j = 0; j < cols; ++j)
if (decimal) {
if (numerics::abs(A(i,j) - B(i,j)) > static_cast<T>(tol)){
return false;
}
} else {
if (A(i,j) != B(i,j)){
return false;
}
}
return true;
}
// -------------- Parallel ----------------
template <typename T>
bool matequal_omp(const utils::Matrix<T>& A, const utils::Matrix<T>& B, double tol = 1e-9) {
if (A.rows() != B.rows() || A.cols() != B.cols()) {
return false;
}
bool decimal = std::is_floating_point<T>::value;
bool eq = true;
const uint64_t rows=A.rows(), cols=A.cols();
#pragma omp parallel for collapse(2) schedule(static) reduction(&&:eq)
for (uint64_t i = 0; i < rows; ++i)
for (uint64_t j = 0; j < cols; ++j)
if (decimal) {
eq = eq && (numerics::abs(A(i,j) - B(i,j)) <= static_cast<T>(tol));
} else {
eq = eq && (A(i,j) == B(i,j));
}
return eq;
}
// -------------- Auto OpenMP ----------------
template <typename T>
bool matequal_auto(const utils::Matrix<T>& A, const utils::Matrix<T>& B, double tol = 1e-9) {
if (A.rows() != B.rows() || A.cols() != B.cols()) {
return false;
}
uint64_t work = A.rows() * A.cols();
#ifdef _OPENMP
bool can_parallel = omp_config::omp_parallel_allowed();
uint64_t threads = static_cast<uint64_t>(omp_get_max_threads());
#else
bool can_parallel = false;
uint64_t threads = 1;
#endif
if (can_parallel || work > threads * 4ull) {
return matequal_omp(A,B,tol);
}
else{
// Safe fallback
return matequal(A,B,tol);
}
}
} // namespace numerics
include/numerics/matmul.h
+5 -114
View File
@@ -1,124 +1,15 @@
#ifndef _matmul_n_
#define _matmul_n_
#pragma once
#include "./utils/matrix.h"
#include "./core/omp_config.h"
#include "detail/matmul_serial.h"
namespace numerics{
// ---------------- Serial baseline ----------------
// ---------------- matmul ----------------
template <typename T>
utils::Matrix<T> matmul(const utils::Matrix<T>& A, const utils::Matrix<T>& B){
if(A.cols() != B.rows()){
throw std::runtime_error("matmul: dimension mismatch");
}
const uint64_t m = A.rows();
const uint64_t n = A.cols(); // also B.rows()
const uint64_t p = B.cols();
T tmp;
utils::Matrix<T> C(m, p, T{0});
for (uint64_t i = 0; i < m; ++i){
for (uint64_t j = 0; j < n; ++j){
tmp = A(i,j);
for (uint64_t k = 0; k < p; ++k){
C(i,k) += tmp * B(j,k);
}
}
}
return C;
}
// ---------------- Rows-only OpenMP ----------------
template <typename T>
utils::Matrix<T> matmul_rows_omp(const utils::Matrix<T>& A,
const utils::Matrix<T>& B) {
if (A.cols() != B.rows()) throw std::runtime_error("matmul_rows_omp: dim mismatch");
const uint64_t m=A.rows(), n=A.cols(), p=B.cols();
utils::Matrix<T> C(m, p, T{0});
#pragma omp parallel for schedule(static)
for (uint64_t i=0;i<m;++i) {
for (uint64_t j=0;j<p;++j) {
T acc=T{0};
for (uint64_t k=0;k<n;++k) {
acc += A(i,k)*B(k,j);
}
C(i,j)=acc;
}
inline utils::Matrix<T> matmul(const utils::Matrix<T>& A, const utils::Matrix<T>& B){
return detail::matmul_serial(A, B);
}
return C;
}
// -------------- Collapse(2) OpenMP ----------------
template <typename T>
utils::Matrix<T> matmul_collapse_omp(const utils::Matrix<T>& A,
const utils::Matrix<T>& B) {
if (A.cols() != B.rows()) throw std::runtime_error("matmul_collapse_omp: dim mismatch");
const uint64_t m=A.rows(), n=A.cols(), p=B.cols();
utils::Matrix<T> C(m, p, T{0});
#pragma omp parallel for collapse(2) schedule(static)
for (uint64_t i=0;i<m;++i) {
for (uint64_t j=0;j<p;++j) {
T acc=T{0};
for (uint64_t k=0;k<n;++k){
acc += A(i,k)*B(k,j);
}
C(i,j)=acc;
}
}
return C;
}
// -------------------- Auto selector ---------------------
template <typename T>
utils::Matrix<T> matmul_auto(const utils::Matrix<T>& A,
const utils::Matrix<T>& B) {
const uint64_t m=A.rows(), p=B.cols();
const uint64_t work = m * p;
#ifdef _OPENMP
bool can_parallel = omp_config::omp_parallel_allowed();
uint64_t threads = static_cast<uint64_t>(omp_get_max_threads());
#else
bool can_parallel = false;
uint64_t threads = 1;
#endif
// Tiny problems: serial is cheapest.
if (!can_parallel || work < threads*4ull) {
return matmul(A,B);
}
// Plenty of (i,j) work → collapse(2) is a great default.
else if (work >= 8ull * threads) {
return matmul_collapse_omp(A,B);
}
// Many rows and very few columns → rows-only cheaper overhead.
else if (m >= static_cast<uint64_t>(threads) && p <= 4) {
return matmul_rows_omp(A,B);
}
else{
// Safe fallback
return matmul(A,B);
}
}
} // namespace numerics
#endif // _matmul_n_
include/numerics/matrandom.h
-55
View File
@@ -1,55 +0,0 @@
#ifndef _matrandom_n_
#define _matrandom_n_
#include "./utils/vector.h"
#include "./utils/matrix.h"
#include "./core/omp_config.h"
namespace numerics{
template <typename T>
void inplace_matrandom_add(utils::Matrix<T>& A, const T lower, const T higher) {
const uint64_t rows = A.rows();
const uint64_t cols = A.cols();
utils::Matrix<T> B;
B.random(rows,cols, lower, higher);
numerics::inplace_matadd_mat(A, B);
}
template <typename T>
void inplace_matrandom_mul(utils::Matrix<T>& A, const T lower, const T higher) {
const uint64_t rows = A.rows();
const uint64_t cols = A.cols();
utils::Matrix<T> B;
B.random(rows,cols, lower, higher);
for (uint64_t i = 0; i < rows; ++i){
for (uint64_t j = 0; j < cols; ++j){
A(i,j) *= B(i,j);
}
}
}
template <typename T>
utils::Matrix<T> matrandom_add(const utils::Matrix<T>& A, const T lower, const T higher) {
const uint64_t rows = A.rows();
const uint64_t cols = A.cols();
utils::Matrix<T> B = A;
numerics::inplace_matadd_mat(B, lower, higher);
return B;
}
} // namespace numerics
#endif // _matrandom_n_
include/numerics/matvec.h
-156
View File
@@ -1,156 +0,0 @@
#ifndef _matvec_n_
#define _matvec_n_
#include "./utils/matrix.h"
#include "./core/omp_config.h"
namespace numerics{
// =================================================
// y = A * x (Matrix–Vector product)
// =================================================
template <typename T>
utils::Vector<T> matvec(const utils::Matrix<T>& A, const utils::Vector<T>& x) {
if (A.cols() != x.size()) {
throw std::runtime_error("matvec: dimension mismatch");
}
const uint64_t m = A.rows();
const uint64_t n = A.cols();
utils::Vector<T> y(m, T{0});
for (uint64_t i = 0; i < m; ++i) {
for (uint64_t j = 0; j < n; ++j) {
y[i] += A(i, j) * x[j];
}
}
return y;
}
// -------------- Collapse(2) OpenMP ----------------
template <typename T>
utils::Vector<T> matvec_omp(const utils::Matrix<T>& A, const utils::Vector<T>& x) {
if (A.cols() != x.size()) {
throw std::runtime_error("matvec: dimension mismatch");
}
const uint64_t m = A.rows();
const uint64_t n = A.cols();
utils::Vector<T> y(m, T{0}); // <-- y has length m (rows)
const T* xptr = x.data();
const T* Aptr = A.data(); // row-major: A(i,j) == Aptr[i*n + j]
// Each row i is an independent dot product: y[i] = dot(A[i,*], x)
#pragma omp parallel for schedule(static)
for (uint64_t i = 0; i < m; ++i) {
const T* row = Aptr + i * n; // contiguous row i
T acc = T{0};
#pragma omp simd reduction(+:acc)
for (uint64_t j = 0; j < n; ++j) {
acc += row[j] * xptr[j];
}
y[i] = acc;
}
return y;
}
// -------------- Auto OpenMP ----------------
template <typename T>
utils::Vector<T> matvec_auto(const utils::Matrix<T>& A,
const utils::Vector<T>& x) {
uint64_t work = A.rows() * A.cols();
bool can_parallel = omp_config::omp_parallel_allowed();
#ifdef _OPENMP
int threads = omp_get_max_threads();
#else
int threads = 1;
#endif
if (can_parallel || work > static_cast<uint64_t>(threads) * 4ull) {
return matvec_omp(A,x);
}
else{
// Safe fallback
return matvec(A,x);
}
}
// =================================================
// y = x * A (Vector–Matrix product)
// =================================================
template <typename T>
utils::Vector<T> vecmat(const utils::Vector<T>& x, const utils::Matrix<T>& A) {
if (x.size() != A.rows()) {
throw std::runtime_error("vecmat: dimension mismatch");
}
const uint64_t m = A.rows();
const uint64_t n = A.cols();
utils::Vector<T> y(n, T{0});
for (uint64_t j = 0; j < n; ++j) {
for (uint64_t i = 0; i < m; ++i) {
y[j] += x[i] * A(i, j);
}
}
return y;
}
// -------------- Collapse(2) OpenMP ----------------
template <typename T>
utils::Vector<T> vecmat_omp(const utils::Vector<T>& x, const utils::Matrix<T>& A) {
if (x.size() != A.rows()) {
throw std::runtime_error("vecmat: dimension mismatch");
}
const uint64_t m = A.rows();
const uint64_t n = A.cols();
utils::Vector<T> y(n, T{0});
#pragma omp parallel for schedule(static)
for (uint64_t j = 0; j < n; ++j) {
T acc = T{0};
for (uint64_t i = 0; i < m; ++i) {
acc += x[i] * A(i, j);
}
y[j] = acc;
}
return y;
}
// -------------- Auto OpenMP ----------------
template <typename T>
utils::Vector<T> vecmat_auto(const utils::Vector<T>& x,
const utils::Matrix<T>& A) {
uint64_t work = A.rows() * A.cols();
bool can_parallel = omp_config::omp_parallel_allowed();
#ifdef _OPENMP
int threads = omp_get_max_threads();
#else
int threads = 1;
#endif
if (can_parallel || work > static_cast<uint64_t>(threads) * 4ull) {
return vecmat_omp(x,A);
}
else{
// Safe fallback
return vecmat(x,A);
}
}
} // namespace numerics
#endif // _matvec_n_
include/numerics/random.h
+97
View File
@@ -0,0 +1,97 @@
#pragma once
#include "./core/omp_config.h"
#include "detail/random_serial.h"
#include "add.h"
#include "mul.h"
namespace numerics{
// ---------------- inplace_random ----------------
template <typename T>
inline void inplace_random(utils::Vector<T>& v, const T low, const T high) {
detail::inplace_random_serial(v, low, high);
}
template <typename T>
inline void inplace_random(utils::Matrix<T>& A, const T low, const T high) {
detail::inplace_random_serial(A, low, high);
}
// ---------------- random ----------------
template <typename T>
inline utils::Vector<T> random_vector(const uint64_t n, const T low, const T high) {
utils::Vector<T> v(n);
inplace_random(v, low, high);
return v;
}
template <typename T>
inline utils::Matrix<T> random_matrix(const uint64_t n, const uint64_t m, const T low, const T high) {
utils::Matrix<T> A(n, m);
inplace_random(A, low, high);
return A;
}
// ---------------- inplace random add----------------
template <typename T>
inline void inplace_random_add(utils::Vector<T>& v, const T low, const T high){
utils::Vector<T> noise = random_vector(v.size(), low, high);
inplace_add(v, noise);
}
template <typename T>
inline void inplace_random_add(utils::Matrix<T>& A, const T low, const T high){
utils::Matrix<T> noise = random_matrix(A.rows(), A.cols(), low, high);
inplace_add(A, noise);
}
// ---------------- random add----------------
template <typename T>
inline utils::Vector<T> random_add(const utils::Vector<T>& v, const T low, const T high){
utils::Vector<T> out = v;
inplace_random_add(out, low, high);
return out;
}
template <typename T>
inline utils::Matrix<T> random_add(const utils::Matrix<T>& A, const T low, const T high){
utils::Matrix<T> out = A;
inplace_random_add(out, low, high);
return out;
}
// ---------------- inplace random mul----------------
template <typename T>
inline void inplace_random_mul(utils::Vector<T>& v, const T low, const T high){
utils::Vector<T> noise = random_vector(v.size(), low, high);
inplace_mul(v, noise);
}
template <typename T>
inline void inplace_random_mul(utils::Matrix<T>& A, const T low, const T high){
utils::Matrix<T> noise = random_matrix(A.rows(), A.cols(), low, high);
inplace_mul(A, noise);
}
// ---------------- random mul----------------
template <typename T>
inline utils::Vector<T> random_mul(const utils::Vector<T>& v, const T low, const T high){
utils::Vector<T> out = v;
inplace_random_mul(out, low, high);
return out;
}
template <typename T>
inline utils::Matrix<T> random_mul(const utils::Matrix<T>& A, const T low, const T high){
utils::Matrix<T> out = A;
inplace_random_mul(out, low, high);
return out;
}
}
include/numerics/vecrandom.h
-48
View File
@@ -1,48 +0,0 @@
#ifndef _vecrandom_n_
#define _vecrandom_n_
#include "./utils/vector.h"
#include "./utils/matrix.h"
#include "./core/omp_config.h"
namespace numerics{
template <typename T>
void inplace_vecrandom_add(utils::Vector<T>& a, const T lower, const T higher) {
const uint64_t N = a.size();
utils::Vector<T> b;
b.random(N, lower, higher);
a.inplace_add(b);
}
template <typename T>
void inplace_vecrandom_mul(utils::Vector<T>& a, const T lower, const T higher) {
const uint64_t N = a.size();
utils::Vector<T> b;
b.random(N, lower, higher);
a.inplace_multiply(b);
}
template <typename T>
utils::Vector<T> vecrandom_add(const utils::Vector<T>& a, const T lower, const T higher) {
const uint64_t N = a.size();
utils::Vector<T> b;
inplace_vecrandom_add(b, lower, higher);
return b;
}
} // namespace numerics
#endif // _vecrandom_n_