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

ready for parralization

Browse Source
This commit is contained in:
Bausager
2025-09-12 22:58:52 +02:00
parent cb825aec40
commit 320436ce98
14 changed files with 920 additions and 294 deletions
Download Patch File Download Diff File Expand all files Collapse all files
include/core/omp_config.h
+34
View File
@@ -0,0 +1,34 @@
#pragma once
#include <omp.h>
// Configure OpenMP behavior at runtime.
inline void omp_configure(int max_active_levels,
bool dynamic_threads,
const std::vector<int>& threads_per_level = {},
bool bind_close = true)
{
// 1) Allow nested parallel regions (levels of teams)
// Example: outer #pragma omp parallel ... { inner #pragma omp parallel ... }
omp_set_max_active_levels(max_active_levels); // 1 = only top-level; 2+ enables nesting
// 2) Let the runtime shrink/grow thread counts if it thinks it should
// (helps avoid oversubscription when you accidentally ask for too many threads)
omp_set_dynamic(dynamic_threads ? 1 : 0);
// 3) Thread binding (keep threads near their cores) is controlled via env vars,
// so here we just *recommend* a good default (see below). You *can* setenv()
// in code, but it’s cleaner to do it outside the program.
(void)bind_close; // documented below in env var section
// 4) Top-level default thread count (inner levels are usually set per region)
if (!threads_per_level.empty()) {
omp_set_num_threads(threads_per_level[0]); // e.g. 16 for the outermost team
// Inner levels:
// - Use num_threads(threads_per_level[L]) on the inner #pragma omp parallel
// - or set OMP_NUM_THREADS="outer,inner,inner2" as an environment variable
}
}
include/numerics/inverse.h
+113
View File
@@ -0,0 +1,113 @@
#ifndef _inverse_n_
#define _inverse_n_
#include "./utils/vector.h"
#include "./utils/matrix.h"
namespace numerics{
template <typename T>
void inplace_inverse(utils::Matrix<T>& A, std::string method = "Gauss-Jordan"){
if (method == "Gauss-Jordan"){
utils::Matrix<T> B(A.rows(),A.cols(), T{0});
uint64_t icol{0}, irow{0}, rows{A.rows()}, cols{A.cols()};
double big, dum, pivinv, temp;
utils::Vi indxc(rows,0), indxr(rows,0), ipiv(rows,0);
//for (uint64_t j = 0; j < N; ++j){ ipiv[j] = 0;}
for (uint64_t i = 0; i < rows; i++){
big = 0.0;
for (uint64_t j = 0; j < rows; j++){
if (ipiv[j] != 1){
for (uint64_t k = 0; k < rows; k++){
if (ipiv[k] == 0){
if (abs(A(j,k)) >= big){
big = abs(A(j,k));
irow = j;
icol = k;
}
}
}
}
}
ipiv[icol]++;
if (irow != icol){
for (uint64_t l = 0; l < rows; l++){ // SWAP
temp = A(irow,l);
A(irow,l) = A(icol,l);
A(icol,l) = temp;
}
for (uint64_t l = 0; l < cols; l++){ // SWAP temp matrix
temp = B(irow,l);
B(irow,l) = B(icol,l);
B(icol,l) = temp;
}
}
indxr[i] = irow;
indxc[i] = icol;
if (A(icol,icol) == 0.0){
throw std::runtime_error("utill:inplace_inverse('Gauss-Jordan' - Singular Matrix");
}
pivinv= 1.0/A(icol,icol);
A(icol,icol)=1.0;
for (uint64_t l = 0; l < rows; l++){
A(icol,l) *= pivinv;
}
for (uint64_t l = 0; l < cols; l++){
B(icol,l) *= pivinv;
}
for (uint64_t ll = 0; ll < rows; ll++){
if (ll != icol){
dum = A(ll,icol);
A(ll,icol) = 0;
for (uint64_t l = 0; l < rows; l++){
A(ll,l) -= A(icol,l)*dum;
}
for (uint64_t l = 0; l < rows; l++){
B(ll,l) -= B(icol,l)*dum;
}
}
}
}
//m = temp_m;
for (int64_t l = rows-1; l >= 0; l--){
if (indxr[l] != indxc[l]){
for (uint64_t k = 0; k < rows; k++){
temp = A(k,indxr[l]);
A(k,indxr[l]) = A(k,indxc[l]);
A(k,indxc[l]) = temp;
}
}
}
}
else{
throw std::runtime_error("numerics::inplace_inverse(" + method + ") - Not implemented yet \r \nImplemented: 'Gauss-Jordan',");
}
}
template <typename T>
utils::Matrix<T> inverse(utils::Matrix<T>& A, std::string method = "Gauss-Jordan"){
utils::Matrix<T> B = A;
inplace_inverse(B, method);
return B;
}
} // namespace numerics
#endif // _inverse_n_
include/numerics/matmul.h
+42
View File
@@ -0,0 +1,42 @@
#ifndef _matmul_n_
#define _matmul_n_
#include "./utils/matrix.h"
namespace numerics{
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, n, T{0});
//#pragma omp parallel for collapse(2) schedule(static)
#pragma omp parallel for
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
#endif // _matmul_n_
include/numerics/matvec.h
+54
View File
@@ -0,0 +1,54 @@
#ifndef _matvec_n_
#define _matvec_n_
#include "./utils/matrix.h"
namespace numerics{
// y = A * x, where A is (m×n) and x is length n and y is length m
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) {
T acc = T{0};
for (uint64_t j = 0; j < n; ++j) {
acc += A(i, j) * x[j];
}
y[i] = acc;
}
return y;
}
// y = x * A, where x is length m and A is (m×n) -> y is length n
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) {
T acc = T{0};
for (uint64_t i = 0; i < m; ++i) {
acc += x[i] * A(i, j);
}
y[j] = acc;
}
return y;
}
} // namespace numerics
#endif // _matvec_n_
include/numerics/numerics.h
+7
View File
@@ -0,0 +1,7 @@
// "./numerics/numerics.h"
#pragma once
#include "./numerics/transpose.h"
#include "./numerics/inverse.h"
#include "./numerics/matmul.h"
#include "./numerics/matvec.h"
include/numerics/transpose.h
+70
View File
@@ -0,0 +1,70 @@
#ifndef _transpose_n_
#define _transpose_n_
#include "./utils/matrix.h"
namespace numerics{
template <typename T>
void inplace_transpose(utils::Matrix<T>& A){
const uint64_t rows = A.rows();
const uint64_t cols = A.cols();
if (rows != cols){
throw std::runtime_error("inplace_transpose only valid for square matrices");
}
for (uint64_t i = 0; i < rows; ++i){
for (uint64_t j = i + 1; j < cols; ++j){
T tmp = A(j,i);
A(j,i) = A(i,j);
A(i,j) = tmp;
//std::swap(A(j,i), A(i,j));
}
}
}
template <typename T>
utils::Matrix<T> transpose(const utils::Matrix<T>& A){
const uint64_t rows = A.rows();
const uint64_t cols = A.cols();
utils::Matrix<T> B(cols, rows, T{});
for (uint64_t i = 0; i < rows; ++i){
for (uint64_t j = 0; j < cols; ++j){
B(j,i) = A(i,j);
}
}
return B;
}
} // namespace numerics
#endif // _transpose_n_
include/utils/Grid1D.h
-39
View File
@@ -1,39 +0,0 @@
#ifndef _grid1d_n_
#define _grid1d_n_
#include "./utils/matrix.h"
namespace utils{
//#######################################
//# Grid1D TYPE #
//#######################################
template <typename T>
struct Grid1D{
utils::Vector<T> grid;
utils::Vector<T> vertices;
utils::Vector<T> vertices_norm;
void create_vertices_norm(){
vertices_norm.fill(vertices.size()*2, 0);
uint64_t k = 0;
for (uint64_t i = 0; i < grid.size(); i++){
for (uint64_t j = 1; j <= 2; j++){
vertices_norm[k] = grid[i] - vertices[i+j];
k++;
}
//vertices_norm[(i*2)+1] = grid[i] - vertices[(i*2)+1];
}
vertices_norm.print();
}
};
typedef Grid1D<int> Gridi;
typedef Grid1D<float> Gridf;
typedef Grid1D<double> Gridd;
}
#endif // _grid1d_n_
include/utils/matrix.h
+146 -217
View File
@@ -3,249 +3,178 @@
#include "./utils/vector.h"
#include <iomanip>
namespace utils{
//#######################################
//# MATRIX TYPE #
//# Backed by utils::Vector<T> #
//#######################################
template <typename T>
struct Matrix{
utils::Vector<T> m;
T& operator[](uint64_t idx) { return m[idx]; } // Makes it able to do matr[1][1]
const T& operator[](uint64_t idx) const { return m[idx]; } // Makes it able to do matr[1][1]
using vector_type = typename decltype(std::declval<T>().v)::value_type;
class Matrix{
public:
Matrix() : rows_(0), cols_(0), data_() {} // Default constructor
// Constructor to initialize matrix with rows × cols and a fill value
Matrix(uint64_t rows, uint64_t cols, typename T::value_type value = {}) {
fill(rows, cols, value);
}
Matrix() = default; // Default constructor
Matrix(uint64_t rows, uint64_t cols, const T& value = T())
: rows_(rows), cols_(cols), data_(rows * cols, value) {}
//# MATRIX: basic properties #
uint64_t rows() const noexcept {return rows_;}
uint64_t cols() const noexcept {return cols_;}
//# MATRIX: element access (fast; unchecked) #
T& operator()(uint64_t i, uint64_t j) { return data_[i * cols_ + j]; }
const T& operator()(uint64_t i, uint64_t j) const { return data_[i * cols_ + j]; }
void fill(uint64_t rows, uint64_t cols, const vector_type num=0){
m.clear();
for (uint64_t i = 0; i < rows; i++){
T temp_vec;
//# MATRIX: data access #
T* data() noexcept { return data_.data(); }
const T* data() const noexcept { return data_.data(); }
for (uint64_t j = 0; j < cols; j++){
temp_vec.v.push_back(num);
}
m.push_back(temp_vec);
}
//# MATRIX: equal operator #
bool operator==(const Matrix<T>& A) const {
if (rows_ != A.rows_ || cols_ != A.cols_) return false;
for (uint64_t i = 0; i < rows_; ++i)
for (uint64_t j = 0; j < cols_; ++j)
if (data_[i*cols_ + j] != A(i,j))
return false;
return true;
}
void fill_RNG(const uint64_t rows, const uint64_t cols, const vector_type min = 0, const vector_type max = 1){
m.clear();
std::mt19937_64 rng{};
rng.seed( std::random_device{}());
for (uint64_t i = 0; i < rows; i++){
T temp_vec;
for (uint64_t j = 0; j < cols; j++){
temp_vec.v.push_back(std::uniform_real_distribution<>{min, max}(rng));
}
m.push_back(temp_vec);
}
bool operator!=(const Matrix<T>& A) const {
return !(*this == A);
}
inline friend std::ostream& operator << (std::ostream& out, const Matrix& mat){
out << "[";
for (uint64_t i = 0; i < mat.m.size(); i++){
out << "[";
for (uint64_t j = 0; j < mat.m[i].v.size(); j++){
if (j % mat.m[i].v.size() == mat.m[i].v.size() -1 && i == mat.m.size()-1){
out << mat.m[i].v[j] << "]";
}
else if ((j % mat.m[i].v.size() == mat.m[i].v.size() -1)){
out << mat.m[i].v[j] << "]," << std::endl;
}
else{
out << mat.m[i].v[j] << ", ";
}
}
}
out << "]";
return out;
}
void print() const{
std::cout << *this << std::endl;
}
void inplace_transpose(){
utils::Vector<T> temp_m = m;
m.clear();
uint64_t rows = temp_m.size();
uint64_t cols = temp_m[0].v.size();
for (uint64_t i = 0; i < cols; i++){
T temp_vec;
for (uint64_t j = 0; j < rows; j++){
temp_vec.v.push_back(temp_m[j].v[i]);
}
m.push_back(temp_vec);
}
}
Matrix<T> transpose()const{
Matrix<T> copy = *this;
copy.inplace_transpose();
return copy;
}
void inplace_inverse(std::string method = "Gauss-Jordan"){
//Matrix<T> temp_m = *this; // Copies the m into temp_m correctly (Before: utils::Vector<T> temp_m = m;)
if (method == "Gauss-Jordan"){
Matrix<T> temp_m(m.v.size(),m[0].v.size(),0);
//std::cout << temp_m.m.v[0].size() << std::endl;
//std::cout << m.v.size() << std::endl;
uint64_t icol,irow,N=m.v.size(),M=temp_m.m.v[0].size();
double big,dum,pivinv;
Vi indxc(N,0),indxr(N,0),ipiv(N,0);
//for (uint64_t j = 0; j < N; ++j){ ipiv[j] = 0;}
for (uint64_t i = 0; i < N; i++){
big=0.0;
for (uint64_t j = 0; j < N; j++){
if (ipiv[j] != 1){
for (uint64_t k = 0; k < N; k++){
if (ipiv[k] == 0){
if (abs(m[j].v[k]) >= big){
big = abs(m[j].v[k]);
irow = j;
icol = k;
}
}
}
}
}
ipiv[icol]++;
if (irow != icol){
for (uint64_t l = 0; l < N; l++){ // SWAP
double temp = m[irow].v[l];
m[irow].v[l] = m[icol].v[l];
m[icol].v[l] = temp;
}
for (uint64_t l = 0; l < M; l++){ // SWAP temp matrix
double temp = temp_m.m[irow].v[l];
temp_m.m[irow].v[l] = temp_m.m[icol].v[l];
temp_m.m[icol].v[l] = temp;
}
}
indxr[i] = irow;
indxc[i] = icol;
if (m[icol].v[icol] == 0.0){
throw std::runtime_error("utill:Matrix.Gauss-Jordan - Singular Matrix");
}
pivinv= 1.0/m[icol].v[icol];
m[icol].v[icol]=1.0;
for (uint64_t l = 0; l < N; l++){
m[icol].v[l] *= pivinv;
}
for (uint64_t l = 0; l < M; l++){
temp_m.m[icol].v[l] *= pivinv;
}
for (uint64_t ll = 0; ll < N; ll++){
if (ll != icol){
dum = m[ll].v[icol];
m[ll].v[icol] = 0;
for (uint64_t l = 0; l < N; l++){
m[ll].v[l] -= m[icol].v[l]*dum;
}
for (uint64_t l = 0; l < N; l++){
temp_m.m[ll].v[l] -= temp_m.m[icol].v[l]*dum;
}
}
}
}
//m = temp_m;
for (int64_t l = N-1; l >= 0; l--){
if (indxr[l] != indxc[l]){
for (uint64_t k = 0; k < N; k++){
double temp = m[k].v[indxr[l]];
m[k].v[indxr[l]] = m[k].v[indxc[l]];
m[k].v[indxc[l]] = temp;
}
}
}
}
else{
throw std::runtime_error("utill:Matrix." + method + " - Not implemented yet \r \nImplemented: 'Gauss-Jordan',");
}
}
Matrix<T> inverse(std::string method = "Gauss-Jordan")const{
Matrix<T> copy = *this;
copy.inplace_inverse(method);
return copy;
}
utils::Vector<vector_type> vecmult(const utils::Vector<vector_type>& Vec)const{
if (m[0].size() != Vec.size()){
throw std::runtime_error("utill:Matrix.vecmult - Dimentions does not fit");
}
// Create a temporary result vector
utils::Vector<vector_type> copy(Vec.size(), 0);
for (uint64_t i = 0; i < m.size(); ++i) {
for (uint64_t j = 0; j < m[0].size(); ++j) {
copy[i] += m[i][j] * Vec[j];
bool nearly_equal(const Matrix<T>& A, T tol = static_cast<T>(1e-9)) const {
if (rows_ != A.rows_ || cols_ != A.cols_) return false;
for (uint64_t i = 0; i < rows_; ++i)
for (uint64_t j = 0; j < cols_; ++j) {
T a = (*this)(i,j);
T b = A(i,j);
if (std::is_floating_point<T>::value) {
if (std::fabs(a - b) > tol) return false;
} else {
if (a != b) return false;
}
}
return true;
}
//# MATRIX: row helpers (copy out) #
// Read whole row as an owning Vector<T>
// utils::Vf v = M.get_row(2);
Vector<T> get_row(const uint64_t row) const {
if (row >= rows_) {
throw std::out_of_range("Matrix::get_row -> row index");
}
utils::Vector<T> result(cols_, T{});
for (uint64_t i = 0; i < cols_; ++i){
result[i] = data_[row * cols_ + i];
}
return result;
}
//# MATRIX: row helpers (copy in) #
// Assign a whole Vector<T> to a row
// M.set_row(2) = v;
void set_row(const uint64_t row, const Vector<T>& vector){
if (row >= rows_) {
throw std::out_of_range("Matrix::set_row -> row index");
}
if (vector.size() != cols_){
throw std::runtime_error("Matrix::set_row -> size mismatch");
}
return copy;
}
for (uint64_t i = 0; i < cols_; ++i){
data_[row * cols_ + i] = vector[i];
}
}
void inplace_matmult(const Matrix<T>& Mat){
//# MATRIX: col helpers (copy out) #
// Read whole col as an owning Vector<T>
// utils::Vf v = M.get_col(2);
Vector<T> get_col(const uint64_t col) const {
if (col >= cols_) {
throw std::out_of_range("Matrix::get_col -> col index");
}
utils::Vector<T> result(rows_, T{});
for (uint64_t i = 0; i < rows_; ++i){
result[i] = data_[i * cols_ + col];
}
return result;
}
if (m.v[0].size() != Mat.m.v.size()){
throw std::runtime_error("utill:Matrix.matmult - Dimentions does not fit");
}
//# MATRIX: col helpers (copy in) #
// Assign a whole Vector<T> to a col
// M.set_col(2) = v;
void set_col(const uint64_t col, const Vector<T>& vector){
if (col >= cols_) {
throw std::out_of_range("Matrix::set_col -> col index");
}
if (vector.size() != rows_){
throw std::runtime_error("Matrix::set_col -> size mismatch");
}
for (uint64_t i = 0; i < rows_; ++i){
data_[i * cols_ + col] = vector[i];
}
}
// Dimensions of the result
uint64_t rows = m.v.size(); // rows in *this
uint64_t cols = Mat.m[0].v.size(); // columns in Mat
uint64_t inner = m.v[0].size(); // shared dimension
void swap_rows(uint64_t a, uint64_t b){
if (a >= rows_ || b >= rows_) {
throw std::out_of_range("Matrix::swap_rows -> row index");
}
if (a == b){
return;
}
for (uint64_t i = 0; i < cols_; ++i){
T tmp = data_[a * cols_ + i];
data_[a * cols_ + i] = data_[b * cols_ + i];
data_[b * cols_ + i] = tmp;
}
}
void swap_cols(uint64_t a, uint64_t b){
if (a >= cols_ || b >= cols_) {
throw std::out_of_range("Matrix::swap_cols -> col index");
}
if (a == b){
return;
}
for (uint64_t i = 0; i < rows_; ++i){
T tmp = data_[i * cols_ + a];
data_[i * cols_ + a] = data_[i * cols_ + b];
data_[i * cols_ + b] = tmp;
}
}
// Create a temporary result matrix
Matrix<T> temp_m(rows, cols, 0);
inline friend std::ostream& operator<<(std::ostream& out, const Matrix& M) {
out << "[";
for (uint64_t i = 0; i < M.rows_; ++i) {
out << "[";
for (uint64_t j = 0; j < M.cols_; ++j) {
out << std::setw(4) << std::setprecision(3) << std::fixed << M(i, j);
if (j + 1 < M.cols_) out << ", ";
}
out << "]";
if (i + 1 < M.rows_) out << ",\n ";
}
out << "]";
return out;
}
// Perform matrix multiplication
for (uint64_t i = 0; i < rows; i++){
for (uint64_t j = 0; j < cols; j++){
for (uint64_t k = 0; k < inner; k++){
temp_m.m[i].v[j] += m[i].v[k] * Mat.m[k].v[j];
}
}
}
*this = temp_m;
}
Matrix<T> matmult(const Matrix<T>& Mat)const{
Matrix<T> copy = *this;
copy.inplace_matmult(Mat);
return copy;
}
void print() const {
std::cout << *this << std::endl;
}
private:
uint64_t rows_, cols_;
std::vector<T> data_;
};
typedef Matrix<Vi> Mi;
typedef Matrix<Vf> Mf;
typedef Matrix<Vd> Md;
};
typedef Matrix<int> Mi;
typedef Matrix<float> Mf;
typedef Matrix<double> Md;
}
#endif // _numerics_n_
#endif // _matrix_n_
include/utils/utils.h
-1
View File
@@ -3,4 +3,3 @@
#include "./utils/vector.h"
#include "./utils/matrix.h"
#include "./utils/Grid1D.h"
include/utils/vector.h
+24 -25
View File
@@ -44,6 +44,9 @@ public:
// vector.size();
uint64_t size() const noexcept { return v.size(); }
void resize(uint64_t new_size, const T& value = T()) {
v.resize(new_size, value);
}
//###########################################
//# VECTOR: == and != #
@@ -303,38 +306,19 @@ public:
Vector<T> result = *this;
result.inplace_power(a);
return result;
}
//################################################
//# VECTOR: Scalar Square #
//################################################
template <typename U, typename = typename std::enable_if<std::is_convertible<U, T>::value>::type>
void inplace_square(const U a){
const uint64_t n = v.size();
for (uint64_t i = 0; i < n; ++i){
v[i] = static_cast<T>(std::sqrt(v[i], a));
}
}
template <typename U, typename = typename std::enable_if<std::is_convertible<U, T>::value>::type>
Vector<T> square(const U a) const{
Vector<T> result = *this;
result.inplace_square(a);
return result;
}
//################################################
//# VECTOR: Vector square #
//################################################
void inplace_square(const Vector<T>& a){
if (a.size() != v.size()){
throw std::runtime_error("utill:Vector.inplace_square -> Dimensions does not fit");
}
uint64_t n = a.size();
void inplace_sqrt(){
uint64_t n = v.size();
for (uint64_t i = 0; i < n; ++i){
v[i] = static_cast<T>(std::sqrt(v[i], a[i]));
v[i] = static_cast<T>(std::sqrt(v[i]));
}
}
Vector<T> square(const Vector<T>& a) const{
Vector<T> sqrt() const{
Vector<T> result = *this;
result.inplace_square(a);
result.inplace_sqrt();
return result;
}
//###################################################
@@ -344,7 +328,7 @@ public:
if (a.size() != v.size()){
throw std::runtime_error("utill:Vector.dot -> Dimensions does not fit");
}
T result;
T result = T{0};
const uint64_t n = v.size();
for (uint64_t i = 0; i < n; ++i){
result += a[i]*v[i];
@@ -368,6 +352,21 @@ public:
T norm() const{
return static_cast<T>(std::sqrt(this->dot(*this)));
}
//############################################
//# VECTOR: Normalize #
//############################################
void inplace_normalize() {
T norm = this->norm();
if (norm == T{0}){
throw std::runtime_error("utils::Vector.normalize -> zero norm");
}
this->inplace_divide(norm);
}
Vector<T> normalize() const{
Vector<T> result = *this;
result.inplace_normalize();
return result;
}
//######################################################
//# VECTOR: Support Functions #
//######################################################