row_to_column.hpp

Go to the documentation of this file.

#ifndef TATAMI_MULT_SPARSE_MATRIX_SPARSE_ROW_ROW_TO_COLUMN_HPP
#define TATAMI_MULT_SPARSE_MATRIX_SPARSE_ROW_ROW_TO_COLUMN_HPP
 
#include <cstddef>
#include <vector>
 
#include "tatami/tatami.hpp"
#include "sanisizer/sanisizer.hpp"
 
#include "../utils.hpp"
#include "../../utils.hpp"
 
namespace tatami_mult {
 
/* See https://github.com/tatami-inc/test-multiplication/tree/master/sparse_row/sparse_matrix
 * for an explanation of the choice of algorithm.
 */
 
struct MultiplySparseRowWithSparseRowMatrixToColumnOutputOptions {
    int num_threads = 1;
};
 
template<typename LeftValue_, typename LeftIndex_, typename RightValue_, typename RightIndex_, typename Output_>
void multiply_sparse_row_with_sparse_row_matrix_to_column_output(
    const tatami::Matrix<LeftValue_, LeftIndex_>& left,
    const tatami::Matrix<RightValue_, RightIndex_>& right,
    Output_* const output,
    const MultiplySparseRowWithSparseRowMatrixToColumnOutputOptions& options
) {
    const auto left_NR = left.nrow();
    const auto common_dim = left.ncol();
    const auto right_NC = right.ncol();
 
    auto right_vbuffers = tatami::create_container_of_Index_size<std::vector<std::vector<RightValue_> > >(common_dim);
    auto right_ibuffers = tatami::create_container_of_Index_size<std::vector<std::vector<RightIndex_> > >(common_dim);
    auto right_ranges = tatami::create_container_of_Index_size<std::vector<tatami::SparseRange<RightValue_, RightIndex_> > >(common_dim);
    populate_sparse_buffers(true, common_dim, right_NC, right, right_vbuffers, right_ibuffers, right_ranges, options.num_threads);
 
    tatami::parallelize([&](int, LeftIndex_ start, LeftIndex_ length) -> void {
        auto ext = tatami::consecutive_extractor<true>(left, true, start, length);
        auto vbuffer = tatami::create_container_of_Index_size<std::vector<LeftValue_> >(common_dim);
        auto ibuffer = tatami::create_container_of_Index_size<std::vector<LeftIndex_> >(common_dim);
 
        // Use a temporary buffer to (i) improve data locality and (ii) avoid false sharing.
        auto tmp_row = tatami::create_container_of_Index_size<std::vector<Output_> >(right_NC);
 
        for (LeftIndex_ lr = 0; lr < length; ++lr) {
            const auto lrange = ext->fetch(vbuffer.data(), ibuffer.data());
 
            for (LeftIndex_ x = 0; x < lrange.number; ++x) {
                const Output_ mult = lrange.value[x];
                const auto rrange = right_ranges[lrange.index[x]];
                for (RightIndex_ y = 0; y < rrange.number; ++y) {
                    tmp_row[rrange.index[y]] += mult * static_cast<Output_>(rrange.value[y]);
                }
            }
 
            // Technically, we could limit the transposition to only those RHS columns that are not empty.
            // This avoids unnecessary accesses to non-contiguous memory that will always be zero.
            // To implement this, we'd need to do another pass through the RHS matrix to find the empty columns and store their IDs.
            // Such extra complexity is probably not worth it; we're not in the hot loop,
            // and accessing non-consecutive IDs may be a pessimization if it prevents compiler optimizations and/or strided prefetching.
            // Ensuring correct zeroing of the output columns also becomes a bit complicated if we skip empty RHS columns.
            for (RightIndex_ rc = 0; rc < right_NC; ++rc) {
                output[sanisizer::nd_offset<std::size_t>(start + lr, left_NR, rc)] = tmp_row[rc];
            }
            std::fill(tmp_row.begin(), tmp_row.end(), 0);
        }
    }, left_NR, options.num_threads);
}
 
}
 
#endif