matrix.hpp

Go to the documentation of this file.

// -------------------------------------------------------//
//
// SHAMROCK code for hydrodynamics
// Copyright (c) 2021-2026 Timothée David--Cléris <tim.shamrock@proton.me>
// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1
// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information
//
// -------------------------------------------------------//
 
#pragma once

 
#include "shambase/assert.hpp"
#include "shambackends/sycl.hpp"
#include "shammath/matrix_op.hpp"
#include <experimental/mdspan>
#include <array>
 
namespace shammath {

    template<class T, int m, int n>
    class mat {
        public:
        std::array<T, m * n> data;

        inline constexpr auto get_mdspan() {
            return std::mdspan<T, std::extents<size_t, m, n>>(data.data());
        }

        inline constexpr auto get_mdspan() const {
            return std::mdspan<const T, std::extents<size_t, m, n>>(data.data());
        }

        inline constexpr T &operator()(int i, int j) { return get_mdspan()(i, j); }

        inline constexpr const T &operator()(int i, int j) const { return get_mdspan()(i, j); }

        bool operator==(const mat<T, m, n> &other) const { return data == other.data; }
 
        // Addition operator for matrices
        inline mat &operator+=(const mat &other) {
#pragma unroll
            for (size_t i = 0; i < m * n; i++) {
                data[i] += other.data[i];
            }
            return *this;
        }

        bool equal_at_precision(const mat<T, m, n> &other, const T precision) const {
            bool res = true;
            for (auto i = 0; i < m; i++) {
                for (auto j = 0; j < n; j++) {
                    if (sham::abs(data[i * n + j] - other.data[i * n + j]) >= precision) {
                        res = false;
                    }
                }
            }
            return res;
        }
    };

    template<class T, int n>
    inline constexpr mat<T, n, n> mat_identity() {
        mat<T, n, n> res{};
        mat_set_identity(res.get_mdspan());
        return res;
    }

    template<class T, int n>
    class vec {
        public:
        std::array<T, n> data;

        inline constexpr auto get_mdspan() {
            return std::mdspan<T, std::extents<size_t, n>>(data.data());
        }

        inline constexpr auto get_mdspan_mat_col() {
            return std::mdspan<T, std::extents<size_t, n, 1>>(data.data());
        }

        inline constexpr auto get_mdspan_mat_row() {
            return std::mdspan<T, std::extents<size_t, 1, n>>(data.data());
        }

        inline constexpr T &operator[](int i) { return get_mdspan()(i); }

        bool operator==(const vec<T, n> &other) const { return data == other.data; }
    };

    template<class T>
    class mat_d {
        public:
        std::vector<T> data;
        int rows;
        int columns;

        mat_d(int rows, int columns) : rows(rows), columns(columns), data(rows * columns) {}

        inline constexpr auto get_mdspan() {
            return std::mdspan<T, std::dextents<size_t, 2>>(data.data(), rows, columns);
        }

        inline constexpr auto get_mdspan() const {
            return std::mdspan<const T, std::dextents<size_t, 2>>(data.data(), rows, columns);
        }

        inline constexpr T &operator()(int i, int j) { return get_mdspan()(i, j); }

        inline constexpr const T &operator()(int i, int j) const { return get_mdspan()(i, j); }

        bool operator==(const mat_d<T> &other) const { return data == other.data; }

        inline mat_d &operator+=(const mat_d &other) {
            for (size_t i = 0; i < get_mdspan().extent(0) * get_mdspan().extent(1); i++) {
                data[i] += other.data[i];
            }
            return *this;
        }

        bool equal_at_precision(const mat_d<T> &other, const T precision) const {
            bool res = true;
            for (auto i = 0; i < rows; i++) {
                for (auto j = 0; j < columns; j++) {
                    if (sham::abs(data[i * columns + j] - other.data[i * columns + j])
                        >= precision) {
                        res = false;
                    }
                }
            }
            return res;
        }
    };

    template<class T>
    class vec_d {
        public:
        std::vector<T> data;
        int size;

        vec_d(int size) : size(size), data(size) {}

        inline constexpr auto get_mdspan() {
            return std::mdspan<T, std::dextents<size_t, 1>>(data.data(), size);
        }

        inline constexpr auto get_mdspan_mat_col() {
            return std::mdspan<T, std::dextents<size_t, 2>>(data.data(), size, 1);
        }

        inline constexpr auto get_mdspan_mat_row() {
            return std::mdspan<T, std::dextents<size_t, 2>>(data.data(), 1, size);
        }

        inline constexpr T &operator[](int i) { return get_mdspan()(i); }

        bool operator==(const vec_d<T> &other) const { return data == other.data; }
    };
} // namespace shammath
 
template<class T, int m, int n>
struct sham::VectorProperties<shammath::mat<T, m, n>> {
    using component_type           = T;
    static constexpr u32 dimension = m * n;
 
    static constexpr bool is_float_based
        = std::is_same<T, f16>::value || std::is_same<T, f32>::value || std::is_same<T, f64>::value;
    static constexpr bool is_uint_based = std::is_same<T, u8>::value || std::is_same<T, u16>::value
                                          || std::is_same<T, u32>::value
                                          || std::is_same<T, u64>::value;
    static constexpr bool is_int_based  = std::is_same<T, i8>::value || std::is_same<T, i16>::value
                                          || std::is_same<T, i32>::value
                                          || std::is_same<T, i64>::value;
    static constexpr bool has_info      = is_float_based || is_int_based || is_uint_based;
 
    static constexpr shammath::mat<T, m, n> get_min() {
        constexpr T min = shambase::get_min<T>();
        return {min};
    }
    static constexpr shammath::mat<T, m, n> get_max() {
        constexpr T max = shambase::get_max<T>();
        return {max};
    }
    static constexpr shammath::mat<T, m, n> get_zero() {
        constexpr T zero = 0;
        return {zero};
    }
};
 
using f32_3x3 = shammath::mat<f32, 3, 3>; 
using f32_4x4 = shammath::mat<f32, 4, 4>; 
using f64_3x3 = shammath::mat<f64, 3, 3>; 
using f64_4x4 = shammath::mat<f64, 4, 4>;