matrix.cpp

Go to the documentation of this file.

/// @file core/matrix.cpp
/// @brief Matrix constructors, GPU lifecycle, and backend dispatch for matrix ops.
///
/// Adding a new backend:
///   1. Add the enumerator to enum class Backend in include/core/policy.hpp
///   2. Create src/core/backends/<name>/ with impl.hpp and matrix.cpp
///   3. Add `case Backend::<name>:` to each switch below
///   4. Register the .cpp in cmake/sources.cmake
 
#include "core/matrix.hpp"
#include "core/parallel/cuda_ops.hpp"
#include <algorithm>
 
#include "backends/seq/impl.hpp"
#include "backends/blas/impl.hpp"
#include "backends/omp/impl.hpp"
#include "backends/gpu/impl.hpp"
#include "backends/simd/impl.hpp"
 
namespace num {
 
Matrix::Matrix(idx rows, idx cols)
    : rows_(rows), cols_(cols), data_(new real[rows * cols]()) {}
 
Matrix::Matrix(idx rows, idx cols, real val)
    : rows_(rows), cols_(cols), data_(new real[rows * cols]) {
    std::fill_n(data_.get(), size(), val);
}
 
Matrix::~Matrix() {
    if (d_data_) cuda::free(d_data_);
}
 
Matrix::Matrix(const Matrix& o)
    : rows_(o.rows_), cols_(o.cols_), data_(new real[o.size()]) {
    std::copy_n(o.data_.get(), size(), data_.get());
}
 
Matrix::Matrix(Matrix&& o) noexcept
    : rows_(o.rows_), cols_(o.cols_), data_(std::move(o.data_)), d_data_(o.d_data_) {
    o.rows_ = o.cols_ = 0;
    o.d_data_ = nullptr;
}
 
Matrix& Matrix::operator=(const Matrix& o) {
    if (this != &o) {
        rows_ = o.rows_;
        cols_ = o.cols_;
        data_.reset(new real[size()]);
        std::copy_n(o.data_.get(), size(), data_.get());
    }
    return *this;
}
 
Matrix& Matrix::operator=(Matrix&& o) noexcept {
    if (this != &o) {
        if (d_data_) cuda::free(d_data_);
        rows_ = o.rows_;
        cols_ = o.cols_;
        data_ = std::move(o.data_);
        d_data_ = o.d_data_;
        o.rows_ = o.cols_ = 0;
        o.d_data_ = nullptr;
    }
    return *this;
}
 
void Matrix::to_gpu() {
    if (!d_data_) {
        d_data_ = cuda::alloc(size());
        cuda::to_device(d_data_, data_.get(), size());
    }
}
 
void Matrix::to_cpu() {
    if (d_data_) {
        cuda::to_host(data_.get(), d_data_, size());
        cuda::free(d_data_);
        d_data_ = nullptr;
    }
}
 
void matmul(const Matrix& A, const Matrix& B, Matrix& C, Backend b) {
    switch (b) {
    case Backend::seq:      backends::seq::matmul(A, B, C);              break;
    case Backend::blocked:  backends::seq::matmul_blocked(A, B, C, 64);  break;
    case Backend::simd:     backends::simd::matmul(A, B, C, 64);         break;
    case Backend::blas:     backends::blas::matmul(A, B, C);             break;
    case Backend::omp:      backends::omp::matmul(A, B, C);              break;
    case Backend::gpu:      backends::gpu::matmul(A, B, C);              break;
    }
}
 
void matvec(const Matrix& A, const Vector& x, Vector& y, Backend b) {
    switch (b) {
    case Backend::seq:      backends::seq::matvec(A, x, y);    break;
    case Backend::blocked:  backends::seq::matvec(A, x, y);    break;
    case Backend::simd:     backends::simd::matvec(A, x, y);   break;
    case Backend::blas:     backends::blas::matvec(A, x, y);   break;
    case Backend::omp:      backends::omp::matvec(A, x, y);    break;
    case Backend::gpu:      backends::gpu::matvec(A, x, y);    break;
    }
}
 
void matadd(real alpha, const Matrix& A, real beta, const Matrix& B, Matrix& C,
            Backend b) {
    switch (b) {
    case Backend::seq:
    case Backend::blocked:
    case Backend::simd:     backends::seq::matadd(alpha, A, beta, B, C);  break;
    case Backend::blas:     backends::blas::matadd(alpha, A, beta, B, C); break;
    case Backend::omp:      backends::omp::matadd(alpha, A, beta, B, C);  break;
    case Backend::gpu:      backends::seq::matadd(alpha, A, beta, B, C);  break;
    }
}
 
void matmul_blocked(const Matrix& A, const Matrix& B, Matrix& C, idx block_size) {
    backends::seq::matmul_blocked(A, B, C, block_size);
}
 
void matmul_register_blocked(const Matrix& A, const Matrix& B, Matrix& C,
                              idx block_size, idx reg_size) {
    backends::seq::matmul_register_blocked(A, B, C, block_size, reg_size);
}
 
void matmul_simd(const Matrix& A, const Matrix& B, Matrix& C, idx block_size) {
    backends::simd::matmul(A, B, C, block_size);
}
 
void matvec_simd(const Matrix& A, const Vector& x, Vector& y) {
    backends::simd::matvec(A, x, y);
}
 
} // namespace num