graduate_code/cpp/matrix.hpp

#pragma once
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <stdexcept>
#include <vector>
namespace graduate {
inline void validate(const std::vector<double>& a,
                     const std::vector<double>& b,
                     std::size_t m, std::size_t k, std::size_t n) {
    if (m==0 || k==0 || n==0 || m>512 || k>512 || n>512)
        throw std::invalid_argument("dimensions must be in 1..512");
    if (a.size()!=m*k || b.size()!=k*n)
        throw std::invalid_argument("shape mismatch");
    for (double x:a) if (!std::isfinite(x))
        throw std::invalid_argument("nonfinite A");
    for (double x:b) if (!std::isfinite(x))
        throw std::invalid_argument("nonfinite B");
}
inline void check_result(const std::vector<double>& c) {
    for (double x:c) if (!std::isfinite(x))
        throw std::overflow_error("nonfinite result");
}
inline std::vector<double> naive(const std::vector<double>& a,
                                 const std::vector<double>& b,
                                 std::size_t m,std::size_t k,std::size_t n) {
    validate(a,b,m,k,n);
    std::vector<double> c(m*n,0.);
    for (std::size_t i=0;i<m;++i)
        for (std::size_t j=0;j<n;++j)
            for (std::size_t p=0;p<k;++p)
                c[i*n+j]+=a[i*k+p]*b[p*n+j];
    check_result(c);
    return c;
}
inline std::vector<double> blocked(const std::vector<double>& a,
                                   const std::vector<double>& b,
                                   std::size_t m,std::size_t k,std::size_t n,
                                   std::size_t tile=24) {
    validate(a,b,m,k,n);
    if (tile==0 || tile>1024)
        throw std::invalid_argument("tile must be in 1..1024");
    std::vector<double> c(m*n,0.);
    for (std::size_t ii=0;ii<m;ii+=tile)
        for (std::size_t pp=0;pp<k;pp+=tile)
            for (std::size_t jj=0;jj<n;jj+=tile)
                for (std::size_t i=ii;i<std::min(ii+tile,m);++i)
                    for (std::size_t p=pp;p<std::min(pp+tile,k);++p) {
                        const double value=a[i*k+p];
                        for (std::size_t j=jj;j<std::min(jj+tile,n);++j)
                            c[i*n+j]+=value*b[p*n+j];
                    }
    check_result(c);
    return c;
}
} // namespace graduate