matrix.cpp

Go to the documentation of this file.

/// @file core/backends/seq/matrix.cpp
/// @brief Sequential and blocked C++ matrix operations
 
#include "core/matrix.hpp"
#include "kernel/raw.hpp"
#include <algorithm>
 
namespace num::backends::seq {
 
void matmul(const Matrix& A, const Matrix& B, Matrix& C) {
    const idx M = A.rows(), K = A.cols(), N = B.cols();
    for (idx i = 0; i < M; ++i)
        for (idx j = 0; j < N; ++j) {
            C(i, j) = 0;
            for (idx k = 0; k < K; ++k)
                C(i, j) += A(i, k) * B(k, j);
        }
}
 
void matvec(const Matrix& A, const Vector& x, Vector& y) {
    kernel::raw::matvec(y.data(), A.data(), x.data(), A.rows(), A.cols());
}
 
void matadd(real alpha, const Matrix& A, real beta, const Matrix& B, Matrix& C) {
    kernel::raw::axpbyz(C.data(), A.data(), B.data(), alpha, beta, A.size());
}
 
void matmul_blocked(const Matrix& A, const Matrix& B, Matrix& C, idx block_size) {
    const idx M = A.rows(), K = A.cols(), N = B.cols();
    std::fill_n(C.data(), M * N, real(0));
 
    for (idx ii = 0; ii < M; ii += block_size) {
        const idx i_end = std::min(ii + block_size, M);
        for (idx jj = 0; jj < N; jj += block_size) {
            const idx j_end = std::min(jj + block_size, N);
            for (idx kk = 0; kk < K; kk += block_size) {
                const idx k_end = std::min(kk + block_size, K);
                for (idx i = ii; i < i_end; ++i) {
                    for (idx k = kk; k < k_end; ++k) {
                        const real a_ik = A(i, k);
                        for (idx j = jj; j < j_end; ++j)
                            C(i, j) += a_ik * B(k, j);
                    }
                }
            }
        }
    }
}
 
void matmul_register_blocked(const Matrix& A,
                             const Matrix& B,
                             Matrix& C,
                             idx block_size,
                             idx reg_size) {
    const idx M = A.rows(), K = A.cols(), N = B.cols();
    std::fill_n(C.data(), M * N, real(0));
 
    for (idx ii = 0; ii < M; ii += block_size) {
        const idx i_lim = std::min(ii + block_size, M);
        for (idx jj = 0; jj < N; jj += block_size) {
            const idx j_lim = std::min(jj + block_size, N);
            for (idx kk = 0; kk < K; kk += block_size) {
                const idx k_lim = std::min(kk + block_size, K);
                for (idx ir = ii; ir < i_lim; ir += reg_size) {
                    const idx ri = std::min(ir + reg_size, i_lim);
                    for (idx jr = jj; jr < j_lim; jr += reg_size) {
                        const idx rj = std::min(jr + reg_size, j_lim);
                        real c[4][4] = {};
                        for (idx i = ir; i < ri; ++i)
                            for (idx j = jr; j < rj; ++j)
                                c[i - ir][j - jr] = C(i, j);
                        for (idx k = kk; k < k_lim; ++k) {
                            for (idx i = ir; i < ri; ++i) {
                                const real a_ik = A(i, k);
                                for (idx j = jr; j < rj; ++j)
                                    c[i - ir][j - jr] += a_ik * B(k, j);
                            }
                        }
                        for (idx i = ir; i < ri; ++i)
                            for (idx j = jr; j < rj; ++j)
                                C(i, j) = c[i - ir][j - jr];
                    }
                }
            }
        }
    }
}
 
} // namespace num::backends::seq