Accept a RHS tatami::Matrix in the public dense_row x dense matrix functions.

This matches the dense row x dense matrix function signature in the dispatcher.

Also cleaned up some of the variable names for clarity, and replace std::min
with sanisizer::min to avoid hard-coding the expected input type.

Author: LTLA

include/tatami_mult/dense_matrix/dense_row/column_to_column.hpp

@@ -58,44 +58,47 @@ struct MultiplyDenseRowWithDenseColumnMatrixToColumnOutputOptions {
 template<std::size_t accumulators_ = 4, typename LeftValue_, typename LeftIndex_, typename RightValue_, typename RightIndex_, typename Output_>
 void multiply_dense_row_with_dense_column_matrix_to_column_output(
     const tatami::Matrix<LeftValue_, LeftIndex_>& left,
-    const std::vector<const RightValue_*>& right,
-    const RightIndex_ num_rhs,
+    const tatami::Matrix<RightValue_, RightIndex_>& right,
     Output_* const output,
     const MultiplyDenseRowWithDenseColumnMatrixToColumnOutputOptions& options
 ) {
-    const auto NR = left.nrow();
-    const auto NC = left.ncol();
+    const auto left_NR = left.nrow();
+    const auto common_dim = left.ncol();
+    const auto right_NC = right.ncol();
+
+    auto right_buffers = tatami::create_container_of_Index_size<std::vector<std::vector<RightValue_> > >(right_NC);
+    auto right_ptrs = tatami::create_container_of_Index_size<std::vector<const RightValue_*> >(right_NC);
+    populate_dense_buffers(false, right_NC, common_dim, right, right_buffers, right_ptrs, options.num_threads);
 
     if (options.primary_block_size == 1) {
         tatami::parallelize([&](int, LeftIndex_ start, LeftIndex_ length) -> void {
             auto lext = tatami::consecutive_extractor<false>(left, true, start, length);
-            auto lbuffer = tatami::create_container_of_Index_size<std::vector<LeftValue_> >(NC);
-            for (LeftIndex_ r = 0; r < length; ++r) {
+            auto lbuffer = tatami::create_container_of_Index_size<std::vector<LeftValue_> >(common_dim);
+            for (LeftIndex_ lr = 0; lr < length; ++lr) {
                 const auto lptr = lext->fetch(lbuffer.data());
-                for (RightIndex_ h = 0; h < num_rhs; ++h) {
-                    output[sanisizer::nd_offset<std::size_t>(start + r, NR, h)] = dense_dot_product<accumulators_>(NC, lptr, right[h], static_cast<Output_>(0));
+                for (RightIndex_ rc = 0; rc < right_NC; ++rc) {
+                    // cast of common_dim to size_t is safe due to tatami's contract.
+                    const auto val = dense_dot_product<accumulators_>(common_dim, lptr, right_ptrs[rc], static_cast<Output_>(0));
+                    output[sanisizer::nd_offset<std::size_t>(start + lr, left_NR, rc)] = val;
                 }
             }
-        }, NR, options.num_threads);
+        }, left_NR, options.num_threads);
         return;
     } 
 
-    const auto primary_block_size = sanisizer::cast<LeftIndex_>(options.primary_block_size);
-    const auto secondary_block_size = sanisizer::cast<LeftIndex_>(options.secondary_block_size);
-
     const bool do_parallel = options.num_threads > 1;
     if (!do_parallel) {
-        std::fill_n(output, sanisizer::product_unsafe<std::size_t>(NR, num_rhs), 0);
+        std::fill_n(output, sanisizer::product_unsafe<std::size_t>(left_NR, right_NC), 0);
     }
 
     tatami::parallelize([&](int, LeftIndex_ start, LeftIndex_ length) -> void {
         auto ext = tatami::consecutive_extractor<false>(left, true, start, length);
 
-        const LeftIndex_ max_block_rows = std::min(length, primary_block_size);
+        const LeftIndex_ max_block_rows = sanisizer::min(length, options.primary_block_size);
         std::vector<std::vector<LeftValue_> > lbuffers;
         lbuffers.reserve(max_block_rows);
         for (LeftIndex_ b = 0; b < max_block_rows; ++b) {
-            lbuffers.emplace_back(tatami::cast_Index_to_container_size<std::vector<LeftValue_> >(NC));
+            lbuffers.emplace_back(tatami::cast_Index_to_container_size<std::vector<LeftValue_> >(common_dim));
         }
         auto lptrs = tatami::create_container_of_Index_size<std::vector<const LeftValue_*> >(max_block_rows);
 
@@ -106,70 +109,71 @@ void multiply_dense_row_with_dense_column_matrix_to_column_output(
             // This aims to mitigate false sharing as we update each block's partial dot products in the loop over the common dimension.
             // There is still some potential for false sharing when we transfer the results to the output buffers,
             // but this is the same as the unblocked case so we won't worry about it.
-            const auto max_block_cols = sanisizer::min(num_rhs, primary_block_size);
-            const auto max_block_rows = std::min(length, primary_block_size);
+            const RightIndex_ max_block_cols = sanisizer::min(right_NC, options.primary_block_size);
+            const LeftIndex_ max_block_rows = sanisizer::min(length, options.primary_block_size);
             tmp_output.emplace(sanisizer::product<I<decltype(tmp_output->size())> >(max_block_cols, max_block_rows));
             optr = tmp_output->data();
         } else {
             optr = output; // no need to add 'start' as this is zero in the single-threaded case.
         }
 
-        LeftIndex_ r = 0;
-        while (r < length) {
-            const LeftIndex_ rnum = std::min<LeftIndex_>(primary_block_size, length - r);
-            for (LeftIndex_ rcounter = 0; rcounter < rnum; ++rcounter) {
-                lptrs[rcounter] = ext->fetch(lbuffers[rcounter].data());
+        LeftIndex_ lr = 0;
+        while (lr < length) {
+            const LeftIndex_ lr_num = sanisizer::min(options.primary_block_size, length - lr);
+            for (LeftIndex_ lr_counter = 0; lr_counter < lr_num; ++lr_counter) {
+                lptrs[lr_counter] = ext->fetch(lbuffers[lr_counter].data());
             }
 
-            RightIndex_ h = 0;
-            while (h < num_rhs) {
-                const RightIndex_ hnum = sanisizer::min(primary_block_size, num_rhs - h);
+            RightIndex_ rc = 0;
+            while (rc < right_NC) {
+                const RightIndex_ rc_num = sanisizer::min(options.primary_block_size, right_NC - rc);
 
-                LeftIndex_ row_offset, col_offset, stride;
+                LeftIndex_ out_row_offset, out_stride;
+                RightIndex_ out_col_offset;
                 if (do_parallel) {
-                    std::fill_n(optr, sanisizer::product_unsafe<std::size_t>(hnum, rnum), 0);
-                    row_offset = 0;
-                    col_offset = 0;
-                    stride = rnum;
+                    std::fill_n(optr, sanisizer::product_unsafe<std::size_t>(rc_num, lr_num), 0);
+                    out_row_offset = 0;
+                    out_col_offset = 0;
+                    out_stride = lr_num;
                 } else {
-                    row_offset = r; // again, no need to add start, as it should be zero if there's only one thread.
-                    col_offset = h;
-                    stride = NR;
+                    out_row_offset = lr; // again, no need to add start, as it should be zero if there's only one thread.
+                    out_col_offset = rc;
+                    out_stride = left_NR;
                 }
 
-                LeftIndex_ c = 0;
-                while (c < NC) {
-                    const LeftIndex_ cnum = std::min<LeftIndex_>(secondary_block_size, NC - c);
-                    for (RightIndex_ hcounter = 0; hcounter < hnum; ++hcounter) {
-                        const auto& rightcol = right[h + hcounter];
-                        for (LeftIndex_ rcounter = 0; rcounter < rnum; ++rcounter) {
-                            auto& dest = optr[sanisizer::nd_offset<std::size_t>(row_offset + rcounter, stride, col_offset + hcounter)]; 
+                LeftIndex_ cd = 0;
+                while (cd < common_dim) {
+                    const LeftIndex_ cd_num = sanisizer::min(options.secondary_block_size, common_dim - cd);
+                    for (RightIndex_ rc_counter = 0; rc_counter < rc_num; ++rc_counter) {
+                        const auto& rightcol = right_ptrs[rc + rc_counter];
+                        for (LeftIndex_ lr_counter = 0; lr_counter < lr_num; ++lr_counter) {
+                            auto& dest = optr[sanisizer::nd_offset<std::size_t>(out_row_offset + lr_counter, out_stride, out_col_offset + rc_counter)]; 
                             dest = dense_dot_product<accumulators_>(
-                                cnum,
-                                rightcol + c,
-                                lptrs[rcounter] + c,
+                                cd_num, // cast to size_t is safe due to tatami's contract.
+                                rightcol + cd,
+                                lptrs[lr_counter] + cd,
                                 dest
                             );
                         }
                     }
-                    c += cnum;
+                    cd += cd_num;
                 }
 
                 if (do_parallel) {
-                    for (RightIndex_ hcounter = 0; hcounter < hnum; ++hcounter) {
+                    for (RightIndex_ rc_counter = 0; rc_counter < rc_num; ++rc_counter) {
                         std::copy_n(
-                            optr + sanisizer::product_unsafe<std::size_t>(rnum, hcounter),
-                            rnum,
-                            output + sanisizer::nd_offset<std::size_t>(start + r, NR, h + hcounter)
+                            optr + sanisizer::product_unsafe<std::size_t>(lr_num, rc_counter),
+                            lr_num,
+                            output + sanisizer::nd_offset<std::size_t>(start + lr, left_NR, rc + rc_counter)
                         );
                     }
                 }
 
-                h += hnum;
+                rc += rc_num;
             }
-            r += rnum;
+            lr += lr_num;
         }
-    }, NR, options.num_threads);
+    }, left_NR, options.num_threads);
 }
 
 }

include/tatami_mult/dense_matrix/dense_row/column_to_row.hpp

@@ -59,44 +59,47 @@ struct MultiplyDenseRowWithDenseColumnMatrixToRowOutputOptions {
 template<std::size_t accumulators_ = 4, typename LeftValue_, typename LeftIndex_, typename RightValue_, typename RightIndex_, typename Output_>
 void multiply_dense_row_with_dense_column_matrix_to_row_output(
     const tatami::Matrix<LeftValue_, LeftIndex_>& left,
-    const std::vector<const RightValue_*>& right,
-    const RightIndex_ num_rhs,
+    const tatami::Matrix<RightValue_, RightIndex_>& right,
     Output_* const output,
     const MultiplyDenseRowWithDenseColumnMatrixToRowOutputOptions& options
 ) {
-    const auto NR = left.nrow();
-    const auto NC = left.ncol();
+    const auto left_NR = left.nrow();
+    const auto common_dim = left.ncol();
+    const auto right_NC = right.ncol();
+
+    auto right_buffers = tatami::create_container_of_Index_size<std::vector<std::vector<RightValue_> > >(right_NC);
+    auto right_ptrs = tatami::create_container_of_Index_size<std::vector<const RightValue_*> >(right_NC);
+    populate_dense_buffers(false, right_NC, common_dim, right, right_buffers, right_ptrs, options.num_threads);
 
     if (options.primary_block_size == 1) {
         tatami::parallelize([&](int, LeftIndex_ start, LeftIndex_ length) -> void {
             auto lext = tatami::consecutive_extractor<false>(left, true, start, length);
-            auto lbuffer = tatami::create_container_of_Index_size<std::vector<LeftValue_> >(NC);
-            for (LeftIndex_ r = 0; r < length; ++r) {
+            auto lbuffer = tatami::create_container_of_Index_size<std::vector<LeftValue_> >(common_dim);
+            for (LeftIndex_ lr = 0; lr < length; ++lr) {
                 const auto lptr = lext->fetch(lbuffer.data());
-                for (RightIndex_ h = 0; h < num_rhs; ++h) {
-                    output[sanisizer::nd_offset<std::size_t>(h, num_rhs, start + r)] = dense_dot_product<accumulators_>(NC, lptr, right[h], static_cast<Output_>(0));
+                for (RightIndex_ rc = 0; rc < right_NC; ++rc) {
+                    // cast of common_dim to size_t is safe due to tatami's contract.
+                    const auto res = dense_dot_product<accumulators_>(common_dim, lptr, right_ptrs[rc], static_cast<Output_>(0));
+                    output[sanisizer::nd_offset<std::size_t>(rc, right_NC, start + lr)] = res;
                 }
             }
-        }, NR, options.num_threads);
+        }, left_NR, options.num_threads);
         return;
     } 
 
-    const auto primary_block_size = sanisizer::cast<LeftIndex_>(options.primary_block_size);
-    const auto secondary_block_size = sanisizer::cast<LeftIndex_>(options.secondary_block_size);
-
     const bool do_parallel = options.num_threads > 1;
     if (!do_parallel) {
-        std::fill_n(output, sanisizer::product_unsafe<std::size_t>(NR, num_rhs), 0);
+        std::fill_n(output, sanisizer::product_unsafe<std::size_t>(left_NR, right_NC), 0);
     }
 
     tatami::parallelize([&](int, LeftIndex_ start, LeftIndex_ length) -> void {
         auto ext = tatami::consecutive_extractor<false>(left, true, start, length);
 
-        const LeftIndex_ max_block_rows = std::min(length, primary_block_size);
+        const LeftIndex_ max_block_rows = sanisizer::min(length, options.primary_block_size);
         std::vector<std::vector<LeftValue_> > lbuffers;
         lbuffers.reserve(max_block_rows);
         for (LeftIndex_ b = 0; b < max_block_rows; ++b) {
-            lbuffers.emplace_back(tatami::cast_Index_to_container_size<std::vector<LeftValue_> >(NC));
+            lbuffers.emplace_back(tatami::cast_Index_to_container_size<std::vector<LeftValue_> >(common_dim));
         }
         auto lptrs = tatami::create_container_of_Index_size<std::vector<const LeftValue_*> >(max_block_rows);
 
@@ -107,70 +110,71 @@ void multiply_dense_row_with_dense_column_matrix_to_row_output(
             // This aims to mitigate false sharing as we update each block's partial dot products in the loop over the common dimension.
             // There is still some potential for false sharing when we transfer the results to the output buffers,
             // but this is the same as the unblocked case so we won't worry about it.
-            const auto max_block_cols = sanisizer::min(num_rhs, primary_block_size);
-            const auto max_block_rows = std::min(length, primary_block_size);
+            const auto max_block_cols = sanisizer::min(right_NC, options.primary_block_size);
+            const auto max_block_rows = sanisizer::min(length, options.primary_block_size);
             tmp_output.emplace(sanisizer::product<I<decltype(tmp_output->size())> >(max_block_cols, max_block_rows));
             optr = tmp_output->data();
         } else {
             optr = output; // no need to add 'start' as this is zero in the single-threaded case.
         }
 
-        LeftIndex_ r = 0;
-        while (r < length) {
-            const LeftIndex_ rnum = std::min<LeftIndex_>(primary_block_size, length - r);
-            for (LeftIndex_ rcounter = 0; rcounter < rnum; ++rcounter) {
-                lptrs[rcounter] = ext->fetch(lbuffers[rcounter].data());
+        LeftIndex_ lr = 0;
+        while (lr < length) {
+            const LeftIndex_ lr_num = sanisizer::min(options.primary_block_size, length - lr);
+            for (LeftIndex_ lr_counter = 0; lr_counter < lr_num; ++lr_counter) {
+                lptrs[lr_counter] = ext->fetch(lbuffers[lr_counter].data());
             }
 
-            RightIndex_ h = 0;
-            while (h < num_rhs) {
-                const RightIndex_ hnum = sanisizer::min(primary_block_size, num_rhs - h);
+            RightIndex_ rc = 0;
+            while (rc < right_NC) {
+                const RightIndex_ rc_num = sanisizer::min(options.primary_block_size, right_NC - rc);
 
-                RightIndex_ row_offset, col_offset, stride;
+                LeftIndex_ out_row_offset;
+                RightIndex_ out_col_offset, out_stride;
                 if (do_parallel) {
-                    std::fill_n(optr, sanisizer::product_unsafe<std::size_t>(hnum, rnum), 0);
-                    row_offset = 0;
-                    col_offset = 0;
-                    stride = hnum;
+                    std::fill_n(optr, sanisizer::product_unsafe<std::size_t>(rc_num, lr_num), 0);
+                    out_row_offset = 0;
+                    out_col_offset = 0;
+                    out_stride = rc_num;
                 } else {
-                    row_offset = r; // start == 0 in serial code, so no need to add start here.
-                    col_offset = h;
-                    stride = num_rhs;
+                    out_row_offset = lr; // start == 0 in serial code, so no need to add start here.
+                    out_col_offset = rc;
+                    out_stride = right_NC;
                 }
 
-                LeftIndex_ c = 0;
-                while (c < NC) {
-                    const LeftIndex_ cnum = std::min<LeftIndex_>(secondary_block_size, NC - c);
-                    for (LeftIndex_ rcounter = 0; rcounter < rnum; ++rcounter) {
-                        const auto leftrow = lptrs[rcounter];
-                        for (RightIndex_ hcounter = 0; hcounter < hnum; ++hcounter) {
-                            auto& dest = optr[sanisizer::nd_offset<std::size_t>(col_offset + hcounter, stride, row_offset + rcounter)]; 
+                LeftIndex_ cd = 0;
+                while (cd < common_dim) {
+                    const LeftIndex_ cd_num = sanisizer::min(options.secondary_block_size, common_dim - cd);
+                    for (LeftIndex_ lr_counter = 0; lr_counter < lr_num; ++lr_counter) {
+                        const auto leftrow = lptrs[lr_counter];
+                        for (RightIndex_ rc_counter = 0; rc_counter < rc_num; ++rc_counter) {
+                            auto& dest = optr[sanisizer::nd_offset<std::size_t>(out_col_offset + rc_counter, out_stride, out_row_offset + lr_counter)]; 
                             dest = dense_dot_product<accumulators_>(
-                                cnum,
-                                right[h + hcounter] + c,
-                                leftrow + c,
+                                cd_num, // cast to size_t is safe due to tatami's contract.
+                                right_ptrs[rc + rc_counter] + cd,
+                                leftrow + cd,
                                 dest
                             );
                         }
                     }
-                    c += cnum;
+                    cd += cd_num;
                 }
 
                 if (do_parallel) {
-                    for (LeftIndex_ rcounter = 0; rcounter < rnum; ++rcounter) {
+                    for (LeftIndex_ lr_counter = 0; lr_counter < lr_num; ++lr_counter) {
                         std::copy_n(
-                            optr + sanisizer::product_unsafe<std::size_t>(hnum, rcounter),
-                            hnum,
-                            output + sanisizer::nd_offset<std::size_t>(h, num_rhs, start + r + rcounter)
+                            optr + sanisizer::product_unsafe<std::size_t>(rc_num, lr_counter),
+                            rc_num,
+                            output + sanisizer::nd_offset<std::size_t>(rc, right_NC, start + lr + lr_counter)
                         );
                     }
                 }
 
-                h += hnum;
+                rc += rc_num;
             }
-            r += rnum;
+            lr += lr_num;
         }
-    }, NR, options.num_threads);
+    }, left_NR, options.num_threads);
 }
 
 }