dense.cpp

Go to the documentation of this file.

/// @file kernel/dense.cpp
/// @brief Implementations for num::kernel::dense (seq_t, par_t, and no-policy
/// ops).
 
#include "kernel/dense.hpp"
#include "kernel/raw.hpp"
#include <stdexcept>
 
namespace num::kernel::dense {
 
void ger(real alpha, const Vector& x, const Vector& y, Matrix& A, seq_t) noexcept {
    raw::ger(A.data(), x.data(), y.data(), alpha, x.size(), y.size());
}
 
void ger(real alpha, const Vector& x, const Vector& y, Matrix& A, par_t) {
#ifdef NUMERICS_HAS_OMP
    const idx m = x.size();
    const idx n = y.size();
    const real* xd = x.data();
    const real* yd = y.data();
    real* Ad = A.data();
    #pragma omp parallel for schedule(static)
    for (idx i = 0; i < m; ++i) {
        const real axi = alpha * xd[i];
        real* row = Ad + (i * n);
        for (idx j = 0; j < n; ++j) {
            row[j] += axi * yd[j];
        }
    }
#else
    ger(alpha, x, y, A, seq_t{});
#endif
}
 
void trsv_lower(const Matrix& L, const Vector& b, Vector& x) {
    const idx n = L.rows();
    if (L.cols() != n || b.size() != n) {
        throw std::invalid_argument("kernel::dense::trsv_lower: dimension mismatch");
    }
    if (x.size() != n) {
        x = Vector(n);
    }
    raw::trsv_lower(x.data(), L.data(), b.data(), n);
}
 
void trsv_upper(const Matrix& U, const Vector& b, Vector& x) {
    const idx n = U.rows();
    if (U.cols() != n || b.size() != n) {
        throw std::invalid_argument("kernel::dense::trsv_upper: dimension mismatch");
    }
    if (x.size() != n) {
        x = Vector(n);
    }
    raw::trsv_upper(x.data(), U.data(), b.data(), n);
}
 
} // namespace num::kernel::dense