core/algebra/matrix.hpp

0001 /*
0002  * SPDX-PackageName: "covfie, a part of the ACTS project"
0003  * SPDX-FileCopyrightText: 2022 CERN
0004  * SPDX-License-Identifier: MPL-2.0
0005  */
0006
0007 #pragma once
0008
0009 #include <cstddef>
0010
0011 #include <covfie/core/array.hpp>
0012 #include <covfie/core/qualifiers.hpp>
0013
0014 namespace covfie::algebra {
0015 template <
0016     std::size_t N,
0017     std::size_t M,
0018     typename T = float,
0019     typename I = std::size_t>
0020 struct matrix {
0021     COVFIE_HOST_DEVICE matrix()
0022     {
0023     }
0024
0025     COVFIE_HOST_DEVICE matrix(array::array<array::array<T, M>, N> l)
0026     {
0027         for (I i = 0; i < l.size(); ++i) {
0028             for (I j = 0; j < l[i].size(); ++j) {
0029                 m_elems[i][j] = l[i][j];
0030             }
0031         }
0032     }
0033
0034     matrix(const matrix<N, M, T, I> &) = default;
0035
0036     matrix & operator=(const matrix<N, M, T, I> &) = default;
0037
0038     COVFIE_HOST_DEVICE T operator()(const I i, const I j) const
0039     {
0040         return m_elems[i][j];
0041     }
0042
0043     COVFIE_HOST_DEVICE T & operator()(const I i, const I j)
0044     {
0045         return m_elems[i][j];
0046     }
0047
0048     template <std::size_t P>
0049     COVFIE_HOST_DEVICE matrix<N, P, T, I> operator*(const matrix<M, P, T, I> & o
0050     ) const
0051     {
0052         matrix<N, P, T, I> r;
0053
0054         for (I i = 0; i < N; ++i) {
0055             for (I j = 0; j < P; ++j) {
0056                 T t = static_cast<T>(0.);
0057
0058                 for (I k = 0; k < M; ++k) {
0059                     t += this->operator()(i, k) * o(k, j);
0060                 }
0061
0062                 r(i, j) = t;
0063             }
0064         }
0065
0066         return r;
0067     }
0068
0069     COVFIE_HOST_DEVICE static matrix<N, M, T, I> identity()
0070     {
0071         matrix<N, M, T, I> result;
0072
0073         for (I i = 0; i < N; ++i) {
0074             for (I j = 0; j < M; ++j) {
0075                 result(i, j) = (i == j) ? static_cast<T>(1) : static_cast<T>(0);
0076             }
0077         }
0078
0079         return result;
0080     }
0081
0082 private:
0083     T m_elems[N][M];
0084 };
0085 }