dense_column.hpp

Go to the documentation of this file.

#ifndef TATAMI_MULT_SINGLE_VECTOR_DENSE_COLUMN_HPP
#define TATAMI_MULT_SINGLE_VECTOR_DENSE_COLUMN_HPP
 
#include <cstddef>
#include <vector>
#include <optional>
#include <algorithm>
 
#include "tatami/tatami.hpp"
#include "sanisizer/sanisizer.hpp"
 
#include "../utils.hpp"
 
namespace tatami_mult {
 
/* See https://github.com/tatami-inc/test-multiplication/tree/master/dense_column/single_vector
 * for an explanation of the choice of algorithm.
 */
 
struct MultiplyDenseColumnWithSingleVectorOptions {
    int num_threads = 1;
};
 
template<typename LeftValue_, typename LeftIndex_, typename RightValue_, typename Output_>
void multiply_dense_column_with_single_vector(
    const tatami::Matrix<LeftValue_, LeftIndex_>& left,
    const RightValue_* const right,
    Output_* const output,
    const MultiplyDenseColumnWithSingleVectorOptions& options
) {
    const auto NR = left.nrow();
    const auto NC = left.ncol();
 
    std::optional<std::vector<std::optional<std::vector<Output_> > > > tmp_results;
    const bool do_parallel = options.num_threads > 1;
    if (do_parallel) {
        tmp_results.emplace(sanisizer::cast<I<decltype(tmp_results->size())> >(options.num_threads - 1));
    }
    std::fill_n(output, NR, 0);
 
    const auto num_used = tatami::parallelize([&](int t, LeftIndex_ start, LeftIndex_ length) -> void {
        auto ext = tatami::consecutive_extractor<false>(left, false, start, length);
        auto buffer = tatami::create_container_of_Index_size<std::vector<LeftValue_> >(NR);
 
        Output_* optr;
        std::optional<std::vector<Output_> > cur_output;
        if (!do_parallel || t == 0) {
            optr = output;
        } else {
            cur_output.emplace(tatami::cast_Index_to_container_size<std::vector<Output_> >(NR));
            optr = cur_output->data();
        }
 
        for (LeftIndex_ c = 0; c < length; ++c) {
            auto ptr = ext->fetch(buffer.data());
            const Output_ mult = right[start + c];
            for (LeftIndex_ r = 0; r < NR; ++r) {
                optr[r] += mult * ptr[r];
            }
        }
 
        if (do_parallel && t > 0) {
            (*tmp_results)[t - 1] = std::move(cur_output);            
        }
    }, NC, options.num_threads);
 
    if (do_parallel) {
        for (int u = 1; u < num_used; ++u) {
            const auto& tmp = *((*tmp_results)[u - 1]);
            for (LeftIndex_ r = 0; r < NR; ++r) {
                output[r] += tmp[r];
            }
        }
    }
}
 
}
 
#endif